Repository: snek5000/snek5000 (language: python)

Issue:
Plot methods in snek5000/output/print_stdout.py
Nek5000 logs are quite complicated (and for me mysterious). I get something like:

```
Starting time loop ...

     DT/DTCFL/DTFS/DTINIT   0.100E-01   0.100E-01   0.000E+00   0.100E-01
Step      1, t= 1.0000000E-02, DT= 1.0000000E-02, C=  0.000 0.0000E+00 0.0000E+00
             Solving for Hmholtz scalars
  Temperature/Passive scalar solution
 New CG1-tolerance (RINIT*epsm) =    1.1153045821067478E-014   0.0000000000000000     
          1  Hmholtz TEMP       1   1.1153E-01   7.3922E-05   1.0000E-07   1.0000E+02   F
          1  Hmholtz TEMP       2   1.7017E-03   7.3922E-05   1.0000E-07   1.0000E+02   F
          1  Hmholtz TEMP       3   3.7470E-05   7.3922E-05   1.0000E-07   1.0000E+02   F
          1  Hmholtz TEMP       4   1.1114E-06   7.3922E-05   1.0000E-07   1.0000E+02   F
          1  Hmholtz TEMP       5   1.9912E-08   7.3922E-05   1.0000E-07   1.0000E+02   F
          1  Hmholtz TEMP       4   1.9912E-08   1.1153E-01   1.0000E-07
          1  Scalars done  1.0000E-02  1.0160E-03
             Solving for fluid
          1  PRES gmres         0   0.0000E+00   0.0000E+00   1.0000E-05   0.0000E+00   8.1685E-05    F
          1  Helmh3 fluid       0   0.0000E+00   0.0000E+00   1.0000E-07
             L1/L2 DIV(V)           0.0000E+00   0.0000E+00
             L1/L2 QTL              0.0000E+00   0.0000E+00
             L1/L2 DIV(V)-QTL       0.0000E+00   0.0000E+00
          1  Fluid done  1.0000E-02  7.0113E-04
Step      2, t= 2.0000000E-02, DT= 1.0000000E-02, C=  0.000 3.4612E-03 3.4612E-03
             Solving for Hmholtz scalars
          2  Hmholtz TEMP       4   6.1961E-09   1.6318E-01   1.0000E-07
          2  Scalars done  2.0000E-02  3.2448E-04
             Solving for fluid
          2  PRES gmres         6   7.4452E-06   2.4753E-03   1.0000E-05   6.4027E-04   9.9262E-04    F
          2  Helmh3 fluid       3   3.3066E-08   4.4121E-05   1.0000E-07
             L1/L2 DIV(V)          -6.4656E-20   1.5455E-03
             L1/L2 QTL              0.0000E+00   0.0000E+00
             L1/L2 DIV(V)-QTL      -6.4656E-20   1.5455E-03
          2  Fluid done  2.0000E-02  1.9669E-03
Step      3, t= 3.0000000E-02, DT= 1.0000000E-02, C=  0.001 6.6836E-03 3.2224E-03
             Solving for Hmholtz scalars
          3  Hmholtz TEMP       4   3.4024E-09   1.9467E-01   1.0000E-07
          3  Scalars done  3.0000E-02  3.1270E-04
             Solving for fluid
          3  PRES gmres         7   4.5838E-06   3.9343E-03   1.0000E-05   6.8044E-04   1.0463E-03    F
          3  Helmh3 fluid       3   1.9763E-08   4.4491E-05   1.0000E-07
             L1/L2 DIV(V)          -1.3031E-19   7.3944E-04
             L1/L2 QTL              0.0000E+00   0.0000E+00
             L1/L2 DIV(V)-QTL      -1.3031E-19   7.3944E-04
          3  Fluid done  3.0000E-02  1.8685E-03
[...]
```

It would be great to also print the clock time for each time step so that we could get figures as `plot_clock_times` in https://foss.heptapod.net/fluiddyn/fluidsim/-/blob/branch/default/fluidsim/base/output/print_stdout.py#L222. It is really very useful.

I guess that other plots representing the data printed by Nek5000 could also be useful.


Discussion / hints:
Surprisingly, I can't find documentation on what is printed by Nek5000 :slightly_frowning_face: 
We could look at https://github.com/Nek5000/Nek5000/blob/master/core/drive2.f
Good find @akhoubani. To be more precise `TTIME_STP` seems to be what you want:

https://github.com/Nek5000/Nek5000/blob/71669e62586c7786a936a29306e4f37a0b9c1490/core/drive2.f#L127-L136

