From 8a29f486dd994209752bc7f97e608107f7644e96 Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Fri, 20 Oct 2023 10:52:42 +0200 Subject: [PATCH] Fixed Neumann-BCs for up to second order for advection and heat equation --- pySDC/helpers/problem_helper.py | 98 ++++++++++++++++--- .../problem_classes/AllenCahn_1D_FD.py | 4 +- .../problem_classes/AllenCahn_2D_FD.py | 2 +- .../GeneralizedFisher_1D_FD_implicit.py | 2 +- .../HeatEquation_ND_FD_CuPy.py | 2 +- .../implementations/problem_classes/Quench.py | 34 ++----- .../problem_classes/generic_ND_FD.py | 2 +- pySDC/projects/Resilience/quench.py | 2 +- .../tests/test_helpers/test_problem_helper.py | 65 +++++++++++- 9 files changed, 162 insertions(+), 49 deletions(-) diff --git a/pySDC/helpers/problem_helper.py b/pySDC/helpers/problem_helper.py index b8e35642dd..544840aced 100644 --- a/pySDC/helpers/problem_helper.py +++ b/pySDC/helpers/problem_helper.py @@ -77,7 +77,18 @@ def get_finite_difference_stencil(derivative, order, stencil_type=None, steps=No def get_finite_difference_matrix( - derivative, order, stencil_type=None, steps=None, dx=None, size=None, dim=None, bc=None, cupy=False + derivative, + order, + stencil_type=None, + steps=None, + dx=None, + size=None, + dim=None, + bc=None, + cupy=False, + bc_val_left=0.0, + bc_val_right=0.0, + neumann_bc_order=None, ): """ Build FD matrix from stencils, with boundary conditions. @@ -93,9 +104,13 @@ def get_finite_difference_matrix( dim (int): Number of dimensions bc (str): Boundary conditions for both sides cupy (bool): Construct a GPU ready matrix if yes + bc_val_left (float): Value for left boundary condition + bc_val_right (float): Value for right boundary condition + neumann_bc_order (int): Order of accuracy of Neumann BC, optional Returns: Sparse matrix: Finite difference matrix + numpy.ndarray: Vector containing information about the boundary conditions """ if cupy: import cupyx.scipy.sparse as sp @@ -107,17 +122,74 @@ def get_finite_difference_matrix( derivative=derivative, order=order, stencil_type=stencil_type, steps=steps ) + b = np.zeros(size**dim) + if "dirichlet" in bc: A_1d = sp.diags(coeff, steps, shape=(size, size), format='csc') elif "neumann" in bc: - A_1d = sp.diags(coeff, steps, shape=(size, size), format='lil') - - # replace the first and last lines with one-sided stencils for a first derivative of the same order as the rest - bc_coeff, bc_steps = get_finite_difference_stencil(derivative=1, order=order, stencil_type='forward') - A_1d[0, :] = 0.0 - A_1d[0, bc_steps] = bc_coeff * (dx ** (derivative - 1)) - A_1d[-1, :] = 0.0 - A_1d[-1, -bc_steps - 1] = bc_coeff * (dx ** (derivative - 1)) + """ + We will solve only for values within the boundary because the values on the boundary can be determined from the + discretization of the boundary condition (BC). Therefore, we discretize both the BC and the original problem + such that the stencils reach into the boundary with one element. We then proceed to eliminate the values on the + boundary, which modifies the finite difference matrix and yields an inhomogeneity if the BCs are inhomogeneous. + + Keep in mind that centered stencils are often more efficient in terms of sparsity of the resulting matrix. High + order centered discretizations will reach beyond the boundary and hence need to be replaced with lopsided + stencils near the boundary, changing the number of non-zero diagonals. + """ + # check if we need to alter the sparsity structure because of the BC + if steps.min() < -1 or steps.max() > 1: + A_1d = sp.diags(coeff, steps, shape=(size, size), format='lil') + else: + A_1d = sp.diags(coeff, steps, shape=(size, size), format='csc') + + if derivative == 1 and (steps.min() < -1 or steps.max() > 1): + neumann_bc_order = neumann_bc_order if neumann_bc_order else order + 1 + assert neumann_bc_order != order, 'Need different stencils for BC and the rest' + else: + neumann_bc_order = neumann_bc_order if neumann_bc_order else order + + if dim > 1 and (bc_val_left != 0.0 or bc_val_right != 0): + raise NotImplementedError( + f'Non-zero Neumann BCs are only implemented in 1D. You asked for {dim} dimensions.' + ) + + # ---- left boundary ---- + # generate the one-sided stencil to discretize the first derivative at the boundary + bc_coeff_left, bc_steps_left = get_finite_difference_stencil( + derivative=1, order=neumann_bc_order, stencil_type='forward' + ) + + # check if we can just use the available stencil or if we need to generate a new one + if steps.min() == -1: + coeff_left = coeff.copy() + else: # need to generate lopsided stencils + raise NotImplementedError( + 'Neumann BCs on the right are not implemented for your desired stencil. Maybe try a lower order' + ) + + # modify system matrix and inhomogeneity according to BC + b[0] = bc_val_left * (coeff_left[0] / dx**derivative) / (bc_coeff_left[0] / dx) + A_1d[0, : len(bc_coeff_left) - 1] -= coeff_left[0] / bc_coeff_left[0] * bc_coeff_left[1:] + + # ---- right boundary ---- + # generate the one-sided stencil to discretize the first derivative at the boundary + bc_coeff_right, bc_steps_right = get_finite_difference_stencil( + derivative=1, order=neumann_bc_order, stencil_type='backward' + ) + + # check if we can just use the available stencil or if we need to generate a new one + if steps.max() == +1: + coeff_right = coeff.copy() + else: # need to generate lopsided stencils + raise NotImplementedError( + 'Neumann BCs on the right are not implemented for your desired stencil. Maybe try a lower order' + ) + + # modify system matrix and inhomogeneity according to BC + b[-1] = bc_val_right * (coeff_right[-1] / dx**derivative) / (bc_coeff_right[0] / dx) + A_1d[-1, -len(bc_coeff_right) + 1 :] -= coeff_right[-1] / bc_coeff_right[0] * bc_coeff_right[::-1][:-1] + elif bc == 'periodic': A_1d = 0 * sp.eye(size, format='csc') for i in steps: @@ -129,6 +201,7 @@ def get_finite_difference_matrix( else: raise NotImplementedError(f'Boundary conditions \"{bc}\" not implemented.') + # TODO: extend the BCs to higher dimensions if dim == 1: A = A_1d elif dim == 2: @@ -144,7 +217,7 @@ def get_finite_difference_matrix( A /= dx**derivative - return A + return A, b def get_1d_grid(size, bc, left_boundary=0.0, right_boundary=1.0): @@ -165,12 +238,9 @@ def get_1d_grid(size, bc, left_boundary=0.0, right_boundary=1.0): if bc == 'periodic': dx = L / size xvalues = np.array([left_boundary + dx * i for i in range(size)]) - elif "dirichlet" in bc: + elif "dirichlet" in bc or "neumann" in bc: dx = L / (size + 1) xvalues = np.array([left_boundary + dx * (i + 1) for i in range(size)]) - elif "neumann" in bc: - dx = L / (size - 1) - xvalues = np.array([left_boundary + dx * i for i in range(size)]) else: raise NotImplementedError(f'Boundary conditions \"{bc}\" not implemented.') diff --git a/pySDC/implementations/problem_classes/AllenCahn_1D_FD.py b/pySDC/implementations/problem_classes/AllenCahn_1D_FD.py index b8979ee2b8..b288ed0c54 100644 --- a/pySDC/implementations/problem_classes/AllenCahn_1D_FD.py +++ b/pySDC/implementations/problem_classes/AllenCahn_1D_FD.py @@ -98,7 +98,7 @@ def __init__( self.dx = (self.interval[1] - self.interval[0]) / (self.nvars + 1) self.xvalues = np.array([(i + 1 - (self.nvars + 1) / 2) * self.dx for i in range(self.nvars)]) - self.A = problem_helper.get_finite_difference_matrix( + self.A, _ = problem_helper.get_finite_difference_matrix( derivative=2, order=2, type='center', @@ -570,7 +570,7 @@ def __init__( self.dx = (self.interval[1] - self.interval[0]) / self.nvars self.xvalues = np.array([self.interval[0] + i * self.dx for i in range(self.nvars)]) - self.A = problem_helper.get_finite_difference_matrix( + self.A, _ = problem_helper.get_finite_difference_matrix( derivative=2, order=2, type='center', diff --git a/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py b/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py index 6d7aaa02c0..a391dd5bec 100644 --- a/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py +++ b/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py @@ -108,7 +108,7 @@ def __init__( # compute dx and get discretization matrix A self.dx = 1.0 / self.nvars[0] - self.A = problem_helper.get_finite_difference_matrix( + self.A, _ = problem_helper.get_finite_difference_matrix( derivative=2, order=self.order, stencil_type='center', diff --git a/pySDC/implementations/problem_classes/GeneralizedFisher_1D_FD_implicit.py b/pySDC/implementations/problem_classes/GeneralizedFisher_1D_FD_implicit.py index c28dc3334c..c260307a36 100644 --- a/pySDC/implementations/problem_classes/GeneralizedFisher_1D_FD_implicit.py +++ b/pySDC/implementations/problem_classes/GeneralizedFisher_1D_FD_implicit.py @@ -102,7 +102,7 @@ def __init__( # compute dx and get discretization matrix A self.dx = (self.interval[1] - self.interval[0]) / (self.nvars + 1) - self.A = problem_helper.get_finite_difference_matrix( + self.A, _ = problem_helper.get_finite_difference_matrix( derivative=2, order=2, stencil_type='center', diff --git a/pySDC/implementations/problem_classes/HeatEquation_ND_FD_CuPy.py b/pySDC/implementations/problem_classes/HeatEquation_ND_FD_CuPy.py index b8fd1aee3c..00763f7cc8 100644 --- a/pySDC/implementations/problem_classes/HeatEquation_ND_FD_CuPy.py +++ b/pySDC/implementations/problem_classes/HeatEquation_ND_FD_CuPy.py @@ -148,7 +148,7 @@ def __init__( else: raise ProblemError(f'Boundary conditions {self.bc} not implemented.') - self.A = problem_helper.get_finite_difference_matrix( + self.A, _ = problem_helper.get_finite_difference_matrix( derivative=2, order=self.order, stencil_type=self.stencil_type, diff --git a/pySDC/implementations/problem_classes/Quench.py b/pySDC/implementations/problem_classes/Quench.py index a6ffc4e8a3..350499526d 100644 --- a/pySDC/implementations/problem_classes/Quench.py +++ b/pySDC/implementations/problem_classes/Quench.py @@ -89,8 +89,8 @@ def __init__( self, Cv=1000.0, K=1000.0, - u_thresh=1e-2, - u_max=2e-2, + u_thresh=3e-2, + u_max=6e-2, Q_max=1.0, leak_range=(0.45, 0.55), leak_type='linear', @@ -158,7 +158,7 @@ def __init__( # setup finite difference discretization from problem helper self.dx, xvalues = problem_helper.get_1d_grid(size=self.nvars, bc=self.bc) - self.A = problem_helper.get_finite_difference_matrix( + self.A, self.b = problem_helper.get_finite_difference_matrix( derivative=2, order=self.order, stencil_type=self.stencil_type, @@ -179,26 +179,6 @@ def __init__( if not self.direct_solver: self.work_counters['linear'] = WorkCounter() - def apply_BCs(self, f): - """ - Apply boundary conditions to the right hand side. Please supply only the nonlinear part as f! - - Args: - dtype_f: nonlinear (!) part of the right hand side - - Returns: - dtype_f: nonlinear part of the right hand side with boundary conditions - """ - if self.bc == 'neumann-zero': - f[0] = 0.0 - f[-1] = 0.0 - elif self.bc == 'periodic': - pass - else: - raise NotImplementedError(f'BCs \"{self.bc}\" not implemented!') - - return f - def eval_f_non_linear(self, u, t): """ Get the non-linear part of f. @@ -239,7 +219,7 @@ def eval_f_non_linear(self, u, t): me[:] /= self.Cv - return self.apply_BCs(me) + return me def eval_f(self, u, t): """ @@ -258,7 +238,7 @@ def eval_f(self, u, t): The right-hand side of the problem. """ f = self.dtype_f(self.init) - f[:] = self.A.dot(u.flatten()).reshape(self.nvars) + self.eval_f_non_linear(u, t) + f[:] = self.A.dot(u.flatten()).reshape(self.nvars) + self.b + self.eval_f_non_linear(u, t) self.work_counters['rhs']() return f @@ -301,7 +281,7 @@ def get_non_linear_Jacobian(self, u): me[:] /= self.Cv - return sp.diags(self.apply_BCs(me), format='csc') + return sp.diags(me, format='csc') def solve_system(self, rhs, factor, u0, t): r""" @@ -523,7 +503,7 @@ def eval_f(self, u, t): f = self.dtype_f(self.init) f.impl[:] = self.A.dot(u.flatten()).reshape(self.nvars) - f.expl[:] = self.eval_f_non_linear(u, t) + f.expl[:] = self.eval_f_non_linear(u, t) + self.b self.work_counters['rhs']() return f diff --git a/pySDC/implementations/problem_classes/generic_ND_FD.py b/pySDC/implementations/problem_classes/generic_ND_FD.py index ac5663bfa0..2704940557 100644 --- a/pySDC/implementations/problem_classes/generic_ND_FD.py +++ b/pySDC/implementations/problem_classes/generic_ND_FD.py @@ -117,7 +117,7 @@ def __init__( dx, xvalues = problem_helper.get_1d_grid(size=nvars[0], bc=bc, left_boundary=0.0, right_boundary=1.0) - self.A = problem_helper.get_finite_difference_matrix( + self.A, _ = problem_helper.get_finite_difference_matrix( derivative=derivative, order=order, stencil_type=stencil_type, diff --git a/pySDC/projects/Resilience/quench.py b/pySDC/projects/Resilience/quench.py index fec74ab9b2..3d176878da 100644 --- a/pySDC/projects/Resilience/quench.py +++ b/pySDC/projects/Resilience/quench.py @@ -334,7 +334,7 @@ def compare_imex_full(plotting=False, leak_type='linear'): custom_description['step_params'] = {'maxiter': maxiter} custom_description['sweeper_params'] = {'num_nodes': num_nodes} custom_description['convergence_controllers'] = { - Adaptivity: {'e_tol': 1e-6}, + Adaptivity: {'e_tol': 1e-7}, } custom_controller_params = {'logger_level': 15} diff --git a/pySDC/tests/test_helpers/test_problem_helper.py b/pySDC/tests/test_helpers/test_problem_helper.py index 8514244e19..2b22d15d04 100644 --- a/pySDC/tests/test_helpers/test_problem_helper.py +++ b/pySDC/tests/test_helpers/test_problem_helper.py @@ -6,7 +6,7 @@ def fd_stencil_single(derivative, order, stencil_type): """ Make a single tests where we generate a finite difference stencil using the generic framework above and compare to - harscoded stencils that were implemented in a previous version of the code. + hardcoded stencils that were implemented in a previous version of the code. Args: derivative (int): Order of the derivative @@ -133,3 +133,66 @@ def test_fd_stencils(): # test if we get the correct result when we put in steps rather than a stencil_type new_coeff, _ = get_finite_difference_stencil(derivative=2, order=2, steps=steps) assert np.allclose(coeff, new_coeff), f"Error when setting steps yourself! Expected {expect_coeff}, got {coeff}." + + +@pytest.mark.base +@pytest.mark.parametrize('bc_left', [0.0, 7.0]) +@pytest.mark.parametrize('bc_right', [0.0, 9.0]) +@pytest.mark.parametrize('dx', [0.1, 10.0]) +@pytest.mark.parametrize('derivative', [1, 2]) +@pytest.mark.parametrize('order', [2]) +def test_Neumann_bcs(derivative, bc_left, bc_right, dx, order): + from pySDC.helpers.problem_helper import get_finite_difference_matrix + + A, b = get_finite_difference_matrix( + derivative=derivative, + order=order, + stencil_type='center', + bc='neumann', + bc_val_left=bc_left, + bc_val_right=bc_right, + dim=1, + size=4, + dx=dx, + ) + + if derivative == 1: + expect = np.zeros(A.shape[0]) + expect[0] = -2 / (3 * dx) + expect[1] = +2 / (3 * dx) + assert np.allclose( + A.toarray()[0, :], expect + ), f'Error in left boundary, expected {expect} got {A.toarray()[0, :]}!' + expect = np.zeros(A.shape[0]) + expect[-2] = -2 / (3 * dx) + expect[-1] = +2 / (3 * dx) + assert np.allclose( + A.toarray()[-1, :], expect + ), f'Error in right boundary, expected {expect} got {A.toarray()[-1, :]}!' + + assert np.isclose( + b[-1], bc_right / 3.0 + ), f'Error in right boundary value! Expected {bc_right / 3.}, got {b[-1]}' + assert np.isclose(b[0], bc_left / 3.0), f'Error in left boundary value! Expected {bc_left / 3}, got {b[0]}' + + if derivative == 2: + expect = np.zeros(A.shape[0]) + expect[0] = -2 / (3 * dx**2) + expect[1] = +2 / (3 * dx**2) + assert np.allclose( + A.toarray()[0, :], expect + ), f'Error in left boundary, expected {expect} got {A.toarray()[0, :]}!' + assert np.allclose( + A.toarray()[-1, :], expect[::-1] + ), f'Error in right boundary, expected {expect[::-1]} got {A.toarray()[-1, :]}!' + + assert np.isclose( + b[-1], bc_right * 2 / (3 * dx) + ), f'Error in right boundary value! Expected {bc_right * 2 / (3*dx)}, got {b[-1]}' + assert np.isclose( + b[0], -bc_left * 2 / (3 * dx) + ), f'Error in left boundary value! Expected {-bc_left * 2 / (3*dx)}, got {b[0]}' + + +if __name__ == '__main__': + test_Neumann_bcs(2, 0, 1, 3.0 / 3.0, 2)