Switched the method for the scipy reference solution to the implicit BDF method

Author: Thomas Baumann

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

@@ -145,47 +145,47 @@ 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)
 
-        def get_non_linear_Jacobian(u):
-            """
-            Evaluate the non-linear part of the Jacobian only
+        me[:] = Q_max / (u_max - u_thresh)
+        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 +209,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 +251,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 +277,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 +339,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 +366,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

@@ -203,13 +203,15 @@ 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'] = {
@@ -228,12 +230,13 @@ def compare_imex_full(plotting=False):
             custom_controller_params=custom_controller_params,
             imex=imex,
             Tend=5e2,
-            hook_class=LogWork,
+            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 +261,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)