Fixed Neumann-BCs for up to second order for advection and heat equation

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.')
 

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',

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',

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',

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,

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

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,

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}

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)