Merge pull request #283 from brownbaerchen/getdp

Reference solution for Quench problem

Author: pancetta

pySDC/core/Problem.py

@@ -70,13 +70,17 @@ def apply_mass_matrix(self, u):
         """
         raise NotImplementedError('ERROR: if you want a mass matrix, implement apply_mass_matrix(u)')
 
-    def generate_scipy_reference_solution(self, eval_rhs, t, u_init=None, t_init=None):
+    def generate_scipy_reference_solution(self, eval_rhs, t, u_init=None, t_init=None, **kwargs):
         """
         Compute a reference solution using `scipy.solve_ivp` with very small tolerances.
         Keep in mind that scipy needs the solution to be a one dimensional array. If you are solving something higher
         dimensional, you need to make sure the function `eval_rhs` takes a flattened one-dimensional version as an input
         and output, but reshapes to whatever the problem needs for evaluation.
 
+        The keyword arguments will be passed to `scipy.solve_ivp`. You should consider passing `method='BDF'` for stiff
+        problems and to accelerate that you can pass a function that evaluates the Jacobian with arguments `jac(t, u)`
+        as `jac=jac`.
+
         Args:
             eval_rhs (function): Function evaluate the full right hand side. Must have signature `eval_rhs(float: t, numpy.1darray: u)`
             t (float): current time
@@ -94,4 +98,6 @@ def generate_scipy_reference_solution(self, eval_rhs, t, u_init=None, t_init=Non
         t_init = 0 if t_init is None else t_init
 
         u_shape = u_init.shape
-        return solve_ivp(eval_rhs, (t_init, t), u_init.flatten(), rtol=tol, atol=tol).y[:, -1].reshape(u_shape)
+        return (
+            solve_ivp(eval_rhs, (t_init, t), u_init.flatten(), rtol=tol, atol=tol, **kwargs).y[:, -1].reshape(u_shape)
+        )

pySDC/implementations/problem_classes/LeakySuperconductor.py

@@ -40,6 +40,7 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=mesh):
             'u_max': 2e-2,
             'Q_max': 1.0,
             'leak_range': (0.45, 0.55),
+            'leak_type': 'linear',
             'order': 2,
             'stencil_type': 'center',
             'bc': 'neumann-zero',
@@ -116,7 +117,13 @@ def eval_f_non_linear(self, u, t):
         Q_max = self.params.Q_max
         me = self.dtype_u(self.init)
 
-        me[:] = (u - u_thresh) / (u_max - u_thresh) * Q_max
+        if self.params.leak_type == 'linear':
+            me[:] = (u - u_thresh) / (u_max - u_thresh) * Q_max
+        elif self.params.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.params.leak_type} not implemented!')
+
         me[u < u_thresh] = 0
         me[self.leak] = Q_max
         me[u >= u_max] = Q_max
@@ -145,47 +152,53 @@ def eval_f(self, u, t):
         self.work_counters['rhs']()
         return f
 
-    def solve_system(self, rhs, factor, u0, t):
+    def get_non_linear_Jacobian(self, u):
         """
-        Simple Newton solver for (I-factor*f)(u) = rhs
+        Evaluate the non-linear part of the Jacobian only
 
         Args:
-            rhs (dtype_f): right-hand side
-            factor (float): abbrev. for the local stepsize (or any other factor required)
-            u0 (dtype_u): initial guess for the iterative solver
-            t (float): current time (e.g. for time-dependent BCs)
+            u (dtype_u): Current solution
 
         Returns:
-            dtype_u: solution as mesh
+            scipy.sparse.csc: The derivative of the non-linear part of the solution w.r.t. to the solution.
         """
+        u_thresh = self.params.u_thresh
+        u_max = self.params.u_max
+        Q_max = self.params.Q_max
+        me = self.dtype_u(self.init)
+
+        if self.params.leak_type == 'linear':
+            me[:] = Q_max / (u_max - u_thresh)
+        elif self.params.leak_type == 'exponential':
+            me[:] = Q_max * np.exp(u) / (np.exp(u_max) - np.exp(u_thresh))
+        else:
+            raise NotImplementedError(f'Leak type {self.params.leak_type} not implemented!')
 
-        def get_non_linear_Jacobian(u):
-            """
-            Evaluate the non-linear part of the Jacobian only
+        me[u < u_thresh] = 0
+        me[u > u_max] = 0
+        me[self.leak] = 0
 
-            Args:
-                u (dtype_u): Current solution
+        # boundary conditions
+        me[0] = 0.0
+        me[-1] = 0.0
 
-            Returns:
-                scipy.sparse.csc: The derivative of the non-linear part of the solution w.r.t. to the solution.
-            """
-            u_thresh = self.params.u_thresh
-            u_max = self.params.u_max
-            Q_max = self.params.Q_max
-            me = self.dtype_u(self.init)
+        me[:] /= self.params.Cv
 
-            me[:] = Q_max / (u_max - u_thresh)
-            me[u < u_thresh] = 0
-            me[u > u_max] = 0
-            me[self.leak] = 0
+        return sp.diags(me, format='csc')
 
-            # boundary conditions
-            me[0] = 0.0
-            me[-1] = 0.0
+    def solve_system(self, rhs, factor, u0, t):
+        """
+        Simple Newton solver for (I-factor*f)(u) = rhs
 
-            me[:] /= self.params.Cv
+        Args:
+            rhs (dtype_f): right-hand side
+            factor (float): abbrev. for the local stepsize (or any other factor required)
+            u0 (dtype_u): initial guess for the iterative solver
+            t (float): current time (e.g. for time-dependent BCs)
 
-            return sp.diags(me, format='csc')
+        Returns:
+            dtype_u: solution as mesh
+        """
 
         u = self.dtype_u(u0)
         res = np.inf
