Eliminating Nyquist mode for even resolution in RBC

pySDC/helpers/spectral_helper.py

@@ -568,9 +568,21 @@ def itransform(self, u, axis=-1):
     def get_BC(self, kind, **kwargs):
         if kind.lower() == 'integral':
             return self.get_integ_BC_row(**kwargs)
+        elif kind.lower() == 'nyquist':
+            assert (
+                self.N % 2 == 0
+            ), f'Do not eliminate the Nyquist mode with odd resolution as it is fully resolved. You chose {self.N} in this axis'
+            BC = self.xp.zeros(self.N)
+            BC[self.get_Nyquist_mode_index()] = 1
+            return BC
         else:
             return super().get_BC(kind, **kwargs)
 
+    def get_Nyquist_mode_index(self):
+        k = self.get_wavenumbers()
+        Nyquist_mode = min(k)
+        return self.xp.where(k == Nyquist_mode)[0][0]
+
     def get_integ_BC_row(self, **kwargs):
         """
         Only the 0-mode has non-zero integral with FFT basis in periodic BCs
@@ -861,6 +873,7 @@ def add_BC(self, component, equation, axis, kind, v, line=-1, scalar=False, **kw
             )
             N = self.axes[axis].N
             if (N + line) % N in self.xp.arange(N)[self.local_slice[axis]]:
+                slices[axis + 1] -= self.local_slice[axis].start
                 self.BC_rhs_mask[(*slices,)] = True
 
     def setup_BCs(self):
@@ -932,7 +945,10 @@ def put_BCs_in_rhs_hat(self, rhs_hat):
                 )
                 if axis == bc['axis']:
                     _slice = [self.index(bc['equation'])] + slices
-                    rhs_hat[(*_slice,)] = 0
+                    N = self.axes[axis].N
+                    if (N + bc['line']) % N in self.xp.arange(N)[self.local_slice[axis]]:
+                        _slice[axis + 1] -= self.local_slice[axis].start
+                        rhs_hat[(*_slice,)] = 0
 
         return rhs_hat + self.rhs_BCs_hat
 
@@ -963,7 +979,12 @@ def put_BCs_in_rhs(self, rhs):
                 )
                 if axis == bc['axis']:
                     _slice = [self.index(bc['equation'])] + slices
-                    _rhs_hat[(*_slice,)] = bc['v']
+
+                    N = self.axes[axis].N
+                    if (N + bc['line']) % N in self.xp.arange(N)[self.local_slice[axis]]:
+                        _slice[axis + 1] -= self.local_slice[axis].start
+
+                        _rhs_hat[(*_slice,)] = bc['v']
 
             rhs = self.itransform(_rhs_hat, axes=(axis - ndim,))
 

pySDC/implementations/problem_classes/RayleighBenard.py

@@ -128,12 +128,15 @@ def __init__(
             kind='Dirichlet',
             line=-1,
         )
-        self.setup_BCs()
 
+        # eliminate Nyquist mode if needed
         if nx % 2 == 0:
-            self.logger.warning(
-                f'The resolution is x-direction is even at {nx}. Keep in mind that the Nyquist mode is only partially resolved in this case. Consider changing the solution by one.'
-            )
+            Nyquist_mode_index = self.axes[0].get_Nyquist_mode_index()
+            for component in self.components:
+                self.add_BC(
+                    component=component, equation=component, axis=0, kind='Nyquist', line=Nyquist_mode_index, v=0
+                )
+        self.setup_BCs()
 
     def eval_f(self, u, *args, **kwargs):
         f = self.f_init

pySDC/tests/test_problems/test_RayleighBenard.py

@@ -249,6 +249,7 @@ def test_Poisson_problem_v():
     assert np.allclose(u_exact, u)
 
 
+@pytest.mark.mpi4py
 def test_CFL():
     from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard, CFLLimit
     import numpy as np
@@ -271,11 +272,26 @@ def test_CFL():
     assert np.allclose(dt2, 1 / u2[iv])
 
 
+@pytest.mark.mpi4py
+def test_Nyquist_mode_elimination():
+    from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard
+    import numpy as np
+
+    P = RayleighBenard(nx=32, nz=8)
+    u0 = P.u_exact(noise_level=1e-3)
+
+    u = P.solve_system(u0, dt=1e-3)
+
+    Nyquist_mode_index = P.axes[0].get_Nyquist_mode_index()
+    assert np.allclose(u[:, Nyquist_mode_index, :], 0)
+
+
 if __name__ == '__main__':
-    test_eval_f(2**0, 2**2, 'z', True)
+    # test_eval_f(2**0, 2**2, 'z', True)
     # test_Poisson_problem(1, 'T')
     # test_Poisson_problem_v()
     # test_Nusselt_numbers(1)
     # test_buoyancy_computation()
     # test_viscous_dissipation()
     # test_CFL()
+    test_Nyquist_mode_elimination()