|
| 1 | +Exponential SDC for the Monodomain Equation in Cardiac Electrophysiology |
| 2 | +============================================================================== |
| 3 | +This project implements the exponential spectral deferred correction (ESDC) method for the monodomain equation in cardiac electrophysiology. |
| 4 | +The method proposed here is an adaptation of the `ESDC method proposed by T. Buvoli <https://doi.org/10.1137/19M1256166>`_ to the monodomain equation. |
| 5 | +In particular, the implicit-explicit Rush-Larsen method is used as correction scheme. Moreover, not all model components have exponential terms, therefore the resulting method is an hybrid between ESDC and the standard SDC method. |
| 6 | + |
| 7 | +Monodomain equation |
| 8 | +------------------- |
| 9 | +The monodomain equation models the electrical activity in the heart. It is a reaction-diffusion equation coupled with an ordinary differential equation for the ionic model and is given by |
| 10 | + |
| 11 | +.. math:: |
| 12 | + \begin{align} |
| 13 | + \chi (C_m\frac{\partial V}{\partial t}+I_{ion}(V,z_E,z_e, t)) &= \nabla \cdot (D \nabla V) & \quad \text{in } &\Omega \times (0,T), \\ |
| 14 | + \frac{\partial z_E}{\partial t} &= g_E(V,z_E,z_e) & \quad \text{in } &\Omega \times (0,T), \\ |
| 15 | + \frac{\partial z_e}{\partial t} &= \Lambda_e(V)(z_e-z_{e,\infty}(V)) & \quad \text{in } &\Omega \times (0,T), \\ |
| 16 | + \end{align} |
| 17 | + |
| 18 | +plus the boundary conditions, where :math:`V(t,x)\in\mathbb{R}` is the transmembrane potential and :math:`z_E(t,x)\in\mathbb{R}^n`, :math:`z_e(t,x)\in\mathbb{R}^m` are the ionic model state variables. |
| 19 | +The ionic model right-hand side :math:`g_E` is a general nonlinear term, while :math:`\Lambda_e` is a diagonal matrix. The typical range for the number of unknowns :math:`N=1+n+m` is :math:`N\in [4,50]` and depends on the ionic model of choice. |
| 20 | + |
| 21 | +Spatial discretization yields a system of ODEs which can be written in compact form as |
| 22 | + |
| 23 | +.. math:: |
| 24 | + \mathbf y'=f_I(\mathbf y)+f_E(\mathbf y)+f_e(\mathbf y), |
| 25 | +
|
| 26 | +where :math:`\mathbf y(t)\in\mathbb{R}^{M N}` is the vector of unknowns and :math:`M` the number of mesh nodes. |
| 27 | +Concerning the right-hand sides, :math:`f_I` is a linear term for the discrete diffusion, :math:`f_E` is a nonlinear but non-stiff term for :math:`I_{ion},g_E`, and :math:`f_e` is a severely stiff term for :math:`\Lambda_e(V)(z_e-z_{e,\infty}(V))`. |
| 28 | + |
| 29 | +The standard (serial) way of integrating the monodomain equation is by using a splitting method, where :math:`f_I` is integrated implicitly, :math:`f_E` explicitly, and :math:`f_e` using the exponential Euler method (which is inexpensive due to the diagonal structure of :math:`\Lambda_e`). We denote this method as IMEXEXP. |
| 30 | + |
| 31 | +The ESDC method for the monodomain equation |
| 32 | +------------------------------------------- |
| 33 | +A possible way to parallelize the integration of the monodomain equation is by employing the SDC method in combination with the IMEXEXP approach for the correction scheme (preconditioner). |
| 34 | +However, this approach is unstable due to the severe stiffness of :math:`f_e`. |
| 35 | +Therefore we propose a hybrid method, where we employ SDC for the :math:`f_I,f_E` terms and ESDC for the :math:`f_e` term. For the correcttion scheme we still use the IMEXEXP method. |
| 36 | +The resulting method can be seen as a particular case of ESDC and will be denoted by ESDC in the next figures, for simplicity. |
| 37 | + |
| 38 | +Running the code |
| 39 | +---------------- |
| 40 | +Due to their complexity, ionic models are coded in C++ and wrapped to Python. Therefore, before running any example you need to compile the ionic models by running the following command in the root folder: |
| 41 | + |
| 42 | +.. code-block:: |
| 43 | +
|
| 44 | + export IONIC_MODELS_PATH=pySDC/projects/Monodomain/problem_classes/ionicmodels/cpp |
| 45 | + c++ -O3 -Wall -shared -std=c++11 -fPIC -fvisibility=hidden $(python3 -m pybind11 --includes) ${IONIC_MODELS_PATH}/bindings_definitions.cpp -o ${IONIC_MODELS_PATH}/ionicmodels$(python3-config --extension-suffix) |
| 46 | +
|
| 47 | +Then an example can be run: |
| 48 | + |
| 49 | +.. code-block:: |
| 50 | +
|
| 51 | + cd pySDC/projects/Monodomain/run_scripts |
| 52 | + mpirun -n 4 python run_MonodomainODE_cli.py --dt 0.05 --end_time 0.2 --num_nodes 6,3 --domain_name cube_1D --refinements 0 --ionic_model_name TTP --truly_time_parallel --n_time_ranks 4 |
| 53 | +
|
| 54 | +Stability |
| 55 | +--------- |
| 56 | +We display here the stability domain of the ESDC and SDC methods, both with IMEXEXP as correction scheme, applied to the test problem |
| 57 | + |
| 58 | +.. math:: |
| 59 | + y'=\lambda_I y+\lambda_E y+\lambda_e y, |
| 60 | + |
| 61 | +with :math:`\lambda_I,\lambda_E,\lambda_e` representing :math:`f_I,f_E,f_e`, respectively. |
| 62 | +We fix :math:`\lambda_E=-1` and vary the stiff terms :math:`\lambda_I,\lambda_e` only. We see that the ESDC method is stable for all tested values of :math:`\lambda_I,\lambda_e`, while SDC is not. |
| 63 | + |
| 64 | +.. image:: ../../../data/stability_domain_IMEXEXP_EXPRK.png |
| 65 | + :scale: 60 % |
| 66 | +.. image:: ../../../data/stability_domain_IMEXEXP.png |
| 67 | + :scale: 60 % |
| 68 | + |
| 69 | +Convergence |
| 70 | +----------- |
| 71 | +Here we verify convergence of the ESDC method for the monodomain equation. |
| 72 | +We fix the number of collocation nodes to :math:`m=6` and perform a convergence experiment fixing the number of sweeps to either :math:`k=3` or :math:`k=6`. |
| 73 | +We use the ten Tusscher-Panfilov ionic model, which is employed in practical applications. |
| 74 | +We see that we gain one order of accuracy per sweep, as expected. |
| 75 | + |
| 76 | +.. image:: ../../../data/convergence_ESDC_fixed_iter.png |
| 77 | + :scale: 100 % |
| 78 | + |
| 79 | + |
| 80 | +Iterations |
| 81 | +---------- |
| 82 | +Here we consider three methods: |
| 83 | + |
| 84 | +* ESDC: with :math:`m=6` collocation nodes. |
| 85 | +* MLESDC: This is a multilevel version of ESDC with :math:`m=6` collocation nodes on the fine level and :math:`m=3` nodes on the coarse level. |
| 86 | +* PFASST: Combination of the PFASST parallelization method with MLESDC, using 24 processors. |
| 87 | + |
| 88 | +We display the number of iterations required by each method to reach a given tolerance and the residual at convergence. As ionic model we use again the ten Tusscher-Panfilov model. |
| 89 | +We see that PFASST requires a reasonalbly small number of iterations, comparable to the serial counterparts ESDC and MLESDC. |
| 90 | + |
| 91 | +.. image:: ../../../data/niter_VS_time.png |
| 92 | + :scale: 100 % |
| 93 | +.. image:: ../../../data/res_VS_time.png |
| 94 | + :scale: 100 % |
0 commit comments