# Coverage for pySDC/implementations/problem_classes/Quench.py: 77%

## 152 statements

, created at 2024-09-09 14:59 +0000

1import numpy as np

2import scipy.sparse as sp

3from scipy.sparse.linalg import spsolve, gmres

4from scipy.linalg import inv

6from pySDC.core.errors import ProblemError

7from pySDC.core.problem import Problem, WorkCounter

8from pySDC.helpers import problem_helper

9from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh

12# noinspection PyUnusedLocal

13class Quench(Problem):

14 """

15 This is a toy problem [1]_ to emulate a magnet that has been cooled to temperatures where superconductivity is possible.

16 However, there is a leak! Some point in the domain is constantly heated and when this has heated up its environment

17 sufficiently, there will be a runaway effect heating up the entire magnet.

18 This effect has actually lead to huge magnets being destroyed at CERN in the past and hence warrants investigation.

20 The model we use is a 1d heat equation with Neumann-zero boundary conditions, meaning this magnet is totally

21 insulated from its environment except for the leak.

22 We add a non-linear term that heats parts of the domain that exceed a certain temperature threshold as well as the

23 leak itself.

25 The problem is discretised with finite difference in space and treated *fully-implicitly*.

27 Parameters

28 ----------

29 Cv : float, optional

30 Volumetric heat capacity.

31 K : float, optional

32 Thermal conductivity.

33 u_thresh : float, optional

34 Threshold for temperature.

35 u_max : float, optional

36 Maximum temperature.

37 Q_max : float, optional

38 Maximum heat source power density.

39 leak_range : tuple of float

40 Range of the leak.

41 leak_type : str, optional

42 Type of leak, choose between 'linear' or 'exponential'.

43 leak_transition : str, optional

44 Indicates how the heat in the leak propagates, choose between 'step' and 'Gaussian'.

45 order : int, optional

46 Order of the finite difference discretization.

47 stencil_type : str, optional

48 Type of stencil for finite differences.

49 bc : str, optional

50 Type of boundary conditions. Default is 'neumann-zero'.

51 nvars : int, optional

52 Spatial resolution.

53 newton_tol : float, optional

54 Tolerance for Newton to terminate.

55 newton_maxiter : int, optional

56 Maximum number of Newton iterations to be done.

57 lintol : float, optional

58 Tolerance for linear solver to be done.

59 liniter : int, optional

60 Maximum number of linear iterations inside the Newton solver.

61 direct_solver : bool, optional

62 Indicates if a direct solver should be used.

63 inexact_linear_ratio : float, optional

64 Ratio of tolerance of linear solver to the Newton residual, overrides lintol

65 min_lintol : float, optional

66 Minimal tolerance for the linear solver

67 reference_sol_type : str, optional

68 Indicates which method should be used to compute a reference solution.

69 Choose between 'scipy', 'SDC', or 'DIRK'.

71 Attributes

72 ----------

73 A : sparse matrix (CSC)

74 FD discretization matrix of the ND grad operator.

75 Id : sparse matrix (CSC)

76 Identity matrix of the same dimension as A.

77 dx : float

78 Distance between two spatial nodes.

79 xv : np.1darray

80 Spatial grid values.

81 leak : np.1darray of bool

82 Indicates the leak.

84 References

85 ----------

86 .. [1] Thermal thin shell approximation towards finite element quench simulation. E. Schnaubelt, M. Wozniak, S. Schöps.

87 Supercond. Sci. Technol. 36 044004. DOI 10.1088/1361-6668/acbeea

88 """

90 dtype_u = mesh

91 dtype_f = mesh

93 def __init__(

94 self,

95 Cv=1000.0,

96 K=1000.0,

97 u_thresh=3e-2,

98 u_max=6e-2,

99 Q_max=1.0,

100 leak_range=(0.45, 0.55),

101 leak_type='linear',

102 leak_transition='step',

103 order=2,

104 stencil_type='center',

105 bc='neumann-zero',

106 nvars=2**7,

107 newton_tol=1e-8,

108 newton_maxiter=99,

109 lintol=1e-8,

110 liniter=99,

111 direct_solver=True,

112 inexact_linear_ratio=None,

113 min_lintol=1e-12,

114 reference_sol_type='scipy',

115 ):

116 """

117 Initialization routine

119 Args:

120 problem_params (dict): custom parameters for the example

121 dtype_u: mesh data type (will be passed parent class)

122 dtype_f: mesh data type (will be passed parent class)

123 """

124 # invoke super init, passing number of dofs, dtype_u and dtype_f

125 super().__init__(init=(nvars, None, np.dtype('float64')))

