Geometry: fix Tabulation type annotations (np.ndarray), which broke pdoc in the doc CI job

src/Model/Geometry/ProfileXYZ.py

@@ -39,9 +39,9 @@ logger = logging.getLogger()
 
 @dataclass
 class Tabulation:
-    z: np.array([])
-    A: np.array([])
-    L: np.array([])
+    z: np.ndarray
+    A: np.ndarray
+    L: np.ndarray
 
 
 class ProfileXYZ(Profile, SQLSubModel):

src/Solver/RubarBE.py

@@ -33,6 +33,11 @@ from Model.Results.River.River import River, Reach, Profile
 logger = logging.getLogger()
 
 
+class RubarExportError(Exception):
+    """Study data that cannot be exported; the message is shown to the
+    user and no file is written"""
+
+
 class Rubar3(CommandLineSolver):
     _type = "rubar3"
 
@@ -168,21 +173,363 @@ class Rubar3(CommandLineSolver):
                     qlog.put(trad["warning"].format(message=message))
                 return False
 
+        # 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
+
+        for message in hydraulics["warnings"]:
+            logger.warning(message)
+            if qlog is not None:
+                qlog.put(RunSolverTranslate()["warning"].format(
+                    message=message))
+
         self._export_donnee(study, repertory, qlog, name=name)
         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(study, repertory, qlog, name=name)
-        self._export_stricklers(study, repertory, qlog, name=name)
-        self._export_hydro(study, repertory, qlog, name=name)
-        self._export_condav(study, repertory, qlog, name=name)
+        self._export_condin(hydraulics, repertory, qlog, name=name)
+        self._export_stricklers(hydraulics, repertory, qlog, name=name)
+        self._export_boundaries(hydraulics, repertory, qlog, name=name)
         # self._export_abshyd(study, repertory, qlog, name=name)
 
         return True
 
