| Chaste

This tutorial was generated from the file projects/Microvessel/test/tutorials/TestSolvingPdesLiteratePaper.hpp at revision r27265. Note that the code is given in full at the bottom of the page.

Solving PDEs in the Microvessel Project

This tutorial demonstrates methods for solving PDEs in the Microvessel Project. It is noted that the way to set up PDEs differs from that of Cell Based Chaste, although the same solver can be used behind the scenes.

The following is covered:

Solving a linear reaction-diffusion PDE with finite differences.
Solving a linear reaction-diffusion PDE with finite finite elements and discrete sinks and sources.
Solving a non-linear reaction-diffusion PDE with finite differences.
Solving a linear reaction-diffusion PDE using Green’s functions.
Interacting with regular grids and finite element meshes.

The Test

Start by introducing the necessary header files. The first contain functionality for setting up unit tests, smart pointer tools and output management.

#include <cxxtest/TestSuite.h>
#include "AbstractCellBasedWithTimingsTestSuite.hpp"
#include "SmartPointers.hpp"
#include "OutputFileHandler.hpp"
#include "FileFinder.hpp"

Dimensional analysis.

#include "DimensionalChastePoint.hpp"
#include "UnitCollection.hpp"
#include "Owen11Parameters.hpp"
#include "GenericParameters.hpp"
#include "ParameterCollection.hpp"
#include "BaseUnits.hpp"

Geometry tools.

#include "MappableGridGenerator.hpp"
#include "Part.hpp"

Vessel networks.

#include "VesselNetwork.hpp"
#include "VesselNetworkGenerator.hpp"

Grids and PDEs.

#include "DiscreteContinuumMesh.hpp"
#include "VtkMeshWriter.hpp"
#include "FiniteElementSolver.hpp"
#include "FiniteDifferenceSolver.hpp"
#include "DiscreteSource.hpp"
#include "VesselBasedDiscreteSource.hpp"
#include "DiscreteContinuumBoundaryCondition.hpp"
#include "LinearSteadyStateDiffusionReactionPde.hpp"
#include "MichaelisMentenSteadyStateDiffusionReactionPde.hpp"

This should appear last.

#include "PetscSetupAndFinalize.hpp"
class TestSolvingPdesLiteratePaper : public AbstractCellBasedWithTimingsTestSuite
{
public:

Test 1 - Linear Reaction Diffusion With Finite Differences

In the first example we will solve a steady-state linear reaction diffusion PDE with finite differences.

    void TestLinearReactionDiffusionPdeWithFiniteDifferences() throw(Exception)
    {
        MAKE_PTR_ARGS(OutputFileHandler, p_handler, ("TestSolvingPdesLiteratePaper/TestLinearReactionDiffusionPdeWithFiniteDifferences"));

We will work in microns

        units::quantity<unit::length> reference_length(1.0 * unit::microns);
        BaseUnits::Instance()->SetReferenceLengthScale(reference_length);

Set up a simulation domain, which will be a cuboid.

        units::quantity<unit::length> domain_width(100.0 * 1.e-6 * unit::microns);
        units::quantity<unit::length> domain_height(100.0 * 1.e-6 * unit::microns);
        units::quantity<unit::length> domain_depth(20.0 * 1.e-6 * unit::microns);
        boost::shared_ptr<Part<3> > p_domain = Part<3>::Create();
        p_domain->AddCuboid(domain_width, domain_height, domain_depth, DimensionalChastePoint<3>(0.0, 0.0, 0.0));

Make a regular grid on the domain

        boost::shared_ptr<RegularGrid<3> > p_grid = RegularGrid<3>::Create();
        p_grid->GenerateFromPart(p_domain, 10.0*reference_length);

Set up a PDE, we will model oxygen diffusion.

        boost::shared_ptr<LinearSteadyStateDiffusionReactionPde<3> > p_oxygen_pde = LinearSteadyStateDiffusionReactionPde<3>::Create();
        units::quantity<unit::diffusivity> oxygen_diffusivity(1.e-6*unit::metre_squared_per_second);
        p_oxygen_pde->SetIsotropicDiffusionConstant(oxygen_diffusivity);

Add continuum sink term for cells

        units::quantity<unit::rate> oxygen_consumption_rate(1.e-6*unit::per_second);
        p_oxygen_pde->SetContinuumLinearInUTerm(-oxygen_consumption_rate);

Add a Dirichlet boundary condition on the left face of the domain.

        p_domain->GetFacet(DimensionalChastePoint<3>(0.0,
                                                     domain_height/(2.0*reference_length),
                                                     domain_depth/(2.0*reference_length)))->SetLabel("boundary_1");
        boost::shared_ptr<DiscreteContinuumBoundaryCondition<3> > p_left_face_boundary = DiscreteContinuumBoundaryCondition<3>::Create();
        p_left_face_boundary->SetType(BoundaryConditionType::FACET);
        p_left_face_boundary->SetDomain(p_domain);
        p_left_face_boundary->SetValue(10.0*unit::mole_per_metre_cubed);
        p_left_face_boundary->SetLabelName("boundary_1");

Set up the PDE solvers for the oxygen problem

        boost::shared_ptr<FiniteDifferenceSolver<3> > p_oxygen_solver = FiniteDifferenceSolver<3>::Create();
        p_oxygen_solver->SetPde(p_oxygen_pde);
        p_oxygen_solver->SetGrid(p_grid);
        p_oxygen_solver->AddBoundaryCondition(p_left_face_boundary);
        p_oxygen_solver->SetLabel("oxygen");
        p_oxygen_solver->SetFileHandler(p_handler);
        p_oxygen_solver->SetFileName("fd_solution.vti");
        p_oxygen_solver->SetWriteSolution(true);
        p_oxygen_solver->Setup();
        p_oxygen_solver->Solve();
    }
};

