Add configurations for $Ra=10^6$ to Rayleigh-Benard project (#597)

* Added Ra=1e6 configs to RBC project

* Added script for Nek5000 comparison

* linting

Author: brownbaerchen

pySDC/projects/RayleighBenard/RBC3D_configs.py

@@ -101,6 +101,7 @@ def get_initial_condition(self, P, *args, restart_idx=0, **kwargs):
             outfile = FieldsIO.fromFile(self.get_file_name())
 
             t0, solution = outfile.readField(restart_idx)
+            solution = solution[: P.spectral.ncomponents, ...]
 
             u0 = P.u_init
 
@@ -182,6 +183,7 @@ def get_initial_condition(self, P, *args, restart_idx=0, **kwargs):
         filename = ic_config.get_file_name()
         ic_file = FieldsIO.fromFile(filename)
         t0, ics = ic_file.readField(-1)
+        ics = ics[: P.spectral.ncomponents, ...]
         P.logger.info(f'Loaded initial conditions from {filename!r} at t={t0}.')
 
         # interpolate the initial conditions using padded transforms
@@ -265,6 +267,7 @@ def get_description(self, *args, res=-1, dt=-1, **kwargs):
         return desc
 
 
+# --- Ra 1e5 ---
 class RBC3DG4R4SDC22Ra1e5(RBC3DM2K2):
     Tend = 200
     dt = 4e-2
@@ -305,3 +308,36 @@ class RBC3DG4R4EulerRa1e5(RBC3DverificationEuler):
     dt = 8e-2
     res = 32
     converged = 50
+
+
+# --- Ra 1e6 ---
+class RBC3DG4R4SDC44Ra1e6(RBC3DM4K4):
+    Tend = 75
+    dt = 2e-2
+    res = 64
+    # converged = 22
+    ic_config = {'config': RBC3DG4R4SDC34Ra1e5, 'res': 32, 'dt': 0.02}
+
+
+class RBC3DG4R4SDC23Ra1e6(RBC3DM2K3):
+    Tend = 75
+    dt = 1e-2
+    res = 64
+    converged = 22
+    ic_config = {'config': RBC3DG4R4SDC34Ra1e5, 'res': 32, 'dt': 0.02}
+
+
+class RBC3DG4R4RKRa1e6(RBC3DverificationRK):
+    Tend = 75
+    dt = 1e-2
+    res = 64
+    ic_config = {'config': RBC3DG4R4SDC34Ra1e5, 'res': 32, 'dt': 0.02}
+    # converged = 22
+
+
+class RBC3DG4R4EulerRa1e6(RBC3DverificationEuler):
+    Tend = 75
+    dt = 1e-2
+    res = 64
+    ic_config = {'config': RBC3DG4R4SDC34Ra1e5, 'res': 32, 'dt': 0.02}
+    # converged = 22

pySDC/projects/RayleighBenard/analysis_scripts/RBC3D_spectrum.py

@@ -1,27 +1,32 @@
 from pySDC.projects.RayleighBenard.analysis_scripts.process_RBC3D_data import get_pySDC_data
 from pySDC.projects.RayleighBenard.analysis_scripts.plotting_utils import (
-    figsize_by_journal,
     get_plotting_style,
     savefig,
+    figsize,
 )
 import matplotlib.pyplot as plt
 
 
-def plot_spectrum(res, dt, config_name, ax):  # pragma: no cover
+def plot_spectrum(res, dt, config_name, ax, **plotting_params):  # pragma: no cover
     data = get_pySDC_data(res=res, dt=dt, config_name=config_name)
 
     spectrum = data['avg_spectrum']
     k = data['k']
-    ax.loglog(k[spectrum > 1e-16], spectrum[spectrum > 1e-16], **get_plotting_style(config_name), markevery=5)
+    ax.loglog(
+        k[spectrum > 1e-16],
+        spectrum[spectrum > 1e-16],
+        **{**get_plotting_style(config_name), **plotting_params},
+        markevery=5,
+    )
     ax.set_xlabel('$k$')
     ax.set_ylabel(r'$\|\hat{u}_x\|$')
 
 
 def plot_spectra_Ra1e5():  # pragma: no cover
-    fig, ax = plt.subplots(figsize=figsize_by_journal('Nature_CS', 1, 0.6))
+    fig, ax = plt.subplots(figsize=figsize(scale=1, ratio=0.6))
 
-    configs = [f'RBC3DG4R4{name}Ra1e5' for name in ['SDC34', 'SDC23', '', 'Euler', 'RK']]
-    dts = [0.06, 0.06, 0.06, 0.02, 0.04]
+    configs = [f'RBC3DG4R4{name}Ra1e5' for name in ['SDC23', 'SDC44', 'Euler', 'RK']]
+    dts = [0.06, 0.06, 0.02, 0.04]
     res = 32
 
     for config, dt in zip(configs, dts, strict=True):
@@ -31,7 +36,40 @@ def plot_spectra_Ra1e5():  # pragma: no cover
     savefig(fig, 'RBC3D_spectrum_Ra1e5')
 
 
+def plot_spectra_Ra1e6():  # pragma: no cover
+    fig, ax = plt.subplots(figsize=figsize(scale=1, ratio=0.6))
+
+    configs = [f'RBC3DG4R4{name}Ra1e6' for name in ['SDC44', 'SDC23', 'RK']]
+    dts = [0.01, 0.01, 0.01]
+    res = 64
+
+    for config, dt in zip(configs, dts, strict=True):
+        plot_spectrum(res, dt, config, ax)
+
+    ax.legend(frameon=False)
+    savefig(fig, 'RBC3D_spectrum_Ra1e6')
+
+
+def plot_all_spectra():  # pragma: no cover
+    fig, ax = plt.subplots(figsize=figsize(scale=1, ratio=0.6))
+
+    Ras = ['1e5', '1e6']
+    dts = [0.06, 0.01]
+    res = [32, 64]
+    colors = ['tab:blue', 'tab:orange']
+    markers = ['x', 'o']
+
+    for Ra, dt, _res, color, marker in zip(Ras, dts, res, colors, markers, strict=True):
+        config = f'RBC3DG4R4SDC44Ra{Ra}'
+        plot_spectrum(_res, dt, config, ax, label=f'$Ra$={Ra}', color=color, marker=marker)
+
+    ax.legend(frameon=False)
+    savefig(fig, 'RBC3D_all_spectra')
+
+
 if __name__ == '__main__':
-    plot_spectra_Ra1e5()
+    # plot_spectra_Ra1e5()
+    # plot_spectra_Ra1e6()
+    plot_all_spectra()
 
     plt.show()

pySDC/projects/RayleighBenard/analysis_scripts/compare_Nek5000.py

@@ -0,0 +1,167 @@
+"""
+This script plots macroscopic verification of simulation data via comparison with data from https://doi.org/10.5281/zenodo.14205874
+In order to run this script, you need to download the dataset and copy it to `pySDC/projects/RayleighBenard/data/Nek5000`.
+And, you need to have generated the pySDC simulation data, of course.
+"""
+
+import numpy as np
+from scipy import integrate
+import matplotlib.pyplot as plt
+
+ints = {'1e5': 1e5, '1e6': 1e6, '1e7': 1e7}
+
+
+def get_Nek5000_Data(Ra, base_path='data/Nek5000'):  # pragma: no cover
+    assert type(Ra) == str
+
+    # append relative path to base path
+    path = __file__
+    base_path = f'{path[::-1][path[::-1].index('/'):][::-1]}../{base_path}'
+    Pr = 0.7
+
+    if Ra == '1e5':
+        dir_name = '1_1e5'
+        start_time = 3500
+        nelZ = 64
+        nPoly = 5
+    elif Ra == '1e6':
+        dir_name = '2_1e6'
+        start_time = 3500
+        nelZ = 64
+        nPoly = 7
+    elif Ra == '1e7':
+        dir_name = '3_1e7'
+        start_time = 3100
+        nelZ = 64
+        nPoly = 9
+    elif Ra == '1e8':
+        dir_name = '4_1e8'
+        start_time = 3000
+        nelZ = 96
+        nPoly = 7
+    elif Ra == '1e9':
+        dir_name = '5_1e9'
+        start_time = 4000
+        nelZ = 96
+        nPoly = 9
+    elif Ra == '1e10':
+        dir_name = '6_1e10'
+        start_time = 1700
+        nelZ = 200
+        nPoly = 7
+    elif Ra == '1e11':
+        dir_name = '7_1e11'
+        start_time = 260
+        nelZ = 256
+        nPoly = 7
+    else:
+        raise
+
+    path = f'{base_path}/{dir_name}'
+    data = {}
+
+    visc = np.sqrt(Pr / ints[Ra])
+
+    # get averaged data
+    avg = np.load(f'{path}/average.npy')
+    avg_Nu = np.mean(avg[avg[:, 0] > start_time, 3])
+    data['Nu'] = avg_Nu
+    data['std_Nu'] = np.std(avg[avg[:, 0] > start_time, 3])
+
+    Re = np.sqrt(avg[:, 1]) / visc
+    data['Re'] = np.mean(Re[avg[:, 0] > start_time])
+
+    # get profile data
+    profiles = np.load(f'{path}/profile.npy')
+    nzPts = nelZ * nPoly + 1
+    nSnap = int(profiles.shape[0] / nzPts)
+    tVal = profiles[:, 0].reshape((nSnap, nzPts))[:, 0]
+    tInterval = tVal[-1] - tVal[0]
+
+    data['z'] = profiles[:nzPts, 1]
+    data['profile_T'] = profiles[:, 3].reshape((nSnap, nzPts))
+
+    tRMS = profiles[:, 4].reshape((nSnap, nzPts))
+    data['rms_profile_T'] = np.sqrt(integrate.simpson(tRMS**2, tVal, axis=0) / tInterval)
+
+    return data
+
+
+def get_pySDC_data(Ra):
+    from pySDC.projects.RayleighBenard.analysis_scripts.process_RBC3D_data import get_pySDC_data as _get_data
+
+    dts = {'1e5': 0.06, '1e6': 0.02}
+    res = {'1e5': 32, '1e6': 64}
+    return _get_data(config_name=f'RBC3DG4R4SDC34Ra{Ra}', dt=dts[Ra], res=res[Ra])
+
+
+def plot_Nu_scaling(ax):  # pragma: no cover
+
+    # reference values
+    for Ra in ['1e5', '1e6']:
+        dat = get_Nek5000_Data(Ra)
+        ax.errorbar(ints[Ra], dat['Nu'], yerr=dat['std_Nu'], fmt='o', color='black')
+
+    # pySDC values
+    for Ra in ['1e5', '1e6']:
+        dat = get_pySDC_data(Ra)
+        ax.errorbar(ints[Ra], dat['avg_Nu']['V'], yerr=dat['std_Nu']['V'], fmt='.', color='tab:blue')
+
+    ax.errorbar(None, None, fmt='o', color='black', label='Nek5000')
+    ax.errorbar(None, None, fmt='.', color='tab:blue', label='pySDC')
+    ax.legend(frameon=False, loc='lower right')
+
+    ax.set_xscale('log')
+    ax.set_xlabel('$Ra$')
+    ax.set_ylabel('$Nu$')
+
+
+def plot_T_profile(ax):  # pragma: no cover
+    colors = {'1e5': 'tab:blue', '1e6': 'tab:orange'}
+
+    # reference values
+    for Ra in ['1e5', '1e6']:
+        dat = get_Nek5000_Data(Ra)
+        ax.plot(dat['profile_T'].mean(axis=0), dat['z'], color=colors[Ra], label=f'Nek5000 Ra={Ra}')
+
+    # pySDC values
+    for Ra in ['1e5', '1e6']:
+        dat = get_pySDC_data(Ra)
+        ax.scatter(dat['profile_T'], dat['z'], color=colors[Ra], label=f'pySDC Ra={Ra}')
+
+    ax.set_ylabel('$z$')
+    ax.set_xlabel('$T$')
+    ax.legend()
+
+
+def plot_T_rms_profile(ax):  # pragma: no cover
+    colors = {'1e5': 'tab:blue', '1e6': 'tab:orange'}
+
+    # reference values
+    for Ra in ['1e5', '1e6']:
+        dat = get_Nek5000_Data(Ra)
+        ax.plot(dat['rms_profile_T'], dat['z'], color=colors[Ra], label=f'Nek5000 Ra={Ra}')
+
+    # pySDC values
+    for Ra in ['1e5', '1e6']:
+        dat = get_pySDC_data(Ra)
+        ax.scatter(dat['rms_profile_T'], dat['z'], color=colors[Ra], label=f'pySDC Ra={Ra}')
+
+    ax.set_ylabel('$z$')
+    ax.set_xlabel('$T$')
+    ax.legend()
+
+
+def plot_verification():  # pragma: no cover
+    from pySDC.projects.RayleighBenard.analysis_scripts.plotting_utils import savefig, figsize
+
+    fig, axs = plt.subplots(1, 2, figsize=figsize(scale=1, ratio=0.5))
+    plot_Nu_scaling(axs[0])
+    plot_T_profile(axs[1])
+    fig.tight_layout()
+    fig.savefig('plots/verification.pdf')
+
+
+if __name__ == '__main__':
+    plot_verification()
+    plt.show()