126 self._makeAttributeAndRegister(

127 'Cv',

128 'K',

129 'u_thresh',

130 'u_max',

131 'Q_max',

132 'leak_range',

133 'leak_type',

134 'leak_transition',

135 'order',

136 'stencil_type',

137 'bc',

138 'nvars',

139 'direct_solver',

140 'reference_sol_type',

141 localVars=locals(),

143 )

144 self._makeAttributeAndRegister(

145 'newton_tol',

146 'newton_maxiter',

147 'lintol',

148 'liniter',

149 'inexact_linear_ratio',

150 'min_lintol',

151 localVars=locals(),

153 )

155 self._makeAttributeAndRegister(

156 'newton_tol',

157 'newton_maxiter',

158 'lintol',

159 'liniter',

160 'direct_solver',

161 localVars=locals(),

163 )

165 # setup finite difference discretization from problem helper

166 self.dx, xvalues = problem_helper.get_1d_grid(size=self.nvars, bc=self.bc)

168 self.A, self.b = problem_helper.get_finite_difference_matrix(

169 derivative=2,

170 order=self.order,

171 stencil_type=self.stencil_type,

172 dx=self.dx,

173 size=self.nvars,

174 dim=1,

175 bc=self.bc,

176 )

177 self.A *= self.K / self.Cv

179 self.xv = xvalues

180 self.Id = sp.eye(np.prod(self.nvars), format='csc')

182 self.leak = np.logical_and(self.xv > self.leak_range[0], self.xv < self.leak_range[1])

184 self.work_counters['newton'] = WorkCounter()

185 self.work_counters['rhs'] = WorkCounter()

186 if not self.direct_solver:

187 self.work_counters['linear'] = WorkCounter()

189 def eval_f_non_linear(self, u, t):

190 """

191 Get the non-linear part of f.

193 Parameters

194 ----------

195 u : dtype_u

196 Current values of the numerical solution:

197 t : float

198 Current time at which the numerical solution is computed.

200 Returns

201 -------

202 me : dtype_u

203 The non-linear part of the right-hand side.

204 """

205 u_thresh = self.u_thresh

206 u_max = self.u_max

207 Q_max = self.Q_max

208 me = self.dtype_u(self.init)

210 if self.leak_type == 'linear':

211 me[:] = (u - u_thresh) / (u_max - u_thresh) * Q_max

212 elif self.leak_type == 'exponential':

213 me[:] = Q_max * (np.exp(u) - np.exp(u_thresh)) / (np.exp(u_max) - np.exp(u_thresh))

214 else:

215 raise NotImplementedError(f'Leak type \"{self.leak_type}\" not implemented!')

217 me[u < u_thresh] = 0

218 if self.leak_transition == 'step':

219 me[self.leak] = Q_max

220 elif self.leak_transition == 'Gaussian':

221 me[:] = np.max([me, Q_max * np.exp(-((self.xv - 0.5) ** 2) / 3e-2)], axis=0)

222 else:

223 raise NotImplementedError(f'Leak transition \"{self.leak_transition}\" not implemented!')

225 me[u >= u_max] = Q_max

227 me[:] /= self.Cv

229 return me

231 def eval_f(self, u, t):

232 """

233 Evaluate the full right-hand side.

235 Parameters

236 ----------

237 u : dtype_u

238 Current values of the numerical solution.

239 t : float

240 Current time at which the numerical solution is computed.

242 Returns

243 -------

244 f : dtype_f

245 The right-hand side of the problem.

246 """

247 f = self.dtype_f(self.init)

248 f[:] = self.A.dot(u.flatten()).reshape(self.nvars) + self.b + self.eval_f_non_linear(u, t)

250 self.work_counters['rhs']()

251 return f

253 def get_non_linear_Jacobian(self, u):

254 """

255 Evaluate the non-linear part of the Jacobian only.

257 Parameters

258 ----------

259 u : dtype_u

260 Current values of the numerical solution.

262 Returns

263 -------

264 scipy.sparse.csc

265 The derivative of the non-linear part of the solution w.r.t. to the solution.

266 """

267 u_thresh = self.u_thresh

268 u_max = self.u_max

269 Q_max = self.Q_max

270 me = self.dtype_u(self.init)

272 if self.leak_type == 'linear':

273 me[:] = Q_max / (u_max - u_thresh)

274 elif self.leak_type == 'exponential':

275 me[:] = Q_max * np.exp(u) / (np.exp(u_max) - np.exp(u_thresh))

276 else:

277 raise NotImplementedError(f'Leak type {self.leak_type} not implemented!')

279 me[u < u_thresh] = 0

280 if self.leak_transition == 'step':

