Merge pull request #278 from brownbaerchen/resilience_quench

Bugfixes in Quench problem

Author: pancetta

pySDC/helpers/plot_helper.py

@@ -6,6 +6,17 @@
 
 
 def figsize(textwidth, scale, ratio):
+    """
+    Get figsize.
+
+    Args:
+        textwidth (str): Textwdith in your LaTeX file in points
+        scale (float): The width of the figure relative to the textwidth
+        ratio (float): The height of the figure relative to its width
+
+    Returns:
+        list: Width and height of the figure to be passed to matplotlib
+    """
     fig_width_pt = textwidth  # Get this from LaTeX using \the\textwidth
     inches_per_pt = 1.0 / 72.27  # Convert pt to inch
     fig_width = fig_width_pt * inches_per_pt * scale  # width in inches
@@ -14,6 +25,43 @@ def figsize(textwidth, scale, ratio):
     return fig_size
 
 
+def figsize_by_journal(journal, scale, ratio):  # pragma no cover
+    """
+    Get figsize for specific journal. If you supply a text height, we will rescale the figure to fit on the page instead
+    of the parameters supplied.
+
+    Args:
+        journal (str): Name of journal
+        scale (float): The width of the figure relative to the textwidth
+        ratio (float): The height of the figure relative to its width
+
+    Returns:
+        list: Width and height of the figure to be passed to matplotlib
+    """
+    # store text width in points here, get this from LaTeX using \the\textwidth
+    textwidths = {
+        'JSC_beamer': 426.79135,
+        'Springer_Numerical_Algorithms': 338.58778,
+    }
+    # store text height in points here, get this from LaTeX using \the\textheight
+    textheights = {
+        'JSC_beamer': 214.43411,
+    }
+    assert (
+        journal in textwidths.keys()
+    ), f"Textwidth only available for {list(textwidths.keys())}. Please implement one for \"{journal}\"! Get the textwidth using \"\\the\\textwidth\" in your tex file."
+
+    # see if the figure fits on the page or if we need to apply the scaling to the height instead
+    if scale * ratio * textwidths[journal] > textheights.get(journal, 1e9):
+        if textheights[journal] / scale / ratio > textwidths[journal]:
+            raise ValueError(
+                f"We cannot fit figure with scale {scale:.2f} and ratio {ratio:.2f} on the page for journal {journal}!"
+            )
+        return figsize(textheights[journal] / (scale * ratio), 1, ratio)
+
+    return figsize(textwidths[journal], scale, ratio)
+
+
 def setup_mpl(font_size=8, reset=False):
     if reset:
         mpl.rcParams.update(default_mpl_params)

pySDC/implementations/convergence_controller_classes/check_convergence.py

@@ -24,7 +24,35 @@ def setup(self, controller, params, description, **kwargs):
         Returns:
             (dict): The updated params dictionary
         """
-        return {"control_order": +200, **super().setup(controller, params, description, **kwargs)}
+        defaults = {'control_order': +200, 'use_e_tol': 'e_tol' in description['level_params'].keys()}
+
+        return {**defaults, **super().setup(controller, params, description, **kwargs)}
+
+    def dependencies(self, controller, description, **kwargs):
+        """
+        Load the embedded error estimator if needed.
+
+        Args:
+            controller (pySDC.Controller): The controller
+            description (dict): The description object used to instantiate the controller
+
+        Returns:
+            None
+        """
+        super().dependencies(controller, description)
+
+        if self.params.use_e_tol:
+            from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import (
+                EstimateEmbeddedError,
+            )
+            from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+
+            controller.add_convergence_controller(
+                EstimateEmbeddedError.get_implementation("nonMPI" if type(controller) == controller_nonMPI else "MPI"),
+                description=description,
+            )
+
+        return None
 
     @staticmethod
     def check_convergence(S):
@@ -42,10 +70,16 @@ def check_convergence(S):
         L = S.levels[0]
         L.sweep.compute_residual()
 
-        # get residual and check against prescribed tolerance (plus check number of iterations
-        res = L.status.residual
+        # get residual and check against prescribed tolerance (plus check number of iterations)
+        iter_converged = S.status.iter >= S.params.maxiter
+        res_converged = L.status.residual <= L.params.restol
+        e_tol_converged = (
+            L.status.error_embedded_estimate < L.params.e_tol
+            if (L.params.get('e_tol') and L.status.get('error_embedded_estimate'))
+            else False
+        )
         converged = (
-            S.status.iter >= S.params.maxiter or res <= L.params.restol or S.status.force_done
+            iter_converged or res_converged or e_tol_converged or S.status.force_done
         ) and not S.status.force_continue
         if converged is None:
             converged = False

pySDC/implementations/hooks/log_embedded_error_estimate.py

@@ -6,7 +6,7 @@ class LogEmbeddedErrorEstimate(hooks):
     Store the embedded error estimate at the end of each step as "error_embedded_estimate".
     """
 
