Merge pull request #207 from snek5000/doc/check-from-py

Improve tuto cbox

docs/examples/scripts/tuto_cbox.py

@@ -1,13 +1,21 @@
+import builtins
+import os
+import signal
+from time import perf_counter, sleep
+
 import numpy as np
+from scipy.signal import argrelmax
 from snek5000_cbox.solver import Simul
 
+from fluiddyn.util import time_as_str
+
 params = Simul.create_default_params()
 
 aspect_ratio = 1.0
 params.prandtl = 0.71
 
 # The onset of oscillatory flow for aspect ratio 1.0 is at Ra_c = 1.825e8
-params.Ra_side = 1.94e8
+params.Ra_side = 2e8
 
 params.output.sub_directory = "examples_snek/tuto"
 
@@ -47,7 +55,7 @@
 params.nek.velocity.residual_tol = 1e-08
 params.nek.pressure.residual_tol = 1e-05
 
-params.nek.general.end_time = 600
+params.nek.general.end_time = 1000
 params.nek.general.stop_at = "endTime"
 params.nek.general.target_cfl = 2.0
 params.nek.general.time_stepper = "BDF3"
@@ -58,6 +66,87 @@
 
 params.output.history_points.write_interval = 10
 
+
 if __name__ == "__main__":
     sim = Simul(params)
-    sim.make.exec("run_fg", nproc=2)
+
+    # to save the PID of the simulation (a number associated with the process)
+    sim.output.write_snakemake_config(
+        custom_env_vars={"MPIEXEC_FLAGS": "--report-pid PID.txt"}
+    )
+    # run the simulation in the background (non blocking call)
+    sim.make.exec("run", nproc=2)
+
+    pid_file = sim.path_run / "PID.txt"
+
+    # we wait for the simulation to be started and the PID file to be created
+    n = 0
+    while not pid_file.exists():
+        sleep(1)
+        n += 1
+        if n > 60:
+            raise RuntimeError(f"{pid_file} does not exist.")
+
+    with open(pid_file) as file:
+        pid = int(file.read().strip())
+
+    # we know the PID of the simulation so we can control it!
+
+    # a simple way to print in stdout and in a file
+    path_log_py = sim.path_run / f"log_py_{time_as_str()}.txt"
+
+    def print(*args, sep=" ", end="\n", **kwargs):
+        builtins.print(*args, **kwargs)
+        with open(path_log_py, "a") as file:
+            file.write(sep.join(str(arg) for arg in args) + end)
+
+    print(f"{pid = }")
+
+    def check_running():
+        """Check for the existence of a running process"""
+        try:
+            os.kill(pid, 0)
+        except OSError:
+            return False
+        else:
+            return True
+
+    while check_running():
+        sleep(2)
+        t0 = perf_counter()
+        coords, df = sim.output.history_points.load_1point(
+            index_point=12, key="temperature"
+        )
+        t_last = df.time.max()
+        print(
+            f"{time_as_str()}, {t_last = :.2f}: "
+            f"history_points loaded in {perf_counter() - t0:.2f} s"
+        )
+
+        if t_last < 500:
+            continue
+
+        temperature = df.temperature.to_numpy()
+        times = df.time.to_numpy()
+
+        temperature = temperature[times > 500]
+        times = times[times > 500]
+
+        # we know that there are oscillations growing exponentially
+        # we look for the positive maxima of the signal
+        indices_maxima = argrelmax(temperature)
+        temp_maxima = temperature[indices_maxima]
+        temp_maxima = temp_maxima[temp_maxima > 0]
+
+        # similar to a second order derivative
+        diff_maxima = np.diff(np.diff(temp_maxima))
+
+        if any(diff_maxima < 0):
+            # as soon as the growth starts to saturate
+            print(
+                f"Saturation of the instability detected at t = {t_last}\n"
+                "Terminate simulation."
+            )
+
+            os.kill(pid, signal.SIGTERM)
+            break

docs/examples/scripts/tuto_phill.py

@@ -15,7 +15,7 @@
 params.nek.general.time_stepper = "bdf2"
 params.nek.general.write_interval = 10
 
-params.nek.stat.av_step = 3
+params.nek.stat.av_step = 2
 params.nek.stat.io_step = 10
 
 sim = Simul(params)

docs/tuto_cbox.md

@@ -23,13 +23,22 @@ kernelspec:
 This example is based on
 [this study](https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/from-onset-of-unsteadiness-to-chaos-in-a-differentially-heated-square-cavity/617F4CB2C23DD74C3D0CB872AE7C0045).
 The configuration is a square cavity. The control parameters are Prandtl $= 0.71$ and
