Changed default 2D RBC setup

Author: brownbaerchen

pySDC/implementations/problem_classes/RayleighBenard.py

@@ -58,7 +58,8 @@ def __init__(
         BCs=None,
         dealiasing=3 / 2,
         comm=None,
-        Lx=8,
+        Lx=4,
+        Lz=1,
         **kwargs,
     ):
         """
@@ -77,7 +78,7 @@ def __init__(
         BCs = {} if BCs is None else BCs
         BCs = {
             'T_top': 0,
-            'T_bottom': 2,
+            'T_bottom': 1,
             'v_top': 0,
             'v_bottom': 0,
             'u_top': 0,
@@ -101,11 +102,15 @@ def __init__(
             'dealiasing',
             'comm',
             'Lx',
+            'Lz',
             localVars=locals(),
             readOnly=True,
         )
 
-        bases = [{'base': 'fft', 'N': nx, 'x0': 0, 'x1': self.Lx}, {'base': 'ultraspherical', 'N': nz}]
+        bases = [
+            {'base': 'fft', 'N': nx, 'x0': 0, 'x1': self.Lx},
+            {'base': 'ultraspherical', 'N': nz, 'x0': 0, 'x1': self.Lz},
+        ]
         components = ['u', 'v', 'T', 'p']
         super().__init__(bases, components, comm=comm, **kwargs)
 
@@ -370,12 +375,13 @@ def compute_vorticity(self, u):
             u_hat = u.copy()
         else:
             u_hat = self.transform(u)
+
         Dz = self.Dz
         Dx = self.Dx
         iu, iv = self.index(['u', 'v'])
 
         vorticity_hat = self.spectral.u_init_forward
-        vorticity_hat[0] = (Dx * u_hat[iv].flatten() + Dz @ u_hat[iu].flatten()).reshape(u[iu].shape)
+        vorticity_hat[0] = (Dx * u_hat[iv].flatten() + Dz @ u_hat[iu].flatten()).reshape(u_hat[iu].shape)
         return self.itransform(vorticity_hat)[0].real
 
     def getOutputFile(self, fileName):
@@ -418,14 +424,14 @@ def compute_Nusselt_numbers(self, u):
             dict: Nusselt number averaged over the entire volume and horizontally averaged at the top and bottom.
         """
         iv, iT = self.index(['v', 'T'])
-
-        DzT_hat = self.spectral.u_init_forward
+        zAxis = self.spectral.axes[-1]
 
         if self.spectral_space:
             u_hat = u.copy()
         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)
 
         # compute vT with dealiasing
@@ -435,21 +441,25 @@ def compute_Nusselt_numbers(self, u):
         _me[0] = u_pad[iv] * u_pad[iT]
         vT_hat = self.transform(_me, padding=padding)[0]
 
-        # nusselt_hat = (vT_hat - DzT_hat[iT]) / self.nx
+        if not hasattr(self, '_zInt'):
+            self._zInt = zAxis.get_integration_matrix()
+
         nusselt_hat = (vT_hat / self.kappa - DzT_hat) * self.axes[-1].L
 
-        integral_z = self.xp.sum(nusselt_hat * self.spectral.axes[1].get_BC(kind='integral'), axis=-1).real
-        integral_V = (
-            integral_z[0] * self.axes[0].L
-        )  # only the first Fourier mode has non-zero integral with periodic BCs
-        Nusselt_V = self.comm.bcast(integral_V / self.spectral.V, root=0)
+        # get coefficients for evaluation on the boundary
+        top = zAxis.get_BC(kind='Dirichlet', x=1)
+        bot = zAxis.get_BC(kind='Dirichlet', x=-1)
 
-        Nusselt_t = self.comm.bcast(
-            self.xp.sum(nusselt_hat * self.spectral.axes[1].get_BC(kind='Dirichlet', x=1), axis=-1).real[0], root=0
-        )
-        Nusselt_b = self.comm.bcast(
-            self.xp.sum(nusselt_hat * self.spectral.axes[1].get_BC(kind='Dirichlet', x=-1), axis=-1).real[0], root=0
-        )
+        integral_V = 0
+        if self.comm.rank == 0:
+
+            integral_z = (self._zInt @ nusselt_hat[0, 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)
 
         return {
             'V': Nusselt_V,

pySDC/tests/test_hooks/test_log_to_file.py

@@ -126,7 +126,10 @@ def test_logging(tmpdir, use_pickle, ODE=True):
 
     for us, uf in zip(u, u_file):
         assert us[0] == uf[0], 'time does not match'
-        assert np.allclose(us[1], uf[1]), 'solution does not match'
+        if ODE:
+            assert np.allclose(us[1], uf[1]), 'solution does not match'
+        else:
+            assert np.allclose(us[1], uf[1][:4]), 'solution does not match'
 
 
 @pytest.mark.base

pySDC/tests/test_problems/test_RayleighBenard.py

@@ -110,18 +110,34 @@ def test_Nusselt_numbers(v, nx=6, nz=4):
     from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard
 
     BCs = {
-        'v_top': v,
-        'v_bottom': v,
+        'v_top': 0,
+        'v_bottom': 0,
     }
 
-    P = RayleighBenard(nx=nx, nz=nz, BCs=BCs, dealiasing=1.5)
+    P = RayleighBenard(nx=nx, nz=nz, BCs=BCs, dealiasing=1.0, Rayleigh=1)
+    prob = P
+    xp = prob.xp
+    iv, iT = prob.index(['v', 'T'])
 
     u = P.u_exact(noise_level=0)
 
-    nusselt = P.compute_Nusselt_numbers(u)
-    expect = {'V': 1 + v, 't': 1, 'b': +1 + 2 * v, 'b_no_v': 1, 't_no_v': 1}
-    for key in nusselt.keys():
-        assert np.isclose(nusselt[key], expect[key]), key
+    u = prob.u_init
+    u[iT, ...] = 3 * prob.Z**2 + 1
+    u[iv] = v * (1 + xp.sin(prob.X / prob.axes[0].L * 2 * xp.pi))
+    Nu = prob.compute_Nusselt_numbers(u)
+    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):
@@ -318,7 +334,7 @@ def test_apply_BCs():
     # test_Poisson_problem(1, 'T')
     # test_Poisson_problem_v()
     # test_apply_BCs()
-    test_Nusselt_numbers(3.14)
+    test_Nusselt_numbers(1)
     # test_buoyancy_computation()
     # test_viscous_dissipation()
     # test_CFL()