281 me[self.leak] = 0

282 elif self.leak_transition == 'Gaussian':

283 me[self.leak] = 0

284 me[self.leak][u[self.leak] > Q_max * np.exp(-((self.xv[self.leak] - 0.5) ** 2) / 3e-2)] = 1

285 else:

286 raise NotImplementedError(f'Leak transition \"{self.leak_transition}\" not implemented!')

287 me[u > u_max] = 0

289 me[:] /= self.Cv

291 return sp.diags(me, format='csc')

293 def solve_system(self, rhs, factor, u0, t):

294 r"""

295 Simple Newton solver for :math:(I - factor \cdot f)(\vec{u}) = \vec{rhs}.

297 Parameters

298 ----------

299 rhs : dtype_f

300 Right-hand side.

301 factor : float

302 Abbrev. for the local stepsize (or any other factor required).

303 u0 : dtype_u

304 Initial guess for the iterative solver.

305 t : float

306 Current time (e.g. for time-dependent BCs).

308 Returns

309 -------

310 u : dtype_u

311 The solution as mesh.

312 """

313 u = self.dtype_u(u0)

314 res = np.inf

315 delta = self.dtype_u(self.init, val=0.0)

317 # construct a preconditioner for the space solver

318 if not self.direct_solver:

319 M = inv((self.Id - factor * self.A).toarray())

320 zero = self.dtype_u(self.init, val=0.0)

322 for n in range(0, self.newton_maxiter):

323 # assemble G such that G(u) = 0 at the solution of the step

324 G = u - factor * self.eval_f(u, t) - rhs

325 self.work_counters[

326 'rhs'

327 ].decrement() # Work regarding construction of the Jacobian etc. should count into the Newton iterations only

329 res = np.linalg.norm(G, np.inf)

330 if res <= self.newton_tol and n > 0: # we want to make at least one Newton iteration

331 break

333 if self.inexact_linear_ratio:

334 self.lintol = max([res * self.inexact_linear_ratio, self.min_lintol])

336 # assemble Jacobian J of G

337 J = self.Id - factor * (self.A + self.get_non_linear_Jacobian(u))

339 # solve the linear system

340 if self.direct_solver:

341 delta = spsolve(J, G)

342 else:

343 delta, info = gmres(

344 J,

345 G,

346 x0=zero,

347 M=M,

348 rtol=self.lintol,

349 maxiter=self.liniter,

350 atol=0,

351 callback=self.work_counters['linear'],

352 )

354 if not np.isfinite(delta).all():

355 break

357 # update solution

358 u = u - delta

360 self.work_counters['newton']()

362 return u

364 def u_exact(self, t, u_init=None, t_init=None):

365 r"""

366 Routine to compute the exact solution at time :math:t.

368 Parameters

369 ----------

370 t : float

371 Time of the exact solution.

373 Returns

374 -------

375 me : dtype_u

376 The exact solution.

377 """

379 me = self.dtype_u(self.init, val=0.0)

381 if t > 0:

382 if self.reference_sol_type == 'scipy':

384 def jac(t, u):

385 """

386 Get the Jacobian for the implicit BDF method to use in scipy.solve_ivp

388 Parameters

389 ----------

390 t : float

391 The current time.

392 u : dtype_u

393 Current solution.

395 Returns

396 -------

397 scipy.sparse.csc

398 The derivative of the non-linear part of the solution w.r.t. to the solution.

399 """

400 return self.A + self.get_non_linear_Jacobian(u)

402 def eval_rhs(t, u):

403 """

404 Function to pass to scipy.solve_ivp to evaluate the full right-hand side.

406 Parameters

407 ----------

408 t : float

409 Current time.

410 u : numpy.1darray

411 Current solution.

413 Returns

414 -------

415 numpy.1darray

416 Right-hand side.

417 """

418 return self.eval_f(u.reshape(self.init[0]), t).flatten()

420 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac)

422 elif self.reference_sol_type in ['DIRK', 'SDC']:

423 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI

424 from pySDC.implementations.hooks.log_solution import LogSolution

425 from pySDC.helpers.stats_helper import get_sorted

427 description = {}

428 description['problem_class'] = Quench

429 description['problem_params'] = {

430 'newton_tol': 1e-10,

431 'newton_maxiter': 99,

432 'nvars': 2**10,

433 **self.params,

434 }

436 if self.reference_sol_type == 'DIRK':

437 from pySDC.implementations.sweeper_classes.Runge_Kutta import DIRK43

440 description['sweeper_class'] = DIRK43

441 description['sweeper_params'] = {}

442 description['step_params'] = {'maxiter': 1}

