Fixed problem tests for 2D RBC

Author: brownbaerchen

pySDC/implementations/problem_classes/RayleighBenard.py

@@ -431,8 +431,7 @@ def compute_Nusselt_numbers(self, u):
         else:
             u_hat = self.transform(u)
 
-        DzT_hat = self.spectral.u_init_forward
-        DzT_hat[iT] = (self.Dz @ u_hat[iT].flatten()).reshape(DzT_hat[iT].shape)
+        DzT_hat = (self.Dz @ u_hat[iT].flatten()).reshape(u_hat[iT].shape)
 
         # compute vT with dealiasing
         padding = (self.dealiasing, self.dealiasing)
@@ -453,13 +452,13 @@ def compute_Nusselt_numbers(self, u):
         integral_V = 0
         if self.comm.rank == 0:
 
-            integral_z = (self._zInt @ nusselt_hat[0, 0]).real
+            integral_z = (self._zInt @ nusselt_hat[0]).real
             integral_z[0] = zAxis.get_integration_constant(integral_z, axis=-1)
             integral_V = ((top - bot) * integral_z).sum() * self.axes[0].L / self.nx
 
         Nusselt_V = self.comm.bcast(integral_V / self.spectral.V, root=0)
-        Nusselt_t = self.comm.bcast(self.xp.sum(nusselt_hat.real[0, 0] * top, axis=-1) / self.nx, root=0)
-        Nusselt_b = self.comm.bcast(self.xp.sum(nusselt_hat.real[0, 0] * bot, axis=-1) / self.nx, root=0)
+        Nusselt_t = self.comm.bcast(self.xp.sum(nusselt_hat.real[0] * top, axis=-1) / self.nx, root=0)
+        Nusselt_b = self.comm.bcast(self.xp.sum(nusselt_hat.real[0] * bot, axis=-1) / self.nx, root=0)
 
         return {
             'V': Nusselt_V,
@@ -536,8 +535,8 @@ def compute_max_step_size(P, u):
         grid_spacing_x = P.X[1, 0] - P.X[0, 0]
 
         cell_wallz = P.xp.zeros(P.nz + 1)
-        cell_wallz[0] = 1
-        cell_wallz[-1] = -1
+        cell_wallz[0] = P.Lz
+        cell_wallz[-1] = 0
         cell_wallz[1:-1] = (P.Z[0, :-1] + P.Z[0, 1:]) / 2
         grid_spacing_z = cell_wallz[:-1] - cell_wallz[1:]
 

pySDC/tests/test_problems/test_RayleighBenard.py

@@ -105,7 +105,7 @@ def test_vorticity(nx, nz, direction):
 
 @pytest.mark.mpi4py
 @pytest.mark.parametrize('v', [0, 3.14])
-def test_Nusselt_numbers(v, nx=6, nz=4):
+def test_Nusselt_numbers(v, nx=1, nz=10):
     import numpy as np
     from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard
 
@@ -114,12 +114,22 @@ def test_Nusselt_numbers(v, nx=6, nz=4):
         'v_bottom': 0,
     }
 
-    P = RayleighBenard(nx=nx, nz=nz, BCs=BCs, dealiasing=1.0, Rayleigh=1)
+    P = RayleighBenard(nx=nx, nz=nz, BCs=BCs, dealiasing=1.0, Rayleigh=1, spectral_space=False)
     prob = P
     xp = prob.xp
     iv, iT = prob.index(['v', 'T'])
 
-    u = P.u_exact(noise_level=0)
+    u = prob.u_init
+    u[iT, ...] = prob.Z
+    Nu = prob.compute_Nusselt_numbers(u)
+    for key, expect in zip(['t', 'b', 'V'], [-1, -1, -1]):
+        assert xp.isclose(Nu[key], expect), f'Expected Nu_{key}={expect}, but got {Nu[key]}'
+
+    u = prob.u_init
+    u[iT, ...] = 3 * prob.Z**2 + 1
+    Nu = prob.compute_Nusselt_numbers(u)
+    for key, expect in zip(['t', 'b', 'V'], [-6, 0, -3]):
+        assert xp.isclose(Nu[key], expect), f'Expected Nu_{key}={expect}, but got {Nu[key]}'
 
     u = prob.u_init
     u[iT, ...] = 3 * prob.Z**2 + 1
