RubarBE: export sediment layers and bed limits in geomac-i, and refuse invalid layers or markers

src/Solver/RubarBE.py

@@ -176,6 +176,7 @@ class Rubar3(CommandLineSolver):
         # Check the hydraulic data before writing any file
         try:
             hydraulics = self._hydraulic_data(study, values)
+            geometry = self._geometry_data(hydraulics["edge"], trased)
         except RubarExportError as error:
             self._cancel(qlog, str(error))
             return False
@@ -190,7 +191,7 @@ class Rubar3(CommandLineSolver):
         if trased:
             if not self._export_ts(study, repertory, qlog, name=name):
                 return False
-        self._export_geomac_i(study, repertory, qlog, name=name)
+        self._export_geomac_i(geometry, repertory, qlog, name=name)
         self._export_mail(study, repertory, qlog, name=name)
         self._export_devers(study, repertory, qlog, name=name)
         self._export_condin(hydraulics, repertory, qlog, name=name)
@@ -252,6 +253,119 @@ class Rubar3(CommandLineSolver):
             "warnings": warnings,
         }
 
+    ############
+    # Geometry #
+    ############
+
+    # Point names giving the limits of the main channel (G, D) and of the
+    # active bed (g, d) in GEOMAC-I (rubar_init.f90). Other names are not
+    # written
+    _point_markers = {"rg": "G", "rd": "D", "g": "g", "d": "d"}
+    # RubarBE never erodes a compartment whose critical stress is above
+    # 999 Pa (rubarbe.f90, EROSEC)
+    _hard_stress = 999.999
+    # CSMAX (variables.f90)
+    _max_compartments = 10
+
+    @staticmethod
+    def _profile_label(profile):
+        name = profile.name if profile.name != "" else "#"
+        return f"{name} (rk {profile.rk:g})"
+
+    def _geometry_data(self, edge, trased):
+        """Check the GEOMAC-I data of the reach and return its sections
+        [(rk, [(marker, y, [(z_top, d, s, tau_c), ...]), ...]), ...].
+        Raise RubarExportError with a user message"""
+        sections = []
+        for profile in edge.reach.profiles:
+            points = profile.points
+            markers = self._markers(edge, profile, points, trased)
+            station = profile.get_station()
+
+            rows = []
+            for i, (point, marker, y) in enumerate(
+                    zip(points, markers, station)):
+                if trased:
+                    compartments = self._compartments(
+                        edge, profile, i + 1, point)
+                else:
+                    # Only the bed level is used by the hydraulics
+                    compartments = [(point.z, 1.0, 1.0, self._hard_stress)]
+                rows.append((marker, y, compartments))
+
+            sections.append((self._section_rk(profile), rows))
+        return sections
+
+    def _markers(self, edge, profile, points, trased):
+        trad = RunSolverTranslate()
+        markers = [self._point_markers.get(p.name.strip(), " ")
+                   for p in points]
+        names = {v: k for k, v in self._point_markers.items()}
+        for marker in "GDgd":
+            if markers.count(marker) > 1:
+                raise RubarExportError(trad["duplicate_point_name"].format(
+                    name=names[marker], profile=self._profile_label(profile),
+                    reach=edge.name))
+
+        # Defaults of RubarBE: first and last points for G and D, G and D
+        # for g and d. Erosion only applies strictly between g and d
+        def index(marker, default):
+            return markers.index(marker) if marker in markers else default
+
+        left = index("G", 0)
+        right = index("D", len(points) - 1)
+        limits = [("main_channel", "rg", "rd", left, right)]
+        if trased:
+            limits.append(("active_bed", "g", "d",
+                           index("g", left), index("d", right)))
+        for bed, lname, rname, i, j in limits:
+            if j - i < 2:
+                raise RubarExportError(trad["bed_limits_invalid"].format(
+                    bed=trad[bed], left=lname, right=rname,
+                    profile=self._profile_label(profile), reach=edge.name))
+        return markers
+
+    def _compartments(self, edge, profile, number, point):
+        """Sediment layers of a point as RubarBE compartments, from the
+        top: (top elevation, D, S, critical stress). A non erodible
+        compartment is added under an erodible last layer, so that the
+        layer thicknesses are kept"""
+        trad = RunSolverTranslate()
+        sl = point.sl
+        if sl is None or sl.is_deleted() or len(sl) == 0:
+            raise RubarExportError(trad["missing_sediment_layers"].format(
+                point=number, profile=self._profile_label(profile),
+                reach=edge.name))
+
+        layers = sl.layers
+        for layer in layers:
+            if not (layer.height >= 0 and layer.d50 > 0 and
+                    layer.sigma >= 1 and layer.critical_constraint >= 0):
+                raise RubarExportError(trad["invalid_sediment_layer"].format(
+                    layer=layer.name, sl=sl.name,
+                    profile=self._profile_label(profile), reach=edge.name))
+
+        top = point.z
+        compartments = []
+        for layer in layers:
+            if layer.height == 0:
+                continue
+            compartments.append((top, layer.d50, layer.sigma,
+                                 layer.critical_constraint))
+            top -= layer.height
+
+        if not compartments or compartments[-1][3] <= 999:
+            last = layers[-1]
+            compartments.append((top, last.d50, last.sigma,
+                                 self._hard_stress))
+
+        if len(compartments) > self._max_compartments:
+            raise RubarExportError(trad["too_many_compartments"].format(
+                sl=sl.name, count=len(compartments),
+                maximum=self._max_compartments,
+                profile=self._profile_label(profile), reach=edge.name))
+        return compartments
+
     @staticmethod
     def _seconds(value):
         if ":" in value:
