CEPS
Implementation of a problem: example

This page will explain briefly what is going on in CEPS when running a simulation. This will be illustrated with the modonomain problem, for which the source files are

and for the abstract classes:

Initialization of monodomain problem

Let's take a closer look at this line of code:

// in src/applications/monodomain.cpp
p.setupWithInputFile(inputFile,g);

There are actually a lot of things happening when this method is called.

CardiacProblem::setupWithInputFile() is a method that will extract all needed data in a input file by associating keywords and values to actual meaningful parameters for the problem p and the geometry g. After reading the parameters, all structures will be set accordingly.

The source code of this method looks like this:

// in src/cardiac/problem/CardiacProblem.cpp : readParameters
// First read generic parameters of PDE
AbstractPdeProblem::readParameters(inputFile);
// Then cardiac specific parameters
m_cardiacParameters->readFromInputFile(inputFile,this);

These two call will fill the instances of PdeParameters and CardiacParameters to read their respective contents. PdeParameters contains parameters that are common to all PDE problem (start time, end time, time step, output controls, etc). CardiacParameters contains parameters specific to ... cardiac applications ! (Incredible!)

Then in the call to initializeWithParameters :

// in src/cardiac/problem/CardiacProblem.cpp : initializeWithParameters
this->initializeTissue(); // finite elements are created in initDomain
// we call setPdeTimes before initializePdeSolver because
// in the constructor of MonodomainSolver calls getPdeTimeStep
this->setPdeTimes(m_pdeParameters->getPdeStartTime(),
m_pdeParameters->getPdeEndTime(),
m_pdeParameters->getPdeTimeStep());
this->initializePdeSolver();
this->initializeStimuli();
// default ode time step is pde time step. Be aware that this is dangerous!
this->setOdeTimeStep(m_pdeParameters->getPdeTimeStep());
this->initializeIonicModels();
this->initializeBoundaryConditions();

All main structures are set using the previously read parameters. Not that all of these methods are virtual, and can thus be redefined in derived classes. Let's go into details for the geometrical domain.

Initialization of simulation domain

// in pde/problem/AbstractPdeProblem.cpp : initializeDomain
for (uint_t d=1;d<=3;d++)
for (const auto& e : m_parameters->getMeshesOfDim(d))
g->setMeshOfDim(d,e);
[...]
uint32_t uint_t
we make sure ints are 32 bits
Definition: CepsTypes.hpp:101

The Geometry g, which was previously instanciated in the main function, holds all simplices and points of the domain. Then, each linear, surfacic and volumic mesh given in the input file are added to this domain.

// in pde/problem/AbstractPdeProblem.cpp : initializeDomain
[...]
g->setPartitioningMethod(m_parameters->getPartitioningMethod());
g->computePartitioning();
[...]

Here the partitioning method is retrieved from the read parameters, and partitioning occurs. After partitioning, the finite elements are created using the geometry:

// in pde/problem/AbstractPdeProblem.cpp : initializeDomain
[...]
this->m_fe = new FiniteElements(g,m_parameters->getFEType());
[...]
Holds all finite elements corresponding to each geometrical element.
Definition: FiniteElements.hpp:54

The FiniteElements structure holds any additional node to those of the geometry (for example, when P2 elements are used). Each finite element is then just a collection of pointers to a geometrical element, additional nodes, base functions and quadrature formula.

In a nutshell

We just showed how the domain is initialized. Cardiac tissue, ionic models, numerical solver and more are initialized in the same way.

Let's wrap this up.

  • An input file specifies simulation parameters, using key-value pairs
  • In the main file, the user instanciates the monodomain problem class
  • The problem then:
    • instanciates several parameter classes to parse and store these parameters
    • initializes its data structures with those stored parameters

This closes the part on the initialization of the monodomain problem. The next section will focus on running the simulation.

Running the simulation

Let's go back to the main program. After the problem is set up, the virtual method CardiacProblem::run() is called

// in src/applications/Monodomain.cpp
p.run();

This method is responsible for:

  • setting the initial condition
  • asking the numerical solver to solve the equation
// in src/cardiac/problem/CardiacProblem.cpp : run
--> DHVecPtr solution = this->getInitialCondition();
this->m_solver->solve(solution);
std::shared_ptr< DistributedHaloVector > DHVecPtr
Typedef for pointer on Distributed Halo Vector.
Definition: DistributedHaloVector.hpp:152