@@ -128,17 +138,6 @@ def test_Nusselt_numbers(v, nx=6, nz=4):
     for key, expect in zip(['t', 'b', 'V'], [prob.Lz * (3 + 1) * v - 6, v, v * (1 + 1) - 3]):
         assert xp.isclose(Nu[key], expect), f'Expected Nu_{key}={expect}, but got {Nu[key]}'
 
-    return None
-
-    # for key, expect in zip(['t', 'b', 'V'], [prob.Lz * (3 + 1) * w - 6, w, w * (1 + 1) - 3]):
-    #     assert xp.isclose(Nu[key], expect), f'Expected Nu_{key}={expect}, but got {Nu[key]}'
-
-    # nusselt = P.compute_Nusselt_numbers(u)
-    # expect = {'V':v, 't': 0, 'b': + 2 * v}
-    # print(nusselt, expect)
-    # for key in nusselt.keys():
-    #     assert np.isclose(nusselt[key], expect[key]), key
-
 
 def test_viscous_dissipation(nx=2**5 + 1, nz=2**3 + 1):
     import numpy as np
@@ -198,7 +197,7 @@ def test_Poisson_problems(nx, component):
         'T_bottom': 0,
     }
     P = RayleighBenard(
-        nx=nx, nz=6, BCs=BCs, Rayleigh=(max([abs(BCs['T_top'] - BCs['T_bottom']), np.finfo(float).eps]) * 2**3)
+        nx=nx, nz=6, BCs=BCs, Rayleigh=(max([abs(BCs['T_top'] - BCs['T_bottom']), np.finfo(float).eps]) * 2**3), Lz=2
     )
     rhs = P.u_init
 
@@ -237,7 +236,7 @@ def test_Poisson_problem_v():
         'v_top': 0,
         'v_bottom': 0,
         'T_top': 0,
-        'T_bottom': 2,
+        'T_bottom': 1,
     }
     P = RayleighBenard(nx=2, nz=2**3, BCs=BCs, Rayleigh=1.0)
     iv = P.index('v')
@@ -272,7 +271,7 @@ def test_CFL():
     from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard, CFLLimit
     import numpy as np
 
-    P = RayleighBenard(nx=5, nz=2, spectral_space=False)
+    P = RayleighBenard(nx=5, nz=2, spectral_space=False, Lz=2)
     iu, iv = P.index(['u', 'v'])
 
     u = P.u_init
@@ -331,11 +330,11 @@ def test_apply_BCs():
 
 if __name__ == '__main__':
     # test_eval_f(2**0, 2**2, 'z', True)
-    # test_Poisson_problem(1, 'T')
+    # test_Poisson_problems(1, 'T')
     # test_Poisson_problem_v()
     # test_apply_BCs()
-    test_Nusselt_numbers(1)
+    # test_Nusselt_numbers(1)
     # test_buoyancy_computation()
     # test_viscous_dissipation()
-    # test_CFL()
+    test_CFL()
     # test_Nyquist_mode_elimination()

pySDC/tests/test_problems/test_RayleighBenard3D.py

@@ -294,7 +294,7 @@ def test_heterogeneous_implementation():
 
 @pytest.mark.mpi4py
 @pytest.mark.parametrize('w', [0, 1, 3.14])
-def test_Nusselt_number_computation(w, N=4):
+def test_Nusselt_number_computation(w, N=6):
     from pySDC.implementations.problem_classes.RayleighBenard3D import RayleighBenard3D
 
     prob = RayleighBenard3D(nx=N, ny=N, nz=N, dealiasing=1.0, spectral_space=False)