CEPS
|
CEPS uses Finite Elements Methods (FEM) to discretize and solve the Partial Differential Equations from the cardiac model. We won't go much into details about FEMs in this page, so it is advised for newcomers to read introduction courses about FEM, in particular on how the variational formulation of a PDE is obtained.
The computational domain is given as input in the form of a geometrical mesh, in which each region of the domain is split into basic shapes. CEPS current implementation only works with meshes formed with geometrical simplices, i.e.
CEPS will take into account attributes that are given on mesh nodes, which are identifiers associated to points when building the mesh in order to constitute specific regions (e.g. stimulation zones, tissue degradation, ionic models selection, etc.)
Possible mesh file formats are .mesh (generated by e.g. gmsh),
.ele
.node .edge sets of files generated by TetGen, and
.vtk legacy files.
Finite Elements are constituted of four ingredients:
Geometrical data with the nodes defining the shape of the element (a line, a square, a tetrahedron, etc.) , and the jacobian of the transformation from a reference element to the element itself:
Boundary conditions can be given on different sets of elements for Neumann and Robin BC or nodes for Dirichlet BCs. Neumann and Robin BC are taken into account when computing the integrals on boundary elements by replacing the value of by the appropriate value. Dirichlet BCs are imposed by zeroing the rows of the discretization matrix corresponding to Dirichlet nodes (with 1 on the diagonal). It is also possible to zero the column corresponding to these nodes to preserve the matrix symmetry (if any), but this comes to the expense of extra matrix copy operations that may slow down the execution more than necessary.
Several classes of ODE solvers are implemented for ionic models but also for more generic functions. We consider in this section that we want to solve the ODE system
where can be a vectorial function (e.g. ionic variables on all PDE degrees of freedom). Some solvers are written based on the assumption that
can be split into a linear part
and a non-linear part
:
In the following represents the approximation of
and
the time step between
and
.
Note: ODE time steps are constant troughout the simulation. No adaptive time stepping methods are currently implemented.
Currently implemented methods
A bunch of explicit Runge-Kutta (RK) methods is coded in CEPS, from order 1 (which is Explicit Euler) to order 5. The evaluation of given
is based on the computation of
increments
at various times between
and
. Then, a linear combination of those increments is used to compute
. The algorithm is the following:
where the coefficients and
can be written using the Butcher formalism:
The method is explicit if is a lower triangular matrix with 0-diagonal.
The implementation follows this formalism, so it is rather easy to add RK methods. Here are the keywords for the currently implemented methods:
Input file Keyword | Order | Method |
---|---|---|
RK 1 | 1 | Explicit Euler |
RK Heun2 | 2 | Second order Heun method |
RK 2 | 2 | Usual Predictor-Corrector two steps method |
RK 3 | 3 | Kutta's third order method |
RK Heun3 | 3 | Heun's third order method |
RK Ralston3 | 3 | Ralston's third order method [1] |
RK SSP3 | 3 | Strong Stability Preserving Third order method [2] |
RK 4 | 4 | Classic RK 4th order |
RK 38 | 4 | Also a classic, the 3/8 method |
RK Ralston4 | 4 | Ralston's 4th order method, minimizes truncation error [1] |
RK Nystrom5 | 5 | 5th order, REF NEEDED [3] |
RK Fehlberg5 | 5 | Fehlberg's 5th order from the RK-45 Fehlberg scheme. [4] |
These simple multistep solvers are based on Taylor expansions on several steps. For example, for a method of order 3, given values and derivatives at times and
, we have
To achieve order 3, we find a linear combination of the three expansions that nullifies the terms in an
. This leads to
SBDF solvers are implemented for order 1 to 4.
This scheme simply solves
Note that there will trouble when .
Rush-Larsen (RL) solvers are exponential integrators using the decomposition of the evolution function in a linear and a non-linear part, and adapted to the ionic models. Please refer to [5] for explanations on the derivation of the scheme and a complete analysis.
FIXME: valid only for diagonal stabilizers ?
The general formula is
where
and the coefficients and
are computed using data from multiple previous steps. If we note
and
, coefficients are given depending on the order. RL schemes are implemented for order 1 to 4:
Order 1:
Order 2:
Order 3:
Order 4:
Another kind of exponential integrators... FIXME write here.
[6]
EAB schemes are implemented for order 1 to 4.
Details and analysis of the method can be found in [7]. Here we extract the algorithm, following the formalism of Runge-Kutta solvers introduced above.
The RK coefficients now depend on , and are calculated using the same
functions as in Rush-Larsen and EAB methods. Refer to section 5 of [7] to find the coefficients for order 2, 3 and 4.
ERK schemes are implemented for order 1 to 4.
[1] Ralston, Anthony. "Runge-Kutta Methods with Minimum Error Bounds". Math. Compu.
[2] Dale E. Durran, “Numerical Methods for Fluid Dynamics”, Springer. Second Edition.
[4] Ernst Hairer, Syvert Nørsett, and Gerhard Wanner. Solving Ordinary Differential Equations I: Nonstiff Problems, second edition, Springer-Verlag, Berlin. ISBN 3-540-56670-8.
[5] Yves Coudière, Charlie Douanla Lontsi, and Charles Pierre. Rush-Larsen time-stepping methods of high order for stiff problems in cardiac electrophysiology. arXiv
[6] Yves Coudière, Charlie Douanla Lontsi, and Charles Pierre. Exponential Adams–Bashforth integrators for stiff ODEs, application to cardiac electrophysiology. Math Comput Simulat
[7] Explicit Exponential Runge–Kutta Methods for Semilinear Parabolic Problems Marlis Hochbruck and Alexander Ostermann SIAM Journal on Numerical Analysis