CEPS  24.01
Cardiac ElectroPhysiology Simulator
Automated convergence tests

In addition to PDE solvers for cardiac applications, CEPS provides a tool that can be used to measure the numerical convergence rate of the solvers, in both spatial and time dimensions.

Numerical error is computed either with an anamytic solution or with respect to a reference solution $\bar{u}$ that is computed beforehand, usually on a very fine mesh with a very small time-step. Then several numerical solutions $u_i$ are computed on coarser meshes and with larger time steps.

The errors are computed as a combination of the following norms:

  • $L_\infty([0,T])$, $L_1([0,T])$, or value at final time $T$,
  • $L_\infty(\Omega)$, $L_1(\Omega)$, $L_2(\Omega)$, where $\Omega$ is the whole computational domain. Both absolute and relative errors are output for each combination of norms.

If required, the difference between the approximation and the reference solution is interpolated in space and time on the coarse mesh. The interpolation is performed by linear interpolation in space (vtk routines), and P3 Lagrange interpolation in time. All the computed errors are written in a text file as final result of the study.

Files for the reference solution solution are put in a subdirectory named cv_ref under the main output directory. The currently computed numerical solution will be put in a subdirectory named cv_num. For cardiac problems, copies of the activation maps and time profiles at probe points are placed in the cv_num directory, suffixed by the value of time step and mesh diameter. This avoids the overwriting of these post-processed files.

The required parameters for the text input file are listed in this section.

Here are some examples of convergence results for small test problems

Laplacian problem, Dirichlet BC.

Corresponding code and parameters can be found in tests/pde/convergence/TestCvLaplacian.hpp.

CepsString meshBase = CepsString(CEPS_ABSOLUTE_PATH)
+ "/tests/data/meshes/convergence/square.";
p.set("convergence 2d meshes", meshBase + "0.1.mesh "
+ meshBase + "0.05.mesh "
+ meshBase + "0.025.mesh "
+ meshBase + "0.0125.mesh "
+ meshBase + "0.00625.mesh ");
p.set("linear solver type","CG");
p.set("linear solver absolute tolerance","1.e-20");
p.set("linear solver relative tolerance","1.e-20");
p.set("linear solver max iterations","10000");
p.set("convergence report file","./cvLaplacian.md");
ConvergenceStudy<SineProductLaplacianProblem> cs(&p,false);
cs.run();
CepsString
std::basic_string< CepsChar > CepsString
C++ format string.
Definition: CepsTypes.hpp:128
ConvergenceStudy
Generic class that regroups common elements of convergence studies.
Definition: ConvergenceStudy.hpp:56

Error: L2 in space (relative to analytic solution)

hmax Number of cells Error (P1) Slope (P1) Error (P2) Slope (P2)
1.41e-01 242 2.29e-02 - 1.10e-04 -
6.85e-02 1054 5.33e-03 2.01e+00 9.73e-06 3.35e+00
3.59e-02 4262 1.26e-03 2.23e+00 1.16e-06 3.29e+00
1.90e-02 16794 3.20e-04 2.17e+00 1.43e-07 3.31e+00
9.60e-03 67836 7.85e-05 2.05e+00 1.78e-08 3.04e+00

Convergence order (linear regression) : 2.13e+00 (P1), 3.26e+00 (P2)

Heat problem, non-zero Dirichlet BC, pure diffusion (no source term).

Corresponding code and parameters can be found in tests/pde/convergence/TestCvHeat.hpp.

p.set("convergence 2d meshes", meshBase + "0.1.mesh "
+ meshBase + "0.05.mesh "
+ meshBase + "0.025.mesh "
+ meshBase + "0.0125.mesh "
+ meshBase + "0.00625.mesh ");
p.set("convergence time steps","0.02 0.01 0.005 0.002 0.001");
p.set("linear solver type","CG");
p.set("linear solver absolute tolerance","1.e-20");
p.set("linear solver relative tolerance","1.e-20");
p.set("linear solver max iterations","10000");
p.set("spatial discretization","FE 1");
p.set("pde solver","SBDF 2");
p.set("convergence report file","./cvHeatNonZeroDir.md");
ConvergenceStudy<NonZeroDirichletHeatProblem> cs(&p,false);
cs.run();

Here, analysis of convergence rate is more complicated as the error in time may dominate the error of discretization in space and vice versa.

Error : L2 in space, L1 in time (relative to analytic solution)

hmax (Ncells) \ dt 0.02 0.01 0.005 0.002 0.001 Slope (2 dt points) Slope (3 dt points) Slope (4 dt points) Slope (5 dt points)
0.14 (242) 7.53e-02 3.41e-02 2.51e-02 2.31e-02 2.29e-02 1.1 0.79 0.49 0.36
0.068 (1054) 6.20e-02 1.84e-02 8.41e-03 5.79e-03 5.46e-03 1.7 1.4 1 0.78
0.036 (4262) 5.89e-02 1.48e-02 4.53e-03 1.77e-03 1.39e-03 2 1.9 1.5 1.3
0.019 (16794) 5.81e-02 1.40e-02 3.63e-03 8.39e-04 4.50e-04 2.1 2 1.8 1.6
0.0096 (67836) 5.80e-02 1.37e-02 3.40e-03 6.01e-04 2.08e-04 2.1 2 2 1.9

|||||||||||||| | Slope (2 hmax points) | 0.27 | 0.85 | 1.5 | 1.9 | 2 | | Slope (3 hmax points) | 0.18 | 0.61 | 1.3 | 1.9 | 2 | | Slope (4 hmax points) | 0.13 | 0.44 | 0.97 | 1.7 | 2 | | Slope (5 hmax points) | 0.089 | 0.32 | 0.73 | 1.4 | 1.8 |

Pseudo wave front propagation (i.e Heat pb with source term that mimics depol-repolarization front)

p.set("convergence 2d meshes", meshBase + "0.1.mesh "
+ meshBase + "0.05.mesh "
+ meshBase + "0.025.mesh "
+ meshBase + "0.0125.mesh "
+ meshBase + "0.00625.mesh ");
p.set("convergence time steps","1.e-1 5.e-2 2.e-2 1.e-2 5.e-3");
p.set("linear solver type","CG");
p.set("linear solver absolute tolerance","1.e-20");
p.set("linear solver relative tolerance","1.e-20");
p.set("linear solver max iterations","10000");
p.set("spatial discretization","FE 1");
p.set("pde solver","SBDF 2");
p.set("convergence report file","./cvWaveFront.md");
p.set("pde start time","0.");
p.set("pde end time","1");
p.set("wave front normal","1. 1. 0.");
p.set("wave front velocity","0.5");
p.set("wave front start depol","0.5 0.5 0");
p.set("wave front start repol","-0.5 -0.5 0");
p.set("wave front width depol","0.5");
p.set("wave front width repol","0.5");
ConvergenceStudy<PseudoWaveFrontProblem> cs(&p,false);
cs.run();

Error : L2 in space, L1 in time (relative to analytic solution)

hmax (Ncells) \ dt 0.1 0.05 0.02 0.01 0.005 Slope (2 dt points) Slope (3 dt points) Slope (4 dt points) Slope (5 dt points)
0.14 (242) 5.40e-04 1.27e-03 1.71e-03 1.78e-03 1.79e-03 1.2 0.69 0.49 0.36
0.068 (1054) 1.66e-03 1.33e-04 3.76e-04 4.42e-04 4.60e-04 3.6 0.82 0.34 0.17
0.036 (4262) 2.01e-03 4.37e-04 3.50e-05 9.64e-05 1.13e-04 2.2 2.5 1.5 0.99
0.019 (16794) 2.10e-03 5.24e-04 6.16e-05 8.54e-06 2.39e-05 2 2.2 2.4 1.7
0.0096 (67836) 2.12e-03 5.45e-04 8.30e-05 1.56e-05 2.07e-06 2 2 2.1 2.3

|||||||||||||| | Slope (2 hmax points) | 1.6 | 3.1 | 2.1 | 1.9 | 1.9 | | Slope (3 hmax points) | 0.97 | 0.83 | 2.8 | 2.1 | 2 | | Slope (4 hmax points) | 0.65 | 0.25 | 1.9 | 2.6 | 2.2 | | Slope (5 hmax points) | 0.45 | 0.065 | 1.2 | 2 | 2.5 |

The error on the "diagonal" is surprisingly low!