Structure

The CARP examples framework encompasses a number of repositories, the purpose of which can be initially confusing. This section aims to explain the structure.

A typical installation might look like:

software
├── carputils
├── carp-examples
│   ├── benchmarks
│   └── devtests
├── carp-meshes
└── carp-reference
    ├── benchmarks
    └── devtests

carputils

A repository containing most of the Python functionality used in setting up and postprocessing CARP simulations. The bin subdirectory should be added to your PATH environment variable, and the top level directory should be added to your PYTHONPATH. On first run of an experiment, or when using the automatic clone script, the settings.yaml file will be generated, which must be tailored to your setup.

carp-examples

This is a normal directory which contains 2 (or more) Python packages whose modules define experiments. Each run.py script in these packages runs an experiment. In fact, as long as carputils is set up correctly, you can move these run scripts anywhere on your system (along with any input files they depend on) and they will work, but in order to run regression tests they must be inside these packages and carp-examples must be added to the PYTHONPATH.

carp-meshes

This is an SVN repository containing meshes used by examples. As it is very large, you may instead use the --get-mesh flag with run scripts to get only their needed meshes. Note that most run scripts generate meshes programmatically, in which case this is not needed. The location of this repository is set in carputils/settings.yaml with MESH_DIR.

carp-reference

These directories contain the reference solutions for regression tests. The repositories under carp-reference must have the same name as those under carp-examples. These repositories are only needed for running regression tests. The location of carp-reference is set by the REGRESSION_REF settings in carputils/settings.yaml

Concepts

The following concepts are used in the carputils framework:

sim

A single CARP execution - consecutive, related runs are considered to be multiple sims

job

One or more related sims launched by a single execution of a carputils run script

example

A carputils run script, capable of launching a single job per execution, though different jobs can be launched by changing the command line arguments

test

A fixed (according to the command line arguments) run configuration of an example plus several checks to determine the accuracy of the output

check

A method to determine the accuracy of the output of a test - most commonly by comparing with a previous reference solution but also may be by comparison with an analytic solution and other custom methods

Generate HTML Page

This documentation is based on Sphinx and can be generated using the Makefile in the carputils/doc folder.

cd carputils/doc
make clean
make html

The html webpages are found then in the folder carputils/doc/build/html. To view, point your favorite browser to

firefox ${CARPUTILS_DIR}/doc/build/html/index.html

Anatomy of an Example

A CARP run script will contain at a minimum:

from carputils import tools

@tools.carpexample()
def run(args, job):
    # Generate CARP command line and run example
    pass # Placeholder

if __name__ == '__main__':
    run()

The script defines a run function which takes as its arguments a namespace object containing the command line parameters (generated by the standard library argparse module) and a Job object. The run function should use the carputils.tools.carpexample() function decorator, as shown, which generates the above arguments and completes some important pre- and post-processing steps.

Many of the tasks of run scripts are common between examples, for example for setting up the correct CARP command line options for a particular solver. carputils therefore performs many of these common tasks automatically or by providing convenience functions. Below is a more complete example:

"""
Description of the example.
"""

EXAMPLE_DESCRIPTIVE_NAME = 'Example of Feature X'
EXAMPLE_AUTHOR = 'Joe Bloggs <joe.bloggs@example.com>'

from datetime import date
from carputils import tools

def parser():
    # Generate the standard command line parser
    parser = tools.standard_parser()
    group  = parser.add_argument_group('experiment specific options')
    # Add an extra argument
    group.add_argument('--experiment',
                        default='bidomain',
                        choices=['bidomain', 'monodomain'])
    return parser

def jobID(args):
    today = date.today()
    ID = '{}_{}_{}_np{}'.format(today.isoformat(), args.experiment,
                                args.flv, args.np)
    return ID

@tools.carpexample(parser, jobID)
def run(args, job):

    # Generate general command line
    cmd = tools.carp_cmd('example.par')

    # Set output directory
    cmd += ['-simID', job.ID]

    # Add some example-specific command line options
    cmd += ['-bidomain', int(args.experiment == 'bidomain'),
            '-meshname', '/path/to/mesh']

    # Run example
    job.carp(cmd)

