wiki:WorkshopSessionD2

Creating a project in Eclipse

We recommend you use the Eclipse development environment (at least for this workshop). Open eclipse (type eclipse in the command line), and if it prompts you for a workbench location, enter

/home/scratch/workbench/

(Then choose the 'Go to workbench' icon).

Then, outside of Eclipse, move the Chaste code into the workbench folder:

#This will uncompress the release from the CD (in the local drive) 
#and then move it to the correct location:
cd /home/scratch/workbench/
tar xvfz /media/Chaste\ workshop\ 2009/chaste_release_1.tgz 
mv release_1 Chaste

Then, go to Eclipse and do

Project -> untick 'Build Automatically' (if it is ticked)
File -> New -> Project

Double-click General and then double-click Project. Type in 'Chaste' as the project name, so that the location should be the location of the code.

Then do

Window -> Open Perspective -> Other -> C/C++

Finally, do

File -> New -> Convert to a C/C++ Make Project 

Double-click Finish.

Looking at documentation

Tutorial examples for this release 1 are available at: https://chaste.comlab.ox.ac.uk/chaste/tutorials/release_1 and release 1.1 at:. https://chaste.comlab.ox.ac.uk/chaste/tutorials/release_1.1

Documentation generated from the code by Doxygen is available at: https://chaste.comlab.ox.ac.uk/chaste/docs/release_1 or https://chaste.comlab.ox.ac.uk/chaste/docs/release_1.1

Writing a simple test

Create a new test file in heart/test, by double-clicking on heart in the 'C/C++ Projects' side pane, and right-clicking on test -> New -> File. Name your new file Test[Something].hpp.

Read the tutorial on writing tests. Write a simple test that, for example, calculates the sum from n=1 to n=N of 1/n2, for some large N, and uses TS_ASSERT_DELTA to verify that the answer is pi2/6, to within some tolerance. (Note 1: Recall that in C or C++, if n>1 is an int/unsigned 1/(n*n) will be zero - to get the expected answer, you must write 1.0/(n*n). Note 2: pi is M_PI).

Now create a make target for the new file, by going to the Make Targets side pane (on the right-hand side) and right-clicking on Chaste -> Add Make Target. Choose a 'Target Name', delete 'all' from 'Make Target' and in 'Build command' type

scons test_suite=heart/test/Test[Something].hpp

Double-click the new make target to run the newly-written test, and make sure it compiles and passes. Now open a web browser and open the following (local, automatically-generated) webpage

file:///[CHASTE_FULL_PATH]/testoutput/[MACHINE].default/index.html

This is just a file on the computer. For example, it might be something like,

[CHASTE_FULL_PATH] = /home/scratch/workbench/Chaste
[MACHINE] = msc31.ecs.ox.ac.uk

Bookmark this page in the browser. This page displays the test results, using green or red bars to indicate whether the test passed on not. You can view the output that was printed to the screen for a particular test, by clicking on the link for that test.

Solving ODEs and PDEs (optional)

Read the tutorials on solving ODEs and linear PDEs. How much time you spend in these will depend on whether you have already spent time installing Chaste, your familiarity with C++ and in particular abstract classes and inheritence, and what you aim to use Chaste for - it may be sensible to skip these tutorials or only look at them briefly. However, they will give an understanding of the 'core' functionality in Chaste and also introduce the components used in cardiac problems - monodomain and bidomains problems are composed of linear PDEs (with a nonlinear source term) coupled to a set of ODEs.

Run the tutorials, then use them as a basis for the following exercises (either use the test file you already created above, or create another test file)

