Once the ceps executable is generated, you can run it from command line:

$> ./ceps myInputFile [-v] [-p]

A namespace for all utility methods.

CEPS uses a text file as primary input of the program. This file sets most of the simulation parameters. Input files commonly have the .in extension, but it can be any or none.

After a word about the syntax of this file, you will find information on parameters that are used for cardiac simulations. Several input files can be found in the example directory.

Syntax

Parameters are given by lines containing the colon : character:

a keyword : a value

Keywords are used by the application to identify parameters. If any mandatory keyword is missing, the program prints an error message about the missing key and stops. If the parameter has a built-in default value, a warning message is displayed when run with option -w.

Keywords are not case sensitive, and all white spaces are removed from them.

Values can be given on several lines if the backslash \ character is used:

another keyword : a value \
                  on \
                  several lines

Several entries can be given for the same keyword. In this case a single parameter is created by concatenating all the entries related to the same keyword, using commas , as separators of the new parameter string:

key : set of params 1, set of params 2
# is equivalent to
key : set of params 1
key : set of params 2

Run controls

Keyword	Value format	Default value	Description
disable output	YES/NO	NO	Set this entry to YES to disable all outputs. Some empty files be still be created.
output file name	file name prefix	./result	All output files will be prefixed with the given value. The prefix can contain a path to an output directory, which will be created if needed. See section Outputs.
output format	One of `MUSICARDIO`, `VTK_LEGACY`, `PARAVIEW_SERIES` (or `PSERIES`), `CEPS`	`PSERIES`	`MUSIC`, `VTK_LEGACY` : writes vtk legacy files with an XML collection file that can be opened in MUSICARDIO. `PSERIES` : pvtu files with `.series` collection file `CEPS` : custom format that avoids duplication of mesh data (cf Paraview Plugin)
output global indices	YES/NO	NO	adds a integer data array to the input with the global indices of points, as determined after partitioning.
output period	scalar real value	0.1	For PDE problems that include time, give the periodicity of spatial outputs (in simulation units).
output initial guess	YES/NO	NO	For static PDE problems, write the initial value given to the solver.
probe points	comma separated list of real triplets	–	Optional parameter. The values of unknowns at given points over time will be written in a separate output file
output analytic solution	YES/NO	NO	if the PDE has an analytic solution, or a reference solution (for convergence tests), then the solution is also written
output source terms	YES/NO	NO	writes all the regular source terms of the equation. If the problem has source terms of the form $\nabla\cdot(k\nabla f)$ , then f is also written.

Several run parameters are specific to cardiac problems, see Activation tracking section.

Geometry

Keyword	Value format	Default value	Description
3d mesh, 2d mesh, 1d mesh	Comma separated list of files	–	At least one `Xd mesh` entry must be given. Specify meshes that constitute the computational domain with these entries, respectively for volumic, surfacic and linear elements of the geometry. See this section for additional information on meshes. For Tetgen files, only the `.ele` file must be specified.
geometry scale	scalar real value	1.	Multiplies coordinates of all points of all meshes by this scalar.
coupled nodes	file name	–	Optional parameter. Give a list of coupled nodes to make links between two meshes. Please refer to this section to learn how to format the given file.

PDE definition

Keyword	Value format	Default value	Description
problem type	One of: - dirichlet anode cathode - flux anode cathode - monodomain - CL monodomain - bidomain - bilayer monodomain - pacemaker bidomain - pacemaker poisson	–	Mandatory entry.
pde start time	scalar real value	0.	For time based PDEs: starting time
pde end time	scalar real value	1.	For time based PDEs: end time

Solver options