I think plotting number of pressure iterations (`iter`) and divergence (`divex`) is also a useful diagnostic - something's wrong if it goes up. There are different pressure solvers, but the default `gmres` should do for now:

https://github.com/Nek5000/Nek5000/blob/3e3eb9ac9b0099ddb1f385be04e3784cde4a58a8/core/gmres.f#L537-L540

Fix patch:
diff --git a/src/snek5000/output/print_stdout.py b/src/snek5000/output/print_stdout.py
index f1ccde60..806ff9b6 100644
--- a/src/snek5000/output/print_stdout.py
+++ b/src/snek5000/output/print_stdout.py
@@ -3,12 +3,14 @@
 """
 import re
 from pathlib import Path
+from warnings import warn
 
-import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
 
+from fluiddyn.util.util import modification_date
 from snek5000 import mpi
 from snek5000.log import logger
-from snek5000.solvers.base import Simul
 
 # Uses awkupy
 # %load_ext iawk
@@ -19,96 +21,172 @@ class PrintStdOut:
 
     _tag = "print_stdout"
 
-    def __init__(self, output=None, file=None):
+    def __init__(self, output=None, path_file=None):
         self.output = output
-        self._file = file
+        self._path_file = path_file
+        self.data = None
+        self._data_modif_time = None
 
     @property
     def file(self):
+        warn(
+            "The property `file` is deprecated. Use `path_file` instead.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        return self.path_file
+
+    @property
+    def path_file(self):
         output = self.output
-        if output and not self._file:
+        if output and not self._path_file:
             logger.info("Searching for a log file...")
             path_run = Path(output.path_run)
             logfiles = sorted(path_run.glob("*.log"))
             if logfiles:
                 try:
-                    self._file = next(
+                    self._path_file = next(
                         file
                         for file in logfiles
                         if file.name == f"{output.name_solver}.log"
                     )
                 except StopIteration:
-                    self._file = logfiles[-1]
+                    self._path_file = logfiles[-1]
             else:
                 logger.info(f"Cannot find a .log to parse in {path_run}.")
-                self._file = path_run / f"{output.name_solver}.log"
+                self._path_file = path_run / f"{output.name_solver}.log"
 
-        return self._file
+        return self._path_file.resolve()
 
-    @file.setter
-    def file(self, path_log_file):
-        self._file = path_log_file
+    @path_file.setter
+    def path_file(self, path_log_file):
+        self._path_file = path_log_file
 
     @property
     def text(self):
-        with open(self.file) as fp:
+        with open(self.path_file) as fp:
             return fp.read()
 
-    @property
-    def path_run(self):
-        """Parse path_run from log file"""
-        par_file = re.findall("Reading.*par", self.text)[-1].split()[-1]
-        path = Path(par_file).parent
-        return path
+    def load(self, pressure_solver="gmres"):
+        """Load time data from the log file"""
+        # Parse text starting with Step
+        # https://regex101.com/r/enFOAg/1
 
-    @property
-    def params(self):
-        return Simul.load_params_from_file(path_xml=self.path_run / "params_simul.xml")
-
-    @property
-    def text_steps(self):
-        """Parse text starting with Step
+        path_file_time = modification_date(self.path_file)
+        if self.data is not None and path_file_time <= self._data_modif_time:
+            return self.data
 