Code

The full code is given below

File name `TestSolvingPdesLiteratePaper.hpp`

#include <cxxtest/TestSuite.h>
#include "AbstractCellBasedWithTimingsTestSuite.hpp"
#include "SmartPointers.hpp"
#include "OutputFileHandler.hpp"
#include "FileFinder.hpp"
#include "DimensionalChastePoint.hpp"
#include "UnitCollection.hpp"
#include "Owen11Parameters.hpp"
#include "GenericParameters.hpp"
#include "ParameterCollection.hpp"
#include "BaseUnits.hpp"
#include "MappableGridGenerator.hpp"
#include "Part.hpp"
#include "VesselNetwork.hpp"
#include "VesselNetworkGenerator.hpp"
#include "DiscreteContinuumMesh.hpp"
#include "VtkMeshWriter.hpp"
#include "FiniteElementSolver.hpp"
#include "FiniteDifferenceSolver.hpp"
#include "DiscreteSource.hpp"
#include "VesselBasedDiscreteSource.hpp"
#include "DiscreteContinuumBoundaryCondition.hpp"
#include "LinearSteadyStateDiffusionReactionPde.hpp"
#include "MichaelisMentenSteadyStateDiffusionReactionPde.hpp"
#include "PetscSetupAndFinalize.hpp"
class TestSolvingPdesLiteratePaper : public AbstractCellBasedWithTimingsTestSuite
{
public:
    void TestLinearReactionDiffusionPdeWithFiniteDifferences() throw(Exception)
    {
        MAKE_PTR_ARGS(OutputFileHandler, p_handler, ("TestSolvingPdesLiteratePaper/TestLinearReactionDiffusionPdeWithFiniteDifferences"));
        units::quantity<unit::length> reference_length(1.0 * unit::microns);
        BaseUnits::Instance()->SetReferenceLengthScale(reference_length);
        units::quantity<unit::length> domain_width(100.0 * 1.e-6 * unit::microns);
        units::quantity<unit::length> domain_height(100.0 * 1.e-6 * unit::microns);
        units::quantity<unit::length> domain_depth(20.0 * 1.e-6 * unit::microns);
        boost::shared_ptr<Part<3> > p_domain = Part<3>::Create();
        p_domain->AddCuboid(domain_width, domain_height, domain_depth, DimensionalChastePoint<3>(0.0, 0.0, 0.0));
        boost::shared_ptr<RegularGrid<3> > p_grid = RegularGrid<3>::Create();
        p_grid->GenerateFromPart(p_domain, 10.0*reference_length);
        boost::shared_ptr<LinearSteadyStateDiffusionReactionPde<3> > p_oxygen_pde = LinearSteadyStateDiffusionReactionPde<3>::Create();
        units::quantity<unit::diffusivity> oxygen_diffusivity(1.e-6*unit::metre_squared_per_second);
        p_oxygen_pde->SetIsotropicDiffusionConstant(oxygen_diffusivity);
        units::quantity<unit::rate> oxygen_consumption_rate(1.e-6*unit::per_second);
        p_oxygen_pde->SetContinuumLinearInUTerm(-oxygen_consumption_rate);
        p_domain->GetFacet(DimensionalChastePoint<3>(0.0,
                                                     domain_height/(2.0*reference_length),
                                                     domain_depth/(2.0*reference_length)))->SetLabel("boundary_1");
        boost::shared_ptr<DiscreteContinuumBoundaryCondition<3> > p_left_face_boundary = DiscreteContinuumBoundaryCondition<3>::Create();
        p_left_face_boundary->SetType(BoundaryConditionType::FACET);
        p_left_face_boundary->SetDomain(p_domain);
        p_left_face_boundary->SetValue(10.0*unit::mole_per_metre_cubed);
        p_left_face_boundary->SetLabelName("boundary_1");
        boost::shared_ptr<FiniteDifferenceSolver<3> > p_oxygen_solver = FiniteDifferenceSolver<3>::Create();
        p_oxygen_solver->SetPde(p_oxygen_pde);
        p_oxygen_solver->SetGrid(p_grid);
        p_oxygen_solver->AddBoundaryCondition(p_left_face_boundary);
        p_oxygen_solver->SetLabel("oxygen");
        p_oxygen_solver->SetFileHandler(p_handler);
        p_oxygen_solver->SetFileName("fd_solution.vti");
        p_oxygen_solver->SetWriteSolution(true);
        p_oxygen_solver->Setup();
        p_oxygen_solver->Solve();
    }
};

Solving PDEs in the Microvessel Project#

The Test#

Test 1 - Linear Reaction Diffusion With Finite Differences#

Code#

File name TestSolvingPdesLiteratePaper.hpp#

Solving PDEs in the Microvessel Project

The Test

Test 1 - Linear Reaction Diffusion With Finite Differences

Code

File name `TestSolvingPdesLiteratePaper.hpp`