Keyword	Value format	Default value	Description
spatial discretization	One of: - FE 1 - FE 2	FE 1	For now, only Lagrange finite elements of order 1 and 2 are available.
pde time step	scalar real value	0.1	For time based PDEs: time step in milliseconds. (Runge-Kutta methods span over several time steps.)
pde solver	One of: - FBE - SBDF 2, SBDF 3, SBDF 4 - CN	FBE	Select the time numerical scheme to solve the PDE. See the Numerics page for more details.
use mass lumping	YES/NO	NO	Lumps Finite Elements mass matrix (put sum of coefficients of each row in the diagonal)
linear solver type	One of CG, BICGSTAB, GMRES	GMRES	Conjugate gradient, BICGSTAB or GMRES solver
linear solver absolute tolerance	real number	PETSc defaults	Stopping criterion of linear solver
linear solver relatice tolerance	real number	PETSc defaults	Stopping criterion of linear solver
linear solver divergence tolerance	real number	PETSc defaults	Stopping criterion of linear solver
linear solver max iterations	integer	PETSc defaults	Stopping criterion of linear solver

Convergence tools parameters

Those parameters are required on top of other regular parameters to run a convergence study for a PDE problem. Time parameters are not required for static problems. See Automated convergence tests for additional info.

Keyword	Value format	Default value	Description
convergence 3d meshes, convergence 2d meshes, convergence 1d meshes	list of meshes (no commas)	–	Replaces `Nd mesh` entry. Error will be computed for each of these meshes. Only meshes in one file can be used for convergence tests
convergence time steps	scalars	–	time steps for which error will be computed, size can be different than list of meshes
reference 3d mesh, reference 2d mesh, reference 1d mesh	file name	–	If no analytic solution is available, solution computed on this mesh with `reference time step` will be used as reference for error computation.
reference time step	scalar	–	Time step used to compute reference solution if no analytic expression is available.
reference use existing	YES/NO	NO	If the problem requires a reference solution, you may not recompute it. Caution : reference solution will be searched where it is written by a run with reference computation (`outputDir/cvRef`).
convergence output file	file name	–	Name of file in which errors will be written.
convergence ignore 0D unknowns	YES/NO	NO	Activate this option to discard 0D unknowns in error computation.
convergence closest point tolerance	real value	1.e-8	Maximum search distance when looking for closest point in coarse meshes. This should be raised if the reference mesh has boundaries with high curvature.

In addition, a report on convergence error can be automatically generated using the following option. This report is a markdown file that appears on this documentation on each problem page.

Keyword	Value format	Default value	Description
convergence report file	.md file name	<>	optional If provided, the convergence report will be generated with the given file name.

Parameters of Anode Cathode problems (Dirichlet or Neumann)

Keyword	Value format	Default value	Description
conductivity	function string	CONSTANT 1	See Define new functions with text input
anode attributes	list of integers	–	Mandatory. Region of the boundary defining the anode.
cathode attributes	list of integers	–	Mandatory. Region of the boundary defining the cathode.
robin coefficient	comma separated list of couples <real,attribute>	0. -1	To use a Robin boundary condition instead of homogeneous Neumann on non-electrode boundaries, you can specify the alpha coefficient in $\partial_n u = \alpha u$ .

Parameters of cardiac problems

Tissue parameters

Keyword	Value format	Default value	Description
Am or membrane surface	STRING	CONSTANT 1000	Set the value of the membrane surface per unit volume, in $\si{\per\centi\metre}$ . See section Define new functions with text input to get information on how to format this entry.
Cm or membrane capacitance	STRING	–	Set the value of the membrane capacitance, in $\si{\micro\farad\per\centi\metre\squared}$ . See section Define new functions with text input to get information on how to format this entry.
sigma i or intracellular conductivity	STRING	CONSTANT 0.4 0.2 0.2	Set the value of the diagonalized intracellular conductivity, in $\si{\milli\siemens\per\centi\metre}$ . See section Define new functions with text input to get information on how to format this entry.
sigma e or extracellular conductivity	STRING	CONSTANT 0.4 0.2 0.2	Set the value of the diagonalized extracellular conductivity, in $\si{\milli\siemens\per\centi\metre}$ . See section Define new functions with text input to get information on how to format this entry.
fibers	Either: -<fileName> [arrayName [arrayName [arrayName]]] - Comma-separated list of 1 2 or 3 physical coefficient options strings, as described in this section	–	If fibers should be read from a file, it can be the same as the mesh file. See Fibers. Array names can be given for longitudinal, and transverse directions of fibers. Alternatively, fiber orientation for each direction can be given as a function of space
volume fraction	fileName arrayName	–	Volume fraction file and data array name. See Volume Fraction
volume fraction map intracellular	fileName	–	Conductivity interpolation map based on volume fraction. See Volume Fraction for file formatting
volume fraction map extracellular	fileName	–	Conductivity interpolation map based on volume fraction. See Volume Fraction for file formatting
output current	YES/NO	NO	adds a scalar field to the spatial output files, which contains the value of $\partial_t v_m$ . It accounts for the stimulation and the ionic current.