The first part is simple. In our case the initial condition is, for each degree of freedom, the resting membrane potential. DHVecPtr is an alias for std::unique_ptr<DistributedHaloVector>. It is a C++11 pointer on an array of data that is distributed on all the CPUs, including data on the overlapping domains between local domains of each CPU.

The second part is much more interesting, as plenty of stuff happens.

// in src/cardiac/problem/CardiacProblem.cpp : run
DHVecPtr solution = this->getInitialCondition();
--> this->m_solver->solve(solution);

m_solver is an instance of a CardiacSolver class. Its role is described in the next section.

Numerical solver

In CEPS, a solver is a class where implementation of the PDE time discretization scheme takes places. Thus, it regroups:

  • a linear system of type $Ax=b$
  • assembly tools to build this linear system following the chosen time discretization
  • time iterations

 Global view of what is happening

As seen before, m_solver is an instance of a CEPS PDE solver. In this particular case, we are using a MonodomainSolver, which is derived from CardiacSolver. There is very little that is different from a cardiac PDE problem to another. So most of the solving methods are to be found in the CardiacSolver class.

Let us consider that we have chosen the FBE solver type in the input file: the forward-backward (semi-implicit) Euler time scheme. The linear system to solve is written:

\[ \left(M + \dfrac{\delta t}{A_\mathrm{m}C_\mathrm{m}} K[\sigma]\right) U^{n+1} = M U^n + M \dfrac{\delta t}{C_\mathrm{m}} \left(I_\mathrm{app} + I_\mathrm{ion}\right). \]

There are three important virtual methods in each solver class. We are going to take a look at each of them. The first one is the CardiacSolver::solve() method. It is the method called by derived cardiac classes to solve the equation. Its organization is a bit complicated, because it can handle many different solving methods, in particular numerical schemes involving data at several time steps.

The two other essential methods are setupLinearSystem() and updateAndSolveLinearSystem(). The first method is intended to be called as few times as possible, as it assembles the left hand side of the linear system, which is generally a discretization matrixi. In the case of FBE, it is called once. It is called several times when using multi-step methods.

The updateAndSolveLinearSystem() regroups the computation of the right-hand side of the linear system, and the resolution of the linear system. Unfortunately, the two operations have to be grouped in a single method in order to ensure compatibility with different methods (I am looking at you, Runge-Kutta!).

Let's go into details.

Setting up the linear system

The purpose of the setupLinearSystem() method is to setup all data structures that need to be initialized only once for the duration of the entire simulation. In our case, the left-hand side of the linear system is not time dependent: indeed, recall that the left-hand side matrix is of the form $ (M + \delta t K[\sigma]) $, where $ \delta t $ does not change throughout the simulation.

