# RubarBE.py -- Pamhyr # Copyright (C) 2024-2025 INRAE # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # -*- coding: utf-8 -*- import os import logging import math from decimal import Decimal, InvalidOperation, localcontext import numpy as np from tools import timer, trace, old_pamhyr_date_to_timestamp, logger_exception from Solver.CommandLine import CommandLineSolver from View.RunSolver.Translate import RunSolverTranslate from Model.Results.Results import Results 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" def __init__(self, name): super(Rubar3, self).__init__(name) self._type = "rubar3" self._cmd_input = "" self._cmd_solver = "@path @args @input" self._cmd_output = "" @classmethod def default_parameters(cls): # Keep historical rubarbe_* keys for existing Rubar3 studies. lst = [ ("rubarbe_cfl", "0.50000E+00"), ("rubarbe_condam", "1"), ("rubarbe_condav", "3"), ("rubarbe_regime", "0"), ("rubarbe_iodev", "n"), ("rubarbe_iodebord", "o"), ("rubarbe_iostockage", "o"), ("rubarbe_iopdt", "y"), ("rubarbe_iovis", "n"), ("rubarbe_rep", "n"), ("rubarbe_tinit", "000:00:00:00"), ("rubarbe_tmax", "001:00:00:00"), # 1 day ("rubarbe_tiopdt", "000:00:00:00"), ("rubarbe_dt", "5.0"), ("rubarbe_ts", "999:99:99:00"), ("rubarbe_dtsauv", "00:00:05:00"), ("rubarbe_psave", "00:00:05:00"), ("rubarbe_fdeb1", "1"), ("rubarbe_fdeb2", "10"), ("rubarbe_fdeb3", "100"), ("rubarbe_tf_1", "y"), ("rubarbe_tf_2", "y"), ("rubarbe_tf_3", "y"), ("rubarbe_tf_4", "y"), ("rubarbe_tf_5", "y"), ("rubarbe_tf_6", "n"), ("rubarbe_trased", "n"), ("rubarbe_optfpc", "0"), # trased parameters ("rubarbe_ros", "2650.0"), ("rubarbe_dm", "0.1"), ("rubarbe_segma", "1.0"), # Sediment parameters ("rubarbe_sediment_ros", "2650.0"), ("rubarbe_sediment_por", "0.4"), ("rubarbe_sediment_dcharg", "0.0"), ("rubarbe_sediment_halfa", "1.0"), ("rubarbe_sediment_mult_1", "1.0"), ("rubarbe_sediment_mult_2", ""), ("rubarbe_sediment_mult_3", ""), ("rubarbe_sediment_mult_4", ""), ("rubarbe_sediment_mult_5", ""), ("rubarbe_sediment_tcadim", "0.047"), ("rubarbe_sediment_visc", "1.3E-06"), ("rubarbe_sediment_opts", "3"), ("rubarbe_sediment_odchar", "0"), ("rubarbe_sediment_unisol", "1"), ("rubarbe_sediment_typdef", "3"), ("rubarbe_sediment_depot", "2"), ("rubarbe_sediment_choixc", "2"), ("rubarbe_sediment_option", "2"), ("rubarbe_sediment_capsol", "1"), ("rubarbe_sediment_bmiu", "0.85"), ("rubarbe_sediment_demix", "0"), ("rubarbe_sediment_defond", "1"), ("rubarbe_sediment_varcons", "1"), ("rubarbe_sediment_dchard", "0.0"), ("rubarbe_sediment_dchars", "0.0"), ] return lst @classmethod def checkers(cls): lst = [ ] return lst ########## # Export # ########## def cmd_args(self, study): lst = super(Rubar3, self).cmd_args(study) return lst def input_param(self): name = self._study.name return f"{name}" def output_param(self): name = self._study.name return f"profil.{name}" def log_file(self): name = self._study.name return f"geomac-i.{name}" def export(self, study, repertory, qlog=None): self._study = study name = study.name.replace(" ", "_") values = self._parameter_values(study) if (values.get("rubarbe_regime") == "1" and values.get("rubarbe_condam") == "0"): trad = RunSolverTranslate() message = trad["incompatible_upstream_boundary"] logger.warning(message) if qlog is not None: qlog.put(trad["warning"].format(message=message)) return False trased = self._parameter_flag(study, "trased") == "O" # Check the hydraulic data before writing any file try: hydraulics = self._hydraulic_data(study, values) geometry = self._geometry_data(hydraulics["edge"], trased) solid = self._solid_boundary( study.river, hydraulics["edge"], values, 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(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) self._export_stricklers(hydraulics, repertory, qlog, name=name) self._export_boundaries(hydraulics, repertory, qlog, name=name) if solid is not None: self._export_hydros(solid, repertory, qlog, name=name) if trased: self._export_condins(geometry, hydraulics["tinit"], 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", "SL": "solid discharge Qs(t)"} # NCLMAX (variables.f90): points of all the solid hydrographs _max_solid_points = 1000 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)) left, right, active_left, active_right = self._bed_limits(markers) limits = [("main_channel", "rg", "rd", left, right)] if trased: limits.append(("active_bed", "g", "d", active_left, active_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 @staticmethod def _bed_limits(markers): """Indices of the G, D, g and d points of a section. 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(markers) - 1) return left, right, index("g", left), index("d", right) 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 def _solid_boundary(self, river, edge, values, trased): """Solid boundary condition written in HYDROS. RubarBE reads it at the upstream node when TRASED is enabled and CONDAM=1 (CONDAM=2 is refused), and has no downstream solid condition (rubar_init.f90, inits). Without sediment transport the solid conditions are ignored. Return the condition or None""" if not trased: return None trad = RunSolverTranslate() condam = values.get("rubarbe_condam", "") solid = [ bc for bc in river.boundary_condition.get_tab("solid") if not bc.is_deleted() and bc.bctype == "SL" ] for bc in solid: if bc.node == edge.node2: raise RubarExportError( trad["solid_boundary_downstream"].format( name=bc.name, node=edge.node2.name)) expected = ("SL",) if condam == "1" else () selected = self._select_boundaries( [bc for bc in solid if bc.node == edge.node1], expected, edge.node1, {"SL": f"TRASED=O, CONDAM={condam}"}, f"CONDAM={condam}" ) bc = selected.get("SL") if bc is None: return None if len(bc) > self._max_solid_points: raise RubarExportError(trad["solid_boundary_too_long"].format( name=bc.name, node=edge.node1.name, count=len(bc), maximum=self._max_solid_points)) # RubarBE stops on a zero D or S if not (bc.d50 > 0 and bc.sigma >= 1): raise RubarExportError(trad["solid_boundary_invalid"].format( name=bc.name, node=edge.node1.name)) return bc @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 = ( ("cfl", "real"), ("condam", "integer"), ("condav", "integer"), ("regime", "integer"), ("iodev", "flags"), ("iopdt", "flag"), ("iovis", "flag"), ("rep", "flag"), ("tinit", "time"), ("tmax", "time"), ("tiopdt", "time"), ("dt", "real"), ("ts", "time"), ("dtsauv", "time"), ("psave", "time"), ("fdeb1", "real"), ("fdeb2", "real"), ("fdeb3", "real"), ("tf_1", "flag"), ("tf_2", "flag"), ("tf_3", "flag"), ("tf_4", "flag"), ("tf_5", "flag"), ("tf_6", "flag"), ("trased", "flag"), ) _ts_fields = ( ("ros", "real"), ("por", "real"), ("dcharg", "real"), ("halfa", "real"), ("mult_1", "multipliers"), ("visc", "visc"), ("opts", "law"), ("odchar", "integer"), ("unisol", "integer"), ("typdef", "integer"), ("depot", "depot"), ("choixc", "integer"), ("option", "options"), ("capsol", "integer"), ("bmiu", "real"), ("demix", "integer"), ("defond", "integer"), ("varcons", "integer"), ) # rubar_init.f90 (inits): laws whose critical stress may come from the # Shields diagram, and laws without critical stress _shields_laws = {1, 10, 11, 12, 14, 25, 26} _no_threshold_laws = {2, 3, 8, 13, 18, 19, 22, 23, 24} _brownlie_law = 16 @classmethod def tcadim_required(cls, opts, option): """True if RubarBE reads the dimensionless critical stress TCADIM for transport law OPTS and options OPTION (notice §3.10)""" try: law = int(opts) except ValueError: return False shields = len(option) == 3 and option[1] != "1" if shields and law in cls._shields_laws: return False if law in cls._no_threshold_laws or law == cls._brownlie_law: return False return True def _parameter_values(self, study): return {p.name: str(p.value).strip() for p in study.river.get_params(self.type).parameters} @staticmethod def _required(values, key): value = values.get(key, "") if not value: raise ValueError(f"Missing Rubar parameter: {key}") return value @staticmethod def _flag(value, key, mixed=False): result = value.upper() if result == "Y": result = "O" if result not in (("O", "N", "M") if mixed else ("O", "N")): raise ValueError(f"{key}: expected O/Y or N" + (" or M" if mixed else "")) return result def _parameter_flag(self, study, suffix): key = "rubarbe_" + suffix return self._flag(self._required(self._parameter_values(study), key), key) @staticmethod def _integer(value, key, width=1): if not value.isascii() or not value.isdigit() or len(value) > width: raise ValueError(f"{key}: expected an integer in {width} columns") return value.rjust(width) @staticmethod def _real(value, key, scientific=False): try: number = Decimal(value.replace("d", "E").replace("D", "E")) except InvalidOperation: raise ValueError(f"{key}: invalid real value {value!r}") from None if not number.is_finite(): raise ValueError(f"{key}: expected a finite real value") if scientific: # DONNEE: E12.5. Reject overflow instead of shifting columns. floating = float(number) if not math.isfinite(floating) or (number and floating == 0): raise ValueError(f"{key}: real value outside E12.5 range") result = f"{floating:12.5E}" width = 12 else: # TS: F10.0 input accepts an explicit decimal point/exponent. # Preserve the value exactly; do not round porosity etc. to int. with localcontext() as context: context.prec = max(28, len(number.as_tuple().digits)) number = number.normalize() result = str(number).upper() if len(result) > 10: result = format(number, "E") width = 10 if len(result) > width: raise ValueError(f"{key}: value {value!r} exceeds {width} columns") return result.rjust(width) @staticmethod def _write_parameters(repertory, filename, lines): # Validate all records before opening/truncating the destination. path = os.path.join(repertory, filename) with open(path, "w", encoding="ascii") as f: f.writelines(lines) def _export_donnee(self, study, repertory, qlog=None, name="0"): values = self._parameter_values(study) fields = list(self._donnee_fields) # 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") lines = [] for suffix, kind in fields: key = "rubarbe_" + suffix value = self._required(values, key) if kind == "integer": result = self._integer(value, key) elif kind in ("flag", "flags"): result = self._flag(value, key, mixed=suffix == "iopdt") if kind == "flags": # Missing optional flags mean O (notice section 3.1). # Write both explicitly so STOCKAGE cannot shift left. for extra in ("iodebord", "iostockage"): extra_key = "rubarbe_" + extra result += self._flag(values.get(extra_key) or "O", extra_key) else: if kind == "time" and ":" in value: parts = value.split(":") if len(parts) != 4 or not all(p.isdigit() for p in parts): raise ValueError( f"{key}: expected days:hours:minutes:seconds") value = str(old_pamhyr_date_to_timestamp(value)) result = self._real(value, key, scientific=True) lines.append(f"{key:<50}{result}\n") if qlog is not None: qlog.put("Export DONNEE file") self._write_parameters(repertory, f"donnee.{name}", lines) def _export_ts(self, study, repertory, qlog=None, name="0"): trad = RunSolverTranslate() values = self._parameter_values(study) prefix = "rubarbe_sediment_" fields = list(self._ts_fields) demix = self._required(values, prefix + "demix") if demix not in ("0", "1"): raise ValueError(f"{prefix}demix: expected 0 or 1") if demix == "1": fields += [("dchard", "real"), ("dchars", "real")] lines = [] for suffix, kind in fields: key = prefix + suffix value = self._required(values, key) if kind == "real": result = self._real(value, key) elif kind == "visc": # Notice annex 1: TCADIM then VISC, read as (A50, A20) then # two F10.0. A single value would be read as TCADIM, with # VISC = 0 (rubar_init.f90, inits). if Decimal(self._real(value, key)) <= 0: raise ValueError(f"{key}: expected a positive viscosity") tcadim_key = prefix + "tcadim" tcadim = values.get(tcadim_key, "") if tcadim: result = self._real(tcadim, tcadim_key) elif self.tcadim_required(values.get(prefix + "opts", ""), values.get(prefix + "option", "")): message = trad["missing_tcadim"].format(key=tcadim_key) logger.warning(message) if qlog is not None: qlog.put(trad["warning"].format(message=message)) return False else: result = " " * 10 result += self._real(value, key) elif kind == "multipliers": result = self._real(value, key) extras = [values.get(prefix + f"mult_{i}", "") for i in range(2, 6)] # Preserve empty interior F10.0 fields; omit trailing fields # to retain the solver's defaults for optional coefficients. while extras and not extras[-1]: extras.pop() for i, extra in enumerate(extras, 2): result += (self._real(extra, prefix + f"mult_{i}") if extra else " " * 10) elif kind == "options": if len(value) != 3 or any(c not in "12" for c in value): message = trad["invalid_options"].format(key=key) logger.warning(message) if qlog is not None: qlog.put(trad["warning"].format(message=message)) return False result = value elif kind == "law": # Notice section 3.10 / 8.7: 4--7 are not listed. if value not in {str(i) for i in [1, 2, 3, *range(8, 31)]}: message = trad["unsupported_transport_law"].format( key=key, value=value) logger.warning(message) if qlog is not None: qlog.put(trad["warning"].format(message=message)) return False result = self._integer(value, key, width=2) else: result = self._integer(value, key, width=10 if kind == "depot" else 1) lines.append(f"{key:<50}{result}\n") if qlog is not None: qlog.put("Export TS file") self._write_parameters(repertory, f"ts.{name}", lines) return True 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") 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" ) self._write_parameters(repertory, f"geomac-i.{name}", lines) @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") 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]]) 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") self._write_parameters(repertory, f"mail.{name}", lines) def _export_abshyd(self, study, repertory, qlog, name="0"): if qlog is not None: qlog.put("Export ABSHYD file") with open( os.path.join( repertory, f"abshyd.{name}" ), "w+" ) as f: reach_ind = 1 for edge in study.river.enable_edges(): 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: qlog.put("Export DEVERS file") with open( os.path.join( repertory, f"devers.{name}" ), "w+" ) as f: reach_ind = 1 for edge in study.river.enable_edges(): 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} {0.0:15.3f} {0.0:15.3f} {1.0:>15.3f}\n") ind += 1 reach_ind += 1 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") 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) 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") sections = hydraulics["initial"] cells = ([sections[0]] + [((h0 + h1) / 2, (q0 + q1) / 2) for (h0, q0), (h1, q1) in zip(sections, sections[1:])] + [sections[-1]]) 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) 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"] def records(bc, columns): return [columns(d0, d1) for d0, d1 in bc.data] 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")) 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 _export_hydros(self, bc, repertory, qlog, name="0"): """HYDROS: number of points, then Qs, t, D, S (notice annex 1). Qs is a solid discharge (kg/s) or a concentration (kg/m3) according to UNISOL, D and S are those of the condition""" if qlog is not None: qlog.put("Export HYDROS file") d50, sigma = float(bc.d50), float(bc.sigma) lines = [f"{len(bc)}\n"] + [ f"{float(qs)!r} {float(t)!r} {d50!r} {sigma!r}\n" for t, qs in bc.data ] self._write_parameters(repertory, f"hydros.{name}", lines) def _export_condins(self, geometry, tinit, repertory, qlog, name="0"): """CONDINS: TINIT, then for each section (intermaille) its number, Qs, D and S of the active layer (notice annex 1, rubar_init.f90, inits). Pamhyr2 holds no initial solid condition: Qs = 0 (empty active layer) and D, S are the mean of the top compartments of the points strictly between g and d, as advised by the notice §3.5""" if qlog is not None: qlog.put("Export CONDINS file") lines = [f"{tinit!r}\n"] for i, (_, rows) in enumerate(geometry, 1): markers = [marker for marker, _, _ in rows] _, _, left, right = self._bed_limits(markers) tops = [compartments[0] for _, _, compartments in rows[left + 1:right]] # Rounded to 10 significant digits, so that a mean of equal # values is written as these values d = float(f"{sum(c[1] for c in tops) / len(tops):.10g}") s = float(f"{sum(c[2] for c in tops) / len(tops):.10g}") lines.append(f"{i} 0.0 {d!r} {s!r}\n") self._write_parameters(repertory, f"condins.{name}", lines) def read_profil(self, study, fname, results, qlog=None, name="0"): logger.info(f"read: profil.{name}") with open(fname, "r") as f: reachs = [] for i, edge in enumerate(study.river.enable_edges()): reach = edge.reach # Add results reach to reach list res_reach = results.river.add(i) reachs.append((res_reach, len(reach) + 1)) results.set("nb_reach", f"{len(reachs)}") results.set("nb_profile", f"{sum([r[1] - 1 for r in reachs])}") def read_data_line(f): line = f.readline().split() if len(line) > 3: return tuple(map(float, line)) else: return (0.0, 0.0, 0.0, 0.0) def set_and_compute_limites(reach, ind, z, q, s): reach.set(ind, timestamp, "Z", z) reach.set(ind, timestamp, "Q", q) reach.set(ind, timestamp, "V", s) profile = reach.profile(ind) limits = profile.geometry.get_water_limits(z) reach.set( ind, timestamp, "water_limits", limits ) ts = set() timestamp = next(filter( lambda p: p.name == 'rubarbe_tinit', study.river.get_params(self._type).parameters )).value if timestamp.count(':') == 3: timestamp = old_pamhyr_date_to_timestamp(timestamp) timestamp = float(timestamp) ts.add(timestamp) # add initial condition for r, edge in enumerate(study.river.enable_edges()): reach = edge.reach ics = study.river.initial_conditions.get(edge) q = ics.get_discharge() z = ics.get_elevation() k = 0 for i, j in zip(z, q): v = reach.profiles[k].speed(i, j) set_and_compute_limites(reachs[r][0], k, i, j, v) k += 1 # start read end = False while True: line = f.readline() if line == "": break timestamp = float(line) ts.add(timestamp) logger.info(f"read: profil.{name}: timestamp = {timestamp}") for reach, lm in reachs: # First profile h, s, q, z = read_data_line(f) set_and_compute_limites(reach, 0, z+h, q, s) # For each profile ind = 1 ph, ps, pq, pz = read_data_line(f) while ind < lm - 2: h, s, q, z = read_data_line(f) set_and_compute_limites( reach, ind, (pz + z + ph + h) / 2, (pq + q) / 2, (ps + s) / 2 ) ph = h ps = s pq = q pz = z ind += 1 # Last profile h, s, q, z = read_data_line(f) set_and_compute_limites(reach, ind, z+h, q, s) # Tables and plots use timestamp rows and global profile columns. timestamps = sorted(ts) profiles = [] for reach, _ in reachs: offset = len(profiles) profiles.extend(reach.profiles) reach.set_global_index(range(offset, len(profiles))) results.set("timestamps", timestamps) results.compute_timestamp_index() results.set("table", { key: np.array([ [profile.get_ts_key(t, key) for profile in profiles] for t in timestamps ], dtype=float) for key in ("Z", "Q", "V") }) @timer def write_bin(self, study, fname, results, qlog=None, name="0"): logger.info(f"write_bin: Start writing '{fname}' ...") with open(fname, "wb") as f: def newline(j): return np.asarray([j], dtype=np.int32).tofile(f) def endline(j): return np.asarray([j], dtype=np.int32).tofile(f) def write_int(i): newline(len(i)*4) np.array(i, dtype=np.int32).tofile(f) endline(len(i)*4) def write_float(i): newline(len(i)*4) np.array(i, dtype=np.float32).tofile(f) endline(len(i)*4) def write_float64(i): newline(len(i)*8) np.array(i, dtype=np.float64).tofile(f) endline(len(i)*8) def write_float64(i): newline(len(i)*8) np.array(i, dtype=np.float64).tofile(f) endline(len(i)*8) def write_data(npts, t, a, val): newline(npts*4 + 13) np.array(npts, dtype=np.int32).tofile(f) np.array(t, dtype=np.float64).tofile(f) np.array(bytearray(a.encode()), dtype=np.byte).tofile(f) np.array(val, dtype=np.float32).tofile(f) endline(npts*4 + 13) ts_list = sorted(results.get("timestamps")) profiles = [] for r in results.river.reachs: profiles += r.profiles # Meta data (1st line) nb_reach = len(results.river) nb_profile = len(profiles) write_int([nb_reach, nb_profile, "82"]) # Reach information (2nd line) is1 = [] is2 = [] ltmp = [] i = 1 for r in results.river.reachs: j = i+len(r)-1 ltmp += [i, j] is1.append(i) is2.append(j) i += i+len(r) write_int(ltmp) # X (3rd line) rk = [] for r in results.river.reachs: rk += r.geometry.get_rk() write_float(rk) # Z and Y (4th line) # ltmp = [] # for r in results.river.reachs: # for p in r.prifiles: # ltpm.append(p.rk) write_float(3*nb_profile*[0.0]) # Data for timestamp in ts_list: q = list( map( lambda p: p.get_ts_key(timestamp, "Q"), profiles ) ) write_data(nb_profile, timestamp, "Q", q) z = list( map( lambda p: p.get_ts_key(timestamp, "Z"), profiles ) ) write_data(nb_profile, timestamp, "Z", z) logger.info(f"write_bin: ... end with {len(ts_list)} timestamps") @timer def results(self, study, repertory, qlog=None, name=None): results = Results( study=study, solver=self, repertory=repertory, name=name, ) if name is None: name = study.name.replace(" ", "_") results_file = f"profil.{name}" fname = os.path.join(repertory, results_file) if not os.path.isfile(fname): logger.warning(f"Result file {results_file} does not exist") return None try: self.read_profil(study, fname, results, qlog, name=name) except Exception as e: logger.error(f"Failed to read results") logger_exception(e) return None return results class RubarBE(Rubar3): _type = "rubarbe" def __init__(self, name): super(RubarBE, self).__init__(name) self._type = "rubarbe" self._cmd_input = "" self._cmd_solver = "@path @args @input" self._cmd_output = "" @classmethod def default_parameters(cls): # New BE studies evolve the bed. Keep Rubar3 defaults unchanged. overrides = { "rubarbe_trased": "y", "rubarbe_sediment_opts": "1", # Meyer-Peter and Mueller "rubarbe_sediment_option": "212", # no Ikeda, given Shields, total } return [(key, overrides.get(key, value)) for key, value in super().default_parameters()]