diff --git a/src/Solver/RubarBE.py b/src/Solver/RubarBE.py index 939df3eb..ed697e23 100644 --- a/src/Solver/RubarBE.py +++ b/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,6 +173,19 @@ 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) + 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): @@ -175,14 +193,230 @@ class Rubar3(CommandLineSolver): self._export_geomac_i(study, 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, HAMONT 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, warnings + ) + + 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, + } + + @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, warnings): + """Liquid boundary conditions used for CONDAM, REGIME and CONDAV. + Return {'hydro': bc, 'hamont': 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)) + + 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_types = list(self._upstream_types[condam]) + if regime == "1": + upstream_types.append("TZ") + upstream = self._select_boundaries( + at(edge.node1), upstream_types, edge.node1, + {"TD": f"CONDAM={condam}", "TZ": "REGIME=1"}, + f"CONDAM={condam}, REGIME={regime}" + ) + downstream = self._select_boundaries( + at(edge.node2), self._downstream_types[condav], edge.node2, + {"ZD": f"CONDAV={condav}", "TZ": f"CONDAV={condav}"}, + f"CONDAV={condav}" + ) + + # rubar_init.f90: HAMONT is read when REGIME=1, but the 2021 + # code never uses it; CONDAM=2 reads no upstream level at all + if regime == "1" or condam == "2": + warnings.append(trad["upstream_level_unused"]) + + return { + "hydro": upstream.get("TD"), + "hamont": upstream.get("TZ"), + "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 +550,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") @@ -474,32 +714,32 @@ class Rubar3(CommandLineSolver): 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 +785,69 @@ 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), HAMONT (Z, 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""" + 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["hamont"] is not None: + files.append(("hamont", bcs["hamont"], + lambda t, z: f"{float(z):15.6f} {float(t):15.6f}\n", + lambda n: f"{n:>6}\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}") diff --git a/src/View/RunSolver/Translate.py b/src/View/RunSolver/Translate.py index 29ce3b97..dce8b5eb 100644 --- a/src/View/RunSolver/Translate.py +++ b/src/View/RunSolver/Translate.py @@ -56,6 +56,70 @@ 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_unused"] = _translate( + "Solver", + "The upstream water level (REGIME=1 or CONDAM=2) is not used " + "by the current RubarBE version: the computation ignores it." + ) + 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["export_failed"] = _translate( "Solver", "Export failed; task cancelled" ) diff --git a/tests_cases/RubarBE/Seal_run1/Seal_run1.pamhyr b/tests_cases/RubarBE/Seal_run1/Seal_run1.pamhyr index bb0fb548..c1120616 100644 Binary files a/tests_cases/RubarBE/Seal_run1/Seal_run1.pamhyr and b/tests_cases/RubarBE/Seal_run1/Seal_run1.pamhyr differ