Quench resilience

pySDC/implementations/problem_classes/LeakySuperconductor.py → pySDC/implementations/problem_classes/Quench.py

@@ -9,7 +9,7 @@
 
 
 # noinspection PyUnusedLocal
-class LeakySuperconductor(ptype):
+class Quench(ptype):
     """
     This is a toy problem to emulate a magnet that has been cooled to temperatures where superconductivity is possible.
     However, there is a leak! Some point in the domain is constantly heated and when this has heated up its environment
@@ -34,6 +34,7 @@ def __init__(
         Q_max=1.0,
         leak_range=(0.45, 0.55),
         leak_type='linear',
+        leak_transition='step',
         order=2,
         stencil_type='center',
         bc='neumann-zero',
@@ -43,6 +44,7 @@ def __init__(
         lintol=1e-8,
         liniter=99,
         direct_solver=True,
+        reference_sol_type='scipy',
     ):
         """
         Initialization routine
@@ -62,6 +64,7 @@ def __init__(
             'Q_max',
             'leak_range',
             'leak_type',
+            'leak_transition',
             'order',
             'stencil_type',
             'bc',
@@ -71,6 +74,7 @@ def __init__(
             'lintol',
             'liniter',
             'direct_solver',
+            'reference_sol_type',
             localVars=locals(),
             readOnly=True,
         )
@@ -130,15 +134,17 @@ def eval_f_non_linear(self, u, t):
         elif self.leak_type == 'exponential':
             me[:] = Q_max * (np.exp(u) - np.exp(u_thresh)) / (np.exp(u_max) - np.exp(u_thresh))
         else:
-            raise NotImplementedError(f'Leak type {self.leak_type} not implemented!')
+            raise NotImplementedError(f'Leak type \"{self.leak_type}\" not implemented!')
 
         me[u < u_thresh] = 0
-        me[self.leak] = Q_max
-        me[u >= u_max] = Q_max
+        if self.leak_transition == 'step':
+            me[self.leak] = Q_max
+        elif self.leak_transition == 'Gaussian':
+            me[:] = np.max([me, Q_max * np.exp(-((self.xv - 0.5) ** 2) / 3e-2)], axis=0)
+        else:
+            raise NotImplementedError(f'Leak transition \"{self.leak_transition}\" not implemented!')
 
-        # boundary conditions
-        me[0] = 0.0
-        me[-1] = 0.0
+        me[u >= u_max] = Q_max
 
         me[:] /= self.Cv
 
@@ -183,12 +189,14 @@ def get_non_linear_Jacobian(self, u):
             raise NotImplementedError(f'Leak type {self.leak_type} not implemented!')
 
         me[u < u_thresh] = 0
+        if self.leak_transition == 'step':
+            me[self.leak] = 0
+        elif self.leak_transition == 'Gaussian':
+            me[self.leak] = 0
+            me[self.leak][u[self.leak] > Q_max * np.exp(-((self.xv[self.leak] - 0.5) ** 2) / 3e-2)] = 1
+        else:
+            raise NotImplementedError(f'Leak transition \"{self.leak_transition}\" not implemented!')
         me[u > u_max] = 0
-        me[self.leak] = 0
-
-        # boundary conditions
-        me[0] = 0.0
-        me[-1] = 0.0
 
         me[:] /= self.Cv
 
@@ -270,38 +278,92 @@ def u_exact(self, t, u_init=None, t_init=None):
         me = self.dtype_u(self.init, val=0.0)
 
         if t > 0:
+            if self.reference_sol_type == 'scipy':
+
+                def jac(t, u):
+                    """
+                    Get the Jacobian for the implicit BDF method to use in `scipy.solve_ivp`
+
+                    Args:
+                        t (float): The current time
+                        u (dtype_u): Current solution
+
+                    Returns:
+                        scipy.sparse.csc: The derivative of the non-linear part of the solution w.r.t. to the solution.
+                    """
+                    return self.A + self.get_non_linear_Jacobian(u)
+
+                def eval_rhs(t, u):
+                    """
+                    Function to pass to `scipy.solve_ivp` to evaluate the full RHS
+
+                    Args:
+                        t (float): Current time
+                        u (numpy.1darray): Current solution
+
+                    Returns:
+                        (numpy.1darray): RHS
+                    """
+                    return self.eval_f(u.reshape(self.init[0]), t).flatten()
+
+                me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac)
+
+            elif self.reference_sol_type in ['DIRK', 'SDC']:
+                from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+                from pySDC.implementations.hooks.log_solution import LogSolution
+                from pySDC.helpers.stats_helper import get_sorted
+
+                description = {}
+                description['problem_class'] = Quench
+                description['problem_params'] = {
+                    'newton_tol': 1e-10,
+                    'newton_iter': 99,
+                    'nvars': 2**10,
+                    **self.params,
+                }
+
+                if self.reference_sol_type == 'DIRK':
+                    from pySDC.implementations.sweeper_classes.Runge_Kutta import DIRK34
+                    from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK
+
+                    description['sweeper_class'] = DIRK34
+                    description['sweeper_params'] = {}
+                    description['step_params'] = {'maxiter': 1}
+                    description['level_params'] = {'dt': 1e-4}
+                    description['convergence_controllers'] = {AdaptivityRK: {'e_tol': 1e-9, 'update_order': 4}}
+                elif self.reference_sol_type == 'SDC':
+                    from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+
+                    description['sweeper_class'] = generic_implicit
+                    description['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'quad_type': 'RADAU-RIGHT'}
+                    description['step_params'] = {'maxiter': 99}
+                    description['level_params'] = {'dt': 0.5, 'restol': 1e-10}
+
+                controller_params = {'hook_class': LogSolution, 'mssdc_jac': False, 'logger_level': 99}
+
+                controller = controller_nonMPI(
+                    description=description, controller_params=controller_params, num_procs=1
+                )
 
-            def jac(t, u):
-                """
-                Get the Jacobian for the implicit BDF method to use in `scipy.odeint`
-
-                Args:
-                    t (float): The current time
-                    u (dtype_u): Current solution
-
-                Returns:
-                    scipy.sparse.csc: The derivative of the non-linear part of the solution w.r.t. to the solution.
-                """
-                return self.A + self.get_non_linear_Jacobian(u)
+                uend, stats = controller.run(
+                    u0=u_init if u_init is not None else self.u_exact(t=0.0),
+                    t0=t_init if t_init is not None else 0,
+                    Tend=t,
+                )
 
-            def eval_rhs(t, u):
-                """
-                Function to pass to `scipy.solve_ivp` to evaluate the full RHS
+                u_last = get_sorted(stats, type='u', recomputed=False)[-1]
 
-                Args:
-                    t (float): Current time
-                    u (numpy.1darray): Current solution
+                if abs(u_last[0] - t) > 1e-2:
+                    self.logger.warning(
+                        f'Time difference between reference solution and requested time is {abs(u_last[0]-t):.2e}!'
+                    )
 
-                Returns:
-                    (numpy.1darray): RHS
-                """
-                return self.eval_f(u.reshape(self.init[0]), t).flatten()
+                me[:] = u_last[1]
 
-            me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac)
         return me
 
 