443 description['level_params'] = {'dt': 1e-4}

444 description['convergence_controllers'] = {AdaptivityRK: {'e_tol': 1e-9, 'update_order': 4}}

445 elif self.reference_sol_type == 'SDC':

446 from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit

448 description['sweeper_class'] = generic_implicit

450 description['step_params'] = {'maxiter': 99}

451 description['level_params'] = {'dt': 0.5, 'restol': 1e-10}

453 controller_params = {'hook_class': LogSolution, 'mssdc_jac': False, 'logger_level': 99}

455 controller = controller_nonMPI(

456 description=description, controller_params=controller_params, num_procs=1

457 )

459 uend, stats = controller.run(

460 u0=u_init if u_init is not None else self.u_exact(t=0.0),

461 t0=t_init if t_init is not None else 0,

462 Tend=t,

463 )

465 u_last = get_sorted(stats, type='u', recomputed=False)[-1]

467 if abs(u_last[0] - t) > 1e-2:

468 self.logger.warning(

469 f'Time difference between reference solution and requested time is {abs(u_last[0]-t):.2e}!'

470 )

472 me[:] = u_last[1]

474 return me

477class QuenchIMEX(Quench):

478 """

479 This is a toy problem [1]_ to emulate a magnet that has been cooled to temperatures where superconductivity is possible.

480 However, there is a leak! Some point in the domain is constantly heated and when this has heated up its environment

481 sufficiently, there will be a runaway effect heating up the entire magnet.

482 This effect has actually lead to huge magnets being destroyed at CERN in the past and hence warrants investigation.

484 The model we use is a 1d heat equation with Neumann-zero boundary conditions, meaning this magnet is totally

485 insulated from its environment except for the leak.

486 We add a non-linear term that heats parts of the domain that exceed a certain temperature threshold as well as the

487 leak itself.

489 The problem is discretised with finite difference in space and treated *semi-implicitly*.

490 """

492 dtype_f = imex_mesh

494 def eval_f(self, u, t):

495 """

496 Routine to evaluate the right-hand side of the problem.

498 Parameters

499 ----------

500 u : dtype_u

501 Current values of the numerical solution.

502 t : float

503 Current time of the numerical solution is computed.

505 Returns

506 -------

507 f : dtype_f

508 The right-hand side of the problem.

509 """

511 f = self.dtype_f(self.init)

512 f.impl[:] = self.A.dot(u.flatten()).reshape(self.nvars)

513 f.expl[:] = self.eval_f_non_linear(u, t) + self.b

515 self.work_counters['rhs']()

516 return f

518 def solve_system(self, rhs, factor, u0, t):

519 r"""

520 Simple linear solver for :math:(I - factor \cdot f_{expl})(\vec{u}) = \vec{rhs}.

522 Parameters

523 ----------

524 rhs : dtype_f

525 Right-hand side for the linear system.

526 factor : float

527 Abbrev. for the local stepsize (or any other factor required).

528 u0 : dtype_u

529 Initial guess for the iterative solver.

530 t : float

531 Current time (e.g. for time-dependent BCs).

533 Returns

534 -------

535 me : dtype_u

536 The solution as mesh.

537 """

539 me = self.dtype_u(self.init)

540 me[:] = spsolve(self.Id - factor * self.A, rhs.flatten()).reshape(self.nvars)

541 return me

543 def u_exact(self, t, u_init=None, t_init=None):

544 r"""

545 Routine to compute the exact solution at time :math:t.

547 Parameters

548 ----------

549 t : float

550 Time of the exact solution.

552 Returns

553 -------

554 me : dtype_u

555 The exact solution.

556 """

557 me = self.dtype_u(self.init, val=0.0)

559 if t == 0:

560 me[:] = super().u_exact(t, u_init, t_init)

562 if t > 0:

564 def jac(t, u):

565 """

566 Get the Jacobian for the implicit BDF method to use in scipy.solve_ivp.

568 Parameters

569 ----------

570 t : float

571 Current time.

572 u : dtype_u

573 Current solution.

575 Returns

576 -------

577 scipy.sparse.csc

578 The derivative of the non-linear part of the solution w.r.t. to the solution.

579 """

580 return self.A

582 def eval_rhs(t, u):

583 """

584 Function to pass to scipy.solve_ivp to evaluate the full right-hand side.

586 Parameters

587 ----------

588 t : float

589 Current time

590 u : numpy.1darray

591 Current solution

593 Returns

594 -------

595 numpy.1darray

596 The right-hand side.

597 """

598 f = self.eval_f(u.reshape(self.init[0]), t)

599 return (f.impl + f.expl).flatten()

601 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac)

602 return me