Coverage for pySDC/projects/Monodomain/problem_classes/MonodomainODE.py: 95%
186 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1from pathlib import Path
2import logging
3import numpy as np
4from pySDC.core.Problem import ptype
5from pySDC.implementations.datatype_classes.mesh import mesh
6from pySDC.projects.Monodomain.datatype_classes.my_mesh import imexexp_mesh
7from pySDC.projects.Monodomain.problem_classes.space_discretizazions.Parabolic_DCT import Parabolic_DCT
8import pySDC.projects.Monodomain.problem_classes.ionicmodels.cpp as ionicmodels
11class MonodomainODE(ptype):
12 """
13 A class for the discretization of the Monodomain equation. The Monodomain equation is a parabolic PDE composed of
14 a reaction-diffusion equation coupled with an ODE system. The unknowns are the potential V and the ionic model variables (g_1,...,g_N).
15 The reaction-diffusion equation is discretized in another class, where any spatial discretization can be used.
16 The ODE system is the ionic model, which doesn't need spatial discretization, being a system of ODEs.
20 Attributes:
21 -----------
22 parabolic: The parabolic problem class used to discretize the reaction-diffusion equation
23 ionic_model: The ionic model used to discretize the ODE system. This is a wrapper around the actual ionic model, which is written in C++.
24 size: The number of variables in the ionic model
25 vector_type: The type of vector used to store a single unknown (e.g. V). This data type depends on spatial discretization, hence on the parabolic class.
26 dtype_u: The type of vector used to store all the unknowns (V,g_1,...,g_N). This is a vector of vector_type.
27 dtype_f: The type of vector used to store the right-hand side of the ODE system stemming from the monodomain equation. This is a vector of vector_type.
28 output_folder: The folder where the solution is written to file
29 t0: The initial simulation time. This is 0.0 by default but can be changed by the user in order to skip the initial stimulus.
30 Tend: The duration of the simulation.
31 output_V_only: If True, only the potential V is written to file. If False, all the ionic model variables are written to file.
32 read_init_val: If True, the initial value is read from file. If False, the initial value is at equilibrium.
33 """
35 dtype_u = mesh
36 dtype_f = mesh
38 def __init__(self, **problem_params):
39 self.logger = logging.getLogger("step")
41 self.parabolic = Parabolic_DCT(**problem_params)
43 self.define_ionic_model(problem_params["ionic_model_name"])
45 self.init = ((self.size, *self.parabolic.init[0]), self.parabolic.init[1], self.parabolic.init[2])
47 # invoke super init
48 super(MonodomainODE, self).__init__(self.init)
49 # store all problem params dictionary values as attributes
50 self._makeAttributeAndRegister(*problem_params.keys(), localVars=problem_params, readOnly=True)
52 self.define_stimulus()
54 # initial and end time
55 self.t0 = 0.0
56 self.Tend = 50.0 if self.end_time < 0.0 else self.end_time
58 # init output stuff
59 self.output_folder = (
60 Path(self.output_root)
61 / Path(self.parabolic.domain_name)
62 / Path(self.parabolic.mesh_name)
63 / Path(self.ionic_model_name)
64 )
65 self.parabolic.init_output(self.output_folder)
67 def init_exp_extruded(self, new_dim_shape):
68 # The info needed to initialize a new vector of size (M,N) where M is the number of variables in the
69 # ionic model with exponential terms and N is the number of dofs in the mesh.
70 # The vector is further extruded to additional dimensions with shape new_dim_shape.
71 return ((*new_dim_shape, len(self.rhs_exp_indeces), self.init[0][1]), self.init[1], self.init[2])
73 def write_solution(self, uh, t):
74 # write solution to file, only the potential V=uh[0], not the ionic model variables
75 self.parabolic.write_solution(uh[0], t)
77 def write_reference_solution(self, uh, all=False):
78 # write solution to file, only the potential V=uh[0] or all variables if all=True
79 self.parabolic.write_reference_solution(uh, list(range(self.size)) if all else [0])
81 def read_reference_solution(self, uh, ref_file_name, all=False):
82 # read solution from file, only the potential V=uh[0] or all variables if all=True
83 # returns true if read was successful, false else
84 return self.parabolic.read_reference_solution(uh, list(range(self.size)) if all else [0], ref_file_name)
86 def initial_value(self):
87 # Create initial value. Every variable is constant in space
88 u0 = self.dtype_u(self.init)
89 init_vals = self.ionic_model.initial_values()
90 for i in range(self.size):
91 u0[i][:] = init_vals[i]
93 # overwrite the initial value with solution from file if desired
94 if self.read_init_val:
95 read_ok = self.read_reference_solution(u0, self.init_val_name, True)
96 assert read_ok, "ERROR: Could not read initial value from file."
98 return u0
100 def compute_errors(self, uh):
101 """
102 Compute L2 error of uh[0] (potential V)
103 Args:
104 uh (VectorOfVectors): solution as vector of vectors
106 Returns:
107 computed (bool): if error computation was successful
108 error (float): L2 error
109 rel_error (float): relative L2 error
110 """
111 ref_sol_V = self.dtype_u(init=self.init, val=0.0)
112 read_ok = self.read_reference_solution(ref_sol_V, self.ref_sol, False)
113 if read_ok:
114 error_L2, rel_error_L2 = self.parabolic.compute_errors(uh[0], ref_sol_V[0])
116 print(f"L2-errors: {error_L2}")
117 print(f"Relative L2-errors: {rel_error_L2}")
119 return True, error_L2, rel_error_L2
120 else:
121 return False, 0.0, 0.0
123 def define_ionic_model(self, ionic_model_name):
124 self.scale_Iion = 0.01 # used to convert currents in uA/cm^2 to uA/mm^2
125 # scale_im is applied to the rhs of the ionic model, so that the rhs is in units of mV/ms
126 self.scale_im = self.scale_Iion / self.parabolic.Cm
128 if ionic_model_name in ["HodgkinHuxley", "HH"]:
129 self.ionic_model = ionicmodels.HodgkinHuxley(self.scale_im)
130 elif ionic_model_name in ["Courtemanche1998", "CRN"]:
131 self.ionic_model = ionicmodels.Courtemanche1998(self.scale_im)
132 elif ionic_model_name in ["TenTusscher2006_epi", "TTP"]:
133 self.ionic_model = ionicmodels.TenTusscher2006_epi(self.scale_im)
134 elif ionic_model_name in ["TTP_S", "TTP_SMOOTH"]:
135 self.ionic_model = ionicmodels.TenTusscher2006_epi_smooth(self.scale_im)
136 elif ionic_model_name in ["BiStable", "BS"]:
137 self.ionic_model = ionicmodels.BiStable(self.scale_im)
138 else:
139 raise Exception("Unknown ionic model.")
141 self.size = self.ionic_model.size
143 def define_stimulus(self):
145 stim_dur = 2.0
146 if "cuboid" in self.parabolic.domain_name:
147 self.stim_protocol = [[0.0, stim_dur]] # list of stimuli times and stimuli durations
148 self.stim_intensities = [50.0] # list of stimuli intensities
149 self.stim_centers = [[0.0, 0.0, 0.0]] # list of stimuli centers
150 r = 1.5
151 self.stim_radii = [[r, r, r]] * len(
152 self.stim_protocol
153 ) # list of stimuli radii in the three directions (x,y,z)
154 elif "cube" in self.parabolic.domain_name:
155 self.stim_protocol = [[0.0, 2.0], [1000.0, 10.0]]
156 self.stim_intensities = [50.0, 80.0]
157 centers = [[0.0, 50.0, 50.0], [58.5, 0.0, 50.0]]
158 self.stim_centers = [centers[i] for i in range(len(self.stim_protocol))]
159 self.stim_radii = [[1.0, 50.0, 50.0], [1.5, 60.0, 50.0]]
160 else:
161 raise Exception("Unknown domain name.")
163 self.stim_protocol = np.array(self.stim_protocol)
164 self.stim_protocol[:, 0] -= self.init_time # shift stimulus times by the initial time
166 # index of the last stimulus applied. The value -1 means no stimulus has been applied yet.
167 self.last_stim_index = -1
169 def eval_f(self, u, t, fh=None):
170 if fh is None:
171 fh = self.dtype_f(init=self.init, val=0.0)
173 # eval ionic model rhs on u and put result in fh. All indices of the vector of vector fh must be computed (list(range(self.size))
174 self.eval_expr(self.ionic_model.f, u, fh, list(range(self.size)), False)
175 # apply stimulus
176 fh[0] += self.Istim(t)
178 # apply diffusion
179 self.parabolic.add_disc_laplacian(u[0], fh[0])
181 return fh
183 def Istim(self, t):
184 tol = 1e-8
185 for i, (stim_time, stim_dur) in enumerate(self.stim_protocol):
186 # Look for which stimulus to apply at the current time t by checking the stimulus protocol:
187 # Check if t is in the interval [stim_time, stim_time+stim_dur] with a tolerance tol
188 # and apply the corresponding stimulus
189 if (t + stim_dur * tol >= stim_time) and (t + stim_dur * tol < stim_time + stim_dur):
190 # if the stimulus is not the same as the last one applied, update the last_stim_index and the space_stim vector
191 if i != self.last_stim_index:
192 self.last_stim_index = i
193 # get the vector of zeros and ones defining the stimulus region
194 self.space_stim = self.parabolic.stim_region(self.stim_centers[i], self.stim_radii[i])
195 # scale by the stimulus intensity and apply the change of units
196 self.space_stim *= self.scale_im * self.stim_intensities[i]
197 return self.space_stim
199 return self.parabolic.zero_stim_vec
201 def eval_expr(self, expr, u, fh, indeces, zero_untouched_indeces=True):
202 # evaluate the expression expr on u and put the result in fh
203 # Here expr is a wrapper on a C++ function that evaluates the rhs of the ionic model (or part of it)
204 if expr is not None:
205 expr(u, fh)
207 # indeces is a list of integers indicating which variables are modified by the expression expr.
208 # This information is known a priori. Here we use it to zero the variables that are not modified by expr (if zero_untouched_indeces is True)
209 if zero_untouched_indeces:
210 non_indeces = [i for i in range(self.size) if i not in indeces]
211 for i in non_indeces:
212 fh[i][:] = 0.0
215class MultiscaleMonodomainODE(MonodomainODE):
216 """
217 The multiscale version of the MonodomainODE problem. This class is used to solve the monodomain equation with a multirate solver.
218 The main difference with respect to the MonodomainODE class is that the right-hand side of the ODE system is split into three parts:
219 - impl: The discrete Laplacian. This is a stiff term threated implicitly by time integrators.
220 - expl: The non stiff term of the ionic models, threated explicitly by time integrators.
221 - exp: The very stiff but diagonal terms of the ionic models, threated exponentially by time integrators.
222 """
224 dtype_f = imexexp_mesh
226 def __init__(self, **problem_params):
227 super(MultiscaleMonodomainODE, self).__init__(**problem_params)
229 self.define_splittings()
231 self.constant_lambda_and_phi = False
233 def define_splittings(self):
234 """
235 This function defines the splittings used in the problem.
236 The im_* variables are meant for internal use, the rhs_* for external use (i.e. in the sweeper).
237 The *_args and *_indeces are list of integers.
238 The *_args are list of variables that are needed to evaluate a function plus the variables that are modified by the function.
239 The *_indeces are the list of variables that are modified by the function (subset of args).
240 Example: for f(x_0,x_1,x_2,x_3,x_4)=f(x_0,x_2,x_4)=(y_0,y_1,0,0,y_4) we have
241 f_args=[0,1,2,4]=([0,2,4] union [0,1,4]) since f needs x_0,x_2,x_4 and y_0,y_1,y_4 are effective outputs of the function (others are zero).
242 f_indeces=[0,1,4] since only y_0,y_1,y_4 are outputs of the function, y_2,y_3 are zero
244 The ionic model has many variables (say M) and each variable has the same number of dofs as the mesh (say N).
245 Therefore the problem has size N*M and quickly becomes very large. Thanks to args and indeces we can:
246 - avoid to copy the whole vector M*N of variables when we only need a subset, for instance 2*N
247 - avoid unnecessary operations on the whole vector, for instance update only the variables that are effective outputs of a function (indeces),
248 and so on.
250 Yeah, it's a bit a mess, but helpful.
251 """
252 # define nonstiff term (explicit part)
253 # the wrapper to c++ expression that evaluates the nonstiff part of the ionic model
254 self.im_f_nonstiff = self.ionic_model.f_expl
255 # the args of f_expl
256 self.im_nonstiff_args = self.ionic_model.f_expl_args
257 # the indeces of f_expl
258 self.im_nonstiff_indeces = self.ionic_model.f_expl_indeces
260 # define stiff term (implicit part)
261 self.im_f_stiff = None # no stiff part coming from ionic model to be threated implicitly. Indeed, all the stiff terms are diagonal and are threated exponentially.
262 self.im_stiff_args = []
263 self.im_stiff_indeces = []
265 # define exp term (eponential part)
266 # the exponential term is defined by f_exp(u)= lmbda(u)*(u-yinf(u)), hence we only need to define lmbda and yinf
267 # the wrapper to c++ expression that evaluates lmbda(u)
268 self.im_lmbda_exp = self.ionic_model.lmbda_exp
269 # the wrapper to c++ expression that evaluates lmbda(u) and yinf(u)
270 self.im_lmbda_yinf_exp = self.ionic_model.lmbda_yinf_exp
271 # the args of lmbda and yinf (they are the same)
272 self.im_exp_args = self.ionic_model.f_exp_args
273 # the indeces of lmbda and yinf
274 self.im_exp_indeces = self.ionic_model.f_exp_indeces
276 # the spectral radius of the jacobian of non stiff term. We use a bound
277 self.rho_nonstiff_cte = self.ionic_model.rho_f_expl()
279 self.rhs_stiff_args = self.im_stiff_args
280 self.rhs_stiff_indeces = self.im_stiff_indeces
281 # Add the potential V index 0 to the rhs_stiff_args and rhs_stiff_indeces.
282 # Indeed V is used to compute the Laplacian and is affected by the Laplacian, which is the implicit part of the problem.
283 if 0 not in self.rhs_stiff_args:
284 self.rhs_stiff_args = [0] + self.rhs_stiff_args
285 if 0 not in self.rhs_stiff_indeces:
286 self.rhs_stiff_indeces = [0] + self.rhs_stiff_indeces
288 self.rhs_nonstiff_args = self.im_nonstiff_args
289 self.rhs_nonstiff_indeces = self.im_nonstiff_indeces
290 # Add the potential V index 0 to the rhs_nonstiff_indeces. Indeed V is affected by the stimulus, which is a non stiff term.
291 if 0 not in self.rhs_nonstiff_indeces:
292 self.rhs_nonstiff_indeces = [0] + self.rhs_nonstiff_indeces
294 self.im_non_exp_indeces = [i for i in range(self.size) if i not in self.im_exp_indeces]
296 self.rhs_exp_args = self.im_exp_args
297 self.rhs_exp_indeces = self.im_exp_indeces
299 self.rhs_non_exp_indeces = self.im_non_exp_indeces
301 # a vector of ones, useful
302 self.one = self.dtype_u(init=self.init, val=1.0)
304 # some space to store lmbda and yinf
305 self.lmbda = self.dtype_u(init=self.init, val=0.0)
306 self.yinf = self.dtype_u(init=self.init, val=0.0)
308 def solve_system(self, rhs, factor, u0, t, u_sol=None):
309 """
310 Solve the system u_sol[0] = (M-factor*A)^{-1} * M * rhs[0]
311 and sets u_sol[i] = rhs[i] for i>0 (as if A=0 for i>0)
313 Arguments:
314 rhs (dtype_u): right-hand side
315 factor (float): factor multiplying the Laplacian
316 u0 (dtype_u): initial guess
317 t (float): current time
318 u_sol (dtype_u, optional): some space to store the solution. If None, a new space is allocated. Can be the same as rhs.
319 """
320 if u_sol is None:
321 u_sol = self.dtype_u(init=self.init, val=0.0)
323 self.parabolic.solve_system(rhs[0], factor, u0[0], t, u_sol[0])
325 if rhs is not u_sol:
326 for i in range(1, self.size):
327 u_sol[i][:] = rhs[i][:]
329 return u_sol
331 def eval_f(self, u, t, eval_impl=True, eval_expl=True, eval_exp=True, fh=None, zero_untouched_indeces=True):
332 """
333 Evaluates the right-hand side terms.
335 Arguments:
336 u (dtype_u): the current solution
337 t (float): the current time
338 eval_impl (bool, optional): if True, evaluates the implicit part of the right-hand side. Default is True.
339 eval_expl (bool, optional): if True, evaluates the explicit part of the right-hand side. Default is True.
340 eval_exp (bool, optional): if True, evaluates the exponential part of the right-hand side. Default is True.
341 fh (dtype_f, optional): space to store the right-hand side. If None, a new space is allocated. Default is None.
342 zero_untouched_indeces (bool, optional): if True, the variables that are not modified by the right-hand side are zeroed. Default is True.
343 """
345 if fh is None:
346 fh = self.dtype_f(init=self.init, val=0.0)
348 if eval_expl:
349 fh.expl = self.eval_f_nonstiff(u, t, fh.expl, zero_untouched_indeces)
351 if eval_impl:
352 fh.impl = self.eval_f_stiff(u, t, fh.impl, zero_untouched_indeces)
354 if eval_exp:
355 fh.exp = self.eval_f_exp(u, t, fh.exp, zero_untouched_indeces)
357 return fh
359 def eval_f_nonstiff(self, u, t, fh_nonstiff, zero_untouched_indeces=True):
360 # eval ionic model nonstiff terms
361 self.eval_expr(self.im_f_nonstiff, u, fh_nonstiff, self.im_nonstiff_indeces, zero_untouched_indeces)
363 if not zero_untouched_indeces and 0 not in self.im_nonstiff_indeces:
364 fh_nonstiff[0][:] = 0.0
366 # apply stimulus
367 fh_nonstiff[0] += self.Istim(t)
369 return fh_nonstiff
371 def eval_f_stiff(self, u, t, fh_stiff, zero_untouched_indeces=True):
372 # eval ionic model stiff terms
373 self.eval_expr(self.im_f_stiff, u, fh_stiff, self.im_stiff_indeces, zero_untouched_indeces)
375 if not zero_untouched_indeces and 0 not in self.im_stiff_indeces:
376 fh_stiff[0][:] = 0.0
378 # apply diffusion
379 self.parabolic.add_disc_laplacian(u[0], fh_stiff[0])
381 return fh_stiff
383 def eval_f_exp(self, u, t, fh_exp, zero_untouched_indeces=True):
384 # eval ionic model exp terms f_exp(u)= lmbda(u)*(u-yinf(u)
385 self.eval_lmbda_yinf_exp(u, self.lmbda, self.yinf)
386 for i in self.im_exp_indeces:
387 fh_exp[i][:] = self.lmbda[i] * (u[i] - self.yinf[i])
389 if zero_untouched_indeces:
390 fh_exp[self.im_non_exp_indeces] = 0.0
392 return fh_exp
394 def lmbda_eval(self, u, t, lmbda=None):
395 if lmbda is None:
396 lmbda = self.dtype_u(init=self.init, val=0.0)
398 self.eval_lmbda_exp(u, lmbda)
400 lmbda[self.im_non_exp_indeces] = 0.0
402 return lmbda
404 def eval_lmbda_yinf_exp(self, u, lmbda, yinf):
405 self.im_lmbda_yinf_exp(u, lmbda, yinf)
407 def eval_lmbda_exp(self, u, lmbda):
408 self.im_lmbda_exp(u, lmbda)