Coverage for pySDC/implementations/problem_classes/Quench.py: 76%
152 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +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 ptype, WorkCounter
8from pySDC.helpers import problem_helper
9from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh
12# noinspection PyUnusedLocal
13class Quench(ptype):
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(),
142 readOnly=True,
143 )
144 self._makeAttributeAndRegister(
145 'newton_tol',
146 'newton_maxiter',
147 'lintol',
148 'liniter',
149 'inexact_linear_ratio',
150 'min_lintol',
151 localVars=locals(),
152 readOnly=False,
153 )
155 self._makeAttributeAndRegister(
156 'newton_tol',
157 'newton_maxiter',
158 'lintol',
159 'liniter',
160 'direct_solver',
161 localVars=locals(),
162 readOnly=False,
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 tol=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
438 from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK
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
449 description['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'quad_type': 'RADAU-RIGHT'}
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