Changes between Version 2 and Version 3 of PaperTutorials/Frontiers2014/MonodomainConvergence


Ignore:
Timestamp:
Dec 15, 2014, 3:48:50 PM (6 years ago)
Author:
tutorial
Comment:

Auto tutorial submission based on revision r23235

Legend:

Unmodified
Added
Removed
Modified
  • PaperTutorials/Frontiers2014/MonodomainConvergence

    v2 v3  
    1 This tutorial is automatically generated from the file projects/Frontiers2014/test/TestMonodomainConvergenceLiteratePaper.hpp at revision r22835.
     1This tutorial is automatically generated from the file projects/Frontiers2014/test/TestMonodomainConvergenceLiteratePaper.hpp at revision r23235.
    22Note that the code is given in full at the bottom of the page.
    33
     
    88On this wiki page we describe in detail the code that is used to run this example from the paper.
    99
    10 We run a 1D monodomain simulation on a highly refined PDE time step with a handful of the models.
     10This test can be used to run a 1D monodomain simulation with a handful of the models,
     11taking the space and PDE time step to use from command line arguments.  It stimulates the fibre at
     12the left hand end, and saves the voltage trace at the right hand end to file for viewing/comparing
     13at different step settings.  It also checks that the excitation wave has propagated along the fibre,
     14and calculates some summary properties of the action potential, again writing these to file.
    1115
    1216== Code overview ==
     
    1822#include <cxxtest/TestSuite.h>
    1923
     24#include <vector>
     25#include <string>
     26#include <sstream>
    2027#include <boost/assign/list_of.hpp>
    2128#include <boost/foreach.hpp>
     
    3037#include "AbstractCvodeCell.hpp"
    3138
    32 // This header is needed to allow us to run in parallel
     39// Other Chaste headers
     40#include "FileFinder.hpp"
     41#include "OutputFileHandler.hpp"
     42#include "CommandLineArguments.hpp"
     43
     44// This header is needed to allow us to prevent being run in parallel
    3345#include "PetscSetupAndFinalize.hpp"
    3446
     
    4153        EXIT_IF_PARALLEL;
    4254
     55}}}
     56This command line argument says whether to reset the CVODE solver fully at each PDE time step.
     57{{{
     58#!cpp
    4359        bool reset_cvode = false;
    4460        if (CommandLineArguments::Instance()->OptionExists("--reset"))
     
    4864
    4965}}}
    50 This test was run with the following values inserted here:
    51 -> 0.001 ms (much finer than normal)
    52 -> 0.01 ms (typically a fine timestep)
    53 -> 0.1 ms (about average)
    54 -> 1 ms (coarse)
     66This test was run with the following values for PDE time step inserted here:
     67- 0.001 ms (much finer than normal)
     68- 0.01 ms (typically a fine timestep)
     69- 0.1 ms (about average)
     70- 1 ms (coarse)
    5571
    5672{{{
     
    6682        }
    6783
    68         /**
    69          * This test was run with the following values inserted here:
    70          * -> 0.01 cm (fine)
    71          * -> 0.001 cm (much finer than normal)
    72          */
     84}}}
     85This test was run with the following values for mesh spacing inserted here:
     86- 0.01 cm (fine)
     87- 0.001 cm (much finer than normal)
     88
     89{{{
     90#!cpp
    7391        double h;
    7492        if (CommandLineArguments::Instance()->OptionExists("--spacestep"))
     
    8199        }
    82100
    83         // A list of models that we want to do tissue simulations with.
     101}}}
     102Next we loop over the following list of models that we want to do tissue simulations with.
     103{{{
     104#!cpp
    84105        std::vector<std::string> models_to_use = boost::assign::list_of("luo_rudy_1991")
    85                                                  ("beeler_reuter_model_1977")
    86                                                  ("nygren_atrial_model_1998")
    87                                                  ("ten_tusscher_model_2004_epi")
    88                                                  ("grandi_pasqualini_bers_2010_ss")
    89                                                  ("shannon_wang_puglisi_weber_bers_2004")
    90                                                  ("iyer_model_2007");
     106                                                                       ("beeler_reuter_model_1977")
     107                                                                       ("nygren_atrial_model_1998")
     108                                                                       ("ten_tusscher_model_2004_epi")
     109                                                                       ("grandi_pasqualini_bers_2010_ss")
     110                                                                       ("shannon_wang_puglisi_weber_bers_2004")
     111                                                                       ("iyer_model_2007");
    91112
    92113        // Model file locations
     
    98119        {
    99120}}}
    100 Find the CellML file for this model.
     121Find the CellML file for this model, and set up a cell factory that uses this model for the whole mesh.
    101122{{{
    102123#!cpp
     
    107128            if (reset_cvode)
    108129            {
    109                 output_folder_stream<< "_Reset";
     130                output_folder_stream << "_Reset";
    110131            }
    111132            std::string output_folder = output_folder_stream.str();
    112133
    113             // Use a different constructor
    114134            OutputFileHandler base_model_handler(output_folder);
    115135            DynamicModelCellFactory cell_factory(model_to_use,
     
    120140
    121141}}}
    122 We will auto-generate a mesh this time, and pass it in, rather than
    123 provide a mesh file name. This is how to generate a cuboid mesh with
    124 a given spatial stepsize h.
     142We will auto-generate a mesh this time, and pass it in, rather than provide a mesh file name.
     143This is how to generate a cuboid mesh with a given spatial stepsize ''h''.
    125144
    126145Using a `DistributedTetrahedralMesh` is faster than `TetrahedralMesh` when running on multiple processes.
     
    132151#!cpp
    133152            DistributedTetrahedralMesh<1,1> mesh;
    134 
    135 }}}
    136 We need to change the printing time steps if going to less than 0.1 ms output, but note
     153            mesh.ConstructRegularSlabMesh(h, 1 /*length*/);
     154            HeartConfig::Instance()->SetOutputUsingOriginalNodeOrdering(true);
     155
     156}}}
     157We need to change the printing time steps if going to be less than 0.1 ms output, but note
    137158that the AP calculations will necessarily be affected by this. We could do them all
    138159at 1ms output, then they would be more consistent, but far less accurate!
     
    149170            }
    150171
    151             mesh.ConstructRegularSlabMesh(h, 1 /*length*/);
    152             HeartConfig::Instance()->SetOutputUsingOriginalNodeOrdering(true);
    153 
    154 }}}
    155 Set the simulation duration, etc, and create an instance of the cell factory.
     172}}}
     173Set the simulation duration, etc.
    156174One thing that should be noted for monodomain problems, the ''intracellular
    157175conductivity'' is used as the monodomain effective conductivity (not a
    158176harmonic mean of intra and extracellular conductivities). So if you want to
    159 alter the monodomain conductivity call
    160 `HeartConfig::Instance()->SetIntracellularConductivities()`
     177alter the monodomain conductivity call `HeartConfig::Instance()->SetIntracellularConductivities()`.
    161178
    162179{{{
     
    168185
    169186}}}
    170 Now we declare the problem class
     187Now we declare the problem class.
    171188{{{
    172189#!cpp
     
    175192}}}
    176193If a mesh-file-name hasn't been set using `HeartConfig`, we have to pass in
    177 a mesh using the `SetMesh` method (must be called before `Initialise`).
     194a mesh using the `SetMesh` method (which must be called before `Initialise`).
    178195{{{
    179196#!cpp
     
    183200`SetWriteInfo` is a useful method that means that the min/max voltage is
    184201printed as the simulation runs (useful for verifying that cells are stimulated
    185 and the wave propagating, for example) (although note scons does buffer output
    186 before printing to screen)
     202and the wave is propagating, for example). Note that `scons` does buffer output
     203before printing to screen, so don't worry if you don't see any output for a while!
    187204{{{
    188205#!cpp
     
    190207
    191208}}}
    192 Finally, call `Initialise` and `Solve` as before
     209Finally, call `Initialise` to finish the problem setup.
    193210{{{
    194211#!cpp
     
    198215            {
    199216}}}
    200 The cells should now be set up...
    201 We have to hack in to call a method on each one
    202 {{{
    203 #!cpp
     217The cells should now be set up.
     218We have to hack in to call a method on each one.
     219{{{
     220#!cpp
     221                // TODO: use FinaliseCellCreation instead!
    204222                DistributedVectorFactory* p_factory = mesh.GetDistributedVectorFactory();
    205223                Vec monodomain_vec = p_factory->CreateVec();
     
    225243            }
    226244
    227             OutputFileHandler handler(output_folder, false);
    228 
    229 }}}
    230 Repository data location
    231 {{{
    232 #!cpp
    233             FileFinder this_file(__FILE__);
    234             FileFinder repo_data("data/reference_traces", this_file);
    235 
    236 }}}
    237 Read some of the outputted data back in, and evaluate AP properties at the last node,
    238 as per the single cell stuff.
     245}}}
     246Having simulated the system, we now read some of the results data (which gets written to disk) back in,
     247and evaluate AP properties at the last node, as per the single cell simulations.
    239248
    240249{{{
     
    242251            Hdf5DataReader data_reader = monodomain_problem.GetDataReader();
    243252            std::vector<double> times = data_reader.GetUnlimitedDimensionValues();
    244 
    245253            std::vector<double> last_node = data_reader.GetVariableOverTime("V", mesh.GetNumNodes()-1u);
    246254
    247             // Output the raw AP data
     255}}}
     256First we write the raw AP data to file, in the same folder as the full results.  To do so we need to get a
     257handler for that folder, but pass `false` as the second argument to avoid removing the existing results!
     258
     259{{{
     260#!cpp
     261            OutputFileHandler handler(output_folder, false);
    248262            std::stringstream file_suffix;
    249263            file_suffix << "_tissue_pde_" << pde_timestep << "_h_" << h;
     
    259273            p_file->close();
    260274
    261             // Now copy it into the repository.
     275}}}
     276Now we copy this file into the repository, finding the data folder relative to this test file.
     277{{{
     278#!cpp
     279            FileFinder this_file(__FILE__);
     280            FileFinder repo_data("data/reference_traces", this_file);
    262281            FileFinder ref_data = handler.FindFile(model + file_suffix.str() + ".dat");
    263282            ref_data.CopyTo(repo_data);
    264283
    265284}}}
    266 Check that the solution looks like an action potential.
     285Next, check that the solution looks like an action potential, and save summary statistics to file.
    267286{{{
    268287#!cpp
     
    312331#include <cxxtest/TestSuite.h>
    313332
     333#include <vector>
     334#include <string>
     335#include <sstream>
    314336#include <boost/assign/list_of.hpp>
    315337#include <boost/foreach.hpp>
     
    324346#include "AbstractCvodeCell.hpp"
    325347
    326 // This header is needed to allow us to run in parallel
     348// Other Chaste headers
     349#include "FileFinder.hpp"
     350#include "OutputFileHandler.hpp"
     351#include "CommandLineArguments.hpp"
     352
     353// This header is needed to allow us to prevent being run in parallel
    327354#include "PetscSetupAndFinalize.hpp"
    328355
     
    351378        }
    352379
    353         /**
    354          * This test was run with the following values inserted here:
    355          * -> 0.01 cm (fine)
    356          * -> 0.001 cm (much finer than normal)
    357          */
    358380        double h;
    359381        if (CommandLineArguments::Instance()->OptionExists("--spacestep"))
     
    366388        }
    367389
    368         // A list of models that we want to do tissue simulations with.
    369390        std::vector<std::string> models_to_use = boost::assign::list_of("luo_rudy_1991")
    370                                                  ("beeler_reuter_model_1977")
    371                                                  ("nygren_atrial_model_1998")
    372                                                  ("ten_tusscher_model_2004_epi")
    373                                                  ("grandi_pasqualini_bers_2010_ss")
    374                                                  ("shannon_wang_puglisi_weber_bers_2004")
    375                                                  ("iyer_model_2007");
     391                                                                       ("beeler_reuter_model_1977")
     392                                                                       ("nygren_atrial_model_1998")
     393                                                                       ("ten_tusscher_model_2004_epi")
     394                                                                       ("grandi_pasqualini_bers_2010_ss")
     395                                                                       ("shannon_wang_puglisi_weber_bers_2004")
     396                                                                       ("iyer_model_2007");
    376397
    377398        // Model file locations
     
    388409            if (reset_cvode)
    389410            {
    390                 output_folder_stream<< "_Reset";
     411                output_folder_stream << "_Reset";
    391412            }
    392413            std::string output_folder = output_folder_stream.str();
    393414
    394             // Use a different constructor
    395415            OutputFileHandler base_model_handler(output_folder);
    396416            DynamicModelCellFactory cell_factory(model_to_use,
     
    401421
    402422            DistributedTetrahedralMesh<1,1> mesh;
    403 
    404             if (pde_timestep < 0.1)
    405             {
    406                 HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(pde_timestep, pde_timestep, 0.1);
    407             }
    408             else
    409             {
    410                 HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(pde_timestep, pde_timestep, pde_timestep);
    411             }
    412 
    413423            mesh.ConstructRegularSlabMesh(h, 1 /*length*/);
    414424            HeartConfig::Instance()->SetOutputUsingOriginalNodeOrdering(true);
     425
     426            if (pde_timestep < 0.1)
     427            {
     428                HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(pde_timestep, pde_timestep, 0.1);
     429            }
     430            else
     431            {
     432                HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(pde_timestep, pde_timestep, pde_timestep);
     433            }
    415434
    416435            HeartConfig::Instance()->SetSimulationDuration(500); //ms
     
    429448            if (reset_cvode)
    430449            {
     450                // TODO: use FinaliseCellCreation instead!
    431451                DistributedVectorFactory* p_factory = mesh.GetDistributedVectorFactory();
    432452                Vec monodomain_vec = p_factory->CreateVec();
     
    452472            }
    453473
    454             OutputFileHandler handler(output_folder, false);
    455 
    456             FileFinder this_file(__FILE__);
    457             FileFinder repo_data("data/reference_traces", this_file);
    458 
    459474            Hdf5DataReader data_reader = monodomain_problem.GetDataReader();
    460475            std::vector<double> times = data_reader.GetUnlimitedDimensionValues();
    461 
    462476            std::vector<double> last_node = data_reader.GetVariableOverTime("V", mesh.GetNumNodes()-1u);
    463477
    464             // Output the raw AP data
     478            OutputFileHandler handler(output_folder, false);
    465479            std::stringstream file_suffix;
    466480            file_suffix << "_tissue_pde_" << pde_timestep << "_h_" << h;
     
    476490            p_file->close();
    477491
    478             // Now copy it into the repository.
     492            FileFinder this_file(__FILE__);
     493            FileFinder repo_data("data/reference_traces", this_file);
    479494            FileFinder ref_data = handler.FindFile(model + file_suffix.str() + ".dat");
    480495            ref_data.CopyTo(repo_data);