RubarBE: export save abscissas (abshyd) and times (tnprof) from two solver parameters, and read results with tnprof save

src/Solver/RubarBE.py

@@ -17,6 +17,7 @@
 # -*- coding: utf-8 -*-
 
 import os
+import re
 import logging
 import math
 from decimal import Decimal, InvalidOperation, localcontext
@@ -76,11 +77,14 @@ class Rubar3(CommandLineSolver):
             ("rubarbe_fdeb2", "10"),
             ("rubarbe_fdeb3", "100"),
             ("rubarbe_tf_1", "y"),
-            ("rubarbe_tf_2", "y"),
+            # TF(2) and TF(4) need the abscissas and times below
+            ("rubarbe_tf_2", "n"),
             ("rubarbe_tf_3", "y"),
-            ("rubarbe_tf_4", "y"),
+            ("rubarbe_tf_4", "n"),
             ("rubarbe_tf_5", "y"),
             ("rubarbe_tf_6", "n"),
+            ("rubarbe_abshyd", ""),
+            ("rubarbe_tnprof", ""),
             ("rubarbe_trased", "n"),
             ("rubarbe_optfpc", "0"),
             # trased parameters
@@ -167,6 +171,7 @@ class Rubar3(CommandLineSolver):
             geometry = self._geometry_data(hydraulics["edge"], trased)
             solid = self._solid_boundary(
                 study.river, hydraulics["edge"], values, trased)
+            saves = self._save_data(values, hydraulics)
         except RubarExportError as error:
             self._cancel(qlog, str(error))
             return False
@@ -192,7 +197,7 @@ class Rubar3(CommandLineSolver):
         if trased:
             self._export_condins(geometry, hydraulics["tinit"], repertory,
                                  qlog, name=name)
-        # self._export_abshyd(study, repertory, qlog, name=name)
+        self._export_saves(saves, repertory, qlog, name=name)
 
         return True
 
@@ -497,6 +502,99 @@ class Rubar3(CommandLineSolver):
                 name=bc.name, node=edge.node1.name))
         return bc
 
