Coverage for pySDC/implementations/problem_classes/Quench.py: 77%
151 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-29 12:50 +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 1D heat equation with a nonlinear heat source, modelling a magnet quench, fully implicit with Newton.
17 This is a toy problem [1]_ to emulate a magnet that has been cooled to temperatures where superconductivity is possible.
18 However, there is a leak! Some point in the domain is constantly heated and when this has heated up its environment
19 sufficiently, there will be a runaway effect heating up the entire magnet.
20 This effect has actually lead to huge magnets being destroyed at CERN in the past and hence warrants investigation.
22 The model we use is a 1d heat equation with Neumann-zero boundary conditions, meaning this magnet is totally
23 insulated from its environment except for the leak.
24 We add a non-linear term that heats parts of the domain that exceed a certain temperature threshold as well as the
25 leak itself.
27 The problem is discretised with finite difference in space and treated *fully-implicitly*.
29 Parameters
30 ----------
31 Cv : float, optional
32 Volumetric heat capacity.
33 K : float, optional
34 Thermal conductivity.
35 u_thresh : float, optional
36 Threshold for temperature.
37 u_max : float, optional
38 Maximum temperature.
39 Q_max : float, optional
40 Maximum heat source power density.
41 leak_range : tuple of float
42 Range of the leak.
43 leak_type : str, optional
44 Type of leak, choose between ``'linear'`` or ``'exponential'``.
45 leak_transition : str, optional
46 Indicates how the heat in the leak propagates, choose between ``'step'`` and ``'Gaussian'``.
47 order : int, optional
48 Order of the finite difference discretization.
49 stencil_type : str, optional
50 Type of stencil for finite differences.
51 bc : str, optional
52 Type of boundary conditions. Default is ``'neumann-zero'``.
53 nvars : int, optional
54 Spatial resolution.
55 newton_tol : float, optional
56 Tolerance for Newton to terminate.
57 newton_maxiter : int, optional
58 Maximum number of Newton iterations to be done.
59 lintol : float, optional
60 Tolerance for linear solver to be done.
61 liniter : int, optional
62 Maximum number of linear iterations inside the Newton solver.
63 direct_solver : bool, optional
64 Indicates if a direct solver should be used.
65 inexact_linear_ratio : float, optional
66 Ratio of tolerance of linear solver to the Newton residual, overrides `lintol`
67 min_lintol : float, optional
68 Minimal tolerance for the linear solver
69 reference_sol_type : str, optional
70 Indicates which method should be used to compute a reference solution.
71 Choose between ``'scipy'``, ``'SDC'``, or ``'DIRK'``.
73 Attributes
74 ----------
75 A : sparse matrix (CSC)
76 FD discretization matrix of the ND grad operator.
77 Id : sparse matrix (CSC)
78 Identity matrix of the same dimension as A.
79 dx : float
80 Distance between two spatial nodes.
81 xv : np.1darray
82 Spatial grid values.
83 leak : np.1darray of bool
84 Indicates the leak.
86 References
87 ----------
88 .. [1] Thermal thin shell approximation towards finite element quench simulation. E. Schnaubelt, M. Wozniak, S. Schöps.
89 Supercond. Sci. Technol. 36 044004. DOI 10.1088/1361-6668/acbeea
90 """
92 dtype_u = mesh
93 dtype_f = mesh
95 def __init__(
96 self,
97 Cv=1000.0,
98 K=1000.0,
99 u_thresh=3e-2,
100 u_max=6e-2,
101 Q_max=1.0,
102 leak_range=(0.45, 0.55),
103 leak_type='linear',
104 leak_transition='step',
105 order=2,
106 stencil_type='center',
107 bc='neumann-zero',
108 nvars=2**7,
109 newton_tol=1e-8,
110 newton_maxiter=99,
111 lintol=1e-8,
112 liniter=99,
113 direct_solver=True,
114 inexact_linear_ratio=None,
115 min_lintol=1e-12,
116 reference_sol_type='scipy',
117 ):
118 """Initialization routine"""
119 # invoke super init, passing number of dofs, dtype_u and dtype_f
120 super().__init__(init=(nvars, None, np.dtype('float64')))
121 self._makeAttributeAndRegister(
122 'Cv',
123 'K',
124 'u_thresh',
125 'u_max',
126 'Q_max',
127 'leak_range',
128 'leak_type',
129 'leak_transition',
130 'order',
131 'stencil_type',
132 'bc',
133 'nvars',
134 'direct_solver',
135 'reference_sol_type',
136 localVars=locals(),
137 readOnly=True,
138 )
139 self._makeAttributeAndRegister(
140 'newton_tol',
141 'newton_maxiter',
142 'lintol',
143 'liniter',
144 'inexact_linear_ratio',
145 'min_lintol',
146 localVars=locals(),
147 readOnly=False,
148 )
150 # setup finite difference discretization from problem helper
151 self.dx, xvalues = problem_helper.get_1d_grid(size=self.nvars, bc=self.bc)
153 self.A, self.b = problem_helper.get_finite_difference_matrix(
154 derivative=2,
155 order=self.order,
156 stencil_type=self.stencil_type,
157 dx=self.dx,
158 size=self.nvars,
159 dim=1,
160 bc=self.bc,
161 )
162 self.A *= self.K / self.Cv
164 self.xv = xvalues
165 self.Id = sp.eye(np.prod(self.nvars), format='csc')
167 self.leak = np.logical_and(self.xv > self.leak_range[0], self.xv < self.leak_range[1])
169 self.work_counters['newton'] = WorkCounter()
170 self.work_counters['rhs'] = WorkCounter()
171 if not self.direct_solver:
172 self.work_counters['linear'] = WorkCounter()
174 def eval_f_non_linear(self, u, t):
175 """
176 Get the non-linear part of f.
178 Parameters
179 ----------
180 u : dtype_u
181 Current values of the numerical solution:
182 t : float
183 Current time at which the numerical solution is computed.
185 Returns
186 -------
187 me : dtype_u
188 The non-linear part of the right-hand side.
189 """
190 u_thresh = self.u_thresh
191 u_max = self.u_max
192 Q_max = self.Q_max
193 me = self.dtype_u(self.init)
195 if self.leak_type == 'linear':
196 me[:] = (u - u_thresh) / (u_max - u_thresh) * Q_max
197 elif self.leak_type == 'exponential':
198 me[:] = Q_max * (np.exp(u) - np.exp(u_thresh)) / (np.exp(u_max) - np.exp(u_thresh))
199 else:
200 raise NotImplementedError(f'Leak type \"{self.leak_type}\" not implemented!')
202 me[u < u_thresh] = 0
203 if self.leak_transition == 'step':
204 me[self.leak] = Q_max
205 elif self.leak_transition == 'Gaussian':
206 me[:] = np.max([me, Q_max * np.exp(-((self.xv - 0.5) ** 2) / 3e-2)], axis=0)
207 else:
208 raise NotImplementedError(f'Leak transition \"{self.leak_transition}\" not implemented!')
210 me[u >= u_max] = Q_max
212 me[:] /= self.Cv
214 return me
216 def eval_f(self, u, t):
217 """
218 Evaluate the full right-hand side.
220 Parameters
221 ----------
222 u : dtype_u
223 Current values of the numerical solution.
224 t : float
225 Current time at which the numerical solution is computed.
227 Returns
228 -------
229 f : dtype_f
230 The right-hand side of the problem.
231 """
232 f = self.dtype_f(self.init)
233 f[:] = self.A.dot(u.flatten()).reshape(self.nvars) + self.b + self.eval_f_non_linear(u, t)
235 self.work_counters['rhs']()
236 return f
238 def get_non_linear_Jacobian(self, u):
239 """
240 Evaluate the non-linear part of the Jacobian only.
242 Parameters
243 ----------
244 u : dtype_u
245 Current values of the numerical solution.
247 Returns
248 -------
249 scipy.sparse.csc
250 The derivative of the non-linear part of the solution w.r.t. to the solution.
251 """
252 u_thresh = self.u_thresh
253 u_max = self.u_max
254 Q_max = self.Q_max
255 me = self.dtype_u(self.init)
257 if self.leak_type == 'linear':
258 me[:] = Q_max / (u_max - u_thresh)
259 elif self.leak_type == 'exponential':
260 me[:] = Q_max * np.exp(u) / (np.exp(u_max) - np.exp(u_thresh))
261 else:
262 raise NotImplementedError(f'Leak type {self.leak_type} not implemented!')
264 me[u < u_thresh] = 0
265 if self.leak_transition == 'step':
266 me[self.leak] = 0
267 elif self.leak_transition == 'Gaussian':
268 me[self.leak] = 0
269 me[self.leak][u[self.leak] > Q_max * np.exp(-((self.xv[self.leak] - 0.5) ** 2) / 3e-2)] = 1
270 else:
271 raise NotImplementedError(f'Leak transition \"{self.leak_transition}\" not implemented!')
272 me[u > u_max] = 0
274 me[:] /= self.Cv
276 return sp.diags(me, format='csc')
278 def solve_system(self, rhs, factor, u0, t):
279 r"""
280 Simple Newton solver for :math:`(I - factor \cdot f)(\vec{u}) = \vec{rhs}`.
282 Parameters
283 ----------
284 rhs : dtype_f
285 Right-hand side.
286 factor : float
287 Abbrev. for the local stepsize (or any other factor required).
288 u0 : dtype_u
289 Initial guess for the iterative solver.
290 t : float
291 Current time (e.g. for time-dependent BCs).
293 Returns
294 -------
295 u : dtype_u
296 The solution as mesh.
297 """
298 u = self.dtype_u(u0)
299 res = np.inf
300 delta = self.dtype_u(self.init, val=0.0)
302 # construct a preconditioner for the space solver
303 if not self.direct_solver:
304 M = inv((self.Id - factor * self.A).toarray())
305 zero = self.dtype_u(self.init, val=0.0)
307 for n in range(0, self.newton_maxiter):
308 # assemble G such that G(u) = 0 at the solution of the step
309 G = u - factor * self.eval_f(u, t) - rhs
310 self.work_counters[
311 'rhs'
312 ].decrement() # Work regarding construction of the Jacobian etc. should count into the Newton iterations only
314 res = np.linalg.norm(G, np.inf)
315 if res <= self.newton_tol and n > 0: # we want to make at least one Newton iteration
316 break
318 if self.inexact_linear_ratio:
319 self.lintol = max([res * self.inexact_linear_ratio, self.min_lintol])
321 # assemble Jacobian J of G
322 J = self.Id - factor * (self.A + self.get_non_linear_Jacobian(u))
324 # solve the linear system
325 if self.direct_solver:
326 delta = spsolve(J, G)
327 else:
328 delta, info = gmres(
329 J,
330 G,
331 x0=zero,
332 M=M,
333 rtol=self.lintol,
334 maxiter=self.liniter,
335 atol=0,
336 callback=self.work_counters['linear'],
337 )
339 if not np.isfinite(delta).all():
340 break
342 # update solution
343 u = u - delta
345 self.work_counters['newton']()
347 return u
349 def u_exact(self, t, u_init=None, t_init=None):
350 r"""
351 Routine to compute the exact solution at time :math:`t`.
353 Parameters
354 ----------
355 t : float
356 Time of the exact solution.
357 u_init : dtype_u, optional
358 Initial conditions for getting the exact solution.
359 t_init : float, optional
360 The starting time.
362 Returns
363 -------
364 me : dtype_u
365 The exact solution.
366 """
368 me = self.dtype_u(self.init, val=0.0)
370 if t > 0:
371 if self.reference_sol_type == 'scipy':
373 def jac(t, u):
374 """
375 Get the Jacobian for the implicit BDF method to use in `scipy.solve_ivp`
377 Parameters
378 ----------
379 t : float
380 The current time.
381 u : dtype_u
382 Current solution.
384 Returns
385 -------
386 scipy.sparse.csc
387 The derivative of the non-linear part of the solution w.r.t. to the solution.
388 """
389 return self.A + self.get_non_linear_Jacobian(u)
391 def eval_rhs(t, u):
392 """
393 Function to pass to `scipy.solve_ivp` to evaluate the full right-hand side.
395 Parameters
396 ----------
397 t : float
398 Current time.
399 u : numpy.1darray
400 Current solution.
402 Returns
403 -------
404 numpy.1darray
405 Right-hand side.
406 """
407 return self.eval_f(u.reshape(self.init[0]), t).flatten()
409 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac)
411 elif self.reference_sol_type in ['DIRK', 'SDC']:
412 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
413 from pySDC.implementations.hooks.log_solution import LogSolution
414 from pySDC.helpers.stats_helper import get_sorted
416 description = {}
417 description['problem_class'] = Quench
418 description['problem_params'] = {
419 'newton_tol': 1e-10,
420 'newton_maxiter': 99,
421 'nvars': 2**10,
422 **self.params,
423 }
425 if self.reference_sol_type == 'DIRK':
426 from pySDC.implementations.sweeper_classes.Runge_Kutta import DIRK43
427 from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK
429 description['sweeper_class'] = DIRK43
430 description['sweeper_params'] = {}
431 description['step_params'] = {'maxiter': 1}
432 description['level_params'] = {'dt': 1e-4}
433 description['convergence_controllers'] = {AdaptivityRK: {'e_tol': 1e-9, 'update_order': 4}}
434 elif self.reference_sol_type == 'SDC':
435 from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
437 description['sweeper_class'] = generic_implicit
438 description['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'quad_type': 'RADAU-RIGHT'}
439 description['step_params'] = {'maxiter': 99}
440 description['level_params'] = {'dt': 0.5, 'restol': 1e-10}
442 controller_params = {'hook_class': LogSolution, 'mssdc_jac': False, 'logger_level': 99}
444 controller = controller_nonMPI(
445 description=description, controller_params=controller_params, num_procs=1
446 )
448 uend, stats = controller.run(
449 u0=u_init if u_init is not None else self.u_exact(t=0.0),
450 t0=t_init if t_init is not None else 0,
451 Tend=t,
452 )
454 u_last = get_sorted(stats, type='u', recomputed=False)[-1]
456 if abs(u_last[0] - t) > 1e-2:
457 self.logger.warning(
458 f'Time difference between reference solution and requested time is {abs(u_last[0]-t):.2e}!'
459 )
461 me[:] = u_last[1]
463 return me
466class QuenchIMEX(Quench):
467 """
468 1D heat equation with a nonlinear heat source, modelling a magnet quench, IMEX with diffusion implicit.
470 This is a toy problem [1]_ to emulate a magnet that has been cooled to temperatures where superconductivity is possible.
471 However, there is a leak! Some point in the domain is constantly heated and when this has heated up its environment
472 sufficiently, there will be a runaway effect heating up the entire magnet.
473 This effect has actually lead to huge magnets being destroyed at CERN in the past and hence warrants investigation.
475 The model we use is a 1d heat equation with Neumann-zero boundary conditions, meaning this magnet is totally
476 insulated from its environment except for the leak.
477 We add a non-linear term that heats parts of the domain that exceed a certain temperature threshold as well as the
478 leak itself.
480 The problem is discretised with finite difference in space and treated *semi-implicitly*.
481 """
483 dtype_f = imex_mesh
485 def eval_f(self, u, t):
486 """
487 Routine to evaluate the right-hand side of the problem.
489 Parameters
490 ----------
491 u : dtype_u
492 Current values of the numerical solution.
493 t : float
494 Current time of the numerical solution is computed.
496 Returns
497 -------
498 f : dtype_f
499 The right-hand side of the problem.
500 """
502 f = self.dtype_f(self.init)
503 f.impl[:] = self.A.dot(u.flatten()).reshape(self.nvars)
504 f.expl[:] = self.eval_f_non_linear(u, t) + self.b
506 self.work_counters['rhs']()
507 return f
509 def solve_system(self, rhs, factor, u0, t):
510 r"""
511 Simple linear solver for :math:`(I - factor \cdot f_{expl})(\vec{u}) = \vec{rhs}`.
513 Parameters
514 ----------
515 rhs : dtype_f
516 Right-hand side for the linear system.
517 factor : float
518 Abbrev. for the local stepsize (or any other factor required).
519 u0 : dtype_u
520 Initial guess for the iterative solver.
521 t : float
522 Current time (e.g. for time-dependent BCs).
524 Returns
525 -------
526 me : dtype_u
527 The solution as mesh.
528 """
530 me = self.dtype_u(self.init)
531 me[:] = spsolve(self.Id - factor * self.A, rhs.flatten()).reshape(self.nvars)
532 return me
534 def u_exact(self, t, u_init=None, t_init=None):
535 r"""
536 Routine to compute the exact solution at time :math:`t`.
538 Parameters
539 ----------
540 t : float
541 Time of the exact solution.
542 u_init : dtype_u, optional
543 Initial conditions for getting the exact solution.
544 t_init : float, optional
545 The starting time.
547 Returns
548 -------
549 me : dtype_u
550 The exact solution.
551 """
552 me = self.dtype_u(self.init, val=0.0)
554 if t == 0:
555 me[:] = super().u_exact(t, u_init, t_init)
557 if t > 0:
559 def jac(t, u):
560 """
561 Get the Jacobian for the implicit BDF method to use in `scipy.solve_ivp`.
563 Parameters
564 ----------
565 t : float
566 Current time.
567 u : dtype_u
568 Current solution.
570 Returns
571 -------
572 scipy.sparse.csc
573 The derivative of the non-linear part of the solution w.r.t. to the solution.
574 """
575 return self.A
577 def eval_rhs(t, u):
578 """
579 Function to pass to `scipy.solve_ivp` to evaluate the full right-hand side.
581 Parameters
582 ----------
583 t : float
584 Current time
585 u : numpy.1darray
586 Current solution
588 Returns
589 -------
590 numpy.1darray
591 The right-hand side.
592 """
593 f = self.eval_f(u.reshape(self.init[0]), t)
594 return (f.impl + f.expl).flatten()
596 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac)
597 return me