Compare commits

...
4 Commits
3 changed files with 418 additions and 27 deletions

No files matched your search

+1
View File
@@ -44,6 +44,7 @@ __old__
mage8/
adists/
rubarbe/
rubar/
linux/
windows/
+307 -24
View File
@@ -170,8 +170,13 @@ class Rubar3(CommandLineSolver):
hydraulics = self._hydraulic_data(study, values)
geometry = self._geometry_data(hydraulics["edge"], trased)
solid = self._solid_boundary(
study.river, hydraulics["edge"], values, trased)
study.river, hydraulics["edge"], values, trased,
hydraulics["warnings"])
saves = self._save_data(values, hydraulics)
laterals = self._lateral_data(
study.river, hydraulics, values, hydraulics["warnings"])
structures = self._structure_data(
study.river, hydraulics, hydraulics["warnings"])
except RubarExportError as error:
self._cancel(qlog, str(error))
return False
@@ -189,11 +194,15 @@ class Rubar3(CommandLineSolver):
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)
if laterals is not None:
self._export_hydevs(laterals, trased, repertory, qlog, name=name)
self._export_ouvrag(structures, 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)
self._export_hydros(solid, hydraulics["tinit"], repertory, qlog,
name=name)
if trased:
self._export_condins(geometry, hydraulics["tinit"], repertory,
qlog, name=name)
@@ -244,7 +253,7 @@ class Rubar3(CommandLineSolver):
reach=edge.name))
warnings = []
boundaries = self._check_boundaries(river, edge, values)
boundaries = self._check_boundaries(river, edge, values, warnings)
return {
"edge": edge,
@@ -381,7 +390,7 @@ class Rubar3(CommandLineSolver):
return float(old_pamhyr_date_to_timestamp(value))
return float(value)
def _check_boundaries(self, river, edge, values):
def _check_boundaries(self, river, edge, values, warnings):
"""Liquid boundary conditions used for CONDAM, REGIME and CONDAV.
Return {'hydro': bc, 'condav': bc} (or None)"""
trad = RunSolverTranslate()
@@ -417,12 +426,12 @@ class Rubar3(CommandLineSolver):
upstream = self._select_boundaries(
at(edge.node1), self._upstream_types[condam], edge.node1,
{"TD": f"CONDAM={condam}"}, f"CONDAM={condam}"
{"TD": f"CONDAM={condam}"}, f"CONDAM={condam}", warnings
)
downstream = self._select_boundaries(
at(edge.node2), self._downstream_types[condav], edge.node2,
{"ZD": f"CONDAV={condav}", "TZ": f"CONDAV={condav}"},
f"CONDAV={condav}"
f"CONDAV={condav}", warnings
)
return {
@@ -430,9 +439,27 @@ class Rubar3(CommandLineSolver):
"condav": downstream.get("ZD", downstream.get("TZ")),
}
def _select_boundaries(self, bcs, expected, node, required_by, setting):
@staticmethod
def _ignore_empty(bcs, node, warnings, meaningful=()):
"""A boundary condition without value carries no data: it is
ignored with a warning, except for the types where an empty
condition has a meaning. Return (kept, ignored)"""
trad = RunSolverTranslate()
kept, ignored = [], []
for bc in bcs:
if len(bc) == 0 and bc.bctype not in meaningful:
warnings.append(trad["boundary_empty_ignored"].format(
name=bc.name, node=node.name))
ignored.append(bc)
else:
kept.append(bc)
return kept, ignored
def _select_boundaries(self, bcs, expected, node, required_by, setting,
warnings, meaningful=()):
trad = RunSolverTranslate()
selected = {}
bcs, ignored = self._ignore_empty(bcs, node, warnings, meaningful)
for bc in bcs:
kind = self._boundary_kinds[bc.bctype]
@@ -443,9 +470,6 @@ class Rubar3(CommandLineSolver):
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:])):
@@ -455,18 +479,26 @@ class Rubar3(CommandLineSolver):
for bctype in expected:
if bctype not in selected:
raise RubarExportError(trad["boundary_missing"].format(
message = trad["boundary_missing"].format(
setting=required_by[bctype],
kind=self._boundary_kinds[bctype], node=node.name))
kind=self._boundary_kinds[bctype], node=node.name)
# The required condition may be the ignored empty one
empty = [bc.name for bc in ignored if bc.bctype == bctype]
if empty:
message += " " + trad["boundary_missing_empty"].format(
names=", ".join(empty))
raise RubarExportError(message)
return selected
def _solid_boundary(self, river, edge, values, trased):
def _solid_boundary(self, river, edge, values, trased, warnings):
"""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"""
ignored. At the upstream node, an empty solid condition means a
supply at equilibrium in Pamhyr2 (Mage 'default'), written as a
negative Qs (_export_hydros). Return the condition or None"""
if not trased:
return None
trad = RunSolverTranslate()
@@ -476,17 +508,19 @@ class Rubar3(CommandLineSolver):
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))
downstream, _ = self._ignore_empty(
[bc for bc in solid if bc.node == edge.node2], edge.node2,
warnings)
if downstream:
raise RubarExportError(
trad["solid_boundary_downstream"].format(
name=downstream[0].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}"
f"CONDAM={condam}", warnings, meaningful=("SL",)
)
bc = selected.get("SL")
if bc is None:
@@ -502,6 +536,211 @@ class Rubar3(CommandLineSolver):
name=bc.name, node=edge.node1.name))
return bc
############
# Laterals #
############
# Pamhyr2 lateral contribution types; evaporation is a withdrawal
_lateral_signs = {"LC": 1.0, "RA": 1.0, "EV": -1.0}
# QMU, TMU (variables.f90, LMAX*20): RubarBE copies the series of an
# inflow in each cell of its zone, without checking this size
_max_lateral_points = 3000 * 20
def _lateral_data(self, river, hydraulics, values, warnings):
"""Lateral inflows of HYDEVS, read when IODEV=O (rubar_init.f90,
init2). Return [(first cell, last cell, [(t, q), ...]), ...]
sorted from upstream, or None when IODEV=N.
A contribution from section b to section e (1-based) covers the
cells b+1 to e, i.e. from section b to section e (RubarBE divides
the discharge by this length); a contribution at a single section
b covers the cell downstream of it (upstream of it for the last
section). The first and last cells of the reach are never used,
as RubarBE requires. With IODEV=N, the contributions are ignored
with a warning. Otherwise, solid and suspension contributions have
no RubarBE equivalent (their sediment cannot be linked to an
inflow): they are refused"""
trad = RunSolverTranslate()
edge = hydraulics["edge"]
profiles = hydraulics["profiles"]
iodev = values.get("rubarbe_iodev", "n").lower() in ("y", "o")
lst = river.lateral_contribution
def on_reach(tab):
# A contribution without value carries no data: ignored
kept = []
for lc in lst.get_tab(tab):
if lc.is_deleted() or lc.reach != edge:
continue
if len(lc) == 0:
warnings.append(trad["lateral_empty_ignored"].format(
name=lc.name, reach=edge.name))
else:
kept.append(lc)
return kept
lcs = on_reach("liquid")
sediment = on_reach("solid") + on_reach("suspenssion")
# Not read by RubarBE: the user is warned, the computation runs
# without them
if not iodev:
for lc in lcs + sediment:
warnings.append(trad["lateral_iodev_disabled"].format(
name=lc.name, reach=edge.name))
return None
if sediment:
raise RubarExportError(trad["lateral_solid_unsupported"].format(
name=sediment[0].name, reach=edge.name))
index = {id(p): i for i, p in enumerate(profiles, 1)}
count = len(profiles)
blocks = []
for lc in lcs:
if lc.lctype not in self._lateral_signs:
raise RubarExportError(trad["lateral_type_unsupported"].format(
name=lc.name, reach=edge.name))
if (id(lc.begin_section) not in index or
id(lc.end_section) not in index):
raise RubarExportError(trad["lateral_sections"].format(
name=lc.name, reach=edge.name))
b, e = sorted((index[id(lc.begin_section)],
index[id(lc.end_section)]))
first, last = (b + 1, e) if b < e else ((min(b + 1, count),) * 2)
data = [(float(d[0]), float(d[1])) for d in lc.data]
times = [t for t, _ in data]
if any(t1 < t0 for t0, t1 in zip(times, times[1:])):
raise RubarExportError(trad["lateral_not_sorted"].format(
name=lc.name, reach=edge.name))
sign = self._lateral_signs[lc.lctype]
blocks.append((first, last, [(t, sign * q) for t, q in data],
lc.name))
blocks.sort(key=lambda block: block[0])
for (_, last, _, name0), (first, _, _, name1) in zip(
blocks, blocks[1:]):
if first <= last:
raise RubarExportError(trad["lateral_overlap"].format(
first=name0, second=name1, reach=edge.name))
points = sum((last - first + 1) * len(data)
for first, last, data, _ in blocks)
if points > self._max_lateral_points:
raise RubarExportError(trad["lateral_too_long"].format(
count=points, maximum=self._max_lateral_points))
return [block[:3] for block in blocks]
##############
# Structures #
##############
# Pamhyr2 basic structures (those of Mage) with a RubarBE equivalent
# (notice, annex 5): weir and rectangular orifice as a weir-orifice D,
# circular orifice as O
_structure_types = {"S1": "D", "OR": "D", "OC": "O"}
# Loading elevation of a weir that never loads (default of Pamhyr2)
_no_loading = 9999.999
# nou1Dmax, noe1Dmax (variables.f90)
_max_structures = 50
_max_elementary = 10
def _structure_data(self, river, hydraulics, warnings):
"""Structures of OUVRAG (rubar_init.f90, LOUVR1D). A Pamhyr2
structure sits on a section: RubarBE structure on this single
interface, with reference -2 (the discharge is that of the
structure only), as for Mage. Disabled structures are skipped, as
for Mage. Return [(rk, [(type, L, zd, zm, coefficient), ...]), ...]
sorted from upstream. Raise RubarExportError with a user message"""
trad = RunSolverTranslate()
edge = hydraulics["edge"]
sections = {id(p) for p in hydraulics["profiles"]}
structures = {}
for hs in river.hydraulic_structures.lst:
if (hs.is_deleted() or not hs.enabled
or hs.input_reach is None or hs.input_section is None
or hs.input_reach != edge):
continue
if id(hs.input_section) not in sections:
raise RubarExportError(trad["structure_section"].format(
name=hs.name, reach=edge.name))
# Basic structures without RubarBE equivalent are not taken
# into account: the user is warned
lines = []
for bhs in hs.basic_structures:
if not bhs.enabled or bhs.is_deleted():
continue
if bhs.type not in self._structure_types:
warnings.append(trad["structure_type_ignored"].format(
name=self._basic_label(hs, bhs), type=bhs.type))
continue
lines.append(self._elementary_structure(hs, bhs))
if not lines:
warnings.append(trad["structure_empty_ignored"].format(
name=hs.name))
continue
if len(lines) > self._max_elementary:
raise RubarExportError(trad["structure_too_many"].format(
name=hs.name, count=len(lines),
maximum=self._max_elementary))
rk = self._section_rk(hs.input_section)
if rk in structures:
raise RubarExportError(trad["structure_same_section"].format(
first=structures[rk][1], second=hs.name, rk=f"{rk:g}"))
structures[rk] = (lines, hs.name)
if len(structures) > self._max_structures:
raise RubarExportError(trad["structure_too_many_structures"]
.format(count=len(structures),
maximum=self._max_structures))
return [(rk, structures[rk][0]) for rk in sorted(structures)]
@staticmethod
def _basic_label(hs, bhs):
return f"{bhs.name or bhs.type} ({hs.name})"
def _elementary_structure(self, hs, bhs):
"""RubarBE line of a basic structure of a type of _structure_types:
(type, L, zd, zm, coefficient) for D (zm: loading elevation),
(type, L, zd, D, coefficient) for O (L: pipe length, none in
Pamhyr2: 0, a thin wall orifice)"""
trad = RunSolverTranslate()
label = self._basic_label(hs, bhs)
values = [float(v.value) for v in bhs.parameters]
if bhs.type == "S1":
width, elevation, coefficient = values
line = ("D", width, elevation, self._no_loading, coefficient)
elif bhs.type == "OR":
width, elevation, loading, coefficient, maximal = values
if maximal < self._no_loading:
raise RubarExportError(trad["structure_maximal_loading"]
.format(name=label))
line = ("D", width, elevation, loading, coefficient)
else:
diameter, elevation, siltation, coefficient = values
if siltation != 0:
raise RubarExportError(trad["structure_siltation"].format(
name=label))
line = ("O", 0.0, elevation, diameter, coefficient)
kind, size, zd, zm, coefficient = line
if (size <= 0 if kind == "D" else zm <= 0) or coefficient <= 0:
raise RubarExportError(trad["structure_invalid"].format(
name=label))
if kind == "D" and zm <= zd:
raise RubarExportError(trad["structure_loading_below"].format(
name=label))
# Fixed columns (A1, A9, 4F10.0)
if any(len(f"{v:10.3f}") > 10 for v in line[1:]):
raise RubarExportError(trad["structure_format"].format(
name=label))
return line
#########
# Saves #
#########
@@ -978,6 +1217,9 @@ class Rubar3(CommandLineSolver):
self._write_parameters(repertory, f"tnprof.{name}", lines)
def _export_devers(self, study, repertory, qlog, name="0"):
"""DEVERS, required by RubarBE when IODEV=O: per cell, weir
level, coefficient and exponent (notice annex 1). Pamhyr2 has no
lateral weir data: a zero coefficient everywhere (no overflow)"""
if qlog is not None:
qlog.put("Export DEVERS file")
@@ -999,6 +1241,42 @@ class Rubar3(CommandLineSolver):
reach_ind += 1
def _export_hydevs(self, laterals, trased, repertory, qlog, name="0"):
"""HYDEVS: per inflow, first and last cells and number of points,
then (Q, t) per point (rubar_init.f90, init2). With sediment
transport, D and S on the first line and (Q, t, Qs, ts) per point:
Pamhyr2 gives no sediment for a liquid inflow, so Qs = 0 (and
D = S = 0, the sediment of the main channel). Written, possibly
empty, so that a previous HYDEVS of the directory is not used"""
if qlog is not None:
qlog.put("Export HYDEVS file")
lines = []
for first, last, data in laterals:
if trased:
lines.append(f"{first} {last} {len(data)} 0.0 0.0\n")
lines += [f"{q!r} {t!r} 0.0 {t!r}\n" for t, q in data]
else:
lines.append(f"{first} {last} {len(data)}\n")
lines += [f"{q!r} {t!r}\n" for t, q in data]
self._write_parameters(repertory, f"hydevs.{name}", lines)
def _export_ouvrag(self, structures, repertory, qlog, name="0"):
"""OUVRAG (notice annex 1, rubar_init.f90, LOUVR1D): per
structure, upstream abscissa, reference and reach, downstream
abscissa, number of basic structures and reach, then one line per
basic structure in (A1, A9, 4F10.0). Written, possibly empty, so
that a previous OUVRAG of the directory is not used"""
if qlog is not None:
qlog.put("Export OUVRAG file")
lines = []
for rk, elementary in structures:
lines.append(f"{rk:11.3f} {-2:4d} {1:4d}\n")
lines.append(f"{rk:11.3f} {len(elementary):4d} {1:4d}\n")
for kind, size, zd, zm, coefficient in elementary:
lines.append(f"{kind}{'':9}{size:10.3f}{zd:10.3f}"
f"{zm:10.3f}{coefficient:10.3f}\n")
self._write_parameters(repertory, f"ouvrag.{name}", lines)
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)"""
@@ -1060,16 +1338,21 @@ class Rubar3(CommandLineSolver):
lines = [header(len(bc))] + records(bc, columns)
self._write_parameters(repertory, f"{filename}.{name}", lines)
def _export_hydros(self, bc, repertory, qlog, name="0"):
def _export_hydros(self, bc, tinit, 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"""
according to UNISOL, D and S are those of the condition. An
empty condition (supply at equilibrium) is written as a single
negative Qs: RubarBE then takes the transport capacity of the
upstream section at each time step (notice §3.6; rubarbe.f90,
'si Hydros.etude contient un Qs<0')"""
if qlog is not None:
qlog.put("Export HYDROS file")
d50, sigma = float(bc.d50), float(bc.sigma)
lines = [f"{len(bc)}\n"] + [
data = [(t, qs) for t, qs in bc.data] or [(tinit, -1.0)]
lines = [f"{len(data)}\n"] + [
f"{float(qs)!r} {float(t)!r} {d50!r} {sigma!r}\n"
for t, qs in bc.data
for t, qs in data
]
self._write_parameters(repertory, f"hydros.{name}", lines)
+110 -3
View File
@@ -129,10 +129,117 @@ class RunSolverTranslate(PamhyrTranslate):
"Export cancelled: several {kind} boundary conditions at node "
"{node}."
)
self._dict["boundary_empty"] = _translate(
self._dict["structure_section"] = _translate(
"Solver",
"Export cancelled: boundary condition {name} at node {node} "
"has no value."
"Export cancelled: the section of hydraulic structure {name} is "
"not a profile of reach {reach}."
)
self._dict["structure_empty_ignored"] = _translate(
"Solver",
"Hydraulic structure {name} has no enabled basic structure "
"that RubarBE takes into account: it is ignored."
)
self._dict["structure_too_many"] = _translate(
"Solver",
"Export cancelled: hydraulic structure {name} has {count} "
"enabled basic structures, RubarBE accepts at most {maximum}."
)
self._dict["structure_too_many_structures"] = _translate(
"Solver",
"Export cancelled: {count} hydraulic structures, RubarBE accepts "
"at most {maximum}."
)
self._dict["structure_same_section"] = _translate(
"Solver",
"Export cancelled: hydraulic structures {first} and {second} are "
"on the same section (rk {rk}). Put their basic structures in a "
"single hydraulic structure."
)
self._dict["structure_type_ignored"] = _translate(
"Solver",
"Basic structure {name} of type {type} has no RubarBE "
"equivalent: it is not taken into account. RubarBE accepts "
"discharge weirs (S1), rectangular orifices (OR) and circular "
"orifices (OC)."
)
self._dict["structure_maximal_loading"] = _translate(
"Solver",
"Export cancelled: rectangular orifice {name}: RubarBE has no "
"maximal loading elevation. Leave it at 9999.999."
)
self._dict["structure_siltation"] = _translate(
"Solver",
"Export cancelled: circular orifice {name}: RubarBE has no "
"siltation height. Set it to 0."
)
self._dict["structure_invalid"] = _translate(
"Solver",
"Export cancelled: basic structure {name} must have a width (or "
"diameter) > 0 and a discharge coefficient > 0."
)
self._dict["structure_loading_below"] = _translate(
"Solver",
"Export cancelled: basic structure {name}: the loading elevation "
"must be above the weir elevation."
)
self._dict["structure_format"] = _translate(
"Solver",
"Export cancelled: a value of basic structure {name} does not fit "
"the 10 columns of RubarBE (3 decimals)."
)
self._dict["lateral_empty_ignored"] = _translate(
"Solver",
"Lateral contribution {name} on reach {reach} has no value: it "
"is ignored."
)
self._dict["lateral_solid_unsupported"] = _translate(
"Solver",
"Export cancelled: solid or suspension lateral contribution "
"{name} on reach {reach}: RubarBE only reads the sediment of a "
"liquid inflow, which Pamhyr2 cannot link to it. Remove it."
)
self._dict["lateral_iodev_disabled"] = _translate(
"Solver",
"Lateral contribution {name} on reach {reach} is not taken into "
"account: lateral overflows (rubarbe_iodev) are disabled in the "
"solver parameters."
)
self._dict["lateral_type_unsupported"] = _translate(
"Solver",
"Export cancelled: lateral contribution {name} on reach {reach} "
"has no type. Choose a type (lateral sources, rain or "
"evaporation)."
)
self._dict["lateral_sections"] = _translate(
"Solver",
"Export cancelled: the sections of lateral contribution {name} "
"are not profiles of reach {reach}."
)
self._dict["lateral_not_sorted"] = _translate(
"Solver",
"Export cancelled: the times of lateral contribution {name} on "
"reach {reach} must be increasing."
)
self._dict["lateral_overlap"] = _translate(
"Solver",
"Export cancelled: lateral contributions {first} and {second} "
"on reach {reach} overlap, which RubarBE does not accept. Merge "
"them or move one of them."
)
self._dict["lateral_too_long"] = _translate(
"Solver",
"Export cancelled: the lateral contributions give {count} values "
"(number of points times number of cells of their zone), "
"RubarBE accepts at most {maximum}."
)
self._dict["boundary_empty_ignored"] = _translate(
"Solver",
"Boundary condition {name} at node {node} has no value: it is "
"ignored."
)
self._dict["boundary_missing_empty"] = _translate(
"Solver",
"Empty boundary condition(s) ignored: {names}."
)
self._dict["boundary_not_sorted"] = _translate(
"Solver",