parouby/Pamhyr2
1
0 Fork 0
mirror of https://gitlab.com/pamhyr/pamhyr2.git synced 2026-10-29 19:10:42 +01:00
Code Issues Packages Projects Releases Wiki Activity
Files
4e767d86d6734ee7cc8eb8dcfb5674e9fb71630b
Pamhyr2/src/Solver/RubarBE.py
T

1218 lines
47 KiB
Python
Raw Blame History

# 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 <https://www.gnu.org/licenses/>.
# -*- 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)
# 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))
# Defaults of RubarBE: first and last points for G and D, G and D
# for g and d. Erosion only applies strictly between g and d
def index(marker, default):
return markers.index(marker) if marker in markers else default
left = index("G", 0)
right = index("D", len(points) - 1)
limits = [("main_channel", "rg", "rd", left, right)]
if trased:
limits.append(("active_bed", "g", "d",
index("g", left), index("d", right)))
for bed, lname, rname, i, j in limits:
if j - i < 2:
raise RubarExportError(trad["bed_limits_invalid"].format(
bed=trad[bed], left=lname, right=rname,
profile=self._profile_label(profile), reach=edge.name))
return markers
def _compartments(self, edge, profile, number, point):
"""Sediment layers of a point as RubarBE compartments, from the
top: (top elevation, D, S, critical stress). A non erodible
compartment is added under an erodible last layer, so that the
layer thicknesses are kept"""
trad = RunSolverTranslate()
sl = point.sl
if sl is None or sl.is_deleted() or len(sl) == 0:
raise RubarExportError(trad["missing_sediment_layers"].format(
point=number, profile=self._profile_label(profile),
reach=edge.name))
layers = sl.layers
for layer in layers:
if not (layer.height >= 0 and layer.d50 > 0 and
layer.sigma >= 1 and layer.critical_constraint >= 0):
raise RubarExportError(trad["invalid_sediment_layer"].format(
layer=layer.name, sl=sl.name,
profile=self._profile_label(profile), reach=edge.name))
top = point.z
compartments = []
for layer in layers:
if layer.height == 0:
continue
compartments.append((top, layer.d50, layer.sigma,
layer.critical_constraint))
top -= layer.height
if not compartments or compartments[-1][3] <= 999:
last = layers[-1]
compartments.append((top, last.d50, last.sigma,
self._hard_stress))
if len(compartments) > self._max_compartments:
raise RubarExportError(trad["too_many_compartments"].format(
sl=sl.name, count=len(compartments),
maximum=self._max_compartments,
profile=self._profile_label(profile), reach=edge.name))
return compartments
@staticmethod
def _seconds(value):
if ":" in value:
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 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()]
Reference in new issue View Git Blame Copy Permalink