+    @staticmethod
+    def _cancel(qlog, message):
+        logger.warning(message)
+        if qlog is not None:
+            qlog.put(RunSolverTranslate()["warning"].format(message=message))
+
+    ##############
+    # Hydraulics #
+    ##############
+
+    # Notice §3.1: liquid boundary conditions read for each CONDAM, CONDAV
+    _upstream_types = {"0": (), "1": ("TD",), "2": ()}
+    _downstream_types = {"0": (), "1": ("ZD",), "2": (), "3": ("TZ",),
+                         "5": (), "7": ()}
+    _boundary_kinds = {"TD": "hydrograph Q(t)", "TZ": "water level Z(t)",
+                       "ZD": "rating curve Q(Z)",
+                       "PC": "punctual contribution"}
+
+    def _hydraulic_data(self, study, values):
+        """Check the single reach data needed by DONNEE, MAIL, CONDIN,
+        FROT, HYDRO and CONDAV, and return them. Raise
+        RubarExportError with a user message"""
+        trad = RunSolverTranslate()
+        river = study.river
+
+        edges = river.enable_edges()
+        if len(edges) != 1:
+            raise RubarExportError(trad["multireach_unsupported"].format(
+                count=len(edges)))
+        edge = edges[0]
+        profiles = edge.reach.profiles
+        if len(profiles) < 2:
+            raise RubarExportError(trad["not_enough_profiles"].format(
+                reach=edge.name))
+        rks = [p.rk for p in profiles]
+        if any(b <= a for a, b in zip(rks, rks[1:])):
+            raise RubarExportError(trad["rk_not_increasing"].format(
+                reach=edge.name))
+
+        warnings = []
+        boundaries = self._check_boundaries(river, edge, values)
+
+        return {
+            "edge": edge,
+            "profiles": profiles,
+            "tinit": self._seconds(values.get("rubarbe_tinit", "0")),
+            "frictions": self._interpolated_frictions(edge),
+            "initial": self._initial_conditions(river, edge),
+            "boundaries": boundaries,
+            "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:
+            return float(old_pamhyr_date_to_timestamp(value))
+        return float(value)
+
+    def _check_boundaries(self, river, edge, values):
+        """Liquid boundary conditions used for CONDAM, REGIME and CONDAV.
+        Return {'hydro': bc, 'condav': bc} (or None)"""
+        trad = RunSolverTranslate()
+        condam = values.get("rubarbe_condam", "")
+        condav = values.get("rubarbe_condav", "")
+        regime = values.get("rubarbe_regime", "0")
+
+        if condam not in self._upstream_types:
+            raise RubarExportError(trad["boundary_code_unsupported"].format(
+                key="rubarbe_condam", value=condam))
+        if condav not in self._downstream_types:
+            raise RubarExportError(trad["boundary_code_unsupported"].format(
+                key="rubarbe_condav", value=condav))
+
+        # rubar_init.f90 reads HAMONT only when REGIME=1 and the
+        # computation never uses its values (HIMP, TA, NT3); CONDAM=2 reads
+        # no upstream level at all. The upstream level would be silently
+        # ignored: refuse both settings.
+        for key, value, unsupported in (("REGIME", regime, "1"),
+                                        ("CONDAM", condam, "2")):
+            if value == unsupported:
+                raise RubarExportError(
+                    trad["upstream_level_unsupported"].format(
+                        setting=f"{key}={value}"))
+
+        liquid = [
+            bc for bc in river.boundary_condition.get_tab("liquid")
+            if not bc.is_deleted() and bc.bctype in self._boundary_kinds
+        ]
+
+        def at(node):
+            return [bc for bc in liquid if bc.node == node]
+
+        upstream = self._select_boundaries(
+            at(edge.node1), self._upstream_types[condam], edge.node1,
+            {"TD": f"CONDAM={condam}"}, f"CONDAM={condam}"
+        )
+        downstream = self._select_boundaries(
+            at(edge.node2), self._downstream_types[condav], edge.node2,
+            {"ZD": f"CONDAV={condav}", "TZ": f"CONDAV={condav}"},
+            f"CONDAV={condav}"
+        )
+
+        return {
+            "hydro": upstream.get("TD"),
+            "condav": downstream.get("ZD", downstream.get("TZ")),
+        }
+
+    def _select_boundaries(self, bcs, expected, node, required_by, setting):
+        trad = RunSolverTranslate()
+        selected = {}
+
+        for bc in bcs:
+            kind = self._boundary_kinds[bc.bctype]
+            if bc.bctype not in expected:
+                raise RubarExportError(trad["boundary_unexpected"].format(
+                    name=bc.name, kind=kind, node=node.name,
+                    setting=setting))
+            if bc.bctype in selected:
+                raise RubarExportError(trad["boundary_ambiguous"].format(
+                    kind=kind, node=node.name))
+            if len(bc) == 0:
+                raise RubarExportError(trad["boundary_empty"].format(
+                    name=bc.name, node=node.name))
+            times = [d[0] for d in bc.data]
+            if bc.bctype != "ZD" and any(
+                    b < a for a, b in zip(times, times[1:])):
+                raise RubarExportError(trad["boundary_not_sorted"].format(
+                    name=bc.name, node=node.name))
+            selected[bc.bctype] = bc
+
+        for bctype in expected:
+            if bctype not in selected:
+                raise RubarExportError(trad["boundary_missing"].format(
+                    setting=required_by[bctype],
+                    kind=self._boundary_kinds[bctype], node=node.name))
+
+        return selected
+
+    @staticmethod
+    def _interpolated_frictions(edge):
+        """Minor and medium Strickler at each centre of cell. The model
+        holds one Strickler per zone (Friction.end_strickler returns the
+        begin one): constant inside a zone, linearly interpolated between
+        two zones (notice §3.4), constant beyond the first and last ones"""
+        points = []
+        for friction in edge.frictions.frictions:
+            if friction.is_deleted() or not friction.is_full_defined():
+                continue
+            begin = friction.begin_strickler
+            end = friction.end_strickler or begin
+            pts = [(friction.begin_rk, begin.minor, begin.medium),
+                   (friction.end_rk, end.minor, end.medium)]
+            points += sorted(pts, key=lambda p: p[0])
+
+        if len(points) == 0:
+            raise RubarExportError(
+                RunSolverTranslate()["missing_friction"].format(
+                    reach=edge.name))
+
+        points.sort(key=lambda p: p[0])
+
+        def at(rk):
+            if rk <= points[0][0]:
+                return points[0][1:]
+            if rk >= points[-1][0]:
+                return points[-1][1:]
+            for (r0, *k0), (r1, *k1) in zip(points, points[1:]):
+                if r0 <= rk <= r1:
+                    if r1 == r0:
+                        return tuple(k1)
+                    a = (rk - r0) / (r1 - r0)
+                    return tuple(v0 + a * (v1 - v0)
+                                 for v0, v1 in zip(k0, k1))
+
+        return [at(rk) for rk in edge.reach.inter_profiles_rk()]
+
+    @staticmethod
+    def _initial_conditions(river, edge):
+        """(depth, discharge) at each profile, linearly interpolated in rk
+        between the profiles holding an initial condition (notice §3.5)"""
+        profiles = edge.reach.profiles
+        ics = river.initial_conditions.get(edge)
+        known = {}
+        for data in ics.data:
+            profile = data["rk"]
+            if profile is None or profile not in profiles:
+                continue
+            depth = max(data["elevation"] - profile.z_min(), 0.0)
+            known[profile] = (depth, data["discharge"])
+
+        if len(known) == 0:
+            raise RubarExportError(
+                RunSolverTranslate()["missing_initial_conditions"].format(
+                    reach=edge.name))
+
+        points = sorted((p.rk, *v) for p, v in known.items())
+
+        def at(rk):
+            if rk <= points[0][0]:
+                return points[0][1:]
+            if rk >= points[-1][0]:
+                return points[-1][1:]
+            for (r0, h0, q0), (r1, h1, q1) in zip(points, points[1:]):
+                if r0 <= rk <= r1:
+                    a = (rk - r0) / (r1 - r0)
+                    return (h0 + a * (h1 - h0), q0 + a * (q1 - q0))
+
+        result = [known[p] if p in known else at(p.rk) for p in profiles]
+        if any(q <= -500.0 for _, q in result):
+            # CONDIN stores 1000 + Q, read as a discharge above 500
+            raise RubarExportError(
+                RunSolverTranslate()["discharge_out_of_range"].format(
+                    reach=edge.name))
+
+        return result
+
     # Notice Rubar3/RubarBE, annex 1: records are positional, not keyed.
     # Each record starts with an A50 label, including records of integers.
     _donnee_fields = (
@@ -316,10 +663,16 @@ class Rubar3(CommandLineSolver):
     def _export_donnee(self, study, repertory, qlog=None, name="0"):
         values = self._parameter_values(study)
         fields = list(self._donnee_fields)
-        # The diagnostic transport-capacity extension is optional as a block.
-        if values.get("rubarbe_optfpc", ""):
-            fields += [("optfpc", "integer"), ("ros", "real"),
-                       ("dm", "real"), ("segma", "real")]
+        # Optional OPTFPC. ROS, DM and SEGMA follow it only for a friction
+        # computed from sediments without transport (notice §3.1)
+        optfpc = values.get("rubarbe_optfpc", "")
+        if optfpc:
+            fields += [("optfpc", "integer")]
+            trased = self._flag(values.get("rubarbe_trased", "n"),
+                                "rubarbe_trased")
+            if optfpc != "0" and trased == "N":
+                fields += [("ros", "real"), ("dm", "real"),
+                           ("segma", "real")]
         elif any(values.get("rubarbe_" + key, "")
                  for key in ("ros", "dm", "segma")):
             raise ValueError("rubarbe_optfpc: required before ROS/DM/SEGMA")
@@ -425,81 +778,52 @@ 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")
+        self._write_parameters(repertory, f"geomac-i.{name}", lines)
 
-                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
+    @staticmethod
+    def _section_rk(profile):
+        """Abscissa of a section as written in GEOMAC-I (F11.3). MAIL
+        intermailles must match it within 1e-6 (rubar_init.f90, geomac)"""
+        return float(f"{profile.rk:.3f}")
 
     def _export_mail(self, study, repertory, qlog, name="0"):
+        """MAIL: LM, the LM centres of cells, then the LM - 1 intermailles
+        (notice, annex 1). Free format"""
         if qlog is not None:
             qlog.put("Export MAIL file")
 
-        with open(
-                os.path.join(
-                    repertory, f"mail.{name}"
-                ), "w+"
-        ) as f:
-            for edge in study.river.enable_edges():
-                reach = edge.reach
-                lm = len(reach) + 1
-                f.write(f"{lm:>13}\n")
+        lines = []
+        for edge in study.river.enable_edges():
+            sections = [self._section_rk(p) for p in edge.reach.profiles]
+            centres = ([sections[0]] +
+                       [(a + b) / 2 for a, b in zip(sections, sections[1:])] +
+                       [sections[-1]])
 
-                for mails in [reach.inter_profiles_rk(),
-                              reach.get_rk()]:
-                    ind = 0
-                    for mail in mails:
-                        f.write(f"{mail:15.3f}")
+            lines.append(f"{len(centres):>13}\n")
+            for values in (centres, sections):
+                for i in range(0, len(values), 3):
+                    lines.append("".join(f"{v:15.6f}"
+                                         for v in values[i:i + 3]) + "\n")
 
-                        ind += 1
-                        if ind % 3 == 0:
-                            f.write("\n")
-
-                    if ind % 3 != 0:
-                        f.write("\n")
+        self._write_parameters(repertory, f"mail.{name}", lines)
 
     def _export_abshyd(self, study, repertory, qlog, name="0"):
         if qlog is not None:
@@ -545,208 +869,66 @@ class Rubar3(CommandLineSolver):
 
                 reach_ind += 1
 
-    def _export_stricklers(self, study, repertory, qlog, name="0"):
-        self._export_frot(study, repertory, qlog, name=name, version="")
-        self._export_frot(study, repertory, qlog, name=name, version="2")
+    def _export_stricklers(self, hydraulics, repertory, qlog, name="0"):
+        """FROT (minor bed) and FROT2 (medium bed): LM, then one Strickler
+        per cell (notice, annex 1)"""
+        frictions = hydraulics["frictions"]
+        for version, column in (("", 0), ("2", 1)):
+            if qlog is not None:
+                qlog.put(f"Export FROT{version} file")
 
-    def _export_frot(self, study, repertory, qlog, name="0", version=""):
-        if qlog is not None:
-            qlog.put(f"Export FROT{version} file")
+            lines = [f"{len(frictions):>6}\n"]
+            lines += [f"{i:>6} {values[column]:>12.5f}\n"
+                      for i, values in enumerate(frictions, 1)]
+            self._write_parameters(repertory, f"frot{version}.{name}", lines)
 
-        with open(
-                os.path.join(
-                    repertory, f"frot{version}.{name}"
-                ), "w+"
-        ) as f:
-            for edge in study.river.enable_edges():
-                reach = edge.reach
-                lm = len(reach) + 1
-                f.write(f"{lm:>6}\n")
-
-                lst = list(filter(
-                    lambda f: f.is_full_defined(),
-                    edge.frictions.frictions
-                ))
-
-                rk_min = 9999999.9
-                rk_max = -9999999.9
-                coeff_min = -1.0
-                coeff_max = -1.0
-                for s in lst:  # TODO optimise ?
-                    if s.begin_rk > rk_max:
-                        rk_max = s.begin_rk
-                        coeff_max = s.begin_strickler
-                    if s.begin_rk < rk_min:
-                        rk_min = s.begin_rk
-                        coeff_min = s.begin_strickler
-                    if s.end_rk > rk_max:
-                        rk_max = s.end_rk
-                        coeff_max = s.end_strickler
-                    if s.end_rk < rk_min:
-                        rk_min = s.end_rk
-                        coeff_min = s.end_strickler
-
-                def get_stricklers_from_rk(rk, lst):
-
-                    coeff = None
-                    if rk > rk_max:
-                        coeff = coeff_max
-                    elif rk < rk_min:
-                        coeff = coeff_min
-                    else:
-                        for s in lst:
-                            if (rk >= s.begin_rk and rk <= s.end_rk or
-                                    rk <= s.begin_rk and rk >= s.end_rk):
-                                coeff = s.begin_strickler  # TODO: inerpolate
-                                break
-
-                    # TODO interpolation if rk is not in frictons
-
-                    if coeff is None:
-                        logger.error(
-                            "Study frictions are not fully defined"
-                        )
-                        return None
-
-                    return coeff.medium if version == "2" else coeff.minor
-
-                ind = 1
-                for mail in edge.reach.inter_profiles_rk():
-                    coef = get_stricklers_from_rk(mail, lst)
-                    if coef is not None:
-                        f.write(f"{ind:>6} {coef:>12.5f}")
-
-                    ind += 1
-                    f.write("\n")
-
-    def _export_condin(self, study, repertory, qlog, name="0"):
+    def _export_condin(self, hydraulics, repertory, qlog, name="0"):
+        """CONDIN: TINIT LM, then for each cell I, depth and 1000 + Q
+        (notice §3.5 and annex 1). Cell values are the mean of the two
+        adjacent sections; the first and last cells take the values of
+        the first and last sections"""
         if qlog is not None:
             qlog.put("Export CONDIN file")
 
-        with open(
-                os.path.join(
-                    repertory, f"condin.{name}"
-                ), "w+"
-        ) as f:
-            for edge in study.river.enable_edges():
-                reach = edge.reach
+        sections = hydraulics["initial"]
+        cells = ([sections[0]] +
+                 [((h0 + h1) / 2, (q0 + q1) / 2)
+                  for (h0, q0), (h1, q1) in zip(sections, sections[1:])] +
+                 [sections[-1]])
 
-                f.write(f"{0.0}\n")
+        lines = [f"{hydraulics['tinit']!r} {len(cells)}\n"]
+        lines += [f"{i:>5} {h:.6f} {1000.0 + q:.6f}\n"
+                  for i, (h, q) in enumerate(cells, 1)]
+        self._write_parameters(repertory, f"condin.{name}", lines)
 
-                ics = study.river.initial_conditions.get(edge)
-                data = self._export_condin_init_data(ics)
+    def _export_boundaries(self, hydraulics, repertory, qlog, name="0"):
+        """HYDRO (Q, t) and CONDAV ((Q, Z) for a rating curve, (t, Z) for
+        a Z(t) condition), notice annex 1. A file is written only when the
+        matching condition is used. HAMONT is not written: REGIME=1 and
+        CONDAM=2 are refused (see _check_boundaries)"""
+        bcs = hydraulics["boundaries"]
 
-                profiles = reach.profiles
-                first = profiles[0]
-                last = profiles[-1]
+        def records(bc, columns):
+            return [columns(d0, d1) for d0, d1 in bc.data]
 
-                if first not in data or last not in data:
-                    logger.error(
-                        "Study initial condition is not fully defined"
-                    )
-                    return
+        files = []
+        if bcs["hydro"] is not None:
+            files.append(("hydro", bcs["hydro"],
+                          lambda t, q: f"{float(q)!r} {float(t)!r}\n",
+                          lambda n: f"{n}\n"))
+        if bcs["condav"] is not None:
+            bc = bcs["condav"]
+            if bc.bctype == "ZD":
+                def columns(z, q): return f"{float(q)!r} {float(z)!r}\n"
+            else:
+                def columns(t, z): return f"{float(t)!r} {float(z)!r}\n"
+            files.append(("condav", bc, columns, lambda n: f"{n}\n"))
 
-                f_h_s = self._export_condin_profile_height_speed(first, data)
-                l_h_s = self._export_condin_profile_height_speed(last, data)
-
-                # First mail
-                f.write(f"{1:>5} {f_h_s[0]} {f_h_s[1]}\n")
-
-                ind = 2
-                it = iter(profiles)
-                prev = next(it)
-
-                prev_h, prev_s = f_h_s
-                for profile in it:
-                    if profile not in data:
-                        ind += 1
-                        continue
-
-                    cur_h, cur_s = self._export_condin_profile_height_speed(
-                        profile, data
-                    )
-
-                    # Mean of height and speed
-                    h = (prev_h + cur_h) / 2
-                    s = (prev_s + cur_s) / 2
-
-                    f.write(f"{ind:>5} {h} {s}\n")
-
-                    prev_h, prev_s = cur_h, cur_s
-                    ind += 1
-
-                # Last mail
-                f.write(f"{ind:>5} {f_h_s[0]} {f_h_s[1]}\n")
-
-    def _export_condin_init_data(self, ics):
-        data = {}
-
-        for d in ics.data:
-            data[d['rk']] = (
-                d['elevation'],
-                d['discharge'],
-            )
-
-        return data
-
-    def _export_condin_profile_height_speed(self, profile, data):
-        z = data[profile][0]
-        q = data[profile][1]
-
-        height = z - profile.z_min()
-        speed = profile.speed(q, z)
-
-        return height, speed
-
-    def _export_hydro(self, study, repertory, qlog, name="0"):
-        if qlog is not None:
-            qlog.put("Export HYDRO file")
-
-        with open(
-                os.path.join(
-                    repertory, f"hydro.{name}"
-                ), "w+"
-        ) as f:
-            bcs = []
-            for edge in study.river.enable_edges():
-                for bound in study.river.boundary_condition.get_tab("liquid"):
-                    # BC is an hydrogramme
-                    if bound.bctype == "TD" or bound.bctype == "PC":
-                        # BC is on input node of this reach
-                        if bound.node == edge.node1:
-                            bcs.append(bound)
-
-            for bc in bcs:
-                f.write(f"{len(bc)}\n")
-                for d0, d1 in bc.data:
-                    f.write(f"{d1} {d0}\n")
-
-    def _export_condav(self, study, repertory, qlog, name="0"):
-        if qlog is not None:
-            qlog.put("Export CONDAV file")
-
-        with open(
-                os.path.join(
-                    repertory, f"condav.{name}"
-                ), "w+"
-        ) as f:
-            bcs = []
-            for edge in study.river.enable_edges():
-                for bound in study.river.boundary_condition.get_tab("liquid"):
-                    # BC is an hydrogramme
-                    if bound.bctype == "ZD" or bound.bctype == "TZ":
-                        # BC is on input node of this reach
-                        if bound.node == edge.node2:
-                            bcs.append(bound)
-
-            for bc in bcs:
-                f.write(f"{len(bc)}\n")
-                if bc.bctype == "ZD":
-                    for d0, d1 in bc.data:
-                        f.write(f"{d1} {d0}\n")
-                else:
-                    for d0, d1 in bc.data:
-                        f.write(f"{d0} {d1}\n")
+        for filename, bc, columns, header in files:
+            if qlog is not None:
+                qlog.put(f"Export {filename.upper()} file")
+            lines = [header(len(bc))] + records(bc, columns)
+            self._write_parameters(repertory, f"{filename}.{name}", lines)
 
     def read_profil(self, study, fname, results, qlog=None, name="0"):
         logger.info(f"read: profil.{name}")

src/View/RunSolver/Translate.py

@@ -56,6 +56,105 @@ class RunSolverTranslate(PamhyrTranslate):
             "of reach(es): {reaches}. Add a solid boundary condition, or "
             "disable sediment transport in the solver parameters."
         )
+        self._dict["multireach_unsupported"] = _translate(
+            "Solver",
+            "Export cancelled: the Rubar export supports a single reach, "
+            "{count} enabled reaches found."
+        )
+        self._dict["not_enough_profiles"] = _translate(
+            "Solver",
+            "Export cancelled: reach {reach} needs at least two profiles."
+        )
+        self._dict["rk_not_increasing"] = _translate(
+            "Solver",
+            "Export cancelled: the profiles of reach {reach} must have "
+            "strictly increasing rk. Sort the profiles of the geometry."
+        )
+        self._dict["boundary_code_unsupported"] = _translate(
+            "Solver",
+            "Export cancelled: {key}={value} is not supported by the Rubar "
+            "export (single reach)."
+        )
+        self._dict["boundary_missing"] = _translate(
+            "Solver",
+            "Export cancelled: {setting} requires a {kind} boundary "
+            "condition at node {node}."
+        )
+        self._dict["boundary_unexpected"] = _translate(
+            "Solver",
+            "Export cancelled: boundary condition {name} ({kind}) at node "
+            "{node} is not used with {setting}. Remove it or change the "
+            "boundary settings in the solver parameters."
+        )
+        self._dict["boundary_ambiguous"] = _translate(
+            "Solver",
+            "Export cancelled: several {kind} boundary conditions at node "
+            "{node}."
+        )
+        self._dict["boundary_empty"] = _translate(
+            "Solver",
+            "Export cancelled: boundary condition {name} at node {node} "
+            "has no value."
+        )
+        self._dict["boundary_not_sorted"] = _translate(
+            "Solver",
+            "Export cancelled: the times of boundary condition {name} at "
+            "node {node} must be increasing."
+        )
+        self._dict["upstream_level_unsupported"] = _translate(
+            "Solver",
+            "Export cancelled: {setting} imposes an upstream water level, "
+            "which the current RubarBE version reads but does not use in "
+            "the computation. Choose another upstream boundary setting."
+        )
+        self._dict["missing_friction"] = _translate(
+            "Solver",
+            "Export cancelled: no friction defined on reach {reach}."
+        )
+        self._dict["missing_initial_conditions"] = _translate(
+            "Solver",
+            "Export cancelled: no initial condition defined on reach "
+            "{reach}."
+        )
+        self._dict["discharge_out_of_range"] = _translate(
+            "Solver",
+            "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"
         )

src/View/SolverParameters/RubarFields.py

@@ -42,6 +42,16 @@ OPTION_CHOICES = (
 )
 
 
+# Choices read by RubarBE 2021 but never used by its computation: the
+# upstream level of HAMONT (rubar_init.f90). The export refuses them.
+UNSUPPORTED_CHOICES = {
+    ('rubarbe_condam', '2'): 'La version actuelle de RubarBE ne prend pas '
+                             'en compte la cote amont imposée.',
+    ('rubarbe_regime', '1'): 'La version actuelle de RubarBE ne prend pas '
+                             'en compte le tirant d’eau amont imposé.',
+}
+
+
 @dataclass(frozen=True)
 class Field:
     key: str
@@ -89,7 +99,9 @@ def fields():
            'Nécessite les fichiers de reprise du calcul précédent.', YES_NO)
 
     choice('condam', 1, 'Condition à la limite amont',
-           'Le choix doit correspondre aux données de la limite amont.',
+           'Le choix doit correspondre aux données de la limite amont. '
+           'La cote imposée n’est pas prise en compte par la version '
+           'actuelle de RubarBE.',
            {'0': 'Réflexion', '1': 'Débit imposé Q(t)',
             '2': 'Cote imposée Z(t)', '4': 'Confluence'})
     choice('condav', 1, 'Condition à la limite aval',
@@ -99,8 +111,8 @@ def fields():
             '3': 'Cote imposée Z(t)', '4': 'Confluence / défluence',
             '5': 'Régime uniforme', '7': 'Régime critique'})
     choice('regime', 1, 'Imposer aussi le tirant d’eau amont',
-           'Complète le débit imposé ; incompatible avec une limite '
-           'amont de type "Réflexion".',
+           'Complète le débit imposé ; non pris en compte par la version '
+           'actuelle de RubarBE.',
            ON_OFF)
     choice('iodev', 1, 'Déversement latéral',
            'Nécessite les données de déversements latéraux.',

src/View/SolverParameters/RubarForm.py

@@ -10,7 +10,8 @@ from PyQt5.QtWidgets import (
 )
 
 from View.SolverParameters.RubarFields import (
-    FIELDS, GROUPS, OPTION_CHOICES, Field, active, canonical, validate,
+    FIELDS, GROUPS, OPTION_CHOICES, UNSUPPORTED_CHOICES, Field, active,
+    canonical, validate,
 )
 from View.SolverParameters.UndoCommand import SetCommand
 
@@ -293,6 +294,15 @@ class RubarForm(QWidget):
             index = editor.findData(forbidden)
             if index >= 0:
                 editor.model().item(index).setEnabled(not blocked)
+        # Unsupported choices stay visible, so a stored value is still
+        # displayed, but cannot be selected anymore
+        for (key, value), reason in UNSUPPORTED_CHOICES.items():
+            editor = self.editors.get(key)
+            index = editor.findData(value) if editor is not None else -1
+            if index >= 0:
+                item = editor.model().item(index)
+                item.setEnabled(False)
+                item.setToolTip(tr(reason))
         self._refreshing = False
 
     def _update_context(self, values):
@@ -312,6 +322,11 @@ class RubarForm(QWidget):
                 values.get('rubarbe_regime') == '1'):
             warnings.append(tr('Une limite amont réfléchissante ne permet pas '
                                'd’imposer aussi le tirant d’eau.'))
+        for (key, value), reason in UNSUPPORTED_CHOICES.items():
+            if values.get(key) == value:
+                warnings.append(tr(reason) + tr(
+                    ' Le calcul sera refusé : choisissez un autre réglage '
+                    'de la limite amont.'))
         try:
             from tools import old_pamhyr_date_to_timestamp