05. Breaking up the problem

A critical step in building any large, complex project is to break it down into smaller, more manageable tasks. This approach not only structures and streamlines the process but also enhances debugging, optimization, and scalability. Additionally, it facilitates the development of versatile and modular code that can address a variety of problems. In this section we will discuss how we broke down the problem and implemented each step in a modular fashion.

For our PDE solver, the aim is to solve equations of the form

$$\frac{\partial u }{\partial t} = \epsilon \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial x^2} \right) + F(u,v), \\ \frac{\partial v }{\partial t}= \left(\frac{\partial^2 v}{\partial x^2}+ \frac{\partial^2 v}{\partial y^2} \right) + G(u,v).$$

The key steps of achieving this are:

Step 1:

In one spatial dimension $x$, consider the boundary value problem (BVP) with solution $u(x)$:

$$A(x) \frac{\partial^2 u}{\partial x^2}+B(x) \frac{\partial u}{\partial x}+C(x) u = D(x),$$

where $A(x), B(x), C(x)$ and $D(x)$ are functions. We supplemented this equation with Dirichlet or Neumann boundary conditions. We will discuss how we discretised this equation in Section 6.

Step 2:

We then consider the effect of multiplying the second derivative by a small parameter $\epsilon$:

$$\epsilon A(x) \frac{\partial^2 u}{\partial x^2}+B(x) \frac{\partial u}{\partial x}+C(x) u = D(x).$$

The mathematics and code structure will be discussed in Section 7.

Step 3:

We then add in the nonlinearities, by including the non-linear function $F$:

$$\epsilon A(x) \frac{\partial^2 u}{\partial x^2}+B(x) \frac{\partial u}{\partial x}+C(x) u + F(u) = D(x).$$

This non-linear equation is then solved using Newton-iteration in the module Newton_method.f90.

Step 4:

We then add a second equation for $v$:

$$\epsilon A_1(x) \frac{\partial^2 u}{\partial x^2}+B_1(x) \frac{\partial u}{\partial x}+C_1(x) u + F(u,v) = D_1(x),\\ A_2(x) \frac{\partial^2 v}{\partial x^2}+B_2(x) \frac{\partial v}{\partial x}+C_2(x) u + G(u,v) = D_2(x).$$

This is easy to implement with finite differences: if our discretised solution was originally $u = (u(x_1),u(x_2),...,u(x_n))$ then the new discretised solution gets written in the form $u = (u(x_1),v(x_1),u(x_2),v(x_2),...,u(x_n),v(x_n))$.

This is achieved by defining an additional equation in src/equations/definition.f90 and then build both equations in the subroutine build_the_matrix. This highlights the modular nature of the code.

Step 5:

We include a second spatial domain $y$:

$$\epsilon A_1(x,y) (\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}) + C_1(x,y) u + F(u,v) = D_1(x,y), \\ A_2(x,y) (\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}) + C_2(x,y) v + G(u,v) = D_2(x,y).$$

We have dropped the first derivatives here for simplicity.

Importantly, this is again coded in a modular way: the new operators are constructed in src/equations/2D.f90. We built copies of the subroutines equation_setup, build_the_matrix and non_linear_setup and constructed the two-dimensional variants equation_setup2D, build_the_matrix2D and non_linear_setup2D. The logic of both sets of subroutines are similar - except the new use of the Kronecker product in equation_setup2D. As an extension it is possible to combine these subroutines - however, we decided to improve readability by keeping them separate.

Step 6:

The final step is to add in the temporal marching and solve the entire IBVP described above. All the required modules and subroutines have already been built. A new temporal marching matrix is defined in the subroutine src/solve_ibvp.f90: implicit_march.