From 73cb8e30428b01f82728136e6cbc095124cb586b Mon Sep 17 00:00:00 2001 From: Lisa Wimmer <68507897+lisawim@users.noreply.github.com> Date: Fri, 27 Oct 2023 12:07:36 +0200 Subject: [PATCH] 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 --- .../problem_classes/AllenCahn_1D_FD.py | 67 ++--- .../test_problems/test_AllenCahn_1D_FD.py | 245 ++++++++++++++++++ 2 files changed, 268 insertions(+), 44 deletions(-) create mode 100644 pySDC/tests/test_problems/test_AllenCahn_1D_FD.py diff --git a/pySDC/implementations/problem_classes/AllenCahn_1D_FD.py b/pySDC/implementations/problem_classes/AllenCahn_1D_FD.py index b288ed0c54..78bd9575b9 100644 --- a/pySDC/implementations/problem_classes/AllenCahn_1D_FD.py +++ b/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,7 +96,7 @@ 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, @@ -109,10 +104,8 @@ def __init__( ) 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[:] diff --git a/pySDC/tests/test_problems/test_AllenCahn_1D_FD.py b/pySDC/tests/test_problems/test_AllenCahn_1D_FD.py new file mode 100644 index 0000000000..b45f953e01 --- /dev/null +++ b/pySDC/tests/test_problems/test_AllenCahn_1D_FD.py @@ -0,0 +1,245 @@ +import pytest + + +@pytest.mark.base +def test_front_imex_vs_fully_implicit(): + import numpy as np + from pySDC.implementations.problem_classes.AllenCahn_1D_FD import ( + allencahn_front_fullyimplicit, + allencahn_front_semiimplicit, + allencahn_front_finel, + ) + + problem_params = { + 'nvars': 31, + } + + t0 = 0.0 + + # for imex and fully-implicit test if right-hand side does match + imex = allencahn_front_semiimplicit(**problem_params) + full = allencahn_front_fullyimplicit(**problem_params) + finel = allencahn_front_finel(**problem_params) + + u0 = imex.u_exact(t0) + + f_tmp = imex.eval_f(u0, t0) + f_imex = f_tmp.impl + f_tmp.expl + f_full = full.eval_f(u0, t0) + + assert np.allclose( + f_imex, f_full + ), "Evaluation of right-hand side in semi-implicit case and in fully-implicit case do not match!" + + # perform one time step and test the error + dt = 1e-4 + args = { + 'rhs': u0, + 'factor': dt, + 'u0': u0, + 't': t0, + } + + u_imex = imex.solve_system(**args) + u_full = full.solve_system(**args) + u_finel = finel.solve_system(**args) + + u_exact = imex.u_exact(t0 + dt) + + e_imex = abs(u_imex - u_exact) / abs(u_exact) + e_full = abs(u_full - u_exact) / abs(u_exact) + e_finel = abs(u_finel - u_exact) / abs(u_exact) + + assert e_imex < 8e-2, f"Error is too large in semi-implicit case! Got {e_imex}" + assert e_full < 2e-3, f"Error is too large in fully-implicit case! Got {e_full}" + assert e_finel < 2.5e-11, f"Error is too large in case of Finel's trick! Got {e_finel}" + + # check if number of right-hand side evaluations do match + rhs_count_imex = imex.work_counters['rhs'].niter + rhs_count_full = full.work_counters['rhs'].niter + + assert ( + rhs_count_imex == rhs_count_full + ), "Number of right-hand side evaluations for semi-implicit vs fully-implicit do not match!" + + +@pytest.mark.base +def test_periodic_imex_vs_fully_implicit_vs_multi_implicit(): + import numpy as np + from pySDC.implementations.problem_classes.AllenCahn_1D_FD import ( + allencahn_periodic_semiimplicit, + allencahn_periodic_fullyimplicit, + allencahn_periodic_multiimplicit, + ) + + problem_params = { + 'nvars': 32, + } + + t0 = 0.0 + + # for imex and fully-implicit test if right-hand side does match + imex = allencahn_periodic_semiimplicit(**problem_params) + full = allencahn_periodic_fullyimplicit(**problem_params) + multi = allencahn_periodic_multiimplicit(**problem_params) + + u0 = imex.u_exact(t0) + + f_tmp_full = imex.eval_f(u0, t0) + f_imex = f_tmp_full.impl + f_tmp_full.expl + f_full = full.eval_f(u0, t0) + f_tmp_multi = multi.eval_f(u0, t0) + f_multi = f_tmp_multi.comp1 + f_tmp_multi.comp2 + + assert np.allclose( + f_imex, f_full + ), "Evaluation of right-hand side in semi-implicit case and in fully-implicit case do not match!" + assert np.allclose( + f_full, f_multi + ), "Evaluation of right-hand side in fully-implicit case and in multi-implicit case do not match!" + + # perform one time step and check errors + dt = 1e-4 + args = { + 'rhs': u0, + 'factor': dt, + 'u0': u0, + 't': t0, + } + + u_imex = imex.solve_system(**args) + u_full = full.solve_system(**args) + u_multi1 = multi.solve_system_1(**args) + u_multi2 = multi.solve_system_2(**args) + + u_exact = imex.u_exact(t0 + dt) + + e_imex = abs(u_imex - u_exact) / abs(u_exact) + e_full = abs(u_full - u_exact) / abs(u_exact) + e_multi1 = abs(u_multi1 - u_exact) / abs(u_exact) + e_multi2 = abs(u_multi2 - u_exact) / abs(u_exact) + + assert e_imex < 1e-2, f"Error is too large in semi-implicit case! Got {e_imex}" + assert e_full < 1.2e-3, f"Error is too large in fully-implicit case! Got {e_full}" + assert e_multi1 < 1e-2, f"Error is too large in multi-implicit case solving the Laplacian part! Got {e_multi1}" + assert ( + e_multi2 < 1.2e-2 + ), f"Error is too large in multi-implicit case solving the part without the Laplacian! Got {e_multi2}" + + # check if number of right-hand side evaluations do match + rhs_count_imex = imex.work_counters['rhs'].niter + rhs_count_full = full.work_counters['rhs'].niter + rhs_count_multi = multi.work_counters['rhs'].niter + + assert ( + rhs_count_imex == rhs_count_full + ), "Number of right-hand side evaluations for semi-implicit vs fully-implicit do not match!" + assert ( + rhs_count_full == rhs_count_multi + ), "Number of right-hand side evaluations for fully-implicit vs multi-implicit do not match!" + + +@pytest.mark.base +@pytest.mark.parametrize('stop_at_nan', [True, False]) +def test_capture_errors_and_warnings(caplog, stop_at_nan): + """ + Test if errors and warnings are raised correctly. + """ + import numpy as np + from pySDC.core.Errors import ProblemError + from pySDC.implementations.problem_classes.AllenCahn_1D_FD import ( + allencahn_front_fullyimplicit, + allencahn_front_semiimplicit, + allencahn_front_finel, + allencahn_periodic_fullyimplicit, + allencahn_periodic_semiimplicit, + allencahn_periodic_multiimplicit, + ) + + newton_maxiter = 1 + problem_params = { + 'stop_at_nan': stop_at_nan, + 'newton_tol': 1e-13, + 'newton_maxiter': newton_maxiter, + } + + full_front = allencahn_front_fullyimplicit(**problem_params) + imex_front = allencahn_front_semiimplicit(**problem_params) + finel_front = allencahn_front_finel(**problem_params) + + full_periodic = allencahn_periodic_fullyimplicit(**problem_params) + imex_periodic = allencahn_periodic_semiimplicit(**problem_params) + multi_periodic = allencahn_periodic_multiimplicit(**problem_params) + + t0 = 0.0 + dt = 1e-3 + + u0_front = full_front.u_exact(t0) + u0_periodic = full_periodic.u_exact(t0) + + args_front = { + 'rhs': u0_front, + 'factor': np.nan, + 'u0': u0_front, + 't': t0, + } + + args_periodic = { + 'rhs': u0_periodic, + 'factor': np.nan, + 'u0': u0_periodic, + 't': t0, + } + + if stop_at_nan: + # test if ProblemError is raised + with pytest.raises(ProblemError): + full_front.solve_system(**args_front) + imex_front.solve_system(**args_front) + finel_front.solve_system(**args_front) + + with pytest.raises(ProblemError): + full_periodic.solve_system(**args_periodic) + imex_periodic.solve_system(**args_periodic) + multi_periodic.solve_system_2(**args_periodic) + + else: + # test if warnings are raised when nan values arise + full_front.solve_system(**args_front) + assert f'Newton got nan after {newton_maxiter} iterations...' in caplog.text + assert 'Newton did not converge after 1 iterations, error is nan' in caplog.text + assert ( + full_front.work_counters['newton'].niter == newton_maxiter + ), 'Number of Newton iterations in fully-implicit front case does not match with maximum number of iterations!' + caplog.clear() + + finel_front.solve_system(**args_front) + assert f'Newton got nan after {newton_maxiter} iterations...' in caplog.text + assert 'Newton did not converge after 1 iterations, error is nan' in caplog.text + assert ( + finel_front.work_counters['newton'].niter == newton_maxiter + ), "Number of Newton iterations in case of using Finel's trick does not match with maximum number of iterations!" + caplog.clear() + + full_periodic.solve_system(**args_periodic) + assert f'Newton got nan after {newton_maxiter} iterations...' in caplog.text + assert 'Newton did not converge after 1 iterations, error is nan' in caplog.text + assert ( + full_periodic.work_counters['newton'].niter == newton_maxiter + ), 'Number of Newton iterations in fully-implicit periodic case does not match with maximum number of iterations!' + caplog.clear() + + multi_periodic.solve_system_2(**args_periodic) + assert f'Newton got nan after {newton_maxiter} iterations...' in caplog.text + assert 'Newton did not converge after 1 iterations, error is nan' in caplog.text + assert ( + multi_periodic.work_counters['newton'].niter == newton_maxiter + ), 'Number of Newton iterations in multi-implicit periodic case does not match with maximum number of iterations!' + caplog.clear() + + +@pytest.mark.base +def test_solve_system_with_reference(): + """ + Computes a solution using the Newton solver and compares it with a Newton-Krylov solver as reference. + """