Coverage for pySDC/projects/Monodomain/problem_classes/space_discretizazions/Parabolic_DCT.py: 82%
149 statements
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +0000
1import numpy as np
2import scipy as sp
3from pySDC.core.common import RegisterParams
4from pySDC.implementations.datatype_classes.mesh import mesh
5from pathlib import Path
6import os
9class Parabolic_DCT(RegisterParams):
10 """
11 A class for the spatial discretization of the parabolic part of the monodomain equation.
12 Here we discretize the spatial domain with a uniform mesh and use the discrete cosine transform (DCT)
13 to discretize the Laplacian operator. The DCT is a real-to-real type of Fourier transform that is well suited for
14 Neumann boundary conditions.
16 Parameters:
17 -----------
18 problem_params: dict containing the problem parameters
20 Attributes:
21 -----------
22 chi: float
23 Surface-to-volume ratio of the cell membrane
24 Cm: float
25 Membrane capacitance
26 si_l: float
27 Longitudinal intracellular conductivity
28 se_l: float
29 Longitudinal extracellular conductivity
30 si_t: float
31 Transversal intracellular conductivity
32 se_t: float
33 Transversal extracellular conductivity
34 sigma_l: float
35 Longitudinal conductivity
36 sigma_t: float
37 Transversal conductivity
38 diff_l: float
39 Longitudinal diffusion coefficient
40 diff_t: float
41 Transversal diffusion coefficient
42 diff: tuple of floats
43 Tuple containing the diffusion coefficients
44 dom_size: tuple of floats
45 Tuple containing the domain size
46 n_elems: tuple of ints
47 Tuple containing the number of elements in each direction
48 refinements: int
49 Number of refinements with respect to a baseline mesh. Can be both positive (to get finer meshes) and negative (to get coarser meshes).
50 grids: tuple of 1D arrays
51 Tuple containing the grid points in each direction
52 dx: tuple of floats
53 Tuple containing the grid spacings in each direction
54 shape: tuple of ints
55 Tuple containing the number of grid points in each direction. Same as n_elems with with reversed order.
56 n_dofs: int
57 Total number of degrees of freedom
58 dim: int
59 Dimension of the spatial domain. Can be 1, 2 or 3.
60 init: tuple
61 Shape of the mesh, None and data type of the mesh (np.double)
62 mesh_name: str
63 Name of the mesh. Can be cube_ND, cubdoid_ND, cubdoid_ND_smaller, cubdoid_ND_small, cubdoid_ND_large, cubdoid_ND_very_large. Where N=1,2,3.
64 diff_dct: array
65 Array containing the discrete Laplacian operator
66 output_folder: Path
67 Path to the output folder
68 output_file_path: Path
69 Path to the output file
70 output_file: str
71 Name of the output file
72 enable_output: bool
73 If True, the solution is written to file. Else not.
74 t_out: list
75 List containing the output times
76 order: int
77 Order of the spatial discretization. Can be 2 or 4
78 zero_stim_vec: float
79 Used to apply zero stimili.
80 """
82 def __init__(self, **problem_params):
83 self._makeAttributeAndRegister(*problem_params.keys(), localVars=problem_params, readOnly=True)
85 self.define_domain()
86 self.define_coefficients()
87 self.define_diffusion()
88 self.define_stimulus()
90 def __del__(self):
91 if self.enable_output:
92 # Close the output file
93 self.output_file.close()
94 # Save the output times and the grid points
95 with open(
96 self.output_file_path.parent / Path(self.output_file_name + '_txyz').with_suffix(".npy"), 'wb'
97 ) as f:
98 np.save(f, np.array(self.t_out))
99 xyz = self.grids
100 for i in range(self.dim):
101 np.save(f, xyz[i])
103 @property
104 def mesh_name(self):
105 return "ref_" + str(self.refinements)
107 def define_coefficients(self):
108 self.chi = 140.0 # mm^-1
109 self.Cm = 0.01 # uF/mm^2
110 self.si_l = 0.17 # mS/mm
111 self.se_l = 0.62 # mS/mm
112 self.si_t = 0.019 # mS/mm
113 self.se_t = 0.24 # mS/mm
115 if "cube" in self.domain_name:
116 # For this domain we use isotropic conductivities
117 self.si_t = self.si_l
118 self.se_t = self.se_l
120 self.sigma_l = self.si_l * self.se_l / (self.si_l + self.se_l)
121 self.sigma_t = self.si_t * self.se_t / (self.si_t + self.se_t)
122 self.diff_l = self.sigma_l / self.chi / self.Cm
123 self.diff_t = self.sigma_t / self.chi / self.Cm
125 if self.dim == 1:
126 self.diff = (self.diff_l,)
127 elif self.dim == 2:
128 self.diff = (self.diff_l, self.diff_t)
129 else:
130 self.diff = (self.diff_l, self.diff_t, self.diff_t)
132 def define_domain(self):
133 if "cube" in self.domain_name:
134 self.dom_size = (100.0, 100.0, 100.0)
135 self.dim = int(self.domain_name[5])
136 else: # cuboid
137 if "smaller" in self.domain_name:
138 self.dom_size = (10.0, 4.5, 2.0)
139 elif "small" in self.domain_name:
140 self.dom_size = (5.0, 3.0, 1.0)
141 elif "very_large" in self.domain_name:
142 self.dom_size = (280.0, 112.0, 48.0)
143 elif "large" in self.domain_name:
144 self.dom_size = (60.0, 21.0, 9.0)
145 else:
146 self.dom_size = (20.0, 7.0, 3.0)
147 self.dim = int(self.domain_name[7])
149 self.dom_size = self.dom_size[: self.dim]
150 self.n_elems = [int(2 ** np.round(np.log2(5.0 * L * 2**self.refinements))) for L in self.dom_size]
151 self.grids, self.dx = self.get_grids_dx(self.dom_size, self.n_elems)
153 self.shape = tuple(np.flip([x.size for x in self.grids]))
154 self.n_dofs = int(np.prod(self.shape))
155 self.init = ((self.n_dofs,), None, np.dtype('float64'))
157 def define_diffusion(self):
158 N = self.n_elems
159 dx = self.dx
160 dim = len(N)
161 if self.order == 2:
162 diff_dct = self.diff[0] * (2.0 * np.cos(np.pi * np.arange(N[0]) / N[0]) - 2.0) / dx[0] ** 2
163 if dim >= 2:
164 diff_dct = (
165 diff_dct[None, :]
166 + self.diff[1]
167 * np.array((2.0 * np.cos(np.pi * np.arange(N[1]) / N[1]) - 2.0) / dx[1] ** 2)[:, None]
168 )
169 if dim >= 3:
170 diff_dct = (
171 diff_dct[None, :, :]
172 + self.diff[2]
173 * np.array((2.0 * np.cos(np.pi * np.arange(N[2]) / N[2]) - 2.0) / dx[2] ** 2)[:, None, None]
174 )
175 elif self.order == 4:
176 diff_dct = (
177 self.diff[0]
178 * (
179 (-1.0 / 6.0) * np.cos(2.0 * np.pi * np.arange(N[0]) / N[0])
180 + (8.0 / 3.0) * np.cos(np.pi * np.arange(N[0]) / N[0])
181 - 2.5
182 )
183 / dx[0] ** 2
184 )
185 if dim >= 2:
186 diff_dct = (
187 diff_dct[None, :]
188 + self.diff[1]
189 * np.array(
190 (
191 (-1.0 / 6.0) * np.cos(2.0 * np.pi * np.arange(N[1]) / N[1])
192 + (8.0 / 3.0) * np.cos(np.pi * np.arange(N[1]) / N[1])
193 - 2.5
194 )
195 / dx[1] ** 2
196 )[:, None]
197 )
198 if dim >= 3:
199 diff_dct = (
200 diff_dct[None, :, :]
201 + self.diff[2]
202 * np.array(
203 (
204 (-1.0 / 6.0) * np.cos(2.0 * np.pi * np.arange(N[2]) / N[2])
205 + (8.0 / 3.0) * np.cos(np.pi * np.arange(N[2]) / N[2])
206 - 2.5
207 )
208 / dx[2] ** 2
209 )[:, None, None]
210 )
211 else:
212 raise NotImplementedError("Only order 2 and 4 are implemented for Parabolic_DCT.")
214 self.diff_dct = diff_dct
216 def grids_from_x(self, x):
217 dim = len(x)
218 if dim == 1:
219 return (x[0],)
220 elif dim == 2:
221 return (x[0][None, :], x[1][:, None])
222 elif dim == 3:
223 return (x[0][None, None, :], x[1][None, :, None], x[2][:, None, None])
225 def get_grids_dx(self, dom_size, N):
226 # The grid points are the midpoints of the elements, hence x_{n+1/2}=(n+1/2)*dx
227 # This is needed for the DCT.
228 x = [np.linspace(0, dom_size[i], 2 * N[i] + 1) for i in range(len(N))]
229 x = [xi[1::2] for xi in x]
230 dx = [xi[1] - xi[0] for xi in x]
231 return self.grids_from_x(x), dx
233 def define_stimulus(self):
234 self.zero_stim_vec = 0.0
235 # all remaining stimulus parameters are set in MonodomainODE
237 def solve_system(self, rhs, factor, u0, t, u_sol):
238 """
239 Solve the linear system: u_sol = (I - factor * A)^{-1} rhs
241 Arguments:
242 ----------
243 rhs: mesh
244 The right-hand side of the linear system
245 factor: float
246 The factor in the linear system multiplying the Laplacian
247 u0: mesh
248 The initial guess for the solution. Not used here since we use a direct solver.
249 t: float
250 The current time. Not used here since the Laplacian is time-independent.
251 u_sol: mesh
252 The vector to store the solution in.
253 """
254 rhs_hat = sp.fft.dctn(rhs.reshape(self.shape))
255 u_sol_hat = rhs_hat / (1.0 - factor * self.diff_dct)
256 u_sol[:] = sp.fft.idctn(u_sol_hat).ravel()
258 return u_sol
260 def add_disc_laplacian(self, uh, res):
261 """
262 Add the discrete Laplacian operator to res: res += A * uh
263 """
264 res[:] += sp.fft.idctn(self.diff_dct * sp.fft.dctn(uh.reshape(self.shape))).ravel()
266 def init_output(self, output_folder):
267 # Initialize the output parameters and file
268 self.output_folder = output_folder
269 self.output_file_path = self.output_folder / Path(self.output_file_name).with_suffix(".npy")
270 if self.enable_output:
271 if self.output_file_path.is_file():
272 os.remove(self.output_file_path)
273 if not self.output_folder.is_dir():
274 os.makedirs(self.output_folder)
275 self.output_file = open(self.output_file_path, 'wb')
276 self.t_out = []
278 def write_solution(self, uh, t):
279 # Write the solution to file. Meant for visualization, hence we print the time stamp too.
280 # Usually we only write the first component of the solution, hence the electric potential V and not the ionic model state variables.
281 if self.enable_output:
282 np.save(self.output_file, uh.reshape(self.shape))
283 self.t_out.append(t)
285 def write_reference_solution(self, uh, indeces):
286 """
287 Write the solution to file. This is meant to print solutions that will be used as reference solutions
288 and compute errors or as well solutions that will be used as initial values for other simulations.
289 Here we write the model variables listed in indeces.
290 """
291 if self.output_file_path.is_file():
292 os.remove(self.output_file_path)
293 if not self.output_file_path.parent.is_dir():
294 os.makedirs(self.output_file_path.parent)
295 with open(self.output_file_path, 'wb') as file:
296 [np.save(file, uh[i].reshape(self.shape)) for i in indeces]
298 def read_reference_solution(self, uh, indeces, ref_file_name):
299 """
300 Read a reference solution from file. It can be used to set the initial values of a simulation
301 or to compute errors.
302 We read the model variables listed in indeces.
303 """
304 if ref_file_name == "":
305 return False
306 ref_sol_path = Path(self.output_folder) / Path(ref_file_name).with_suffix(".npy")
307 if ref_sol_path.is_file():
308 with open(ref_sol_path, 'rb') as f:
309 for i in indeces:
310 uh[i][:] = np.load(f).ravel()
311 return True
312 else:
313 return False
315 def stim_region(self, stim_center, stim_radius):
316 """
317 Define the region where the stimulus is applied, given the center and the radius of the stimulus.
318 Returns a vector of the same size as the grid, with 1s inside the stimulus region and 0s outside.
319 """
320 grids = self.grids
321 coord_inside_stim_box = []
322 for i in range(len(grids)):
323 coord_inside_stim_box.append(abs(grids[i] - stim_center[i]) < stim_radius[i])
325 inside_stim_box = True
326 for i in range(len(grids)):
327 inside_stim_box = np.logical_and(inside_stim_box, coord_inside_stim_box[i])
329 stim = mesh(self.init)
330 stim[:] = inside_stim_box.ravel().astype(float)
332 return stim
334 def compute_errors(self, uh, ref_sol):
335 # Compute L2 errors with respect to the reference solution
336 error_L2 = np.linalg.norm(uh - ref_sol)
337 sol_norm_L2 = np.linalg.norm(ref_sol)
338 rel_error_L2 = error_L2 / sol_norm_L2
340 return error_L2, rel_error_L2