+    #########
+    # Saves #
+    #########
+
+    # NBHYPR (variables.f90): abscissas of ABSHYD and times of TNPROF.
+    # RubarBE does not check this size
+    _max_saves = 100
+
+    @staticmethod
+    def save_list(text, times=False):
+        """Numbers of a list parameter (ABSHYD abscissas, TNPROF times),
+        separated by spaces or semicolons. Times are in seconds or in
+        days:hours:minutes:seconds. Raise ValueError with the invalid
+        token"""
+        values = []
+        for token in re.split(r"[\s;]+", text.strip()):
+            if token == "":
+                continue
+            try:
+                if times and ":" in token:
+                    parts = token.split(":")
+                    if len(parts) != 4 or not all(
+                            p.isascii() and p.isdigit() for p in parts):
+                        raise ValueError(token)
+                    value = float(sum(
+                        int(p) * factor for p, factor in
+                        zip(parts, (86400, 3600, 60, 1))))
+                else:
+                    value = float(token)
+            except ValueError:
+                raise ValueError(token) from None
+            if not math.isfinite(value):
+                raise ValueError(token)
+            values.append(value)
+        return values
+
+    def _save_data(self, values, hydraulics):
+        """Abscissas of ABSHYD, read when TF(2)=O, and times of TNPROF,
+        always read: they give the saves of PROFIL and PROFILS when
+        TF(4)=O and of the complementary outputs when IOVIS=O
+        (rubar_init.f90, init2; rubarbe.f90, SAVEVISU). Return
+        {'abshyd': [rk] or None, 'tnprof': [t]}"""
+        trad = RunSolverTranslate()
+        lists = {}
+        for key, times in (("abshyd", False), ("tnprof", True)):
+            key = f"rubarbe_{key}"
+            try:
+                lists[key] = self.save_list(values.get(key, ""), times)
+            except ValueError as error:
+                raise RubarExportError(trad["save_list_invalid"].format(
+                    key=key, value=str(error)))
+            if len(lists[key]) > self._max_saves:
+                raise RubarExportError(trad["save_list_too_long"].format(
+                    key=key, count=len(lists[key]),
+                    maximum=self._max_saves))
+
+        def enabled(flag):
+            return values.get(f"rubarbe_{flag}", "n").lower() in ("y", "o")
+
+        for flag, key in (("tf_2", "rubarbe_abshyd"),
+                          ("tf_4", "rubarbe_tnprof")):
+            if enabled(flag) and not lists[key]:
+                raise RubarExportError(trad["save_list_required"].format(
+                    flag=f"rubarbe_{flag}", key=key))
+
+        # RubarBE looks for each abscissa between two centres of cells,
+        # from the first to the last section
+        abshyd = None
+        if enabled("tf_2"):
+            edge = hydraulics["edge"]
+            first = self._section_rk(hydraulics["profiles"][0])
+            last = self._section_rk(hydraulics["profiles"][-1])
+            for rk in lists["rubarbe_abshyd"]:
+                if not first <= round(rk, 3) <= last:
+                    raise RubarExportError(trad["abshyd_outside"].format(
+                        rk=f"{rk:g}", first=f"{first:g}", last=f"{last:g}",
+                        reach=edge.name))
+            abshyd = lists["rubarbe_abshyd"]
+
+        # Times outside the computation do not falsify it: RubarBE saves
+        # the earlier ones at the first time steps and never reaches the
+        # later ones
+        tinit = hydraulics["tinit"]
+        tmax = self._seconds(values.get("rubarbe_tmax", "0"))
+        tnprof = sorted(lists["rubarbe_tnprof"])
+        outside = [t for t in tnprof if not tinit <= t <= tmax]
+        if outside:
+            hydraulics["warnings"].append(trad["tnprof_outside"].format(
+                times=" ".join(f"{t:g}" for t in outside),
+                tinit=f"{tinit:g}", tmax=f"{tmax:g}"))
+
+        return {"abshyd": abshyd, "tnprof": tnprof}
+
     @staticmethod
     def _interpolated_frictions(edge):
         """Minor and medium Strickler at each centre of cell. The model
@@ -856,27 +954,28 @@ class Rubar3(CommandLineSolver):
 
         self._write_parameters(repertory, f"mail.{name}", lines)
 
-    def _export_abshyd(self, study, repertory, qlog, name="0"):
+    def _export_saves(self, saves, repertory, qlog, name="0"):
+        """ABSHYD (number, then I4, 1X, F11.3 per abscissa; single reach,
+        its number is omitted) when TF(2)=O, and TNPROF (number, then
+        index and time per line). TNPROF is always read by RubarBE: an
+        empty list is written as 0, so that a previous TNPROF of the
+        directory is not used"""
+        if saves["abshyd"] is not None:
             if qlog is not None:
                 qlog.put("Export ABSHYD file")
+            lines = [f"{len(saves['abshyd']):>4}\n"] + [
+                f"{i:>4} {rk:>11.3f}\n"
+                for i, rk in enumerate(saves["abshyd"], 1)
+            ]
+            self._write_parameters(repertory, f"abshyd.{name}", lines)
 
-        with open(
+        if qlog is not None:
-                os.path.join(
+            qlog.put("Export TNPROF file")
-                    repertory, f"abshyd.{name}"
+        lines = [f"{len(saves['tnprof']):>12}\n"] + [
-                ), "w+"
+            f"{i:>12} {t:>24.3f}\n"
-        ) as f:
+            for i, t in enumerate(saves["tnprof"], 1)
-            reach_ind = 1
+        ]
-            for edge in study.river.enable_edges():
+        self._write_parameters(repertory, f"tnprof.{name}", lines)
-                reach = edge.reach
-                lm = len(reach) + 1
-                f.write(f"{lm:>13}\n")
-
-                ind = 1
-                for mail in reach.inter_profiles_rk():
-                    f.write(f"{ind:>4} {mail:15.3f} {reach_ind:>4}\n")
-                    ind += 1
-
-                reach_ind += 1
 
     def _export_devers(self, study, repertory, qlog, name="0"):
         if qlog is not None:
@@ -1055,7 +1154,9 @@ class Rubar3(CommandLineSolver):
                 if line == "":
                     break
 
-                timestamp = float(line)
+                # Saves at the times of TNPROF are marked 'tnprof'
+                # (rubarbe.f90, SAVEVISU)
+                timestamp = float(line.split()[0])
                 ts.add(timestamp)
 
                 logger.info(f"read: profil.{name}: timestamp = {timestamp}")

src/View/RunSolver/Translate.py

@@ -60,6 +60,34 @@ class RunSolverTranslate(PamhyrTranslate):
             "Export cancelled: solid boundary condition {name} at node "
             "{node} has {count} points, RubarBE accepts at most {maximum}."
         )
+        self._dict["save_list_invalid"] = _translate(
+            "Solver",
+            "Export cancelled: {key}: '{value}' is not a valid value. "
+            "Separate the values with spaces or semicolons, use a dot as "
+            "decimal separator, and give the times in seconds or "
+            "days:hours:minutes:seconds."
+        )
+        self._dict["save_list_too_long"] = _translate(
+            "Solver",
+            "Export cancelled: {key} has {count} values, RubarBE accepts "
+            "at most {maximum}."
+        )
+        self._dict["save_list_required"] = _translate(
+            "Solver",
+            "Export cancelled: {flag} is enabled but {key} is empty. Fill "
+            "it in the solver parameters (tab Saves), or disable {flag}."
+        )
+        self._dict["abshyd_outside"] = _translate(
+            "Solver",
+            "Export cancelled: rubarbe_abshyd: abscissa {rk} is outside "
+            "reach {reach}, whose sections go from {first} to {last}."
+        )
+        self._dict["tnprof_outside"] = _translate(
+            "Solver",
+            "rubarbe_tnprof: times {times} s are outside the computation, "
+            "from {tinit} s to {tmax} s. RubarBE saves the earlier ones at "
+            "the first time steps and never reaches the later ones."
+        )
         self._dict["solid_boundary_invalid"] = _translate(
             "Solver",
             "Export cancelled: solid boundary condition {name} at node "

src/View/SolverParameters/RubarFields.py

@@ -262,6 +262,17 @@ def fields():
                'Les abscisses ou instants demandés '
                'doivent être définis dans les données correspondantes.',
                YES_NO)
+    add('abshyd', 4, 'Abscisses des hydrogrammes (m)',
+        'Abscisses (PK) où sont enregistrés hydrogrammes et limnigrammes '
+        'au pas de sauvegarde des hydrogrammes, séparées par des espaces '
+        'ou des points-virgules ; 100 au plus. Nécessaires si '
+        '« Hydrogrammes aux abscisses choisies » vaut Oui.',
+        kind='rk_list')
+    add('tnprof', 4, 'Instants de sauvegarde des profils (j:h:min:s ou s)',
+        'Instants séparés par des espaces ou des points-virgules ; 100 au '
+        'plus. Nécessaires si « Profils aux instants choisis » vaut Oui ; '
+        'les sorties complémentaires sont aussi écrites à ces instants.',
+        kind='time_list')
     for i in range(1, 4):
         add(f'fdeb{i}', 4, f'Seuil de débit du front {i} (m³/s)',
             'Débit pour lequel à chaque abscisse, est enregistré le temps '
@@ -298,6 +309,18 @@ def validate(field, value):
         return value
     if field.kind == 'text':
         return value
+    if field.kind in ('rk_list', 'time_list'):
+        try:
+            numbers = Rubar3.save_list(value, field.kind == 'time_list')
+        except ValueError as error:
+            raise ValueError(
+                f'Valeur non reconnue : {error}. Séparez les valeurs par '
+                'des espaces ou des points-virgules, avec un point pour '
+                'les décimales.') from None
+        if len(numbers) > Rubar3._max_saves:
+            raise ValueError(f'{Rubar3._max_saves} valeurs au plus '
+                             '(limite de RubarBE).')
+        return value
     if not value and field.optional:
         return value
     value = value.replace(',', '.').replace('d', 'E').replace('D', 'E')
@@ -364,4 +387,12 @@ def active(field, values):
         return values.get('rubarbe_' + source) in accepted
     if key.startswith('fdeb'):
         return values.get('rubarbe_tf_6', 'n').lower() in ('y', 'o')
+
+    def enabled(flag):
+        return values.get('rubarbe_' + flag, 'n').lower() in ('y', 'o')
+
+    if key == 'abshyd':
+        return enabled('tf_2')
+    if key == 'tnprof':
+        return enabled('tf_4') or enabled('iovis')
     return True