Neumann BCs in Chebychev and ultraspherical bases

pySDC/helpers/spectral_helper.py

@@ -535,6 +535,8 @@ def get_BC(self, kind, **kwargs):
             return self.get_integ_BC_row(**kwargs)
         elif kind.lower() == 'dirichlet':
             return self.get_Dirichlet_BC_row(**kwargs)
+        elif kind.lower() == 'neumann':
+            return self.get_Neumann_BC_row(**kwargs)
         else:
             return super().get_BC(kind)
 
@@ -574,6 +576,27 @@ def get_Dirichlet_BC_row(self, x):
         else:
             raise NotImplementedError(f'Don\'t know how to generate Dirichlet BC\'s at {x=}!')
 
+    def get_Neumann_BC_row(self, x):
+        """
+        Get a row for generating Neumann BCs at x with T polynomials.
+
+        Args:
+            x (float): Position of the boundary condition
+
+        Returns:
+            self.xp.ndarray: Row to put into a matrix
+        """
+        n = self.xp.arange(self.N, dtype='D')
+        nn = n**2
+        if x == -1:
+            me = nn
+            me[1:] *= (-1) ** n[:-1]
+            return me
+        elif x == 1:
+            return nn
+        else:
+            raise NotImplementedError(f'Don\'t know how to generate Neumann BC\'s at {x=}!')
+
     def get_Dirichlet_recombination_matrix(self):
         '''
         Get matrix for Dirichlet recombination, which changes the basis to have sparse boundary conditions.

pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py

@@ -419,5 +419,25 @@ def test_tau_method2D_diffusion(nz, nx, bc_val, plotting=False):
         ), f'Solution is incorrectly transformed back to real space at x={x[i]}'
 
 
-if __name__ == '__main__':
-    test_transform_cupy()
+@pytest.mark.base
+@pytest.mark.parametrize('N', [4, 7, 32])
+@pytest.mark.parametrize('x', [-1, 1])
+@pytest.mark.parametrize('v', [2, 3])
+def test_Neumann_BCs(N, x, v):
+    from pySDC.helpers.spectral_helper import ChebychevHelper
+    import numpy as np
+
+    helper = ChebychevHelper(N)
+
+    BC = helper.get_BC('Neumann', x=x)
+
+    grid = helper.get_1dgrid()
+    u = grid**v
+
+    u_hat = helper.transform(u)
+
+    value_at_BC = np.sum(u_hat * BC, axis=0)
+    if v % 2 == 1:
+        assert np.isclose(value_at_BC, v)
+    else:
+        assert np.isclose(value_at_BC, x * v)