This tutorial was generated from the file projects/CardiacEpVerification/test/TestEpAgainstExactSolutionsLiteratePaper.hpp at revision r20536. Note that the code is given in full at the bottom of the page.

## Testing the cardiac solvers against exact solutions

This is the main code for solving the monodomain, bidomain and bidomain-with-bath model problems.

The following are all standard includes:

Next, the includes for classes that are defined in this project.
`ModelProblemCellModel`

defines the cell model for the model problems:

This file contains several classes which represent the exact solutions of the model problems:

A class for calculating the L2 error (squared) at given time using a computed solution and exact solution:

A class for calculating the second part of the H1 error (squared) at given time using a computed solution and exact solution:

This next file defines three classes which inherit from `MonodomainProblem`

, `BidomainProblem`

and `BidomainWithBathProblem`

but also do error calculations (using the above).

We have to define cell factories to create the cell models as always. This cell factory
creates the custom `ModelProblemCellModel`

for each node, gives each a zero stimulus,
passes in some parameters, and sets up the initial conditions used in the monodomain
and bidomain model problems.

This second cell factory is the same as the above except it passes in the initial conditions defined in the bidomain-with-bath model problem.

The code for running the model problems are defined in ’tests’, as with all code in Chaste (see main documentation).

following function can be mostly ignored as an alternative is used:

A function for computing errors, taking in the output directory of a cardiac problem and a class
saying how to calculate the exact solution. Note: since the output directory must be written for
this function to be used, and with small dt this will lead to a huge datafile (if printing timestep
is dt), we no longer use this method for calculating errors. However the code has been kept as it is
convenient. The errors are now calculated by computing the contributions to the error at the end of
each timestep, as the simulation is progressing. See `MonodomainProblemWithErrorCalculator`

and related classes.

### Monodomain model problem

This function is the code for running the monodomain model problem, which we will walk through.

Define h and dt based on the `parametersScaleFactor`

. Note dt is proportional to h^2^ as required.

Define conductivities in each direction as specified in the paper.

Define the integers m,,1,,, m,,2,,, m,,3,, that go into the function F, where F(x,y,z) = cos(m1*pi*x)*cos(m2*pi*y) cos(m3pi*z).

Set up the mesh to be the unit line/square/cube. Also compute the constant c.

End time is T = 1.

Define an output directory, but see comments below.

Set the timesteps, define the conductivity tensor, the capacitance and the surface-to-volume ratio.

The following cell factory creates cell models for the model problem cell model defined in the paper, and with the required initial values of (V,u,,1,,,u,,2,,,u,,3,,).

This class can be used to get the exact solution for the voltage: (1+t)^1/2^F(x).

Define the monodomain problem class. `MonodomainProblemWithErrorCalculator`

(defined in this project) is just
`MonodomainProblem`

but also does error computing at the end of each timestep.

Set the mesh, **don’t** write output unless in testing mode (since printing_dt = pde_dt, a lot of output would be
written to file (printing dt is small so that error computations can be carried out each dt), initialise
and solve.

Print the errors to screen (the commas and semi-colon are for easy copy & paste into Matlab).

There is also some test code. To check that the code (in particular the error calculators) is doing what it should be, if DIM=1 and the coarsest mesh is being used, we have the following tests.

Check nothing has changed:

Check the two ways of computing the error give the same results:

Finally, check the second way of computing the error: when `ComputeErrors()`

is called with `true`

as the last parameter it
ignores the numerical solution, ie just calculates the norm of the exact solution, which we can also
calculate analytically.

### Bidomain model problem

Next, the main function for solving the bidomain model problem. This is basically the same as the monodomain code, except has an extracellular conductivity, and gets the errors for both voltage and extracellular potential.

Code same as monodomain version:

Define the integers m,,1,,, m,,2,,, m,,3,, that go into the function F, where F(x,y,z) = cos(m1*pi*x)*cos(m2*pi*y) cos(m3pi*z).

Define the constant k:

Code similar to monodomain version:

Classes for returning V=(1+t)^1/2^F(x) and phi_e = -k(1+t)^1/2^F(x) and their derivatives:

Solve and print errors.

Test code similar to monodomain version:

### Bidomain-with-bath model problem

Finally, the bidomain-with-bath-model problem:

All this is as before:

Define the integer m,,1,, that goes into the function F(x,y,z) = cos(m1*pi*x). Note no m2 and m3 as for bath
problem these must be zero.

Set up the domain: x in `[-1,2]`

; y,z in `[0,1]`

. Note the translation at the end.

Set appropriate elements as bath:

As before:

Set up the bath conductivity and the electrodes, I,,E,,=-alpha on x=-1, and I,,E,,=alpha on x=2:

As before:

Use the bath version of the cell factory (different initial conditions):

Solve and output errors:

Similar to before:

### Main test

Finally, we have the public ’tests’, which actually run the simulations. Note that only the first two tests will be run,
as the code is written below, as only those whose name begins with ‘Test’ are run. To run the others, change
`doneTest..`

to `Test..`

.

As integrals over the domain have to be
computed each timestep, virtually all the time is spent in the error calculations, certainly for the larger N (ie smaller h). Also,
dt is proportional to h^2^, so for 3D, the work increases by at least a factor of 32 (2^3^ 2^2^) as N doubles. The N=3 simulations for
monodomain and bidomain take a long time on a desktop. Monodomain will run in parallel, bidomain may sometimes not (see comment in
`CardiacProblemWithErrorCalculatorClasses`

).

## Code

The full code is given below