Added exponential growth to quench problem

Author: Thomas Baumann

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
@@ -160,7 +167,13 @@ def get_non_linear_Jacobian(self, u):
         Q_max = self.params.Q_max
         me = self.dtype_u(self.init)
 
-        me[:] = Q_max / (u_max - u_thresh)
+        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!')
+
         me[u < u_thresh] = 0
         me[u > u_max] = 0
         me[self.leak] = 0

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.
 
@@ -218,6 +218,7 @@ def compare_imex_full(plotting=False):
         '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}
@@ -229,7 +230,7 @@ def compare_imex_full(plotting=False):
             custom_description=custom_description,
             custom_controller_params=custom_controller_params,
             imex=imex,
-            Tend=5e2,
+            Tend=4.3e2,
             hook_class=[LogWork, LogGlobalErrorPostRun],
         )
 

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)