Ionic model parameters

Keyword	Value format	Default value	Description
ionic model solver	STRING	FBE	Check the ODE Solvers section for a list available ODE solvers
ionic model	Comma separated list of groups : <Id> <model> [attributes]	—	Mandatory parameter. See Ionic Models for more information. Id is the identifier of the unknown on which the model is added (problem dependent) model is one of TTP06_Endo, TTP06_MidMyo, TTP06_Epi, MS03, CRN98, BR77 attributes can be given to restrict the model to specific regions

Many options are specific to a given ionic model. They are prefixed with the model name, and have the same name as the variable in the CellML repository they were extracted from. They are regrouped on the Ionic Models page, as there are quite a lot !

Stimulation

Keyword	Value format	Default value	Description
stimulation	Comma-separated list of STRINGS	–	Optional, but cardiac simulations are useless without any stimulation. Each string can contain the elements defining a new function. Additionally, the unknown on which the stimulation is applied can be specified with `UNKNOWN <io>` in the string. If not specified, it will apply on unknown with ID 0. Attributes can also be added after the `ATTRIBUTES` keyword.
boundary stimulation	Comma-separated list of STRINGS	–	Optional. Boundary conditions that can be used for stimulations can be set for cardiac problems with the function syntax defined below. The selection for unknowns is done in the same way as for volumic stimulations. Finally, either the `NEUMANN` or `DIRICHLET` keyword must be found to select the type of BC. Attributes can also be added after the `ATTRIBUTES` keyword. Note that this method sets the $I_{\mathrm{app}}$ volumic term to 0 anywhere.

Activation tracking

Keyword	Value format	Default value	Description
post processing period	scalar real value	Same as PDE time step	Periodicity of activation tracking.
stop at complete activation	YES/NO	NO	Computation will end when the transmembrane voltage has exceeded the activation threshold on all points. See also `activation time data`
stop at complete repolarization	YES/NO	NO	Computation will end when the transmembrane voltage has returned to (almost) its base value on all points. See also `activation time data`
activation time data	Comma separated groups containing: <id> <scalar1> <scalar2> <scalar3>	–	Mandatory for activation tracking. Defines the quantities used to detect depolarization and repolarization. - id is the id of the unknown to track (problem dependant) - scalar1 is the voltage threshold to cross to consider a location depolarized - scalar 2 is the percentage of repolarization necessary to consider a location repolarized - scalar3 is a minimum voltage above which peak detection occurs Usually scalar1 is set just above the resting voltage, scalar3 is a higher value, and scalar2 can be anything between 0 and 1. 0.5 is already computed by default.

Parameters for the bilayer monodomain problem