// in src/cardiac/solver/CardiacSolver.cpp
{
[...]
// create the mass matrix for this problem
massAssembler.assemble(m_mass,nullptr);
[...]
// compute the left hand side matrix: M + alpha*dt*K[sigma]
real_t adt = alpha*m_timeStepper->getTimeStep()/(m_params->getCm()*m_params->getAm());
MonodomainAssembler monoAssembler(m_problem,adt);
monoAssembler.assemble(m_lhs,nullptr);
float real_t
floating point numbers CEPS type, configured at Cmake
Definition: CepsTypes.hpp:134
uint_t getProblemDim() const
Get PDE dimension (number of unknowns per DOF)
Definition: AbstractPdeProblem.cpp:212
DMatPtr m_mass
Mass matrix M.
Definition: CardiacSolver.hpp:157
TimeStepper * m_timeStepper
Link to time management (...)
Definition: CardiacSolver.hpp:150
CardiacProblem * m_problem
Link to math problem.
Definition: CardiacSolver.hpp:148
virtual void setupLinearSystem()
Builds the linear system discretizing the PDE.
Definition: CardiacSolver.cpp:262
DMatPtr m_lhs
Left hand side matrix.
Definition: CardiacSolver.hpp:158
Assembles the mass matrix for a given k-simplexes geometry.
Definition: FEPkMassMatrixAssembler.hpp:38
Matrix building routines for the monodomain equation.
Definition: MonodomainAssembler.hpp:40
real_t getTimeStep() const
PDE Time step.
Definition: TimeStepper.hpp:80

These two code snippets introduce the concept of assemblers. The sole purpose of assemblers is to build distributed matrices and/or vectors, corresponding to the discretization of specific terms of an equation. For example, the FEPkMassMatrixAssembler being used in the code assembles mass matrices for the given finite elements. $ M_{ij} = \int_{\Omega}\varphi_i\varphi_j $ where $ \varphi_i $ is the $ i $-th basis function.

The MonodomainAssembler is used to compute the whole left-hand side matrix $ M + \dfrac{\delta t}{A_\mathrm{m}C_\mathrm{m}} K[\sigma]$.

PDE time iterations

Let's go back to the CardiacSolver::solve() method:

// in src/cardiac/solver/CardiacSolver.cpp
{
[...]
--> while (not p_timeStepper->atEnd())
{
[...]
}
[...]
}
virtual void solve(DHVecPtr solution)
Performs all iterations of the PDE solver.
Definition: CardiacSolver.cpp:151

This while loop uses a CEPS data structure that manages all PDE time parameters, the TimeStepper class. Solving and updating the linear system take place in this loop, as well as output operations, and all pointer shifting needed by multi-step methods.

The updateAndSolveLinearSystem() is a virtual method that is supposed to be called at each time step. This is why we find this method inside the time iterations loop. Its purpose is to assemble or update any structure that changes throughout the simulation.

In our case, each time a new value of $ U $ is computed, the right-hand side must be updated. Recall that: $ rhs = M U^n + M \frac{\delta t}{C_\mathrm{m}} (I_\mathrm{app} + I_\mathrm{ion}) $

Thus, 3 main operations occur:

  • compute ionic currents $ I_\mathrm{ion} $
  • compute any stimuli $ I_\mathrm{app} $
  • compute all vector/matrix, vector/vector operations into rhs

This is exactly what is being done in MonodomainSolver::updateAndSolveLinearSystem(), or, to be exact, in MonodomainSolver::updateAndSolveSBDF() in our case.

// in src/cardiac/solver/CardiacSolver.cpp
{
[...]
auto myocardium = m_problem->getMyocardium();
m_iIonNm[0]->zero();
myocardium->takeOneStep(m_VnM[0],m_problem->getOdeTimeStep());
--> myocardium->computeIonicCurrent(m_iIonNm[0]);
[...]
Myocardium * getMyocardium() const
Pointer on geometry and phyisiological data.
Definition: CardiacProblem.cpp:189
real_t getOdeTimeStep() const
Time stepping of ODE solvers (constant time steps)
Definition: CardiacProblem.cpp:169
std::vector< DVecPtr > m_iIonNm
Ionic current at previous steps.
Definition: CardiacSolver.hpp:161
vector< DHVecPtr > m_VnM
Vector of unknowns at previous steps.
Definition: CardiacSolver.hpp:155
void updateAndSolveSBDF(uint_t nSteps)
Updates matrix and RHS for the SBDF numerical scheme.
Definition: CardiacSolver.cpp:443

New values of ionic currents are computed then stored in m_iIonNm[0] vector (read it " \f$I_\mathrm{ion}^{n-0}\f$ ").

// in src/cardiac/solver/CardiacSolver.cpp
// Check for any applied stimuli, update Iapp if that's the case
m_iAppNm[0].zero();
--> p_p->getStimuliHandler()->applyStimuli(m_iAppNm[0], t);

Similarly, applied stimuli are computed and stored in m_iAppNm[0] vector.

Then, complete right-hand side is computed using operators on matrix and vector classes.

 Output

The MonodomainSolver class embeds an other class that is dedicated to the writing and the process of the solution vectors of cardiac problems. It is called PostProcessOutput, and is instanciated in the constructor of the MonodomainSolver class.

: CardiacSolver(problem), m_problem(problem)
{
[...]
// Initialize outputs
m_writer = new PostProcessOutput(m_params,m_problem,{"Transmembrane potential"});
}
Solves cardiac problems, that all share the same structure.
Definition: CardiacSolver.hpp:54
Defines and solve the monodomain problem, see models page of user doc.
Definition: MonodomainProblem.hpp:35
MonodomainSolver(MonodomainProblem *p)
Constructor.
Definition: MonodomainSolver.hpp:42
Class managing computations from potential to outputs.
Definition: PostProcessOutput.hpp:59

This instance gets its configuration through the CardiacParameters instance we have seen at the start of this page. Here we also define the name of the variable that will appear in the output files.