-    def post_step(self, step, level_number):
+    def post_step(self, step, level_number, appendix=''):
         """
         Record embedded error estimate
 
@@ -27,6 +27,47 @@ def post_step(self, step, level_number):
             level=L.level_index,
             iter=step.status.iter,
             sweep=L.status.sweep,
-            type='error_embedded_estimate',
+            type=f'error_embedded_estimate{appendix}',
+            value=L.status.get('error_embedded_estimate'),
+        )
+
+
+class LogEmbeddedErrorEstimatePostIter(LogEmbeddedErrorEstimate):
+    """
+    Store the embedded error estimate after each iteration as "error_embedded_estimate_post_iteration".
+
+    Because the error estimate is computed after the hook is called, we record the value belonging to the last
+    iteration, which is also why we need to record something after the step, which belongs to the final iteration.
+    """
+
+    def post_iteration(self, step, level_number):
+        """
+        Record embedded error estimate
+
+        Args:
+            step (pySDC.Step.step): the current step
+            level_number (int): the current level number
+
+        Returns:
+            None
+        """
+        # check if the estimate is available at all
+        super().post_iteration(step, level_number)
+
+        L = step.levels[level_number]
+
+        if not L.status.get('error_embedded_estimate'):
+            return None
+
+        self.add_to_stats(
+            process=step.status.slot,
+            time=L.time + L.dt,
+            level=L.level_index,
+            iter=step.status.iter - 1,
+            sweep=L.status.sweep,
+            type='error_embedded_estimate_post_iteration',
             value=L.status.get('error_embedded_estimate'),
         )
+
+    def post_step(self, step, level_number):
+        super().post_step(step, level_number, appendix='_post_iteration')

pySDC/implementations/problem_classes/LeakySuperconductor.py

@@ -1,6 +1,6 @@
 import numpy as np
 import scipy.sparse as sp
-from scipy.sparse.linalg import cg, spsolve
+from scipy.sparse.linalg import spsolve
 
 from pySDC.core.Errors import ParameterError, ProblemError
 from pySDC.core.Problem import ptype
@@ -72,7 +72,7 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=mesh):
             xvalues = np.array([(i + 1) * self.dx for i in range(self.params.nvars)])
         elif self.params.bc == 'neumann-zero':
             self.dx = 1.0 / (self.params.nvars - 1)
-            xvalues = np.array([(i + 1) * self.dx for i in range(self.params.nvars)])
+            xvalues = np.array([i * self.dx for i in range(self.params.nvars)])
         else:
             raise ProblemError(f'Boundary conditions {self.params.bc} not implemented.')
 
@@ -85,13 +85,15 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=mesh):
             dim=1,
             bc=self.params.bc,
         )
-        self.A *= self.params.K
+        self.A *= self.params.K / self.params.Cv
 
         self.xv = xvalues
         self.Id = sp.eye(np.prod(self.params.nvars), format='csc')
 
         self.leak = np.logical_and(self.xv > self.params.leak_range[0], self.xv < self.params.leak_range[1])
 
+        self.total_newton_iter = 0  # store here how many iterations you needed for the Newton solver over the entire run to extract the desired information in the hooks
+
     def eval_f_non_linear(self, u, t):
         """
         Get the non-linear part of f
@@ -117,6 +119,8 @@ def eval_f_non_linear(self, u, t):
         me[0] = 0.0
         me[-1] = 0.0
 
+        me[:] /= self.params.Cv
+
         return me
 
     def eval_f(self, u, t):
@@ -131,7 +135,7 @@ def eval_f(self, u, t):
             dtype_f: The right hand side
         """
         f = self.dtype_f(self.init)
-        f[:] = (self.A.dot(u.flatten()).reshape(self.params.nvars) + self.eval_f_non_linear(u, t)) / self.params.Cv
+        f[:] = self.A.dot(u.flatten()).reshape(self.params.nvars) + self.eval_f_non_linear(u, t)
         return f
 
     def solve_system(self, rhs, factor, u0, t):
@@ -156,7 +160,7 @@ def get_non_linear_Jacobian(u):
                 u (dtype_u): Current solution
 
             Returns:
-                dtype_u: The derivative of the non-linear part of the solution w.r.t. to the solution.
+                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
@@ -172,28 +176,34 @@ def get_non_linear_Jacobian(u):
             me[0] = 0.0
             me[-1] = 0.0
 
-            me = self.dtype_u(self.init)
-            return me
+            me[:] /= self.params.Cv
+
+            return sp.diags(me, format='csc')
 
         u = self.dtype_u(u0)
         res = np.inf
-        for _n in range(0, self.params.newton_iter):
+        for n in range(0, self.params.newton_iter):
             # assemble G such that G(u) = 0 at the solution of the step
             G = u - factor * self.eval_f(u, t) - rhs
 
             res = np.linalg.norm(G, np.inf)
-            if res <= self.params.newton_tol or np.isnan(res):
+            if res <= self.params.newton_tol and n > 0:  # we want to make at least one Newton iteration
                 break
 
             # assemble Jacobian J of G
             J = self.Id - factor * (self.A + get_non_linear_Jacobian(u))
 
             # solve the linear system
-            delta = np.linalg.solve(J, G)
+            delta = spsolve(J, G)
+
+            if not np.isfinite(delta).all():
+                break
 
             # update solution
             u = u - delta
 
+        self.total_newton_iter += n
+
         return u
 
     def u_exact(self, t, u_init=None, t_init=None):
@@ -207,7 +217,7 @@ def u_exact(self, t, u_init=None, t_init=None):
             dtype_u: exact solution
         """
 
-        me = self.dtype_u(self.init, val=0)
+        me = self.dtype_u(self.init, val=0.0)
 
         if t > 0:
 
@@ -245,8 +255,8 @@ def eval_f(self, u, t):
         """
 
         f = self.dtype_f(self.init)
-        f.impl[:] = self.A.dot(u.flatten()).reshape(self.params.nvars) / self.params.Cv
-        f.expl[:] = self.eval_f_non_linear(u, t) / self.params.Cv
+        f.impl[:] = self.A.dot(u.flatten()).reshape(self.params.nvars)
+        f.expl[:] = self.eval_f_non_linear(u, t)
 
         return f