@@ -209,7 +222,7 @@ def get_non_linear_Jacobian(u):
                 break
 
             # assemble Jacobian J of G
-            J = self.Id - factor * (self.A + get_non_linear_Jacobian(u))
+            J = self.Id - factor * (self.A + self.get_non_linear_Jacobian(u))
 
             # solve the linear system
             if self.params.direct_solver:
@@ -251,6 +264,19 @@ def u_exact(self, t, u_init=None, t_init=None):
 
         if t > 0:
 
+            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)
+
             def eval_rhs(t, u):
                 """
                 Function to pass to `scipy.solve_ivp` to evaluate the full RHS
@@ -264,7 +290,7 @@ def eval_rhs(t, u):
                 """
                 return self.eval_f(u.reshape(self.init[0]), t).flatten()
 
-            me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
+            me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac)
         return me
 
 
@@ -326,6 +352,19 @@ def u_exact(self, t, u_init=None, t_init=None):
 
         if t > 0:
 
+            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
+
             def eval_rhs(t, u):
                 """
                 Function to pass to `scipy.solve_ivp` to evaluate the full RHS
@@ -340,5 +379,5 @@ def eval_rhs(t, u):
                 f = self.eval_f(u.reshape(self.init[0]), t)
                 return (f.impl + f.expl).flatten()
 
-            me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
+            me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac)
         return me

pySDC/projects/Resilience/leaky_superconductor.py

@@ -194,7 +194,7 @@ def plot_solution(stats, controller):  # pragma: no cover
     dt_ax.set_ylabel(r'$\Delta t$')
 
 
-def compare_imex_full(plotting=False):
+def compare_imex_full(plotting=False, leak_type='linear'):
     """
     Compare the results of IMEX and fully implicit runs. For IMEX we need to limit the step size in order to achieve convergence, but for fully implicit, adaptivity can handle itself better.
 
@@ -203,19 +203,22 @@ def compare_imex_full(plotting=False):
     """
     from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
     from pySDC.implementations.hooks.log_work import LogWork
+    from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun
 
     maxiter = 5
     num_nodes = 3
     newton_iter_max = 20
 
     res = {}
     rhs = {}
+    error = {}
 
     custom_description = {}
     custom_description['problem_params'] = {
         'newton_tol': 1e-10,
         'newton_iter': newton_iter_max,
         'nvars': 2**9,
+        'leak_type': leak_type,
     }
     custom_description['step_params'] = {'maxiter': maxiter}
     custom_description['sweeper_params'] = {'num_nodes': num_nodes}
@@ -227,13 +230,14 @@ def compare_imex_full(plotting=False):
             custom_description=custom_description,
             custom_controller_params=custom_controller_params,
             imex=imex,
-            Tend=5e2,
-            hook_class=LogWork,
+            Tend=4.3e2,
+            hook_class=[LogWork, LogGlobalErrorPostRun],
         )
 
         res[imex] = get_sorted(stats, type='u')[-1][1]
         newton_iter = [me[1] for me in get_sorted(stats, type='work_newton')]
         rhs[imex] = np.mean([me[1] for me in get_sorted(stats, type='work_rhs')]) // 1
+        error[imex] = get_sorted(stats, type='e_global_post_run')[-1][1]
 
         if imex:
             assert all([me == 0 for me in newton_iter]), "IMEX is not supposed to do Newton iterations!"
@@ -258,6 +262,11 @@ def compare_imex_full(plotting=False):
         rhs[True] == rhs[False]
     ), f"Expected IMEX and fully implicit schemes to take the same number of right hand side evaluations per step, but got {rhs[True]} and {rhs[False]}!"
 
+    assert (
+        error[True] > error[False]
+    ), f"Expected IMEX to be less accurate at the same precision settings than unsplit version, got for IMEX: e={error[True]:.2e} and fully implicit: e={error[False]:.2e}"
+    assert error[True] < 1.1e-4, f'Expected error of IMEX version to be less than 1.1e-4, but got e={error[True]:.2e}!'
+
 
 if __name__ == '__main__':
     compare_imex_full(plotting=True)

pySDC/tests/test_convergence_controllers/test_InterpolateBetweenRestarts.py

@@ -155,6 +155,7 @@ def run_vdp(hook, adaptivity=True):
         'newton_tol': 1e-9,
         'newton_maxiter': 99,
         'u0': np.array([2.0, 0.0]),
+        'crash_at_maxiter': False,
     }
 
     # initialize step parameters

pySDC/tests/test_projects/test_resilience/test_leaky_superconductor.py

@@ -2,10 +2,11 @@
 
 
 @pytest.mark.base
-def test_imex_vs_fully_implicit_leaky_superconductor():
+@pytest.mark.parametrize('leak_type', ['linear', 'exponential'])
+def test_imex_vs_fully_implicit_leaky_superconductor(leak_type):
     """
     Test if the IMEX and fully implicit schemes get the same solution and that the runaway process has started.
     """
     from pySDC.projects.Resilience.leaky_superconductor import compare_imex_full
 
-    compare_imex_full()
+    compare_imex_full(plotting=False, leak_type=leak_type)