Keyword	Value format	Default value	Description
coupling attribute	integer	–	Mandatory for coupling. Attribute of regions of mesh on which will be added to the equation.
coupling strength	real	1.	Coupling strength .
sigma i layer 1 or intracellular conductivity layer 1	STRING	CONSTANT 0.4 0.2 0.2	Conductivity on layer 1. See options for sigma i. Warning Conductivities (intra) must be equal on each layer where coupling occurs.
sigma e layer 1 or extracellular conductivity layer 1	STRING	CONSTANT 0.4 0.2 0.2	Conductivity on layer 1. See options for sigma e. Warning Conductivities (extra) must be equal on each layer where coupling occurs.
sigma i layer 2 or intracellular conductivity layer 2	STRING	CONSTANT 0.4 0.2 0.2	Conductivity on layer 2. See options for sigma i. Warning Conductivities (intra) must be equal on each layer where coupling occurs.
sigma e layer 2 or extracellular conductivity layer 2	STRING	CONSTANT 0.4 0.2 0.2	Conductivity on layer 2. See options for sigma e. Warning Conductivities (extra) must be equal on each layer where coupling occurs.
fibers layer 1	<fileName> [arrayName [arrayName [arrayName]]]	–	Fibers orientation for layer 1. The fibers keyword is ignored.
fibers layer 2	<fileName> [arrayName [arrayName [arrayName]]]	–	Fibers orientation for layer 2. The fibers keyword is ignored.

Parameters for the current-lifted monodomain problem

CL-monodomain problem

Keyword	Value format	Default value	Description
electrodes stimulation	time function options, see below	–	Options describing the evolution of input current .
anode attributes	list of integers	–	List of regions in mesh where from which current flows.
cathode attributes	list of integers	–	List of regions in mesh where into which current flows.
robin coefficient	comma separated list of couples <real,attribute>	0. -1	To use a Robin boundary condition instead of homogeneous Neumann on non-electrode boundaries, you can specify the alpha coefficient in $\partial_n u = \alpha u$ .

Parameters for the extended bidomain problem

extended bidomain problem

Keyword	Value format	Default value	Description
blood attributes, torso attributes, extracardiac attributes	list of integers	–	optional, domains in which the extracellular potential is extended into extracardiac potential. If none is specified, the extension is made in the whole computational domain (except tissue!).
anode attributes	list of integers	–	List of regions in mesh where from which current flows.
cathode attributes	list of integers	–	List of regions in mesh where into which current flows.
robin coefficient	comma separated list of couples <real,attribute>	0. -1	To use a Robin boundary condition instead of homogeneous Neumann on non-electrode boundaries, you can specify the alpha coefficient in $\partial_n u = \alpha u$ .
electrodes stimulation	time function options, see below	–	Options describing the evolution of input current .

Parameters specific to pacemaker problems

Keyword	Value format	Default value	Description
anode attributes	list of integers	–	List of regions in mesh where from which current flows.
cathode attributes	list of integers	–	List of regions in mesh where into which current flows.
robin coefficient	comma separated list of couples <real,attribute>	0. -1	To use a Robin boundary condition instead of homogeneous Neumann on non-electrode boundaries, you can specify the alpha coefficient in $\partial_n u = \alpha u$ .
pacemaker stimulation nspikes	integer	1	Number of consecutive stimulation cycles to simulate.
pacemaker stimulation amplitude	real value	1.0	In mV, initial charge of tank capacitor delivering current into the tissue during pulse phase of stimulation cycle.
pacemaker stimulation duration	real value	1.0	In ms, duration of pulse phase of stimulation cycle.
pacemaker pulse only mode	boolean	false	Simulates only the pulse part of the stimulation cycle.
compute electric field	boolean	false	Computes the electric potential distribution (from which elec field can be computed) generated by the pacemaker as if the whole domain was passive.
pacemaker offset duration	real value	0.5	in ms. Time offset before start of first stim cycle.
pacemaker period duration	real value	666.0	in ms (equiv 90 bpm) stimulation frequency.
pacemaker switch	real value	0.122	in ms. duration between pulse and COD phases of the cycle.
pacemaker ocd duration	real value	12.2	in ms. duration of OCD phase.
pacemaker internal CETS	real value	9.3639	In micro Farad. Tank capacitor (for pulse phase).
pacemaker internal CTS	real value	10.6248	In micro Farad. OCD capacitor (for contact discharge).
pacemaker internal Rgnd	real value	0.018	In kOhm. Device resistance.
pacemaker internal Rpulse	real value	0.007	In kOhm. Device resistance.
pacemaker internal Rwa	real value	0.080	In kOhm. Lead wire resistance (anode).
pacemaker internal Rwc	real value	0.20	In kOhm. Lead wire resistance (cathode).
pacemaker internal Rbig	real value	80.0	In kOhm. Device resistance (sensing and recharge phase).
pacemaker internal Rlittle	real value	0.005	In kOhm. Recharge circuit resistance.
anode resistance	real value	2.0	Two of these three anode parameters are required. R//C contact model resistance (kOhm).
anode capacitance	real value	18.74	Two of these three anode parameters are required. R//C contact model capacitance (micro Farad).
anode time constant	real value	37.48	Two of these three anode parameters are required. R//C contact model characteristic time (RC, ms).
cathode resistance	real value	0.03	Two of these three cathode parameters are required. R//C contact model resistance (kOhm).
cathode capacitance	real value	5.55	Two of these three cathode parameters are required. R//C contact model capacitance (micro Farad).
cathode time constant	real value	0.1665	Two of these three cathode parameters are required. R//C contact model characteristic time (RC, ms).
pacemaker offset time step	real value	0.01	In ms. Time step before any stim.
pacemaker pulse time step	real value	0.01	In ms. Time step during pulse phase.
pacemaker switch time step	real value	0.01	In ms. Time step between pulse and OCD.
pacemaker ocd time step	real value	0.01	In ms. Time step during OCD.
pacemaker waiting time step	real value	0.01	In ms. Time step between two cycles.