Suggested exercises (optional)

  • Solve the ODE dy/dt = f(t,y) for your choice of f, writing the results to a file. Locate and open this file, edit it to remove the header, and load and visualise the results in matlab (or gnuplot).
  • Solve the simple ODE dy/dt = y*cos(t), y(0)=1, until y first goes below 0.5, using state variables. Verify that the state variable has just gone below 0.5. (Remember to have a suitably large end time!). Note: to get the stopping time at the end (no output object is given when state variables are used), you can do
    euler_solver.GetStoppingTime();
    
  • Solve the Poisson equation u_xx + u_yy = 1 on the unit square with boundary conditions u=0, and visualise the results with Matlab (as described in the tutorial). The following code snippets will be useful
    // For returning the identity as the diffusion tensor
    c_matrix<double,2,2> ComputeDiffusionTerm(const ChastePoint<2>& rX)
    {
        return identity_matrix<double>(2);
    }
    
    // for creating a mesh on [0,1]x[0,1]
    TetrahedralMesh<2,2> mesh;
    unsigned num_elem = 10;
    mesh.ConstructRectangularMesh(num_elem,num_elem);   // creates a mesh of size 10cm by 10cm with triangles with sides 1cm
    mesh.Scale(1.0/num_elem, 1.0/num_elem);             // scales the mesh back down to [0,1]x[0,1]
    

Solving Cardiac Problems

First, read the tutorial on solving bidomain problems, and run this tutorial by doing

scons build=GccOpt test_suite=heart/test/TestRunningBidomainSimulationsTutorial.hpp

For running in parallel (with two processors), do

scons build=Parallel test_suite=heart/test/TestRunningBidomainSimulationsTutorial.hpp

For running in optimised mode in parallel (with two processors), do

scons build=GccOpt_2 test_suite=heart/test/TestRunningBidomainSimulationsTutorial.hpp

(See the README_INSTALLATION.txt file for more scons options).

Suggested exercises

  • Run the tutorials and visualise the results using meshalyzer.
  • Copy the first test in the tutorial and convert it to the analogous monodomain problem - you only need to change 'Bi/bi' to 'Mono/mono' and how you deal with the results Vec (which now has the form [V0 .. Vn], not [V0 phi0 .. Vn phin]. Hint we suggest you just use the ReplicatableVector to deal with the solution (fine for these small problems), not DistributedVector). Compare with the bidomain simulation visually.
  • Using the ap_triggered code at the end of the second test in the tutorial (suitably edited), or by calling monodomain_problem.SetWriteInfo() before Solve, determine, approximately, by trial and error, and for this monodomain test, the threshold below which the stimulus magnitude is too small to create an action potential. Hint: in this problem there is no bath so the following line is not needed
    if (mesh.GetNode(i)->GetRegion()==HeartRegionCode::TISSUE)
    
  • Convert the first 2D bidomain test to a 1D problem, stimulated at one end (a larger magnitude of the stimulus will be required). Convert this into a heterogeneous problem where a different cell model is used in a region in centre of the mesh (or, if you prefer, and more realistically, convert the Luo-Rudy cell model (heart/src/odes/LuoRudyIModel1991OdeSystem.hpp and .cpp) so that the parameter plateau_potassium_current_g_Kp can be set to a different value. (To do this: make the parameter non-const and add a Set method, and remember to initialise it in the constructor). The result of this change is a (unexciting) small quantitative difference; you will need to visualise the results, with and without the change, side-by-side in meshalyzer to see it. Note: in the cell factory, you can do something like
    if ( (x>0.4) && (x<0.6) )
    {
        LuoRudyIModel1991OdeSystem* p_luo_rudy_system = new LuoRudyIModel1991OdeSystem(mpSolver, mpZeroStimulus);
        p_luo_rudy_system->SetParam(0.1); // much larger than default value of 0.183
        return new p_luo_system;
    }
    
  • If you get this far, maybe try to reproduce the s1s2 experiment from Day 1, by going through the ChasteParameters.xml for that simulation and using the appropriate calls on HeartConfig and carefully creating a cell factory which gives cells with the required properties in the different regions. You will have to make calls to SetScaleFactorGks (etc) in FaberRudy2000Version3. It may seem like a lot of effort to set up this simulation using the code, compared to running the binary, but recall that the code provides a massive amount of flexibility compared to the binary: for example, using the code you can uses geometries other than cuboids, spatially-varying values of parameters (based on distance from edge, etc), and so on..

Exercise Answers

Last modified 11 years ago Last modified on Aug 6, 2009, 3:51:45 PM