diff --git a/pySDC/implementations/problem_classes/AdvectionDiffusionEquation_1D_FFT.py b/pySDC/implementations/problem_classes/AdvectionDiffusionEquation_1D_FFT.py index 8e164ed88f..beb56d3c79 100644 --- a/pySDC/implementations/problem_classes/AdvectionDiffusionEquation_1D_FFT.py +++ b/pySDC/implementations/problem_classes/AdvectionDiffusionEquation_1D_FFT.py @@ -1,7 +1,7 @@ import numpy as np from pySDC.core.Errors import ParameterError, ProblemError -from pySDC.core.Problem import ptype +from pySDC.core.Problem import ptype, WorkCounter from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh @@ -70,6 +70,8 @@ def __init__(self, nvars=256, c=1.0, freq=-1, nu=0.02, L=1.0): self.ddx = kx * 1j self.lap = -(kx**2) + self.work_counters['rhs'] = WorkCounter() + def eval_f(self, u, t): """ Routine to evaluate the right-hand side of the problem. @@ -94,6 +96,7 @@ def eval_f(self, u, t): f.impl[:] = np.fft.irfft(tmp_impl) f.expl[:] = np.fft.irfft(tmp_expl) + self.work_counters['rhs']() return f def solve_system(self, rhs, factor, u0, t): @@ -153,6 +156,8 @@ def u_exact(self, t): for i in range(self.init[0]): x = self.xvalues[i] - self.c * t + k * self.L me[i] += np.sqrt(t00) / np.sqrt(t00 + t) * np.exp(-(x**2) / (4.0 * self.nu * (t00 + t))) + else: + raise ParameterError("There is no exact solution implemented for negative frequency and negative nu!") return me @@ -177,6 +182,9 @@ class advectiondiffusion1d_implicit(advectiondiffusion1d_imex): This class has the same attributes as the class it inherits from. """ + dtype_u = mesh + dtype_f = mesh + def eval_f(self, u, t): """ Routine to evaluate the right-hand side of the problem. @@ -199,6 +207,7 @@ def eval_f(self, u, t): tmp = self.nu * self.lap * tmp_u - self.c * self.ddx * tmp_u f[:] = np.fft.irfft(tmp) + self.work_counters['rhs'] return f def solve_system(self, rhs, factor, u0, t): diff --git a/pySDC/tests/test_problems/test_AdvectionDiffusionEquation_1D_FFT.py b/pySDC/tests/test_problems/test_AdvectionDiffusionEquation_1D_FFT.py new file mode 100644 index 0000000000..0f5736f143 --- /dev/null +++ b/pySDC/tests/test_problems/test_AdvectionDiffusionEquation_1D_FFT.py @@ -0,0 +1,94 @@ +import pytest + + +@pytest.mark.base +@pytest.mark.parametrize('freq', [-1, 0, 1]) +@pytest.mark.parametrize('nu', [0.02, -0.02]) +def test_imex_vs_implicit(freq, nu): + import numpy as np + from pySDC.core.Errors import ParameterError, ProblemError + from pySDC.implementations.problem_classes.AdvectionDiffusionEquation_1D_FFT import ( + advectiondiffusion1d_imex, + advectiondiffusion1d_implicit, + ) + + problem_params = { + 'nvars': 32, + 'c': 1.0, + 'freq': freq, + 'nu': nu, + 'L': 1.0, + } + + imex = advectiondiffusion1d_imex(**problem_params) + fully_impl = advectiondiffusion1d_implicit(**problem_params) + + t0 = 0.0 + if freq < 0 and nu < 0: + # check if ParameterError is raised correctly for freq < 0 and nu < 0 + with pytest.raises(ParameterError): + imex.u_exact(t0) + else: + # test if evaluations of right-hand side do match + u0 = imex.u_exact(t0) + + tmp = imex.eval_f(u0, t0) + f_imex = tmp.expl + tmp.impl + f_full = fully_impl.eval_f(u0, t0) + + assert np.allclose( + f_imex, f_full + ), 'Evaluation of right-hand side in semi-explicit case and in fully-implicit case do not match!' + + # test if solving one time step satisfies a specific error threshold + dt = 1e-3 + args = { + 'rhs': u0, + 'factor': dt, + 'u0': u0, + 't': t0, + } + + u_ex = imex.u_exact(t0 + dt) + + u_imex = imex.solve_system(**args) + u_full = fully_impl.solve_system(**args) + + e_imex = abs(u_ex - u_imex) + e_full = abs(u_ex - u_full) + + e_tol_imex = { + -1: { + 0.02: 0.011, + }, + 0: { + -0.02: 0.076, + 0.02: 0.055, + }, + 1: { + -0.02: 0.0063, + 0.02: 0.0063, + }, + } + + e_tol_full = { + -1: { + 0.02: 0.00021, + }, + 0: { + -0.02: 0.078, + 0.02: 0.064, + }, + 1: { + -0.02: 2.01e-05, + 0.02: 2e-05, + }, + } + assert e_imex < e_tol_imex[freq][nu], "Error is too large in semi-explicit case!" + assert e_full < e_tol_full[freq][nu], "Error is too large in fully-implicit case!" + + # check if ProblemError is raised correctly in case if nvars % 2 != 0 + problem_params.update({'nvars': 31}) + with pytest.raises(ProblemError): + imex_test = advectiondiffusion1d_imex(**problem_params) + fully_impl_test = advectiondiffusion1d_implicit(**problem_params)