Test for AllenCahn_1D_FD.py (#370)

* Started with test for AC [WIP]

* Added test to capture errors and warnings

* Removed print line

* Tests for multi-implicit periodic case

* Removed lines

* Forgot black, your majesty..

* Removed zero matrix

pySDC/implementations/problem_classes/AllenCahn_1D_FD.py

@@ -3,7 +3,7 @@
 from scipy.sparse.linalg import spsolve
 
 from pySDC.core.Errors import ProblemError
-from pySDC.core.Problem import ptype
+from pySDC.core.Problem import ptype, WorkCounter
 from pySDC.helpers import problem_helper
 from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh, comp2_mesh
 
@@ -53,14 +53,9 @@ class allencahn_front_fullyimplicit(ptype):
         Spatial grid values.
     uext : dtype_u
         Contains additionally the external values of the boundary.
-    newton_itercount : int
-        Counter for iterations in Newton solver.
-    lin_itercount : int
-        Counter for iterations in linear solver.
-    newton_ncalls : int
-        Number of calls of Newton solver.
-    lin_ncalls : int
-        Number of calls of linear solver.
+    work_counters : WorkCounter
+        Counter for statistics. Here, number of Newton calls and number of evaluations
+        of right-hand side are counted.
     """
 
     dtype_u = mesh
@@ -101,18 +96,16 @@ def __init__(
         self.A, _ = problem_helper.get_finite_difference_matrix(
             derivative=2,
             order=2,
-            type='center',
+            stencil_type='center',
             dx=self.dx,
             size=self.nvars + 2,
             dim=1,
             bc='dirichlet-zero',
         )
         self.uext = self.dtype_u((self.init[0] + 2, self.init[1], self.init[2]), val=0.0)
 
-        self.newton_itercount = 0
-        self.lin_itercount = 0
-        self.newton_ncalls = 0
-        self.lin_ncalls = 0
+        self.work_counters['newton'] = WorkCounter()
+        self.work_counters['rhs'] = WorkCounter()
 
     def solve_system(self, rhs, factor, u0, t):
         """
@@ -184,6 +177,7 @@ def solve_system(self, rhs, factor, u0, t):
             # u -= gmres(dg, g, x0=z, tol=self.lin_tol)[0]
             # increase iteration count
             n += 1
+            self.work_counters['newton']()
 
         if np.isnan(res) and self.stop_at_nan:
             raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
@@ -193,9 +187,6 @@ def solve_system(self, rhs, factor, u0, t):
         if n == self.newton_maxiter:
             self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
 
-        self.newton_ncalls += 1
-        self.newton_itercount += n
-
         me = self.dtype_u(self.init)
         me[:] = u[:]
 
@@ -230,6 +221,7 @@ def eval_f(self, u, t):
             - 2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2 * u)
             - 6.0 * self.dw * u * (1.0 - u)
         )
+        self.work_counters['rhs']()
         return f
 
     def u_exact(self, t):
@@ -301,6 +293,7 @@ def eval_f(self, u, t):
         f = self.dtype_f(self.init)
         f.impl[:] = self.A.dot(self.uext)[1:-1]
         f.expl[:] = -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2 * u) - 6.0 * self.dw * u * (1.0 - u)
+        self.work_counters['rhs']()
         return f
 
     def solve_system(self, rhs, factor, u0, t):
@@ -402,7 +395,6 @@ def solve_system(self, rhs, factor, u0, t):
         n = 0
         res = 99
         while n < self.newton_maxiter:
-            # print(n)
             # form the function g(u), such that the solution to the nonlinear problem is a root of g
             self.uext[1:-1] = u[:]
             gprim = 1.0 / self.dx**2 * ((1.0 - a2) / (1.0 - a2 * (2.0 * u - 1.0) ** 2) - 1.0) * (2.0 * u - 1.0)
@@ -432,6 +424,7 @@ def solve_system(self, rhs, factor, u0, t):
             # u -= cg(dg, g, x0=z, tol=self.lin_tol)[0]
             # increase iteration count
             n += 1
+            self.work_counters['newton']()
 
         if np.isnan(res) and self.stop_at_nan:
             raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
@@ -441,9 +434,6 @@ def solve_system(self, rhs, factor, u0, t):
         if n == self.newton_maxiter:
             self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
 
-        self.newton_ncalls += 1
-        self.newton_itercount += n
-
         me = self.dtype_u(self.init)
         me[:] = u[:]
 
@@ -476,6 +466,7 @@ def eval_f(self, u, t):
         gprim = 1.0 / self.dx**2 * ((1.0 - a2) / (1.0 - a2 * (2.0 * u - 1.0) ** 2) - 1) * (2.0 * u - 1.0)
         f = self.dtype_f(self.init)
         f[:] = self.A.dot(self.uext)[1:-1] - 1.0 * gprim - 6.0 * self.dw * u * (1.0 - u)
