# Mage.py -- Pamhyr
# Copyright (C) 2023 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 numpy as np
from tools import timer, trace
from Solver.CommandLine import CommandLineSolver
from Checker.Mage import MageNetworkGraphChecker
from Model.Results.Results import Results
from Model.Results.River.River import River, Reach, Profile
logger = logging.getLogger()
def mage_file_open(filepath, mode):
f = open(filepath, mode)
if "w" in mode:
# Write header
f.write("* This file is generated by PAMHYR, please don't modify\n")
return f
class Mage(CommandLineSolver):
_type = "mage"
def __init__(self, name):
super(Mage, self).__init__(name)
self._type = "mage"
self._cmd_input = ""
self._cmd_solver = "@path @input -o @output"
self._cmd_output = ""
@classmethod
def default_parameters(cls):
lst = super(Mage, cls).default_parameters()
lst += [
("mage_min_timestep", "1.0"),
("mage_timestep_tra", "3600"),
("mage_timestep_bin", "0"),
# ("mage_timestep_melissa", "0"),
("mage_implicitation", "0.70"),
("mage_continuity_discretization", "S"),
("mage_qsj_discretization", "B"),
("mage_stop_criterion_iterations", "R"),
("mage_iteration_type", "0"),
("mage_smooth_coef", "0"),
("mage_cfl_max", "-1."),
("mage_min_height", "0.1"),
("mage_max_niter", "10"),
("mage_timestep_reduction_factor", "2"),
("mage_precision_reduction_factor_Z", "1"),
("mage_precision_reduction_factor_Q", "1"),
("mage_niter_max_precision", "99"),
("mage_niter_before_switch", "99"),
("mage_max_froude", "1.5"),
("mage_diffluence_node_height_balance", "-1"),
("mage_compute_reach_volume_balance", "y"),
("mage_max_reach_volume_balance", "0.001"),
("mage_min_reach_volume_to_check", "1000.0"),
("mage_init_internal", " "),
]
return lst
@classmethod
def checkers(cls):
lst = [
MageNetworkGraphChecker(connectivity=True),
MageNetworkGraphChecker(connectivity=False)
]
return lst
##########
# Export #
##########
def cmd_args(self, study):
lst = super(Mage, self).cmd_args(study)
lst.append("-r")
return lst
def input_param(self):
name = self._study.name
return f"{name}.REP"
def output_param(self):
name = self._study.name
return f"{name}.BIN"
def log_file(self):
name = self._study.name
return f"{name}.TRA"
@timer
def _export_ST(self, study, repertory, qlog, name="0"):
files = []
if qlog is not None:
qlog.put("Export ST file")
os.makedirs(os.path.join(repertory, "net"), exist_ok=True)
# Write header
edges = study.river.edges()
edges = list(
filter(
lambda e: e.is_enable(),
edges
)
)
for edge in edges:
name = edge.name.replace(" ", "_")
if edge._name == "":
name = f"Reach_{edge.id}"
with mage_file_open(
os.path.join(repertory, "net", f"{name}.ST"),
"w+"
) as f:
files.append(str(os.path.join("net", f"{name}.ST")))
cnt_num = 1
for profile in edge.reach.profiles:
self._export_ST_profile_header(
f, files, profile, cnt_num
)
cnt_num += 1
# Points
for point in profile.points:
self._export_ST_point_line(
f, files, point
)
# Profile last line
f.write(f" 999.9990 999.9990 999.9990\n")
return files
def _export_ST_profile_header(self, wfile, files,
profile, cnt):
num = f"{cnt:>6}"
c1 = f"{profile.code1:>6}"
c2 = f"{profile.code2:>6}"
t = f"{len(profile.points):>6}"
kp = f"{profile.kp:>12f}"[0:12]
pname = profile.name
if profile.name == "":
# Generate name from profile id prefixed with
# 'p' (and replace space char with '0' char)
pname = f"p{profile.id:>3}".replace(" ", "0")
name = f"{pname:<19}"
# Generate sediment additional data if available
sediment = ""
if profile.sl is not None:
if not any(filter(lambda f: ".GRA" in f, files)):
files.append(self._gra_file)
# Number of layers
nl = len(profile.sl)
sediment = f" {nl:>3}"
# Layers data
for layer in profile.sl.layers:
sediment += (
f" {layer.height:>10} {layer.d50:>10} " +
f"{layer.sigma:>10} " +
f"{layer.critical_constraint:>10}"
)
# Profile header line
wfile.write(f"{num}{c1}{c2}{t} {kp} {pname} {sediment}\n")
def _export_ST_point_line(self, wfile, files, point):
x = f"{point.x:<12.4f}"[0:12]
y = f"{point.y:<12.4f}"[0:12]
z = f"{point.z:<12.4f}"[0:12]
n = f"{point.name:<3}"
# Generate sediment additional data if available
sediment = ""
prev = point.z
if point.sl is not None:
# Number of layers
nl = len(point.sl)
sediment = f"{nl:>3}"
# Layers data
for layer in point.sl.layers:
prev = round(prev - layer.height, 5)
sediment += (
f" {prev:>10} {layer.d50:>10} " +
f"{layer.sigma:>10} " +
f"{layer.critical_constraint:>10}"
)
# Point line
wfile.write(f"{x} {y} {z} {n} {sediment}\n")
@timer
def _export_BC(self, t, bounds, repertory, qlog, name="0"):
files = []
if len(bounds) == 0:
return files
if qlog is not None:
qlog.put(f"Export {t} file")
with mage_file_open(os.path.join(repertory, f"{name}.{t}"), "w+") as f:
files.append(f"{name}.{t}")
for bound in bounds:
name = f"{bound.node.id:3}".replace(" ", "x")
f.write(f"* {bound.node.name} ({name}) {bound.bctype}\n")
f.write(f"${name}\n")
header = bound.header
f.write(f"*{header[0]:>9}|{header[1]:>10}\n")
for d in bound.data:
v0 = d[0]
v1 = d[1]
if t in ["HYD", "QSO", "LIM"]:
v0 /= 60 # Convert first column to minute
f.write(f"{v0:10}{v1:10}\n")
return files
@timer
def _export_bound_cond(self, study, repertory, qlog, name="0"):
files = []
lst = study.river.boundary_condition
AVA = []
HYD = []
LIM = []
QSO = []
for tab in ["liquid", "solid", "suspenssion"]:
for bound in lst.get_tab(tab):
if bound.node is None:
continue
if bound.bctype == "ZD":
AVA.append(bound)
elif bound.bctype == "TD" or bound.bctype == "PC":
HYD.append(bound)
elif bound.bctype == "TZ":
LIM.append(bound)
elif bound.bctype == "SL":
QSO.append(bound)
files = files + self._export_BC("AVA", AVA, repertory, qlog, name=name)
files = files + self._export_BC("HYD", HYD, repertory, qlog, name=name)
files = files + self._export_BC("LIM", LIM, repertory, qlog, name=name)
files = files + self._export_QSO(QSO, repertory, qlog, name=name)
return files
# @timer
# def _export_LC(self, lateral, repertory, qlog, name="0"):
# files = []
# if qlog is not None:
# qlog.put(f"Export LAT file")
# with mage_file_open(
# os.path.join(repertory, f"{name}.LAT"),
# "w+"
# ) as f:
# files.append(f"{name}.LAT")
# name = f"{lateral.node.id:3}".replace(" ", "x")
# f.write(f"* {lateral.node.name} ({name}) {lateral.bctype}\n")
# f.write(f"${name}\n")
# header = lateral.header
# f.write(f"*{header[0]:>9}|{header[1]:>10}\n")
# for d in lateral.data:
# f.write(f"{d[0]:1{name}.3f}{d[1]:10.3f}\n")
# return files
# @timer
# def _export_lateral_contrib(self, study, repertory, qlog, name="0"):
# files = []
# lst = study.river.lateral_contribution
# for tab in ["liquid", "solid", "suspenssion"]:
# for lateral in lst.get_tab(tab):
# files = files + self._export_LC(lateral, repertory, qlog)
# return files
@timer
def _export_RUG(self, study, repertory, qlog, name="0"):
files = []
if qlog is not None:
qlog.put("Export RUG file")
# Write header
with mage_file_open(os.path.join(repertory, f"{name}.RUG"), "w+") as f:
files.append(f"{name}.RUG")
edges = study.river.edges()
edges = list(
filter(
lambda e: e.is_enable(),
edges
)
)
id = 1
for edge in edges:
frictions = edge.frictions
for friction in frictions.frictions:
num = f"{id:>3}"
bkp = f"{friction.begin_kp:>10.3f}"
ekp = f"{friction.end_kp:>10.3f}"
# if friction.begin_kp != friction.end_kp:
# print("TODO")
strickler = friction.begin_strickler
coef_1 = f"{strickler.minor:>10.3f}"
coef_2 = f"{strickler.medium:>10.3f}"
f.write(f"K{num} {bkp}{ekp}{coef_1}{coef_2}\n")
id += 1
return files
@timer
def _export_INI(self, study, repertory, qlog, name="0"):
files = []
if qlog is not None:
qlog.put("Export INI file")
# Write header
with mage_file_open(os.path.join(repertory, f"{name}.INI"), "w+") as f:
has_ini = False
id = 1
reachs = study.river.enable_edges()
# TODO put real date...
f.write(f"$ date en minutes : 0.00\n")
f.write(f"* IB IS discharge elevation kp\n")
id = 1
for reach in reachs:
cond = study.river.initial_conditions.get(reach)
data = cond.data
if len(data) == 0:
continue
has_ini = True
id_sec = 1
for d in data:
IR = f"{id}"
IS = f"{id_sec}"
discharge = f"{d['discharge']:>10.5f}"
elevation = f"{d['elevation']:>11.6f}"
kp = f"{d['kp']:>9.2f}"
f.write(f"{IR} {IS} {discharge} {elevation} {kp}\n")
id_sec += 1
id += 1
if has_ini:
files.append(f"{name}.INI")
return files
@timer
def _export_CAS(self, study, repertory, qlog, name="0"):
files = []
reservoirs = study.river.reservoir.lst
if len(reservoirs) == 0:
return files
if qlog is not None:
qlog.put("Export CAS file")
with mage_file_open(os.path.join(repertory, f"{name}.CAS"), "w+") as f:
files.append(f"{name}.CAS")
for reservoir in reservoirs:
reservoir.sort()
node = reservoir.node
name = f"{node.id:3}".replace(" ", "x")
f.write(f"* {node.name} ({name}) Reservoir\n")
f.write(f"${name}\n")
f.write(f"*{'Elev(m)':>9}|{'Area(ha)':>10}\n")
for d in reservoir.data:
v0 = d[0]
v1 = d[1]
f.write(f"{v0:>10.3f}{v1:>10.3f}\n")
return files
@timer
def _export_SIN(self, study, repertory, qlog, name="0"):
files = []
sin_dict = {
"ND": "*",
"S1": "D", "S2": "T", "S3": "T",
"OR": "O", "OC": "B", "OV": "F",
"V1": "V", "V2": "W",
"BO": "A",
"UD": "X",
"PO": "P",
}
hydraulic_structures = study.river.hydraulic_structures.lst
if len(hydraulic_structures) == 0:
return files
if qlog is not None:
qlog.put("Export SIN file")
with mage_file_open(os.path.join(repertory, f"{name}.SIN"), "w+") as f:
files.append(f"{name}.SIN")
for hs in hydraulic_structures:
if not hs.input_reach.is_enable():
continue
f.write(
'* ouvrage au pk ' +
f"{hs.input_kp:>12.1f}" + ' ' +
hs.name + '\n'
)
for bhs in hs.basic_structures:
reach_id = study.river.get_edge_id(hs.input_reach) + 1
param_str = ' '.join(
[
f'{p:>10.3f}'
for p in self._export_SIN_parameters(bhs)
]
)
f.write(
f"{sin_dict[bhs._type]} " +
f"{reach_id} {hs.input_kp:>12.3f} {param_str} " +
f"{bhs.name}\n"
)
return files
def _export_SIN_parameters(self, bhs):
res = [9999.999] * 5
if len(bhs) == 5:
res = self._export_SIN_parameters_5(bhs)
elif len(bhs) == 4:
res = self._export_SIN_parameters_4(bhs)
elif len(bhs) == 3:
res = self._export_SIN_parameters_3(bhs)
return res
def _export_SIN_parameters_5(self, bhs):
# S2, OR, V1, V2, UD
return [
bhs._data[0].value,
bhs._data[1].value,
bhs._data[2].value,
bhs._data[3].value,
bhs._data[4].value,
]
def _export_SIN_parameters_4(self, bhs):
# S3, OC
res = [
bhs._data[0].value,
bhs._data[1].value,
bhs._data[2].value,
bhs._data[3].value,
0.0,
]
if bhs._type == "T": # S3
res = [0.0] + res[:-1]
return res
def _export_SIN_parameters_3(self, bhs):
# S1, BO
if bhs._type == "S1":
res = [
bhs._data[0].value,
bhs._data[1].value,
0.0,
bhs._data[2].value,
9999.99,
]
else:
res = [
bhs._data[0].value,
bhs._data[1].value,
bhs._data[2].value,
0.0,
0.0,
]
return res
@timer
def _export_DEV(self, study, repertory, qlog, name="0"):
files = []
if qlog is not None:
qlog.put("Export DEV file")
with mage_file_open(
os.path.join(
repertory, f"{name}.DEV"
), "w+"
) as f:
reachs = study.river.enable_edges()
id = 1
for reach in reachs:
f.write(f"YD{id:3}\n")
f.write(f"YG{id:3}\n")
id += 1
files.append(f"{name}.DEV")
return files
@timer
def _export_REP(self, study, repertory, files, qlog, name="0"):
if qlog is not None:
qlog.put("Export REP file")
# Write header
with mage_file_open(
os.path.join(
repertory, f"{name}.REP"
), "w+"
) as f:
f.write("confirmation=non\n")
for file in files:
EXT = file.split('.')[1]
if EXT not in ["ST", "GRA"]:
f.write(f"{EXT} {file}\n")
f.write("* OUTPUT\n")
f.write(f"TRA {name}.TRA\n")
f.write(f"BIN {name}.BIN\n")
for file in files:
EXT = file.split('.')[1]
if EXT in ["GRA"]:
f.write(f"{EXT} {file}\n")
@timer
def export(self, study, repertory, qlog=None):
self._study = study
name = study.name.replace(" ", "_")
# Define GRA file name
self._gra_file = f"{name}.GRA"
self._bin_file = f"{name}.BIN"
self._export_ST(study, repertory, qlog, name=name)
return True
###########
# RESULTS #
###########
def read_bin(self, study, repertory, results, qlog=None, name="0"):
return
@timer
def results(self, study, repertory, qlog=None, name="0"):
results = Results(
study=study,
solver=self,
repertory=repertory,
name=name,
)
self.read_bin(study, repertory, results, qlog, name=name)
return results
##########
# MAGE 7 #
##########
class Mage7(Mage):
_type = "mage7"
def __init__(self, name):
super(Mage7, self).__init__(name)
self._type = "mage7"
@classmethod
def default_parameters(cls):
lst = super(Mage7, cls).default_parameters()
return lst
##########
# MAGE 8 #
##########
class Mage8(Mage):
_type = "mage8"
def __init__(self, name):
super(Mage8, self).__init__(name)
self._type = "mage8"
@classmethod
def default_parameters(cls):
lst = super(Mage8, cls).default_parameters()
# Insert new parameters at specific position
names = list(map(lambda t: t[0], lst))
i = names.index("mage_precision_reduction_factor_Q")
lst.insert(i+1, ("mage_precision_reduction_factor_r", "1"))
# Mage parameter for sediment module (added in DB 0.0.4)
lst.append(("mage_sediment_masse_volumique", "2650.0"))
lst.append(("mage_sediment_angle_repos", "40.0"))
lst.append(("mage_sediment_porosity", "0.40"))
lst.append(("mage_distance_Han", "0.0"))
lst.append(("mage_distance_chargement_d50", "100.0"))
lst.append(("mage_distance_chargement_sigma", "100.0"))
lst.append(("mage_methode_modification_geometrie", "1"))
lst.append(("mage_shields_critique", "1"))
lst.append(("mage_shields_correction", "1"))
lst.append(("mage_capacite_solide", "1"))
lst.append(("mage_pas_de_temps_charriage", "1"))
lst.append(("mage_facteur_multiplicateur", "1.0"))
return lst
##########
# Export #
##########
def cmd_args(self, study):
lst = super(Mage8, self).cmd_args(study)
if study.river.has_sediment():
lst.append("-c=3")
return lst
@timer
def _export_PAR(self, study, repertory, qlog=None, name="0"):
files = []
if qlog is not None:
qlog.put("Export PAR file")
with mage_file_open(os.path.join(repertory, f"{name}.PAR"), "w+") as f:
files.append(f"{name}.PAR")
params = study.river.get_params(self.type).parameters
for p in params:
name = p.name\
.replace("all_", "")\
.replace("mage_", "")
value = p.value
if name in ["command_line_arguments"]:
continue
if name == "compute_reach_volume_balance":
value = "O" if value.lower() == "y" else "N"
if name == "init_internal":
value = ("p" if value.lower() in ["y", "yes", "true"]
else "")
f.write(f"{name} {value}\n")
return files
@timer
def _export_NET(self, study, repertory, qlog=None, name="0"):
files = []
if qlog is not None:
qlog.put("Export NET file")
with mage_file_open(os.path.join(repertory, f"{name}.NET"), "w+") as f:
files.append(f"{name}.NET")
edges = study.river.edges()
edges = list(
filter(
lambda e: e.is_enable(),
edges
)
)
for e in edges:
name = e.name.replace(" ", "_")
if e._name == "":
name = f"Reach_{e.id}"
id = f"Bief_{e.id+1}"
n1 = f"{e.node1.id:3}".replace(" ", "x")
n2 = f"{e.node2.id:3}".replace(" ", "x")
file = os.path.join("net", name + ".ST")
f.write(f"{id} {n1} {n2} {file}\n")
return files
@timer
def _export_QSO(self, bounds, repertory, qlog, name="0"):
files = []
if len(bounds) == 0:
return files
if qlog is not None:
qlog.put(f"Export QSO file")
with mage_file_open(os.path.join(repertory, f"{name}.QSO"), "w+") as f:
files.append(f"{name}.QSO")
for bound in bounds:
# File header
name = f"{bound.node.id:3}".replace(" ", "x")
f.write(f"* {bound.node.name} ({name}) {bound.bctype}\n")
d50 = bound.d50
sigma = bound.sigma
if len(bound.data) == 0:
f.write(f"${name} {d50} {sigma} default\n")
else:
f.write(f"${name} {d50} {sigma}\n")
# Table header
header = bound.header
f.write(f"*{header[0]:>9}|{header[1]:>10}\n")
# Data
for d in bound.data:
f.write(f"{d[0]:10.3f}{d[1]:10.3f}\n")
return files
@timer
def export(self, study, repertory, qlog=None, name="0"):
self._study = study
name = study.name.replace(" ", "_")
# Define GRA file name
self._gra_file = f"{name}.GRA"
self._bin_file = f"{name}.BIN"
# Generate files
files = []
files = self._export_ST(study, repertory, qlog, name=name)
files = files + self._export_PAR(study, repertory, qlog, name=name)
files = files + self._export_NET(study, repertory, qlog, name=name)
files = files + \
self._export_bound_cond(study, repertory, qlog, name=name)
files = files + self._export_RUG(study, repertory, qlog, name=name)
files = files + self._export_INI(study, repertory, qlog, name=name)
files = files + self._export_SIN(study, repertory, qlog, name=name)
files = files + self._export_CAS(study, repertory, qlog, name=name)
files = files + self._export_DEV(study, repertory, qlog, name=name)
self._export_REP(study, repertory, files, qlog, name=name)
return True
###########
# RESULTS #
###########
@timer
def read_bin(self, study, repertory, results, qlog=None, name="0"):
fname = os.path.join(repertory, f"{name}.BIN")
logger.info(f"read_bin: Start reading '{fname}' ...")
with mage_file_open(fname, "r") as f:
def newline(): return np.fromfile(f, dtype=np.int32, count=1)
def endline(): return np.fromfile(f, dtype=np.int32, count=1)
def read_int(size): return np.fromfile(
f, dtype=np.int32, count=size)
def read_float(size): return np.fromfile(
f, dtype=np.float32, count=size)
def read_float64(size): return np.fromfile(
f, dtype=np.float64, count=size)
# Meta data (1st line)
newline()
data = read_int(3)
nb_reach = data[0]
nb_profile = data[1]
mage_version = data[2]
logger.debug(f"read_bin: nb_reach = {nb_reach}")
logger.debug(f"read_bin: nb_profile = {nb_profile}")
logger.debug(f"read_bin: mage_version = {mage_version}")
if mage_version <= 80:
msg = (
"Read BIN files: " +
f"Possible incompatible mage version '{mage_version}', " +
"please check your solver configuration..."
)
logger.warning(msg)
if qlog is not None:
qlog.put("[WARNING] " + msg)
results.set("solver_version", f"Mage8 ({mage_version})")
results.set("nb_reach", f"{nb_reach}")
results.set("nb_profile", f"{nb_profile}")
endline()
# Reach information (2nd line)
newline()
reachs = []
iprofiles = {}
reach_offset = {}
data = read_int(2*nb_reach)
for i in range(nb_reach):
# Add results reach to reach list
r = results.river.add(i)
reachs.append(r)
# ID of first and last reach profiles
i1 = data[2*i] - 1
i2 = data[2*i+1] - 1
# Add profile id correspondance to reach
key = (i1, i2)
iprofiles[key] = r
# Profile ID offset
reach_offset[r] = i1
logger.debug(f"read_bin: iprofiles = {iprofiles}")
endline()
# X (3rd line)
newline()
_ = read_float(nb_profile)
endline()
# Z and Y (4th line)
newline()
_ = read_float(3*nb_profile)
endline()
# Data
newline()
def ip_to_r(i): return iprofiles[
next(
filter(
lambda k: k[0] <= i <= k[1],
iprofiles
)
)
]
def ip_to_ri(r, i): return i - reach_offset[r]
ts = set()
end = False
while not end:
n = read_int(1)[0]
timestamp = read_float64(1)[0]
key = bytearray(np.fromfile(
f, dtype=np.byte, count=1)).decode()
data = read_float(n)
logger.debug(f"read_bin: timestamp = {timestamp} sec")
ts.add(timestamp)
if key in ["Z", "Q"]:
for i, d in enumerate(data):
# Get reach corresponding to profile ID
reach = ip_to_r(i)
# Get profile id in reach
ri = ip_to_ri(reach, i)
# Set data for profile RI
reach.set(ri, timestamp, key, d)
if key == "Z":
profile = study.river\
.current_reach()\
.reach.profile(ri)
ptX, ptY = profile.get_water_limits(d)
reach.set(ri, timestamp, "ptX", ptX)
reach.set(ri, timestamp, "ptY", ptY)
endline()
end = newline().size <= 0
logger.debug(reachs[0].profiles[0]._data)
results.set("timestamps", ts)
logger.info(f"read_bin: ... end with {len(ts)} timestamp read")
@timer
def read_gra(self, study, repertory, results, qlog=None, name="0"):
if not study.river.has_sediment():
return
fname = os.path.join(repertory, f"{name}.GRA")
logger.info(f"read_gra: Start reading '{fname}' ...")
with mage_file_open(fname, "r") as f:
def newline(): return np.fromfile(f, dtype=np.int32, count=1)
def endline(): return np.fromfile(f, dtype=np.int32, count=1)
def read_int(size): return np.fromfile(
f, dtype=np.int32, count=size)
def read_float(size): return np.fromfile(
f, dtype=np.float32, count=size)
def read_float64(size): return np.fromfile(
f, dtype=np.float64, count=size)
# Meta data (1st line)
newline()
data = read_int(3)
nb_reach = data[0]
nb_profile = data[1]
mage_version = data[2]
logger.debug(f"read_gra: nb_reach = {nb_reach}")
logger.debug(f"read_gra: nb_profile = {nb_profile}")
logger.debug(f"read_gra: mage_version = {mage_version}")
if mage_version < 80:
msg = (
"Read GRA files: " +
f"Possible incompatible mage version '{mage_version}', " +
"please check your solver configuration..."
)
logger.warning(msg)
if qlog is not None:
qlog.put("[WARNING] " + msg)
results.set("gra_solver_version", f"Mage8 ({mage_version})")
results.set("gra_nb_reach", f"{nb_reach}")
results.set("gra_nb_profile", f"{nb_profile}")
endline()
# Reach information (2nd line)
newline()
reachs = []
iprofiles = {}
reach_offset = {}
data = read_int(2*nb_reach)
for i in range(nb_reach):
# Get results reach to reach list
r = results.river.reach(i)
reachs.append(r)
# ID of first and last reach profiles
i1 = data[2*i] - 1
i2 = data[2*i+1] - 1
# Add profile id correspondance to reach
key = (i1, i2)
iprofiles[key] = r
# Profile ID offset
reach_offset[r] = i1
logger.debug(f"read_gra: iprofiles = {iprofiles}")
endline()
# X (3rd line)
newline()
_ = read_float(nb_profile)
endline()
# npts (4th line)
newline()
_ = read_int(nb_profile)
endline()
# Data
def ip_to_r(i): return iprofiles[
next(
filter(
lambda k: k[0] <= i <= k[1],
iprofiles
)
)
]
def ip_to_ri(r, i): return i - reach_offset[r]
ts = set()
end = False
newline()
while not end:
n = read_int(1)[0]
timestamp = read_float64(1)[0]
with_bedload = read_int(1)[0]
logger.debug(f"read_gra: Number of cross section: {n}")
logger.debug(f"read_gra: Timestamp: {timestamp}")
logger.debug(f"read_gra: Type of bedload: {with_bedload}")
endline()
npts = [1] * n
if with_bedload == 1:
newline()
npts = read_int(n)
endline()
sum_npts = sum(npts)
logger.debug(f"read_gra: Number of points: {sum_npts}")
newline()
nsl = read_int(sum_npts)
logger.debug(f"read_gra: Number of sedimentary layers: {nsl}")
endline()
newline()
data = read_float64(sum(nsl) * 3)
endline()
ts.add(timestamp)
i_pts = 0
i_data = 0
# Loop on cross section
for i in range(n):
sec_sl = []
reach = ip_to_r(i)
p_i = ip_to_ri(reach, i)
# Loop on cross section points
for j in range(npts[i]):
sl = []
# Loop on sediment layers
for k in range(nsl[i_pts]):
h = data[i_data]
d50 = data[i_data + 1]
sigma = data[i_data + 2]
sl.append((h, d50, sigma))
i_data += 3
i_pts += 1
sec_sl.append(sl)
reach.set(p_i, timestamp, "sl", sec_sl)
logger.debug(
f"read_gra: data size = {len(data)} ({i_data} readed)"
)
end = newline().size <= 0
results.set("sediment_timestamps", ts)
logger.info(f"read_gra: ... end with {len(ts)} timestamp read")
@timer
def results(self, study, repertory, qlog=None):
self._study = study
name = study.name.replace(" ", "_")
results = super(Mage8, self).results(study, repertory, qlog, name=name)
self.read_gra(study, repertory, results, qlog, name=name)
return results