Define new functions with text input

Functions can be initialized through a single character string. The formatting presented in this section applies to physical coefficients defined on the computational domain, or functions for boundary conditions, source terms, etc. PDE problems use these function depending on their type. The type of values returned by can be either real, vectorial ( $\mathbb{R}^3$ ) or tensorial ( $\mathcal{M}_{3,3}(\mathbb{R})$ ).

The string starts with a keyword that identifies the type of spatial profile the coefficient will have. Options for each keyword will be given afterwards.

CONSTANT : a constant value on the whole domain,
CPIECEWISE : values are piecewise constant
FILE: values defined on points or cells of the mesh are imported from a file,
FUNCTION : a function of time and/or space with compact support
FPIECEWISE : values are defined piecewise, with a previously defined function on each piece
OPERATOR: creates a new functions from two already defined functions and an operator. Only multiplication is supported as of today.

Constant coefficients

After the CONSTANT keyword, one, three or nine scalars must be given depending on the type of return value of $\ff$\f.

Piecewise-constant coefficients

After the CPIECEWISE keyword, one should give a comma separated list of pairs (attribute, value(s)). For each pair, the given value (real, vectorial or tensorial) will be associated to the attribute, which will be used to evaluate the function. Be careful that some functions are evaluated on cells, other on points. Give the attributes and/or build the mesh accordingly.

Coefficients read from file

The following syntax is expected for the FILE keyword:

FILE <fileName> <SCALAR,VECTOR,TENSOR> <arrayName> <CELL,NODE>

Two cases:

1 - if the file is either a .vtk (with POINT_DATA or CELL_DATA) or .msh (with $NodeData or $ElementData) file.

SCALAR VECTOR and TENSOR select the return type of the function.
arrayName is the name of the data array to load from the file.
CELL or NODE selects where the data is defined in the file: on cells or on points.

2 - Any other file name: the file is considered to be a column formatted file, with the 1D variable on the first column, then scalar, vectorial or tensorial data on other columns. No additional required, other than the filename.

Compact support functions

The FUNCTION keyword is restricted to scalar valued functions with compact support, that are defined as a product . Two keywords, SPACE and TIME delimit options in the string that defines the functions and . One or two of those keywords must be present. This way you can defined functions of space only or of time only.

The compact support functions can have the following profiles, plotted as functions of the normalized distance to center of support:

