TL: added script for stability investigation

Author: tlunet

pySDC/playgrounds/dedalus/.gitignore

@@ -11,3 +11,4 @@ scripts/run_*
 scripts/init_*
 scripts/runInit_*
 scripts/initRun_*
+scripts/stab_*

pySDC/playgrounds/dedalus/problems/rbc.py

@@ -1156,13 +1156,13 @@ def checkDNS(spectrum:np.ndarray, kappa:np.ndarray, sRatio:int=4, nThrow:int=0):
 
     .. math::
         \log(s_{tail}) \simeq
-        a\log(\kappa_{tail})^2 + b\log(\kappa_{tail}) + c
+        c_2\log(\kappa_{tail})^2 + c_1\log(\kappa_{tail}) + c_0
 
     where :math:`(\kappa_{tail},s_{tail})` is continuous subset of the
     mapping :math:`(\kappa,s)` for large values of :math:`\kappa`
     (i.e spectrum tail).
     If the quadratic regression produces a convex polynomial
-    (i.e :math:`a > 0`) then the simulation is considered as under-resolved
+    (i.e :math:`c_2 > 0`) then the simulation is considered as under-resolved
     (no DNS).
     Per default, the tail is built considering the
     **last quarter of the spectrum**.
@@ -1187,7 +1187,7 @@ def checkDNS(spectrum:np.ndarray, kappa:np.ndarray, sRatio:int=4, nThrow:int=0):
         Dictionnary containing the results, with keys :
 
         - `DNS` : boolean indicating if the simulation is well resolved
-        - `coeffs` : the :math:`a,b,c` regression coefficients, stored in a tuple
+        - `coeffs` : the :math:`c_2,c_1,c_0` regression coefficients, stored in a tuple
         - `kTail` : the :math:`\kappa_{tail}` values used for the regression
         - `sTail` : the :math:`s_{tail}` values used for the regression
 
@@ -1226,7 +1226,7 @@ def fun(coeffs):
     import matplotlib.pyplot as plt
 
     # dirName = "run_3D_A4_M0.5_R1_Ra1e6"
-    dirName = "run_3D_A4_M1_R2_Ra1e6"
+    dirName = "run_3D_A4_M1_R2_Ra1.5e5"
     # dirName = "run_M4_R2"
     # dirName = "test_M4_R2"
     OutputFiles.VERBOSE = True
@@ -1304,7 +1304,7 @@ def fun(coeffs):
 
         kappa = output.kappa
         plt.figure("spectrum")
-        for name in ["u"]:
+        for name in ["u", "uv", "uh", "b", "p"]:
             vals = spectrum[name]
             check = checkDNS(vals, kappa, sRatio=3)
             a, b, c = check["coeffs"]

pySDC/playgrounds/dedalus/scripts/checkDNS.py

@@ -21,13 +21,18 @@
     "run_3D_A4_M0.5_R1_Ra1e5",
     ]
 
-# simDirs = [
-#     "run_3D_A4_M1_R1_Ra5e4",
-#     "run_3D_A4_M1_R1_Ra1e5",
-#     "run_3D_A4_M1_R1_Ra1.5e5",
-#     "run_3D_A4_M1_R1_Ra2e5",
-#     "run_3D_A4_M1_R1_Ra1e6",
-#     ]
+simDirs = [
+    "run_3D_A4_M1_R1_Ra5e4",
+    "run_3D_A4_M1_R1_Ra1e5",
+    "run_3D_A4_M1_R1_Ra1.5e5",
+    "run_3D_A4_M1_R1_Ra2e5",
+    "run_3D_A4_M1_R1_Ra1e6",
+    ]
+
+simDirs = [
+    "run_3D_A4_M1_R2_Ra1.5e5",
+    "run_3D_A4_M1_R2_Ra1e6",
+    ]
 
 df = pd.DataFrame(
     columns=["Ra", "c_2[u]", "c_2[uv]", "c_2[uh]", "c_2[b]", "c_2[p]"])
@@ -38,10 +43,11 @@
     df.loc[i, "Ra"] = float(dirName.split("_Ra")[-1])
 
     # assert len(output.times) == 61, f"not 61 fields for {dirName}"
-    if len(output.times) == 61:
-        start = 20
-    else:
-        start = 60
+    # if len(output.times) == 61:
+    #     start = 20
+    # else:
+    #     start = 60
+    start = 20
 
     spectrum = output.getSpectrum(
             which="all", zVal="all",

pySDC/playgrounds/dedalus/scripts/stab_rbc.py

@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Sep 30 13:02:50 2025
+
+@author: cpf5546
+"""
+import os
+import matplotlib.pyplot as plt
+
+from pySDC.playgrounds.dedalus.timestepper import SDCIMEX
+from pySDC.playgrounds.dedalus.problems.rbc import RBCProblem3D, OutputFiles
+
+
+dtBase = 2.50e-03
+nSteps = 10
+nCycles = 2
+
+Rayleigh = 1.5e5
+timeScheme = "SDC"
+stabDir = f"stab_A4_M1_R1_{timeScheme}"
+initSol = "init_3D_A4_M1_R1.pySDC"
+
+SDCIMEX.setParameters(
+    nNodes=4, nodeType="LEGENDRE", quadType="RADAU-RIGHT",
+    nSweeps=4, initSweep="COPY", explSweep="PIC", implSweep="MIN-SR-FLEX",
+    )
+
+os.makedirs(stabDir, exist_ok=True)
+
+dtFactors = [f*10**e for e in range(nCycles) for f in [1, 2, 5]]
+fmtSuffix = f":0{nCycles}d"
+
+plt.figure("spectrum")
+for fac in dtFactors:
+    dtRun = dtBase*fac
+    runDir = f"{stabDir}/dt_f" + ("{"+fmtSuffix+"}").format(fac)
+
+    prob = RBCProblem3D.runSimulation(
+        runDir, dtRun*nSteps, dtRun, timeScheme=timeScheme,
+        dtWrite=dtRun*nSteps, initField=initSol,
+        aspectRatio=4, meshRatio=1, resFactor=1,
+        Rayleigh=Rayleigh, Prandtl=0.7,
+        )
+
+    if "tComp" in prob.infos:
+        break
+
+    output = OutputFiles(runDir)
+    spectrum = output.getSpectrum(which="all", start=1)
+
+#     plt.loglog(spectrum["kappa"], spectrum["p"], label="f"+("{"+fmtSuffix+"}").format(fac))
+
+# plt.legend()