# Coverage for pySDC/projects/Monodomain/problem_classes/space_discretizazions/Parabolic_DCT.py: 82%

## 149 statements

, created at 2024-09-09 14:59 +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):

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

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:

311 return True

312 else:

313 return False

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