-        https://regex101.com/r/enFOAg/1
-        """
-        pattern = re.compile(
-            r"""^Step          # For all lines starting with Step
+        pattern_step = r"""^Step          # For all lines starting with Step
                 \s*            # Followed by some whitespaces
-                (\d+)          # 0: Capture timestep which are integers
+                (?P<it>\d+)    # 0: Capture timestep which are integers
                 .*             # Followed by some characters until
                 t=\s*          # t=spaces
-                (\S+)          # 1: Capture t: any non-whitespace char
+                (?P<t>\S+)     # 1: Capture t: any non-whitespace char
                 ,\s*           # comma and spaces
                 DT=\s*         # DT=spaces
-                (\S+)          # 2: Capture DT: any non-whitespace char
+                (?P<dt>\S+)    # 2: Capture DT: any non-whitespace char
                 ,\s*           # comma and spaces
                 C=\s*          # C=spaces
-                (.+)           # 3: Capture C: any char
-    """,
+                (?P<CFL>\S+)   # 3: Capture C: any non-whitespace char
+            """
+        # TODO: add undocumented params.nek.pressure.solver parameter
+        # See Line 445 in Nek5000/core/reader_par.f
+        pattern_pressure = rf"""^\ +
+                (?P<it2>\d+)
+                .*
+                PRES\ {pressure_solver}    # For all lines containing a string like PRES gmres
+                \s+
+                (?P<pres_it>\d+)
+                \s+
+                (?P<pres_div>\S+)
+                \s+
+                (?P<pres_div0>\S+)
+                \s+
+                (?P<pres_tol>\S+)
+                \s+
+                (?P<pres_etime>\S+)
+                \s+
+                (?P<pres_etime1>\S+)
+            """
+
+        expr = re.compile(
+            f"({pattern_step})|({pattern_pressure})",
             re.VERBOSE | re.MULTILINE,
         )
-        return pattern.findall(self.text)
+
+        keys_all = tuple(expr.groupindex)
+        # Divide into two tuples at "it2"
+        idx = keys_all.index("it2")
+        keys_step, keys_pressure = keys_all[:idx], keys_all[idx:]
+
+        matches = tuple(expr.finditer(self.text))
+
+        def make_dict(keys):
+            return {key: [m[key] for m in matches] for key in keys}
+
+        def make_df(keys):
+            return (
+                pd.DataFrame.from_dict(make_dict(keys))
+                .dropna()
+                .astype({key: int if key == "it" else float for key in keys})
+            )
+
+        df_step = make_df(keys_step).set_index("it")
+        df_pressure = (
+            make_df(keys_pressure).rename(columns={"it2": "it"}).set_index("it")
+        )
+
+        self.data = df_step.join(df_pressure)
+        self._data_modif_time = path_file_time
+        return self.data
 
     @property
     def dt(self):
         """Extract time step dt over the course of a simulation."""
-        # Alternate implementations
-        # 1. awkup
-        # dts = %awk -F, -e '/Step/{print $3}' {self.file}
-        #
-        # 2. Basic regex
-        # steps = re.findall("Step.*DT.*", self.text)
-        # try:
-        #     dt_strings = [step.split(',')[2] for step in steps if step]
-        #     dts = [float(s.split("= ")[1]) for s in dt_strings if s]
-        # except IndexError:
-        #     print(steps)
-        dt = [float(field[2]) for field in self.text_steps]
-        return np.array(dt)
+        warn(
+            "The property `dt` is deprecated. Use `load()` instead.",
+            DeprecationWarning,
+        )
+        return self.data.dt
 
     def __call__(self, *args):
         """Print to stdout and log file simultaneously."""
         mpi.printby0(*args)
         if mpi.rank == 0 and self.output._has_to_save:
-            with self.file.open("a") as f:
+            with self.path_file.open("a") as f:
                 f.write(" ".join(str(a) for a in args) + "\n")
+
+    def plot_dt_cfl(self):
+        """Plot the evolution of the time step and the CFL number"""
+
+        df = self.load()
+
+        fig, (ax_top, ax_bot) = plt.subplots(nrows=2, sharex=True)
+
+        ax_bot.plot(df.t, df.dt)
+        ax_bot.set_title("time step")
+        ax_bot.set_xlabel("time")
+
+        ax_top.plot(df.t, df.CFL)
+        ax_top.set_title("CFL number")
+
+        params = self.output.sim.params
+
+        if params.nek.general.variable_dt:
+            ax_top.axhline(params.nek.general.target_cfl, color="r", label="Target CFL")
+            ax_top.legend(loc="lower right")
+
+        fig.tight_layout()
+
+    def plot_nb_iterations(self):
+        """Plot the evolution of the number of iterations for the Poisson solver"""
+
+        df = self.load()
+
+        fig, ax = plt.subplots()
+
+        ax.plot(df.t, df.pres_it)
+        ax.set_ylabel("number of iterations")
+        ax.set_xlabel("time")
+
+        fig.tight_layout()


Test patch:
diff --git a/tests/test_cbox.py b/tests/test_cbox.py
index 806ff782..8b1bcef2 100644
--- a/tests/test_cbox.py
+++ b/tests/test_cbox.py
@@ -46,3 +46,11 @@ def test_phys_fields_get_var_before_load(sim_cbox_executed):
     sim.output.phys_fields.init_reader()
     ux = sim.output.phys_fields.get_var("ux")
     assert isinstance(ux, xr.DataArray)
+
+
+@pytest.mark.slow
+def test_print_stdout(sim_cbox_executed):
+    sim = sim_cbox_executed
+
+    sim.output.print_stdout.plot_nb_iterations()
+    sim.output.print_stdout.plot_dt_cfl()