CEPS

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
src/applications/monodomain.cpp
src/cardiac/problem/MonodomainProblem.cpp
src/cardiac/solver/MonodomainSolver.cpp
and for the abstract classes:
src/cardiac/problem/CardiacProblem.cpp
src/cardiac/solver/CardiacSolver.cpp
pde/problem/AbstractPdeProblem.cpp
Let's take a closer look at this line of code:
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:
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
:
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.
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.
Here the partitioning method is retrieved from the read parameters, and partitioning occurs. After partitioning, the finite elements are created using the geometry:
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.
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.
This closes the part on the initialization of the monodomain problem. The next section will focus on running the simulation.
Let's go back to the main program. After the problem is set up, the virtual method CardiacProblem::run()
is called
This method is responsible for:
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.
m_solver
is an instance of a CardiacSolver class. Its role is described in the next section.
In CEPS, a solver is a class where implementation of the PDE time discretization scheme takes places. Thus, it regroups:
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 forwardbackward (semiimplicit) Euler time scheme. The linear system to solve is written:
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 multistep methods.
The updateAndSolveLinearSystem()
regroups the computation of the righthand 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, RungeKutta!).
Let's go into details.
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 lefthand side of the linear system is not time dependent: indeed, recall that the lefthand side matrix is of the form , where does not change throughout the simulation.
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. where is the th basis function.
The MonodomainAssembler
is used to compute the whole lefthand side matrix .
Let's go back to the CardiacSolver::solve()
method:
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 multistep 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 is computed, the righthand side must be updated. Recall that:
Thus, 3 main operations occur:
This is exactly what is being done in MonodomainSolver::updateAndSolveLinearSystem()
, or, to be exact, in MonodomainSolver::updateAndSolveSBDF()
in our case.
New values of ionic currents are computed then stored in m_iIonNm
[0] vector (read it " \f$I_\mathrm{ion}^{n0}\f$ ").
Similarly, applied stimuli are computed and stored in m_iAppNm
[0] vector.
Then, complete righthand side is computed using operators on matrix and vector classes.
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.
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.