RubarBE: fix hydraulic export (mail, condin, frot, boundary conditions), and prevent user of mistakes instead of run incorrect calculations

This commit is contained in:
Dylan Jeannin committed 2026-10-16 09:55:40 +02:00
1 parent a18accb404
commit ec81a998cb
3 files changed
+381 -216

No files matched your search

+317 -216
View File
@@ -33,6 +33,11 @@ from Model.Results.River.River import River, Reach, Profile
logger = logging.getLogger() 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): class Rubar3(CommandLineSolver):
_type = "rubar3" _type = "rubar3"
@@ -168,6 +173,19 @@ class Rubar3(CommandLineSolver):
qlog.put(trad["warning"].format(message=message)) qlog.put(trad["warning"].format(message=message))
return False 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) self._export_donnee(study, repertory, qlog, name=name)
if trased: if trased:
if not self._export_ts(study, repertory, qlog, name=name): 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_geomac_i(study, repertory, qlog, name=name)
self._export_mail(study, repertory, qlog, name=name) self._export_mail(study, repertory, qlog, name=name)
self._export_devers(study, repertory, qlog, name=name) self._export_devers(study, repertory, qlog, name=name)
self._export_condin(study, repertory, qlog, name=name) self._export_condin(hydraulics, repertory, qlog, name=name)
self._export_stricklers(study, repertory, qlog, name=name) self._export_stricklers(hydraulics, repertory, qlog, name=name)
self._export_hydro(study, repertory, qlog, name=name) self._export_boundaries(hydraulics, repertory, qlog, name=name)
self._export_condav(study, repertory, qlog, name=name)
# self._export_abshyd(study, repertory, qlog, name=name) # self._export_abshyd(study, repertory, qlog, name=name)
return True 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. # Notice Rubar3/RubarBE, annex 1: records are positional, not keyed.
# Each record starts with an A50 label, including records of integers. # Each record starts with an A50 label, including records of integers.
_donnee_fields = ( _donnee_fields = (
@@ -316,10 +550,16 @@ class Rubar3(CommandLineSolver):
def _export_donnee(self, study, repertory, qlog=None, name="0"): def _export_donnee(self, study, repertory, qlog=None, name="0"):
values = self._parameter_values(study) values = self._parameter_values(study)
fields = list(self._donnee_fields) fields = list(self._donnee_fields)
# The diagnostic transport-capacity extension is optional as a block. # Optional OPTFPC. ROS, DM and SEGMA follow it only for a friction
if values.get("rubarbe_optfpc", ""): # computed from sediments without transport (notice §3.1)
fields += [("optfpc", "integer"), ("ros", "real"), optfpc = values.get("rubarbe_optfpc", "")
("dm", "real"), ("segma", "real")] 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, "") elif any(values.get("rubarbe_" + key, "")
for key in ("ros", "dm", "segma")): for key in ("ros", "dm", "segma")):
raise ValueError("rubarbe_optfpc: required before ROS/DM/SEGMA") raise ValueError("rubarbe_optfpc: required before ROS/DM/SEGMA")
@@ -474,32 +714,32 @@ class Rubar3(CommandLineSolver):
ind += 1 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"): 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: if qlog is not None:
qlog.put("Export MAIL file") qlog.put("Export MAIL file")
with open( lines = []
os.path.join(
repertory, f"mail.{name}"
), "w+"
) as f:
for edge in study.river.enable_edges(): for edge in study.river.enable_edges():
reach = edge.reach sections = [self._section_rk(p) for p in edge.reach.profiles]
lm = len(reach) + 1 centres = ([sections[0]] +
f.write(f"{lm:>13}\n") [(a + b) / 2 for a, b in zip(sections, sections[1:])] +
[sections[-1]])
for mails in [reach.inter_profiles_rk(), lines.append(f"{len(centres):>13}\n")
reach.get_rk()]: for values in (centres, sections):
ind = 0 for i in range(0, len(values), 3):
for mail in mails: lines.append("".join(f"{v:15.6f}"
f.write(f"{mail:15.3f}") for v in values[i:i + 3]) + "\n")
ind += 1 self._write_parameters(repertory, f"mail.{name}", lines)
if ind % 3 == 0:
f.write("\n")
if ind % 3 != 0:
f.write("\n")
def _export_abshyd(self, study, repertory, qlog, name="0"): def _export_abshyd(self, study, repertory, qlog, name="0"):
if qlog is not None: if qlog is not None:
@@ -545,208 +785,69 @@ class Rubar3(CommandLineSolver):
reach_ind += 1 reach_ind += 1
def _export_stricklers(self, study, repertory, qlog, name="0"): def _export_stricklers(self, hydraulics, repertory, qlog, name="0"):
self._export_frot(study, repertory, qlog, name=name, version="") """FROT (minor bed) and FROT2 (medium bed): LM, then one Strickler
self._export_frot(study, repertory, qlog, name=name, version="2") per cell (notice, annex 1)"""
frictions = hydraulics["frictions"]
def _export_frot(self, study, repertory, qlog, name="0", version=""): for version, column in (("", 0), ("2", 1)):
if qlog is not None: if qlog is not None:
qlog.put(f"Export FROT{version} file") qlog.put(f"Export FROT{version} file")
with open( lines = [f"{len(frictions):>6}\n"]
os.path.join( lines += [f"{i:>6} {values[column]:>12.5f}\n"
repertory, f"frot{version}.{name}" for i, values in enumerate(frictions, 1)]
), "w+" self._write_parameters(repertory, f"frot{version}.{name}", lines)
) 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( def _export_condin(self, hydraulics, repertory, qlog, name="0"):
lambda f: f.is_full_defined(), """CONDIN: TINIT LM, then for each cell I, depth and 1000 + Q
edge.frictions.frictions (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"""
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"):
if qlog is not None: if qlog is not None:
qlog.put("Export CONDIN file") qlog.put("Export CONDIN file")
with open( sections = hydraulics["initial"]
os.path.join( cells = ([sections[0]] +
repertory, f"condin.{name}" [((h0 + h1) / 2, (q0 + q1) / 2)
), "w+" for (h0, q0), (h1, q1) in zip(sections, sections[1:])] +
) as f: [sections[-1]])
for edge in study.river.enable_edges():
reach = edge.reach
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) def _export_boundaries(self, hydraulics, repertory, qlog, name="0"):
data = self._export_condin_init_data(ics) """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 def records(bc, columns):
first = profiles[0] return [columns(d0, d1) for d0, d1 in bc.data]
last = profiles[-1]
if first not in data or last not in data: files = []
logger.error( if bcs["hydro"] is not None:
"Study initial condition is not fully defined" files.append(("hydro", bcs["hydro"],
) lambda t, q: f"{float(q)!r} {float(t)!r}\n",
return lambda n: f"{n}\n"))
if bcs["hamont"] is not None:
f_h_s = self._export_condin_profile_height_speed(first, data) files.append(("hamont", bcs["hamont"],
l_h_s = self._export_condin_profile_height_speed(last, data) lambda t, z: f"{float(z):15.6f} {float(t):15.6f}\n",
lambda n: f"{n:>6}\n"))
# First mail if bcs["condav"] is not None:
f.write(f"{1:>5} {f_h_s[0]} {f_h_s[1]}\n") bc = bcs["condav"]
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": if bc.bctype == "ZD":
for d0, d1 in bc.data: def columns(z, q): return f"{float(q)!r} {float(z)!r}\n"
f.write(f"{d1} {d0}\n")
else: else:
for d0, d1 in bc.data: def columns(t, z): return f"{float(t)!r} {float(z)!r}\n"
f.write(f"{d0} {d1}\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 read_profil(self, study, fname, results, qlog=None, name="0"): def read_profil(self, study, fname, results, qlog=None, name="0"):
logger.info(f"read: profil.{name}") logger.info(f"read: profil.{name}")
+64
View File
@@ -56,6 +56,70 @@ class RunSolverTranslate(PamhyrTranslate):
"of reach(es): {reaches}. Add a solid boundary condition, or " "of reach(es): {reaches}. Add a solid boundary condition, or "
"disable sediment transport in the solver parameters." "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( self._dict["export_failed"] = _translate(
"Solver", "Export failed; task cancelled" "Solver", "Export failed; task cancelled"
) )
Binary file not shown.