$\begin{tikzpicture}[x=6cm,y=5cm] \draw[line width=0.3pt,blue!80,dotted] (0.8,1) -- (0.8,0); \draw[line width=0.3pt,blue!30] (0.8,0) -- (0.8,-0.12); \draw[line width=0.3pt,blue!30] (1,0) -- (1,-0.12); \draw[line width=0.3pt,<->,>=latex,blue!30] (0.8,-0.1) -- (1,-0.1) node [above,pos=0.5]{$0.2$}; \foreach \x in {0,0.5,1} { \draw (\x,0) node [below]{\x} -- +(0,5pt); } \foreach \y in {0,0.5,1} { \draw (0,\y) node [left]{\y} -- +(5pt,0); } \draw[thick] (0,0) rectangle (1.3,1.1); \node[below left] (xlabel) at (1.3,0) {$d/r$}; \draw[thick,domain=0:1,smooth,variable=\x,blue] plot ({\x},{1.+6*\x^5*(-21+70*\x-90*\x*\x+52.5*\x*\x*\x-35/3.*\x^4)}); \draw[line width=1.25,domain=0:1,smooth,variable=\x,red] plot ({\x},{(1-\x)^6*(1+6*\x+21*\x*\x+56*\x^3+126*\x^4+252*\x^5)}); \draw[line width=1pt,domain=0:0.99,smooth,variable=\x,samples=50,orange] plot ({\x},{exp(1-1/(1-\x*\x))}); \draw[line width=1pt,domain=0:0.8,smooth,variable=\x,dashed,blue] plot ({\x},{1}); \draw[line width=1pt,domain=0.8:1,smooth,variable=\x,dashed,blue] plot ({\x},{0.5*(1+cos(deg(((\x-0.8)/0.2)*3.14159)))}); \draw[line width=1.5pt,blue] (1,0) -- (1.3,0); \draw[line width=1.25pt,red] (1,0) -- (1.3,0); \draw[line width=1pt,orange] (1,0) -- (1.3,0); \draw[line width=1pt,blue,dashed] (1,0) -- (1.3,0); \draw[line width=0.75pt,green] (0,1) -- (1,1) -- (1,0) -- (1.3,0); \draw[line width=1.5pt,blue] (1.5,0.65) -- (1.6,0.65) node[right]{\texttt{C4}}; \draw[line width=1.25pt,red] (1.5,0.55) -- (1.6,0.55) node[right]{\texttt{C5}}; \draw[line width=1pt,orange] (1.5,0.45) -- (1.6,0.45) node[right]{\texttt{CINF}}; \draw[line width=0.75pt,green] (1.5,0.35) -- (1.6,0.35) node[right]{\texttt{CONSTANT}}; \draw[line width=1pt,blue,dashed] (1.5,0.25) -- (1.6,0.25) node[right]{\texttt{SRECTANGLE 0.2}}; \end{tikzpicture}$

The C4 and C5 functions are such that their integral on is equal to 1. The SRECTANGLE function has an extra parameter between 0 and 1 that determines the width of the smooth transition from 1 to 0.

The following options can be given after the SPACE keyword:

Keyword	Value format	Optional - Mandatory - Default value	Description
PROFILE	One of: - CONSTANT - C4 - C5 - CINF - SRECTANGLE + parameter	CONSTANT	See plot above
CENTER	Three real numbers	0. 0. 0.	Position of the center of the support
DIAMETER	Real number	1.	Size of the support, in length units
AMPLITUDE	Real number	1.	Scale factor, maximum value of the function (at center)

The time profile can be periodic. Then the support of the function is a union of smaller supports with fixed duration. The following options can be given after the TIME keyword:

Keyword	Value format	Optional - Mandatory - Default value	Description
START	Real number	0.	Lower bound of support in time
END	Real number	5000.	Upper bound of support in time
DURATION	Real number	1.	Size of smaller periodic supports
PERIOD	Real number	-1.	Periodicity of the time function. Set to negative value to cancel periodicity.
AMPLITUDE	Real number	1.	Scale factor, maximum value of the function (at center of each smaller support)
PROFILE	One of: - CONSTANT - C4 - C5 - CINF - SRECTANGLE + parameter	CONSTANT	See plot above