if __name__ == '__main__':
    run()

def parser():
    # Generate the standard command line parser
    parser = tools.standard_parser()
    group  = parser.add_argument_group('experiment specific options')
    # Add an extra argument
    group.add_argument('--experiment',
                        default='bidomain',
                        choices=['bidomain', 'monodomain'])
    return parser

def jobID(args):
    today = date.today()
    ID = '{}_{}_{}_np{}'.format(today.isoformat(), args.experiment,
                                args.flv, args.np)
    return ID

@tools.carpexample(parser, jobID)
def run(args, job):

# Generate general command line
cmd = tools.carp_cmd('example.par')

# Set output directory
cmd += ['-simID', job.ID]

# Add some example-specific command line options
cmd += ['-bidomain', int(args.experiment == 'bidomain'),
        '-meshname', '/path/to/mesh']

# Run example
job.carp(cmd)

if __name__ == '__main__':
    run()

Documenting Examples

Examples should be well commented, but should also start with a descriptive python docstring to serve as a starting point for new users. This should include:

• A description of the feature(s) the example tests.
• Relevant information about the problem setup (mesh, boundary conditions, etc.)
• The meaning of the command line arguments defined by the example, where this is not obvious.

The docstring must be at the top of the file, before any module imports but after the script interpreter line, and is essentially a comment block started and ended by triple quotes """:

#!/usr/bin/env python

"""
An example testing the ...
"""

EXAMPLE_DESCRIPTIVE_NAME = 'Example of Feature X'
EXAMPLE_AUTHOR = 'Joe Bloggs <joe.bloggs@example.com>'

The docstring should be followed by two variables as above. EXAMPLE_DESCRIPTIVE_NAME should be a short descriptive name of the test to be used as a title in the documentation, e.g. ‘Activation Time Computation’ or ‘Mechanical Unloading’. EXAMPLE_AUTHOR should be the name and email address of the primary author(s) of the test in the format seen above. Multiple authors are specified like:

EXAMPLE_AUTHOR = ('Joe Bloggs <joe.bloggs@example.com>,'
                  'Jane Bloggs <jane.bloggs@example.com> and'
                  'Max Mustermann <max.mustermann@beispiel.at>')

Heading 1
=========

Heading 1.1
-----------

Heading 1.1.1
+++++++++++++

Heading 1.2
-----------

Here is an example of a run function::

    def run(argv, outdir=None):
        for list in argv:
            print argv

.. code-block:: bash

    echo "some message" > file.txt

.. math::

    (a + b)^2 = a^2 + 2ab + b^2

    (a - b)^2 = a^2 - 2ab + b^2

I really like the equation :math:`a^2 + b^2 = c^2`.

.. image:: /images/ring_rigid_bc.png

carputils Mesh Generation

carputils provides extensive functionality for the automatic generation of meshes according to geometric shapes. This currently includes:

• 1D cable (carputils.mesh.Cable)
• 2D grid (carputils.mesh.Grid)
• Block/cuboid (carputils.mesh.Block)
• Ring/cylindrical shell (carputils.mesh.Ring)
• Ellipsoidal shell (carputils.mesh.Ellipsoid)

The carputils.mesh module also provides some utility functions to be used with the above classes. See the API documentation for more information. This page provides a general guide to using the mesh generation code.

Basic Usage

The definition of the mesh geometry is generally done when calling the class constructor. For example, if generating a block mesh of size 2x2x2mm and resolution 0.25mm:

from carputils import mesh
geom = mesh.Block(size=(2,2,2), resolution=0.25)

Note that all dimensional inputs to the meshing classes are in millimetres. Please see individual classes’ documentation, as linked above, for details on the arguments they accept.

To then generate the actual CARP mesh files, there are two approaches depending on the class being used. The block class, which is actually a wrapper for the mesher command line executable, requires use of the mesher_opts method and then for mesher to be called with these arguments. This is done automatically by carputils.mesh.generate(), which also checks if the same mesh already exists and only generates a new one if not:

meshname = mesh.generate(geom) # Returns the meshname of the new mesh

The other classes, which generate the mesh in python directly, generate the mesh files with a call to the generate_carp method:

geom = mesh.Ring(5)
meshname = 'mesh/ring'
geom.generate_carp(meshname)

Note that these classes also support the direct output to VTK with the generate_vtk and generate_vtk_face methods.

Fibre Orientation

Block objects support the assignment of fibre orientations with the set_fibres() method, which takes angles in degrees as arguments:

from carputils import mesh
geom = mesh.Block(size=(2,2,2), resolution=0.25)
geom.set_fibres(0,  # endocardial fibre angle
                0,  # epicardial fibre angle
                90, # endocardial sheet angle
                90) # epicardial sheet angle

Note that the material coordinates in meshes generated by Block are as illustrated below. Specifically, the direction corresponds to the circumferential direction, corresponds to the apical-basal direction and corresponds to the epi- to endocardium direction.

The other 3D mesh classes take the fibre_rule argument to the constructor, which should be a function returning a helix angle in radians when passed a normalised transmural distance (where 0 corresponds to endocardium and 1 corresponds to epicardium):

def rule(t):
    """
    +1 radian endo to -1 radian epi
    """
    return 1 - 2*t

# Or as a one-liner
rule = lambda t: 1 - 2*t

# Then to create a ring with this fibre rule
geom = mesh.Ring(5, fibre_rule=rule)

# Or even more compact
geom = mesh.Ring(5, fibre_rule=lambda t: 1 - 2*t)

The meshing code provides a convenience method for the generation of linear fibre rule in degrees or radians. For example, to generate the common +60 to -60 degree rule:

p60m60 = mesh.linear_fibre_rule(60, -60)
geom = mesh.Ring(5, fibre_rule=p60m60)

The 1D cable mesh always has a constant fibre direction aligned with the cable, while the 2D grid mesh has a constant fibre direction assigned with the fibres argument.

Tagging Regions

Block meshes can be tagged using one of three classes:

• An axis-aligned box (carputils.mesh.BoxRegion)
• A sphere (carputils.mesh.SphereRegion)
• A cylinder (carputils.mesh.CylinderRegion)

See the class documentation for details on the exact arguments they take. They all, however, take the tag argument, which specifies the tag to be assigned to the elements inside the region:

geom = mesh.Block(size=(2,2,2))
reg = mesh.BoxRegion((-1,-1,-1), (1,1,0), tag=10)
geom.add_region(reg)

The other classes have tag regions assigned by passing a python function which returns True/False indicating whether a passed point is inside. For example, to assign a tag region to all points for positive x:

geom = mesh.Ring(5)
reg = lambda coord: coord[0] >= 0.0
geom.add_region(2, reg)

Note that the region tag must be passed separately in this case. The region classes above can also be passed, however note that the tag defined in the region object’s constructor is ignored in favour of that passed to add_region.

Boundary Conditions

An additional convenience function block_boundary_condition() is included to provide for the easy definition of boundary conditions in block meshes. It uses the size and resolution information from the geometry object to compute a box that surrounds only the nodes of interest, and returns a set of CARP command line options. For example, to define a stimulus on the upper-x side of a block:

geom = mesh.Block(size=(2,2,2))
bc = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', False)

The other mesh classes do not have an equivalent method, however the other 3D mesh classes offer functionality for the automatic generation of .surf, .neubc and .vtx files corresponding to relevant surfaces on the geometry. For example, with a ring geometry, you may wish to define Dirichlet boundary conditions on the top and bottom of the Ring and Neumann boundary conditions to the inner face. To generate files for these:

geom = mesh.Ring(5)
meshname = 'mesh/ring'
geom.generate_carp(meshname, faces=['inside', 'top', 'bottom'])

Which will generate files named like mesh/ring_inside.neubc.

geom.generate_carp_rigid_dbc('mesh/ring')

Implementation Details

Tessellation

An important consideration when generating a mesh in a subclass is how the resulting elements are to be tessellated into tetrahedra, which is done for most practical meshes used with CARP. Code to procedurally generate meshes in subclasses of Mesh3D generally use a combination of prisms and hexahedra to fill the geometry, however care must be taken to ensure that the edges of neighbouring elements match after the tessellation into tetrahedra is performed.

Hexahedra are deliberately tesselated in such a way that a grid of them which are all oriented in the same way will join up well when tessellated. There are two types of prisms available, offering the two possible ways for a prism to be tessellated.

When writing new mesh classes, the following nets provide a guide to the way in which these elements are tessellated. Printing these out and constructing the shapes aids considerably in figuring out how the elements mesh together.

etype = 'hexa'
face  = 'front'
surface = 'endo'
self._register_face(etype, face, surface)

Conductivity Regions

Conductivity regions, defined by the class carputils.model.ConductivityRegion, or its alias carputils.model.GRegion, take a list of material IDs, an optional name for the region, and the intra- and extracellular conductivity for the tissue. See the class documentation for the specific parameter names. See also the classmethods for easy definition of passive and isotropic conductivity regions.

Example:

from carputils import model
# Tags 1 and 2, default conductivities
cr1 = model.ConductivityRegion([1, 2], 'myo')
# Tag 3, increased intracellular longitudinal conductivity
cr2 = model.ConductivityRegion([3], 'fast', g_il=0.2)
# Tag 4, isotropic conductivity
cr3 = model.ConductivityRegion.isotropic([4], 'iso', cond=0.1)

Eikonal Regions

Eikonal regions are defined by the class carputils.model.EikonalRegion or its alias carputils.model.EkRegion. Unlike other classes, it accepts only a single ID, which is actually the index of an existing conductivity region. Like conductivity regions, it also has isotropic() and passive() class method alternate constructors.

Examples:

from carputils import model
# Conductivity region 0, default conduction velocity
ek1 = model.EikonalRegion(0, 'myo')
# Conductivity region 2, isotropic conduction velocity
ek2 = model.EikonalRegion.isotropic(2, 'iso', vel=0.3)

Electrical Stimuli

The class carputils.model.Stimulus defines a model stimulus. The user is expected to provide the relevant keyword arguments when using the class directly, for example to define a stimulus with ‘somefile’ as the vertex file:

from carputils import model
stim = model.Stimulus('pacing', vtk_file='somefile')

The keyword arguments (like vtk_file above) correspond with the fields of the CARP command line structure (-stimulus[0].vtk_file in this case).

Additional class methods exist for conveniently setting up a ‘forced foot’ stimulus, whereby the activation times from an Eikonal simulation are used to drive a reaction-Eikonal simulation of electrophysiology. See the class documentation for more information.

Ionic Models

Like the conductivity regions, ionic models also take IDs and name arguments, however each different model takes its own set of optional parameters. The different available models, their allowed arguments and their default parameters can be found in the documentation of the carputils.model.ionic module.

Example:

from carputils import model
imp = model.ionic.TT2IonicModel([1, 2], 'myo', GNa=14)

Active Tension Models

The active tension models are a special case in that they are not models in their own right, but are actually plugins to the ionic models. Active tension models, which only take the model parameters as command line arguments, are assigned to an existing ionic model using its set_stress method:

from carputils import model
# Create an ionic model
imp = model.ionic.TT2IonicModel([1])
# Create an active tension model
act = model.activetension.TanhStress(Tpeak=150)
# Assign the active tension model to the ionic model
imp.set_stress(act)

Available active tension models and their allowed parameters can be found in carputils.model.activetension.

Mechanics Material Models

Like ionic models, mechanics material models take IDs, name and model parameters as their arguments. Available models and their allowed parameters are detailed in carputils.model.mechanics.

Example:

from carputils import model
mat = model.mechanics.NeoHookeanMaterial([1, 2], c=50)

Getting Help

As in any real in vivo or ex vivo experiment there is a large number of parameters that have to be set and tracked. In in silico experiments tend to be even higher since not only the measurement setup must be configured, but also the organ itself. The general purpose design of carpentry renders feasible the configuration of a wide range of experimental conditions. Consequently, the number of supported parameters is large, depending on the specific flavor on the order of a few hundreds of parameters that can be used for experimental descriptions. With such a larger number of parameters it is therefore key to be able to retrieve easily more detailed explanations of meaning and intended use of all parameters.

A list of all available command line parameters is obtained by

carpentry +Help

carpentry [+Default | +F file | +Help topic | +Doc] [+Save file]
Default allocation: Global

Parameters:
         -iout_idx_file RFile
         -eout_idx_file RFile
         -output_level Int
         -num_emission_dbc Int
         -num_illum_dbc Int
         -CV_coupling Short
         -num_cavities Int
         -num_mechanic_nbc Int
         -num_mechanic_dbc Int
         -mesh_statistics Flag
         -save_permute Flag
         -save_mapping Flag
         -mapping_mode Short
         ...
         -rt_lib String
         '-external_imp[Int]' String
         -external_imp '{ num_external_imp x String }'
         -num_external_imp Int

and help for a specific feature is queried by

carpentry +Help stimulus[0].strength

stimulus[0].strength:
        stimulus strength

        type:   Float
        default:(Float)(0.)
        units:  uA/cm^2(2D current), uA/cm^3(3D current), or mV
        Depends on: {
                stimulus[PrMelem1]
        }

More verbose explanation of command line parameters is obtained by using the +Doc parameter which prints the same output as +Help, but provides in addition more details on each individual parameter:

carpentry +Doc

iout_idx_file:
        path to file containing intracellular indices to output

        type:   RFile
        default:(RFile)("")


eout_idx_file:
        path to file containing extracellular indices to output

        type:   RFile
        default:(RFile)("")


output_level:
        Level of terminal output; 0 is minimal output

        type:   Int
        default:(Int)(1)
        min:    (Int)(0)
        max:    (Int)(10)


num_emission_dbc:
        Number of optical Dirichlet boundary conditions

        type:   Int
        default:(Int)(0)
        min:    (Int)(0)
        max:    (Int)(2)
        Changes the allocation of: {
                emission_dbc
                emission_dbc[PrMelem1].bctype
                emission_dbc[PrMelem1].geometry
                emission_dbc[PrMelem1].ctr_def
                emission_dbc[PrMelem1].zd
                emission_dbc[PrMelem1].yd
                emission_dbc[PrMelem1].xd
                emission_dbc[PrMelem1].z0
                emission_dbc[PrMelem1].y0
                emission_dbc[PrMelem1].x0
                emission_dbc[PrMelem1].strength
                emission_dbc[PrMelem1].dump_nodes
                emission_dbc[PrMelem1].vtx_file
                emission_dbc[PrMelem1].name
        }
        Influences the value of: {
                emission_dbc
        }

Details on the executable and the main features it was compiled with are queried by

carpentry -revision

svn revision:             2974
svn path:                 https://carpentry.medunigraz.at/carp-dcse-pt/branches/mechanics/CARP
dependency svn revisions: PT_C=338,elasticity=336,eikonal=111

carpentry +Help experiment

carpentry -meshname slab -tend 100

carpentry +Help bidomain

Post-Processing Mode

Quantities appearing directly as variables in the equations being solved are directly solved for during a simulation. Other quantities which are indirectly derived from the primal variables can also be computed on the fly during a simulation, or, alternatively, subsequent to a simulation in a post-processing steps. On the fly computation is often warranted for derived quantities such as the instant of activation when high temporal resolution is important. For many derived parameters this may not be the case and it may be more convenient then to compute these quantities in a post-processing step from the output data. The computation of such auxiliary quantities such as, for instance, the electrocardiogram is available as postprocessing options. The post-processing mode is triggered by selecting

carpentry -experiment 4 -simID experiment_S1_2_3

In this case, the simulation run is not repeated, but the output of the experiment specified with the simID option is read in and the sought after quantities are computed. Output of data computed in post-processing mode is written to a subdirectory of the original simulation output directory. The data files in the output directory of the simulations are opened for reading, but other than that are left unchanged. A directory name can be passed in with ppID. If no directory name is specified the default name POST_PROC will be used.