@@ -664,54 +778,25 @@ class Rubar3(CommandLineSolver):
         self._write_parameters(repertory, f"ts.{name}", lines)
         return True
 
-    def _export_geomac_i(self, study, repertory, qlog, name="0"):
+    def _export_geomac_i(self, geometry, repertory, qlog, name="0"):
+        """GEOMAC-I: number of sections, then for each section its
+        header and one line per point (A1, 1X, F11.5, then up to 10
+        compartments in F13.5, 2F15.10, F15.5; rubar_init.f90)"""
         if qlog is not None:
             qlog.put("Export GEOMAC-i file")
 
-        with open(
-                os.path.join(
-                    repertory, f"geomac-i.{name}"
-                ), "w+"
-        ) as f:
-            for edge in study.river.enable_edges():
-                reach = edge.reach
-                n_profiles = len(reach)
-                time = 0.0
+        lines = [f"{len(geometry):>4} {0.0:>12.3f}\n"]
+        for i, (rk, rows) in enumerate(geometry, 1):
+            lines.append(f"{i:>4} {rk:>11.3f} {len(rows):>4}\n")
+            for marker, y, compartments in rows:
+                lines.append(
+                    f"{marker} {y:>11.5f}" +
+                    "".join(f"{z:>13.5f}{d:>15.10f}{s:>15.10f}{t:>15.5f}"
+                            for z, d, s, t in compartments) +
+                    "\n"
+                )
 
-                f.write(f"{n_profiles:>5} {time:>11.3f}\n")
-
-                ind = 1
-                for profile in reach.profiles:
-                    rk = profile.rk
-                    n_points = len(profile)
-
-                    f.write(f"{ind:>4} {rk:>11.3f} {n_points:>4}\n")
-
-                    station = profile.get_station()
-                    for i, point in enumerate(profile.points):
-                        label = point.name.lower()
-                        if label != "":
-                            if label[0] == "r":
-                                label = label[1].upper()
-                            else:
-                                label = " "
-                        else:
-                            label = " "
-
-                        y = station[i]
-                        z = point.z
-                        dcs = 1.0
-                        scs = 1.0
-                        tmcs = 0.0
-
-                        f.write(
-                            f"{label} {y:>11.5f}" +
-                            f"{z:>13.5f}{dcs:>15.10f}" +
-                            f"{scs:>15.10f}{tmcs:>15.5f}" +
-                            "\n"
-                        )
-
-                    ind += 1
+        self._write_parameters(repertory, f"geomac-i.{name}", lines)
 
     @staticmethod
     def _section_rk(profile):

src/View/RunSolver/Translate.py

@@ -121,6 +121,40 @@ class RunSolverTranslate(PamhyrTranslate):
             "Export cancelled: initial discharges of reach {reach} must be "
             "greater than -500 m3/s."
         )
+        self._dict["duplicate_point_name"] = _translate(
+            "Solver",
+            "Export cancelled: several points are named '{name}' in "
+            "profile {profile} of reach {reach}."
+        )
+        self._dict["main_channel"] = _translate("Solver", "main channel")
+        self._dict["active_bed"] = _translate("Solver", "active bed")
+        self._dict["bed_limits_invalid"] = _translate(
+            "Solver",
+            "Export cancelled: in profile {profile} of reach {reach}, the "
+            "{bed} limits (points '{left}' and '{right}', or the profile "
+            "ends) must be ordered from left to right with at least one "
+            "point between them."
+        )
+        self._dict["missing_sediment_layers"] = _translate(
+            "Solver",
+            "Export cancelled: point {point} of profile {profile} of reach "
+            "{reach} has no sediment layer, which the sediment transport "
+            "requires."
+        )
+        self._dict["invalid_sediment_layer"] = _translate(
+            "Solver",
+            "Export cancelled: layer '{layer}' of sediment layers '{sl}' "
+            "(profile {profile} of reach {reach}) must have a thickness "
+            ">= 0, D50 > 0 m, sigma >= 1 (1 for a uniform sediment) and a "
+            "critical constraint >= 0 Pa."
+        )
+        self._dict["too_many_compartments"] = _translate(
+            "Solver",
+            "Export cancelled: sediment layers '{sl}' (profile {profile} of "
+            "reach {reach}) give {count} compartments, RubarBE accepts at "
+            "most {maximum}. A non erodible compartment is added under an "
+            "erodible last layer."
+        )
         self._dict["export_failed"] = _translate(
             "Solver", "Export failed; task cancelled"
         )