-Rayleigh $= 1.86 \times 10^{8}$. The mesh size is $64 \times 64$. We want to have $25$
+Rayleigh $= 2 \times 10^{8}$. The mesh size is $64 \times 64$. We want to have $25$
 probes (history points) to record the variable signals. We will use these probe signals
 in monitoring and postprocessing of the simulation. See
 [this example](https://github.com/snek5000/snek5000-cbox/blob/gh-actions/doc/examples/run_side_short.py)
-for the implementation. The simulation will be carried out with the script
-[docs/examples/scripts/tuto_cbox.py](https://github.com/snek5000/snek5000/tree/main/docs/examples/scripts/tuto_cbox.py),
-which contains:
+for the implementation.
+
+The simulation will be carried out with the script
+[docs/examples/scripts/tuto_cbox.py](https://github.com/snek5000/snek5000/tree/main/docs/examples/scripts/tuto_cbox.py).
+Note that this script is more complicated than for the previous tutorial. Here,
+we want to demonstrate that it is possible to check what happen in the
+simulation from Python and to stop the simulation depending on its outputs. We
+know that for moderate Rayleigh number, the side wall convection in a box first
+reach a quasi-steady state from which emerges an oscillatory instability. Here,
+we want to stop the simulation as soon as the linear instability starts to
+saturate, i.e. as soon as the growth of the unstable mode becomes slower than
+exponential.
 
 ```{eval-rst}
 .. literalinclude:: ./examples/scripts/tuto_cbox.py
@@ -55,7 +64,31 @@ print(f"Script executed in {perf_counter() - t_start:.2f} s")
 The script has now been executed. Let's look at its output.
 
 ```{code-cell} ipython3
-lines = process.stdout.split("\n")
+print(process.stdout)
+```
+
+To "load the simulation", i.e. to recreate a simulation object, we now need to extract
+from the output the path of the directory of the simulation:
+
+```{code-cell} ipython3
+path_run = None
+for line in process.stdout.split("\n"):
+    if "path_run: " in line:
+        path_run = line.split("path_run: ")[1].split(" ", 1)[0]
+        break
+if path_run is None:
+    raise RuntimeError
+path_run
+```
+
+We can now read the Nek5000 log file.
+
+```{code-cell} ipython3
+from pathlib import Path
+
+path_log = Path(path_run) / "cbox.log"
+lines = path_log.read_text().split("\n")
+
 index_step2 = 0
 for line in lines:
     if line.startswith("Step      2, t= "):
@@ -73,20 +106,6 @@ for line in lines[::-1]:
 print("\n".join(lines[index_final_step-10:]))
 ```
 
-To "load the simulation", i.e. to recreate a simulation object, we now need to extract
-from the output the path of the directory of the simulation:
-
-```{code-cell} ipython3
-path_run = None
-for line in process.stdout.split("\n"):
-    if "path_run: " in line:
-        path_run = line.split("path_run: ")[1].split(" ", 1)[0]
-        break
-if path_run is None:
-    raise RuntimeError
-path_run
-```
-
 ## Postprocessing
 
 We can load the simulation:
@@ -102,52 +121,60 @@ The command `snek-ipy-load` can be used to start a IPython session and load the
 simulation saved in the current directory.
 ```
 
-Then we are able to plot all the history points for one variable like $u_x$,
+Then we are able to plot all the history points for one variable (here the
+temperature),
 
 ```{code-cell} ipython3
-sim.output.history_points.plot(key='ux');
+sim.output.history_points.plot(key='temperature');
 ```
 
-or just one history point:
-
 ```{code-cell} ipython3
-sim.output.history_points.plot_1point(
-  index_point=0, key='temperature', tmin=400, tmax=800
-);
+ax = sim.output.history_points.plot(key='temperature')
+ax.set_xlim(left=300)
+ax.set_ylim([0.2, 0.36]);
 ```
 
-Also we can load the history points data to compute growth rate:
+or just one history point:
 
 ```{code-cell} ipython3
+ax = sim.output.history_points.plot_1point(
+  index_point=12, key='temperature', tmin=300, tmax=800
+)
+
+coords, df = sim.output.history_points.load()
+
 import numpy as np
 from scipy import stats
 from scipy.signal import argrelmax
 
-coords, df = sim.output.history_points.load()
-```
-
-```{code-cell} ipython3
 df_point = df[df.index_points == 12]
-time = df_point["time"].to_numpy()
-ux = df_point["ux"].to_numpy()
+times = df_point["time"].to_numpy()
+signal = df_point["temperature"].to_numpy()
+
+cond = times > 400
+times = times[cond]
+signal = signal[cond]
+
+indices_maxima = argrelmax(signal)
+times_maxima = times[indices_maxima]
+signal_maxima = signal[indices_maxima]
 
-indx = np.where(time > 450)[0][0]
-time = time[indx:]
-ux = ux[indx:]
-signal = ux
+ax.plot(times_maxima, signal_maxima, "xr")
+```
 
-arg_local_max = argrelmax(signal)
-time_local_max = time[arg_local_max]
-signal_local_max = signal[arg_local_max]
+Also we can also compute an approximation of the growth rate:
 
+```{code-cell} ipython3
 slope, intercept, r_value, p_value, std_err = stats.linregress(
-    time_local_max, np.log(signal_local_max)
+    times_maxima, np.log(abs(signal_maxima))
 )
 
 growth_rate = slope
 print(f"The growth rate is {growth_rate:.2e}")
 ```
 
+## Load the flow field as xarray dataset
+
 There is also the possibility to load to whole field file in
 [xarray dataset](https://docs.xarray.dev/en/stable/index.html)
 

src/snek5000/output/history_points.py

@@ -98,7 +98,6 @@ def load(self):
             nb_data_lines_read = len(self._data)
 
             size_file = self.path_file.stat().st_size
-            print(f"  {size_file = }")
 
             with open(self.path_file) as file:
                 first_line = file.readline()
@@ -114,9 +113,6 @@ def load(self):
                 + nb_data_lines_read * len(line_data)
             )
             nb_chars_not_read = size_file - nb_chars_read
-
-            print(f"  {nb_chars_not_read = }")
-
             if nb_chars_not_read == 0:
                 return self.coords, self._data