Solving Linear Pdes
This tutorial is automatically generated from TestSolvingLinearPdesTutorial.hpp at revision 86da82f2455f. 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
#include <cxxtest/TestSuite.h>
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.
#include "UblasIncludes.hpp"
This is the class that is needed to solve a linear elliptic PDE.
#include "SimpleLinearEllipticSolver.hpp"
This is the class that is needed to solve a linear parabolic PDE.
#include "SimpleLinearParabolicSolver.hpp"
This is a parabolic PDE, one of the PDEs we will solve.
#include "HeatEquationWithSourceTerm.hpp"
We will also solve this PDE.
#include "SimplePoissonEquation.hpp"
This is needed to read mesh datafiles of the ‘Triangles’ format.
#include "TrianglesMeshReader.hpp"
This class represents the mesh internally.
#include "TetrahedralMesh.hpp"
These are used to specify boundary conditions for the PDEs.
#include "BoundaryConditionsContainer.hpp"
#include "ConstBoundaryCondition.hpp"
This class helps us deal with output files.
#include "OutputFileHandler.hpp"
The following header must be included in every test that uses PETSc. Note that it cannot be included in the source code.
#include "PetscSetupAndFinalize.hpp"
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.
class MyPde : public AbstractLinearEllipticPde<2,2>
{
private:
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
).
c_matrix<double,2,2> mDiffusionTensor;
public:
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
MyPde()
{
mDiffusionTensor(0,0) = 2.0;
mDiffusionTensor(0,1) = 0.0;
mDiffusionTensor(1,0) = 0.0;
mDiffusionTensor(1,1) = 1.0;
}
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$.
double ComputeConstantInUSourceTerm(const ChastePoint<2>& rX, Element<2,2>* pElement)
{
return rX[0]*rX[0] + rX[1]*rX[1];
}
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.
double ComputeLinearInUCoeffInSourceTerm(const ChastePoint<2>& rX, Element<2,2>* pElement)
{
return 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.
c_matrix<double,2,2> ComputeDiffusionTerm(const ChastePoint<2>& rX)
{
return mDiffusionTensor;
}
};
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
.
class TestSolvingLinearPdesTutorial : public CxxTest::TestSuite
{
All individual test defined in this test suite must be declared as public.
public:
void TestSolvingEllipticPde()
{
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.
TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_128_elements");
Now declare a tetrahedral mesh with the same dimensions…
TetrahedralMesh<2,2> mesh;
… and construct the mesh using the mesh reader.
mesh.ConstructFromMeshReader(mesh_reader);
Next we instantiate an instance of our PDE we wish to solve.
MyPde pde;
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.
BoundaryConditionsContainer<2,2,1> bcc;
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:
//bcc.DefineZeroDirichletOnMeshBoundary(&mesh);
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.
ConstBoundaryCondition<2>* p_zero_boundary_condition = new ConstBoundaryCondition<2>(0.0);
We then get a boundary node iterator from the mesh…
TetrahedralMesh<2,2>::BoundaryNodeIterator iter = mesh.GetBoundaryNodeIteratorBegin();
…and loop over the boundary nodes, getting the x and y values.
while (iter < mesh.GetBoundaryNodeIteratorEnd())
{
double x = (*iter)->GetPoint()[0];
double y = (*iter)->GetPoint()[1];
If x=0 or y=0…
if ((x==0) || (y==0))
{
…associate the zero boundary condition created above with this boundary node
(*iter
being a pointer to a Node<2>
).
bcc.AddDirichletBoundaryCondition(*iter, p_zero_boundary_condition);
}
iter++;
}
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.
TetrahedralMesh<2,2>::BoundaryElementIterator surf_iter
= mesh.GetBoundaryElementIteratorBegin();
while (surf_iter != mesh.GetBoundaryElementIteratorEnd())
{
Get the x and y values of any node (here, the 0th) in the element.
unsigned node_index = (*surf_iter)->GetNodeGlobalIndex(0);
double x = mesh.GetNode(node_index)->GetPoint()[0];
double y = mesh.GetNode(node_index)->GetPoint()[1];
If x=1 or y=1…
if ((fabs(x-1.0) < 1e-6) || (fabs(y-1.0) < 1e-6))
{
…associate the boundary condition with the surface element.
bcc.AddNeumannBoundaryCondition(*surf_iter, p_zero_boundary_condition);
}
Finally increment the iterator.
surf_iter++;
}
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.
SimpleLinearEllipticSolver<2,2> solver(&mesh, &pde, &bcc);
To solve, just call Solve()
. A PETSc vector is returned.
Vec result = solver.Solve();
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]
.
ReplicatableVector result_repl(result);
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.
OutputFileHandler output_file_handler("TestSolvingLinearPdeTutorial");
Create an out_stream
, which is a stream to a particular file. An out_stream
is a smart pointer to a std::ofstream
.
out_stream p_file = output_file_handler.OpenOutputFile("linear_solution.txt");
Loop over the entries of the solution.
for (unsigned i=0; i<result_repl.GetSize(); i++)
{
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
.
double x = mesh.GetNode(i)->rGetLocation()[0];
double y = mesh.GetNode(i)->rGetLocation()[1];
Get the computed solution at this node from the ReplicatableVector
.
double u = result_repl[i];
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),'.');
(*p_file) << x << " " << y << " " << u << "\n";
}
All PETSc Vec
s should be destroyed when they are no longer needed, or you will have a memory leak.
PetscTools::Destroy(result);
}
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$.
void TestSolvingParabolicPde()
{
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.
TetrahedralMesh<3,3> mesh;
mesh.ConstructRegularSlabMesh(0.1, 1.0, 1.0, 1.0);
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).
HeatEquationWithSourceTerm<3> pde;
Create a new boundary conditions container and specify $u=1.0$ on the boundary.
BoundaryConditionsContainer<3,3,1> bcc;
bcc.DefineConstantDirichletOnMeshBoundary(&mesh, 1.0);
Create an instance of the solver, passing in the mesh, pde and boundary conditions.
SimpleLinearParabolicSolver<3,3> solver(&mesh,&pde,&bcc);
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.
Vec initial_condition = PetscTools::CreateAndSetVec(mesh.GetNumNodes(), 1.0);
solver.SetInitialCondition(initial_condition);
Next define the start time, end time, and timestep, and set them.
double t_start = 0;
double t_end = 1;
double dt = 0.01;
solver.SetTimes(t_start, t_end);
solver.SetTimeStep(dt);
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:
solver.SetOutputDirectoryAndPrefix("ParabolicSolverTutorial","results");
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.)
solver.SetOutputToTxt(true);
solver.SetPrintingTimestepMultiple(10);
Now we can solve the problem. The Vec
that is returned can be passed into a
ReplicatableVector
as before.
Vec solution = solver.Solve();
ReplicatableVector solution_repl(solution);
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.
SimplePoissonEquation<3,3> static_pde;
SimpleLinearEllipticSolver<3,3> static_solver(&mesh, &static_pde, &bcc);
Vec static_solution = static_solver.Solve();
ReplicatableVector static_solution_repl(static_solution);
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.)
for (unsigned i=0; i<static_solution_repl.GetSize(); i++)
{
TS_ASSERT_DELTA( solution_repl[i], static_solution_repl[i], 1e-3);
}
All PETSc vectors should be destroyed when they are no longer needed.
PetscTools::Destroy(initial_condition);
PetscTools::Destroy(solution);
PetscTools::Destroy(static_solution);
}
};
Full code
#include <cxxtest/TestSuite.h>
#include "UblasIncludes.hpp"
#include "SimpleLinearEllipticSolver.hpp"
#include "SimpleLinearParabolicSolver.hpp"
#include "HeatEquationWithSourceTerm.hpp"
#include "SimplePoissonEquation.hpp"
#include "TrianglesMeshReader.hpp"
#include "TetrahedralMesh.hpp"
#include "BoundaryConditionsContainer.hpp"
#include "ConstBoundaryCondition.hpp"
#include "OutputFileHandler.hpp"
#include "PetscSetupAndFinalize.hpp"
class MyPde : public AbstractLinearEllipticPde<2,2>
{
private:
c_matrix<double,2,2> mDiffusionTensor;
public:
MyPde()
{
mDiffusionTensor(0,0) = 2.0;
mDiffusionTensor(0,1) = 0.0;
mDiffusionTensor(1,0) = 0.0;
mDiffusionTensor(1,1) = 1.0;
}
double ComputeConstantInUSourceTerm(const ChastePoint<2>& rX, Element<2,2>* pElement)
{
return rX[0]*rX[0] + rX[1]*rX[1];
}
double ComputeLinearInUCoeffInSourceTerm(const ChastePoint<2>& rX, Element<2,2>* pElement)
{
return 1.0;
}
c_matrix<double,2,2> ComputeDiffusionTerm(const ChastePoint<2>& rX)
{
return mDiffusionTensor;
}
};
class TestSolvingLinearPdesTutorial : public CxxTest::TestSuite
{
public:
void TestSolvingEllipticPde()
{
TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_128_elements");
TetrahedralMesh<2,2> mesh;
mesh.ConstructFromMeshReader(mesh_reader);
MyPde pde;
BoundaryConditionsContainer<2,2,1> bcc;
//bcc.DefineZeroDirichletOnMeshBoundary(&mesh);
ConstBoundaryCondition<2>* p_zero_boundary_condition = new ConstBoundaryCondition<2>(0.0);
TetrahedralMesh<2,2>::BoundaryNodeIterator iter = mesh.GetBoundaryNodeIteratorBegin();
while (iter < mesh.GetBoundaryNodeIteratorEnd())
{
double x = (*iter)->GetPoint()[0];
double y = (*iter)->GetPoint()[1];
if ((x==0) || (y==0))
{
bcc.AddDirichletBoundaryCondition(*iter, p_zero_boundary_condition);
}
iter++;
}
TetrahedralMesh<2,2>::BoundaryElementIterator surf_iter
= mesh.GetBoundaryElementIteratorBegin();
while (surf_iter != mesh.GetBoundaryElementIteratorEnd())
{
unsigned node_index = (*surf_iter)->GetNodeGlobalIndex(0);
double x = mesh.GetNode(node_index)->GetPoint()[0];
double y = mesh.GetNode(node_index)->GetPoint()[1];
if ((fabs(x-1.0) < 1e-6) || (fabs(y-1.0) < 1e-6))
{
bcc.AddNeumannBoundaryCondition(*surf_iter, p_zero_boundary_condition);
}
surf_iter++;
}
SimpleLinearEllipticSolver<2,2> solver(&mesh, &pde, &bcc);
Vec result = solver.Solve();
ReplicatableVector result_repl(result);
OutputFileHandler output_file_handler("TestSolvingLinearPdeTutorial");
out_stream p_file = output_file_handler.OpenOutputFile("linear_solution.txt");
for (unsigned i=0; i<result_repl.GetSize(); i++)
{
double x = mesh.GetNode(i)->rGetLocation()[0];
double y = mesh.GetNode(i)->rGetLocation()[1];
double u = result_repl[i];
(*p_file) << x << " " << y << " " << u << "\n";
}
PetscTools::Destroy(result);
}
void TestSolvingParabolicPde()
{
TetrahedralMesh<3,3> mesh;
mesh.ConstructRegularSlabMesh(0.1, 1.0, 1.0, 1.0);
HeatEquationWithSourceTerm<3> pde;
BoundaryConditionsContainer<3,3,1> bcc;
bcc.DefineConstantDirichletOnMeshBoundary(&mesh, 1.0);
SimpleLinearParabolicSolver<3,3> solver(&mesh,&pde,&bcc);
Vec initial_condition = PetscTools::CreateAndSetVec(mesh.GetNumNodes(), 1.0);
solver.SetInitialCondition(initial_condition);
double t_start = 0;
double t_end = 1;
double dt = 0.01;
solver.SetTimes(t_start, t_end);
solver.SetTimeStep(dt);
solver.SetOutputDirectoryAndPrefix("ParabolicSolverTutorial","results");
solver.SetOutputToTxt(true);
solver.SetPrintingTimestepMultiple(10);
Vec solution = solver.Solve();
ReplicatableVector solution_repl(solution);
SimplePoissonEquation<3,3> static_pde;
SimpleLinearEllipticSolver<3,3> static_solver(&mesh, &static_pde, &bcc);
Vec static_solution = static_solver.Solve();
ReplicatableVector static_solution_repl(static_solution);
for (unsigned i=0; i<static_solution_repl.GetSize(); i++)
{
TS_ASSERT_DELTA( solution_repl[i], static_solution_repl[i], 1e-3);
}
PetscTools::Destroy(initial_condition);
PetscTools::Destroy(solution);
PetscTools::Destroy(static_solution);
}
};