Documentation for Release 2024.1

# Solving Linear Pdes

This tutorial is automatically generated from TestSolvingLinearPdesTutorial.hpp at revision 71ac325b969f. Note that the code is given in full at the bottom of the page.

## Examples showing how to solve linear elliptic and parabolic PDEs

In this tutorial we show how Chaste can be used to solve linear PDEs. The first test
uses the `SimpleLinearEllipticSolver`

to solve a linear elliptic PDE, and the
second test uses the `SimpleLinearParabolicSolver`

to solve a parabolic time-dependent
linear PDE.

The following header files need to be included. First we include the header needed to define this class as a test suite

On some systems (with Boost 1.33.1 and PETSc 2.2) there is a clash between Boost Ublas includes and PETSc. This can be resolved by making sure that Chaste’s interface to the Ublas library is included as early as possible.

This is the class that is needed to solve a linear elliptic PDE.

This is the class that is needed to solve a linear parabolic PDE.

This is a parabolic PDE, one of the PDEs we will solve.

We will also solve this PDE.

This is needed to read mesh datafiles of the ‘Triangles’ format.

This class represents the mesh internally.

These are used to specify boundary conditions for the PDEs.

This class helps us deal with output files.

The following header must be included in every test that uses PETSc. Note that it cannot be included in the source code.

### Test 1: Solving a linear elliptic PDE

Here, we solve the PDE: $\nabla. (D \nabla u) + u + x^2+y^2 = 0$, in 2D, where D is the diffusion tensor (2 0; 0 1) (ie $D_{11}=2$, $D_{12}=D_{21}=0$, $D_{22}=1$), on a square domain, with boundary conditions $u=0$ on $x=0$ or $y=0$, and $(D \nabla u).\mathbf{n} = 0$ on $x=1$ and $y=1$, where $\mathbf{n}$ is the surface normal.

We need to create a class representing the PDE we want to solve, which will be
passed into the solver. The PDE we are solving is of the type
`AbstractLinearEllipticPde`

, which is an abstract class with 3 pure methods
which have to be implemented. The template variables in the following line are both the dimension
of the space.

For efficiency, we will save the diffusion tensor that will be returned by one of the
class’ methods as a member variable. The diffusion tensor which has to be returned
by the `GetDiffusionTensor`

method in PDE classes is of the type
`c_matrix<double,SIZE,SIZE>`

, which is a uBLAS matrix. We use uBLAS vectors
and matrices where small vectors and matrices are needed. Note that uBLAS objects
are only particularly efficient if optimisation is on (`CMAKE_BUILD_TYPE=Release`

).

The constructor just sets up the diffusion tensor. We choose a diffusion tensor which corresponds to twice as much diffusion in the x-direction compared to the y-direction

The first method which has to be implemented returns the constant (not dependent on u) part of the source term, which for our PDE is $x^2 + y^2$.

The second method which has to be implemented returns the coefficient in the linear-in-u part of the source term, which for our PDE is just 1.0.

The third method returns the diffusion tensor $D$. Note that the diffusion tensor should be symmetric and positive definite for a physical, well-posed problem.

Next, we define the test suite (a class). It is sensible to name it the same
as the filename. The class should inherit from `CxxTest::TestSuite`

.

All individual test defined in this test suite **must** be declared as public.

First we declare a mesh reader which reads mesh data files of the ‘Triangle’
format. The path given is relative to the main Chaste directory. As we are in 2d,
the reader will look for three datafiles, [name].nodes, [name].ele and [name].edge.
Note that the first template argument here is the spatial dimension of the
elements in the mesh (`ELEMENT_DIM`

), and the second is the dimension of the nodes,
i.e. the dimension of the space the mesh lives in (`SPACE_DIM`

). Usually
`ELEMENT_DIM`

and `SPACE_DIM`

will be equal.

Now declare a tetrahedral mesh with the same dimensions…

… and construct the mesh using the mesh reader.

Next we instantiate an instance of our PDE we wish to solve.

A set of boundary conditions are stored in a `BoundaryConditionsContainer`

. The
three template arguments are ELEMENT_DIM, SPACE_DIM and PROBLEM_DIM, the latter being
the number of unknowns we are solving for. We have one unknown (ie $u$ is a scalar, not
a vector), so in this case `PROBLEM_DIM`

=1.

Defining the boundary conditions is the only particularly fiddly part of solving PDEs, unless they are very simple, such as $u=0$ on the boundary, which could be done as follows:

We want to specify $u=0$ on $x=0$ and $y=0$. To do this, we first create the boundary condition
object saying what the value of the condition is at any particular point in space. Here
we use the class `ConstBoundaryCondition`

, a subclass of `AbstractBoundaryCondition`

that
yields the same constant value (0.0 here) everywhere it is used.

Note that the object is allocated with `new`

. The `BoundaryConditionsContainer`

object deals
with deleting its associated boundary condition objects. Note too that we could allocate a
separate condition object for each boundary node, but using the same object where possible is
more memory efficient.

We then get a boundary node iterator from the mesh…

…and loop over the boundary nodes, getting the x and y values.

If x=0 or y=0…

…associate the zero boundary condition created above with this boundary node
(`*iter`

being a pointer to a `Node<2>`

).

Now we create Neumann boundary conditions for the *surface elements* on x=1 and y=1. Note that
Dirichlet boundary conditions are defined on nodes, whereas Neumann boundary conditions are
defined on surface elements. Note also that the natural boundary condition statement for this
PDE is $(D \nabla u).\mathbf{n} = g(x)$ (where $\mathbf{n}$ is the outward-facing surface normal), and $g(x)$ is a prescribed
function, *not* something like $\partial u/ \partial n=g(x)$. Hence the boundary condition we are specifying is
$(D \nabla u).\mathbf{n} = 0$.

Note for 1D

If we were solving $2u_{xx}=f(x)$ in 1D, and wanted to specify $\partial u/ \partial x=1$ on the LHS boundary, the Neumann boundary value we have to specify is $-2$, as $n=-1$ (outward facing normal) so $(D \nabla u).n = -2$ when $\partial u/ \partial x=1$.

To define Neumann bcs, we reuse the zero boundary condition object defined above, but apply it at surface elements. We loop over these using another iterator provided by the mesh class.

Get the x and y values of any node (here, the 0th) in the element.

If x=1 or y=1…

…associate the boundary condition with the surface element.

Finally increment the iterator.

Next we define the solver of the PDE.
To solve an `AbstractLinearEllipticPde`

(which is the type of PDE `MyPde`

is),
we use a `SimpleLinearEllipticSolver`

. The solver, again templated over
`ELEMENT_DIM`

and `SPACE_DIM`

, needs to be given (pointers to) the mesh,
pde and boundary conditions.

To solve, just call `Solve()`

. A PETSc vector is returned.

It is a pain to access the individual components of a PETSc vector, even when running only on
one process. A helper class called `ReplicatableVector`

has been created. Create
an instance of one of these, using the PETSc `Vec`

as the data. The $i$th
component of `result`

can now be obtained by simply doing `result_repl[i]`

.

Let us write out the solution to a file. To do this, create an
`OutputFileHandler`

, passing in the directory we want files written to.
This is relative to the directory defined by the `CHASTE_TEST_OUTPUT`

environment
variable - usually `/tmp/$USER/testoutput`

. Note by default the output directory
passed in is emptied by this command. To avoid this, `false`

can be passed in as a second
parameter.

Create an `out_stream`

, which is a stream to a particular file. An `out_stream`

is a smart pointer to a `std::ofstream`

.

Loop over the entries of the solution.

Get the $x$ and $y$-values of the node corresponding to this entry. The method
`GetNode`

on the mesh class returns a pointer to a `Node`

.

Get the computed solution at this node from the `ReplicatableVector`

.

Finally, write $x$, $y$ and $u$ to the output file. The solution could then be
visualised in (eg) matlab, using the commands:
`sol=load('linear_solution.txt'); plot3(sol(:,1),sol(:,2),sol(:,3),'.');`

All PETSc `Vec`

s should be destroyed when they are no longer needed, or you will have a memory leak.

### Test 2: Solving a linear parabolic PDE

Now we solve a parabolic PDE. We choose a simple problem so that the code changes needed from the elliptic case are clearer. We will solve $\frac{\partial u}{\partial t} = \nabla . (\nabla u) + u$, in 3d, with boundary conditions $u=1$ on the boundary, and initial conditions $u=1$.

Create a 10 by 10 by 10 mesh in 3D, this time using the `ConstructRegularSlabMesh`

method
on the mesh. The first parameter is the cartesian space-step and the other three parameters are the width, height and depth of the mesh.

Our PDE object should be a class that is derived from the `AbstractLinearParabolicPde`

.
We could write it ourselves as in the previous test, but since the PDE we want to solve is
so simple, it has already been defined (look it up! - it is located in pde/test/pdes).

Create a new boundary conditions container and specify $u=1.0$ on the boundary.

Create an instance of the solver, passing in the mesh, pde and boundary conditions.

For parabolic problems, initial conditions are also needed. The solver will expect
a PETSc vector, where the $i$-th entry is the initial solution at node $i$, to be passed
in. To create this PETSc `Vec`

, we will use a helper function in the `PetscTools`

class to create a `Vec`

of size num_nodes, with each entry set to `1.0`

. Then we
set the initial condition on the solver.

Next define the start time, end time, and timestep, and set them.

When we call `Solve()`

below we will just get the solution at the final time. If we want
to have intermediate solutions written to file, we do the following. We start by
specifying an output directory and filename prefix for our results file:

When an output directory has been specified, the solver writes output in HDF5 format. To
convert this to another output format, we call the relevant method. Here, we convert
the output to plain text files (VTK or cmgui formats are also possible). We also say how
often to write the data, telling the solver to output results to file every tenth timestep.
The solver will create one file for each variable (in this case there is only one variable)
and for each time, so for example, the file
`results_Variable_0_10`

is the results for u, over all nodes, at the 11th printed time.
Have a look in the output directory after running the test. (For comments on visualising the data in
matlab or octave, see the end of the tutorial UserTutorials/WritingPdeSolvers.)

Now we can solve the problem. The `Vec`

that is returned can be passed into a
`ReplicatableVector`

as before.

Let’s also solve the equivalent static PDE, i.e. set $\partial u/ \partial t=0$, so $0=\nabla . (\nabla u) + u$. This is easy, as the PDE class has already been defined.

We can now compare the solution of the parabolic PDE at $t=1$ with the static solution, to see if the static equilibrium solution was reached in the former. (Ideally we should compute some relative error, but we just compute an absolute error for simplicity.)

All PETSc vectors should be destroyed when they are no longer needed.