+        self.work_counters['rhs']()
         return f
 
 
@@ -523,14 +514,9 @@ class allencahn_periodic_fullyimplicit(ptype):
         Distance between two spatial nodes.
     xvalues : np.1darray
         Spatial grid points.
-    newton_itercount : int
-        Number of iterations for Newton solver.
-    lin_itercount : int
-        Number of iterations for linear solver.
-    newton_ncalls : int
-        Number of calls of Newton solver.
-    lin_ncalls : int
-        Number of calls of linear solver.
+    work_counters : WorkCounter
+        Counter for statistics. Here, number of Newton calls and number of evaluations
+        of right-hand side are counted.
     """
 
     dtype_u = mesh
@@ -573,17 +559,15 @@ def __init__(
         self.A, _ = problem_helper.get_finite_difference_matrix(
             derivative=2,
             order=2,
-            type='center',
+            stencil_type='center',
             dx=self.dx,
             size=self.nvars,
             dim=1,
             bc='periodic',
         )
 
-        self.newton_itercount = 0
-        self.lin_itercount = 0
-        self.newton_ncalls = 0
-        self.lin_ncalls = 0
+        self.work_counters['newton'] = WorkCounter()
+        self.work_counters['rhs'] = WorkCounter()
 
     def solve_system(self, rhs, factor, u0, t):
         """
@@ -616,7 +600,6 @@ def solve_system(self, rhs, factor, u0, t):
         n = 0
         res = 99
         while n < self.newton_maxiter:
-            # print(n)
             # form the function g(u), such that the solution to the nonlinear problem is a root of g
             g = (
                 u
@@ -644,6 +627,7 @@ def solve_system(self, rhs, factor, u0, t):
             # u -= gmres(dg, g, x0=z, tol=self.lin_tol)[0]
             # increase iteration count
             n += 1
+            self.work_counters['newton']()
 
         if np.isnan(res) and self.stop_at_nan:
             raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
@@ -653,9 +637,6 @@ def solve_system(self, rhs, factor, u0, t):
         if n == self.newton_maxiter:
             self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
 
-        self.newton_ncalls += 1
-        self.newton_itercount += n
-
         me = self.dtype_u(self.init)
         me[:] = u[:]
 
@@ -679,6 +660,7 @@ def eval_f(self, u, t):
         """
         f = self.dtype_f(self.init)
         f[:] = self.A.dot(u) - 2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2 * u) - 6.0 * self.dw * u * (1.0 - u)
+        self.work_counters['rhs']()
         return f
 
     def u_exact(self, t):
@@ -735,7 +717,6 @@ def __init__(
         stop_at_nan=True,
     ):
         super().__init__(nvars, dw, eps, newton_maxiter, newton_tol, interval, radius, stop_at_nan)
-        self.A -= sp.eye(self.init) * 0.0 / self.eps**2
 
     def solve_system(self, rhs, factor, u0, t):
         r"""
@@ -785,6 +766,7 @@ def eval_f(self, u, t):
             - 6.0 * self.dw * u * (1.0 - u)
             + 0.0 / self.eps**2 * u
         )
+        self.work_counters['rhs']()
         return f
 
 
@@ -822,7 +804,6 @@ def __init__(
         stop_at_nan=True,
     ):
         super().__init__(nvars, dw, eps, newton_maxiter, newton_tol, interval, radius, stop_at_nan)
-        self.A -= sp.eye(self.init) * 0.0 / self.eps**2
 
     def solve_system_1(self, rhs, factor, u0, t):
         r"""
@@ -872,6 +853,7 @@ def eval_f(self, u, t):
             - 6.0 * self.dw * u * (1.0 - u)
             + 0.0 / self.eps**2 * u
         )
+        self.work_counters['rhs']()
         return f
 
     def solve_system_2(self, rhs, factor, u0, t):
@@ -905,7 +887,6 @@ def solve_system_2(self, rhs, factor, u0, t):
         n = 0
         res = 99
         while n < self.newton_maxiter:
-            # print(n)
             # form the function g(u), such that the solution to the nonlinear problem is a root of g
             g = (
                 u
@@ -932,6 +913,7 @@ def solve_system_2(self, rhs, factor, u0, t):
             # u -= gmres(dg, g, x0=z, tol=self.lin_tol)[0]
             # increase iteration count
             n += 1
+            self.work_counters['newton']()
 
         if np.isnan(res) and self.stop_at_nan:
             raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
@@ -941,9 +923,6 @@ def solve_system_2(self, rhs, factor, u0, t):
         if n == self.newton_maxiter:
             self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
 
-        self.newton_ncalls += 1
-        self.newton_itercount += n
-
         me = self.dtype_u(self.init)
         me[:] = u[:]