-class LeakySuperconductorIMEX(LeakySuperconductor):
+class QuenchIMEX(Quench):
     dtype_f = imex_mesh
 
     def eval_f(self, u, t):
@@ -360,8 +422,7 @@ def u_exact(self, t, u_init=None, t_init=None):
 
             def jac(t, u):
                 """
-                Get the Jacobian for the implicit BDF method to use in `scipy.odeint`
-
+                Get the Jacobian for the implicit BDF method to use in `scipy.solve_ivp`
                 Args:
                     t (float): The current time
                     u (dtype_u): Current solution

pySDC/projects/Resilience/Lorenz.py

@@ -9,6 +9,7 @@
 from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
 from pySDC.core.Errors import ConvergenceError
 from pySDC.projects.Resilience.hook import LogData, hook_collection
+from pySDC.projects.Resilience.strategies import merge_descriptions
 
 
 def run_Lorenz(
@@ -18,7 +19,6 @@ def run_Lorenz(
     hook_class=LogData,
     fault_stuff=None,
     custom_controller_params=None,
-    custom_problem_params=None,
     use_MPI=False,
     **kwargs,
 ):
@@ -32,7 +32,6 @@ def run_Lorenz(
         hook_class (pySDC.Hook): A hook to store data
         fault_stuff (dict): A dictionary with information on how to add faults
         custom_controller_params (dict): Overwrite presets
-        custom_problem_params (dict): Overwrite presets
         use_MPI (bool): Whether or not to use MPI
 
     Returns:
@@ -56,9 +55,6 @@ def run_Lorenz(
         'newton_maxiter': 99,
     }
 
-    if custom_problem_params is not None:
-        problem_params = {**problem_params, **custom_problem_params}
-
     # initialize step parameters
     step_params = dict()
     step_params['maxiter'] = 4
@@ -82,11 +78,7 @@ def run_Lorenz(
     description['step_params'] = step_params
 
     if custom_description is not None:
-        for k in custom_description.keys():
-            if k == 'sweeper_class':
-                description[k] = custom_description[k]
-                continue
-            description[k] = {**description.get(k, {}), **custom_description.get(k, {})}
+        description = merge_descriptions(description, custom_description)
 
     # set time parameters
     t0 = 0.0
@@ -98,18 +90,13 @@ def run_Lorenz(
 
         comm = kwargs.get('comm', MPI.COMM_WORLD)
         controller = controller_MPI(controller_params=controller_params, description=description, comm=comm)
-
-        # get initial values on finest level
         P = controller.S.levels[0].prob
-        uinit = P.u_exact(t0)
     else:
         controller = controller_nonMPI(
             num_procs=num_procs, controller_params=controller_params, description=description
         )
-
-        # get initial values on finest level
         P = controller.MS[0].levels[0].prob
-        uinit = P.u_exact(t0)
+    uinit = P.u_exact(t0)
 
     # insert faults
     if fault_stuff is not None: