diff --git a/src/Solver/RubarBE.py b/src/Solver/RubarBE.py index 52ae75bb..d1012844 100644 --- a/src/Solver/RubarBE.py +++ b/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): diff --git a/src/View/RunSolver/Translate.py b/src/View/RunSolver/Translate.py index 72e663f2..40fd0c7e 100644 --- a/src/View/RunSolver/Translate.py +++ b/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" )