Added spectrum computation to 3D RBC (Parallel-in-Time#581)

* Added spectrum computation to 3D RBC

* Added documentation

Author: brownbaerchen

pySDC/implementations/problem_classes/RayleighBenard3D.py

@@ -28,9 +28,9 @@ class RayleighBenard3D(GenericSpectralLinear):
 
     The domain, vertical boundary conditions and pressure gauge are
 
-        Omega = [0, 8) x (-1, 1)
+        Omega = [0, Lx) x [0, Ly] x (0, Lz)
         T(z=+1) = 0
-        T(z=-1) = 2
+        T(z=-1) = Lz
         u(z=+-1) = v(z=+-1) = 0
         integral over p = 0
 
@@ -53,16 +53,16 @@ class RayleighBenard3D(GenericSpectralLinear):
     def __init__(
         self,
         Prandtl=1,
-        Rayleigh=2e6,
+        Rayleigh=1e6,
         nx=64,
         ny=64,
-        nz=64,
+        nz=32,
         BCs=None,
         dealiasing=1.5,
         comm=None,
         Lz=1,
-        Lx=1,
-        Ly=1,
+        Lx=4,
+        Ly=4,
         useGPU=False,
         **kwargs,
     ):
@@ -375,24 +375,22 @@ def get_vertical_profiles(self, u, components):
 
         return avgs
 
-    def get_Reynolds_number(self, u):
-        raise  # compute RMS velocity incorrectly
-        if self.spectral_space:
-            u = self.itransform(u)
-        else:
-            u = u.copy()
+    def get_frequency_spectrum(self, u):
+        """
+        Compute the frequency spectrum of the velocities in x and y direction in the horizontal plane for every point in
+        z. If the problem is well resolved, the coefficients will decay quickly with the wave number, and the reverse
+        indicates that the resolution is too low.
 
-        indices = list(self.index(['u', 'v', 'w']))
-        vels = u[indices, ...]
-        vels *= vels
-        vels_mean = self.xp.sum(vels, axis=(1, 2, 3)) / self.xp.prod(self.spectral.global_shape[1:])
-        vels_mean = self.comm.allreduce(vels_mean, MPI.SUM)
-        v_RMS = (self.xp.sum(vels_mean)) ** (1 / 2)
+        The returned spectrum has three dimensions. The first is for component (i.e. u or v), the second is for every
+        point in z and the third is the energy in every wave number.
 
-        Re = v_RMS * self.axes[-1].L / self.nu
-        return Re
+        Args:
+            u: The solution you want to compute the spectrum of
 
-    def get_frequency_spectrum(self, u):
+        Returns:
+            RayleighBenard3D.xp.ndarray: wave numbers
+            RayleighBenard3D.xp.ndarray: spectrum
+        """
         xp = self.xp
         indices = slice(0, 2)
 

pySDC/tests/test_problems/test_RayleighBenard3D.py

@@ -11,7 +11,7 @@ def test_eval_f(nx, nz, direction, spectral_space):
     import numpy as np
     from pySDC.implementations.problem_classes.RayleighBenard3D import RayleighBenard3D
 
-    P = RayleighBenard3D(nx=nx, ny=nx, nz=nz, Rayleigh=1, spectral_space=spectral_space)
+    P = RayleighBenard3D(nx=nx, ny=nx, nz=nz, Rayleigh=1, spectral_space=spectral_space, Lx=1, Ly=1, Lz=1)
     iu, iv, iw, ip, iT = P.index(['u', 'v', 'w', 'p', 'T'])
     X, Y, Z = P.X, P.Y, P.Z
     cos, sin = np.cos, np.sin
@@ -297,7 +297,7 @@ def test_heterogeneous_implementation(N=8, useGPU=True):
 def test_Nusselt_number_computation(w, N=4):
     from pySDC.implementations.problem_classes.RayleighBenard3D import RayleighBenard3D
 
-    prob = RayleighBenard3D(nx=N, ny=N, nz=N, dealiasing=1.0, spectral_space=False, Rayleigh=1.0)
+    prob = RayleighBenard3D(nx=N, ny=N, nz=N, dealiasing=1.0, spectral_space=False, Rayleigh=1.0, Prandtl=1.0)
     xp = prob.xp
     iw, iT = prob.index(['w', 'T'])
 
@@ -341,36 +341,6 @@ def test_Nusselt_number_computation(w, N=4):
         assert xp.isclose(Nu[key], w), f'Expected Nu_{key}={w}, but got {Nu[key]} with constant T and perturbed w!'
 
 
-@pytest.mark.mpi4py
-@pytest.mark.mpi(ranks=[1, 4])
-def test_Reynolds_number_computation(mpi_ranks):
-    from pySDC.implementations.problem_classes.RayleighBenard3D import RayleighBenard3D
-
-    N = 4
-    prob = RayleighBenard3D(nx=N, ny=N, nz=N, dealiasing=1.0, spectral_space=False, Rayleigh=1.0)
-    xp = prob.xp
-    iu, iv, iw = prob.index(['u', 'v', 'w'])
-
-    u = prob.u_exact()
-    u[iu] = 1
-    u[iv] = 1
-    u[iw] = 1
-    Re = prob.get_Reynolds_number(u)
-    assert xp.isclose(Re, xp.sqrt(3))
-
-    u = prob.u_exact()
-    u[iv] = (2 * prob.Z) ** (1 / 2)
-    Re = prob.get_Reynolds_number(u)
-    assert xp.isclose(Re, 1)
-
-    u = prob.u_exact()
-    u[iv] = (2 * prob.Z) ** (1 / 2)
-    u[iu] = (2 * prob.Z) ** (1 / 2)
-    u[iw] = (2 * prob.Z) ** (1 / 2)
-    Re = prob.get_Reynolds_number(u)
-    assert xp.isclose(Re, xp.sqrt(3))
-
-
 @pytest.mark.mpi4py
 def test_CFL():
     from pySDC.implementations.problem_classes.RayleighBenard3D import RayleighBenard3D
@@ -452,7 +422,7 @@ def test_spectrum_computation(mpi_ranks):
 
 
 if __name__ == '__main__':
-    # test_eval_f(2**2, 2**1, 'x', False)
+    test_eval_f(2**2, 2**1, 'x', False)
     # test_libraries()
     # test_Poisson_problems(4, 'u')
     # test_Poisson_problem_w()
@@ -461,5 +431,4 @@ def test_spectrum_computation(mpi_ranks):
     # test_heterogeneous_implementation()
     # test_Nusselt_number_computation(N=4, w=3.14)
     test_spectrum_computation(None)
-    # test_Reynolds_number_computation(None)
     # test_CFL()