diff --git a/src/Model/Results/River/River.py b/src/Model/Results/River/River.py index 6494ff62..253e8e44 100644 --- a/src/Model/Results/River/River.py +++ b/src/Model/Results/River/River.py @@ -104,7 +104,10 @@ class Profile(SQLSubModel): return any(map(lambda ts: "sl" in self._data[ts], self._data)) def has_bedload(self): - return any(map(lambda ts: "sl" in self._data[ts], self._data)) + # Sediment layers (Mage) or bed level of the section (RubarBE) + return any(map(lambda ts: ("sl" in self._data[ts] or + "zfd" in self._data[ts]), + self._data)) @classmethod def _db_create(cls, execute, ext=""): diff --git a/src/Solver/RubarBE.py b/src/Solver/RubarBE.py index 8ded5fb2..dc43da93 100644 --- a/src/Solver/RubarBE.py +++ b/src/Solver/RubarBE.py @@ -1113,10 +1113,16 @@ class Rubar3(CommandLineSolver): else: return (0.0, 0.0, 0.0, 0.0) - def set_and_compute_limites(reach, ind, z, q, s): + # Bed level of each section (mean of the CTDF of the + # neighbouring cells): {(reach, ind): {timestamp: level}} + bed = {} + + def set_and_compute_limites(reach, ind, z, q, s, zf=None): reach.set(ind, timestamp, "Z", z) reach.set(ind, timestamp, "Q", q) reach.set(ind, timestamp, "V", s) + if zf is not None: + bed.setdefault((reach, ind), {})[timestamp] = zf profile = reach.profile(ind) limits = profile.geometry.get_water_limits(z) @@ -1164,7 +1170,7 @@ class Rubar3(CommandLineSolver): 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) + set_and_compute_limites(reach, 0, z+h, q, s, z) # For each profile ind = 1 @@ -1175,7 +1181,8 @@ class Rubar3(CommandLineSolver): reach, ind, (pz + z + ph + h) / 2, (pq + q) / 2, - (ps + s) / 2 + (ps + s) / 2, + (pz + z) / 2 ) ph = h @@ -1186,7 +1193,7 @@ class Rubar3(CommandLineSolver): # Last profile h, s, q, z = read_data_line(f) - set_and_compute_limites(reach, ind, z+h, q, s) + set_and_compute_limites(reach, ind, z+h, q, s, z) # Tables and plots use timestamp rows and global profile columns. timestamps = sorted(ts) @@ -1196,6 +1203,21 @@ class Rubar3(CommandLineSolver): profiles.extend(reach.profiles) reach.set_global_index(range(offset, len(profiles))) + # Kept only when the bed moves (sediment transport), so that + # the plots otherwise show the study geometry. The initial + # condition of the study has no bed level: the first one read + # in PROFIL is used + keys = ["Z", "Q", "V"] + moves = any(max(levels.values()) != min(levels.values()) + for levels in bed.values()) + if moves and len(bed) == len(profiles): + for (reach, ind), levels in bed.items(): + first = levels[min(levels)] + for t in timestamps: + first = levels.get(t, first) + reach.set(ind, t, "zfd", first) + keys.append("zfd") + results.set("timestamps", timestamps) results.compute_timestamp_index() results.set("table", { @@ -1203,7 +1225,7 @@ class Rubar3(CommandLineSolver): [profile.get_ts_key(t, key) for profile in profiles] for t in timestamps ], dtype=float) - for key in ("Z", "Q", "V") + for key in keys }) @timer @@ -1321,8 +1343,74 @@ class Rubar3(CommandLineSolver): logger_exception(e) return None + # Without them, the hydraulic results stay readable + try: + self.read_bed_sections(results, repertory, name) + except Exception as e: + logger.warning(f"Failed to read the bed sections") + logger_exception(e) + return results + def read_bed_sections(self, results, repertory, name): + """Bed sections for the cross section plot of the results window: + - 'bed_points' {global index: [bool, ...]}: the points whose level + changes, strictly between g and d (step 4); + - 'bed_sections' {timestamp: {global index: [z, ...]}}: exact + sections saved at the times of TNPROF when IOVIS=O ('geomac-N', + GEOMAC-I format; rubarbe.f90, SAVEVISU). + Between two exact sections, the plot shifts the points of the last + one by the change of the section bed level 'zfd' since then""" + if "zfd" not in results.get("table") or len(results.river) != 1: + return + profiles = results.river.reach(0).profiles + + moving = {} + for profile in profiles: + markers = [self._point_markers.get(p.name.strip(), " ") + for p in profile.geometry.points] + _, _, left, right = self._bed_limits(markers) + moving[profile.global_index] = [ + left < i < right for i in range(len(markers)) + ] + results.set("bed_points", moving) + + timestamps = results.get("timestamps") + sections = {} + pattern = re.compile(rf"geomac-\d+\.{re.escape(name)}") + for fname in sorted(os.listdir(repertory)): + if not pattern.fullmatch(fname): + continue + t, levels = self._read_geomac(os.path.join(repertory, fname)) + if (len(levels) != len(profiles) or + any(len(z) != len(p.geometry.points) + for z, p in zip(levels, profiles))): + logger.warning(f"{fname}: sections differ from the study " + "geometry, ignored") + continue + # Same save as in PROFIL, with fewer digits (F15.3) + timestamp = min(timestamps, key=lambda ts: abs(ts - t)) + if abs(timestamp - t) > 0.01: + continue + sections[timestamp] = { + p.global_index: np.array(z) for p, z in zip(profiles, levels) + } + results.set("bed_sections", sections) + + @staticmethod + def _read_geomac(path): + """Time and level of each point of each section of a GEOMAC file + (GEOMAC-I format: top of the first compartment)""" + with open(path) as f: + lines = iter([line for line in f if line.strip()]) + count, time = next(lines).split()[:2] + levels = [] + for _ in range(int(count)): + points = int(next(lines).split()[2]) + levels.append([float(next(lines)[1:].split()[1]) + for _ in range(points)]) + return float(time), levels + class RubarBE(Rubar3): _type = "rubarbe" diff --git a/src/View/Results/PlotAC.py b/src/View/Results/PlotAC.py index 0487177a..b2e59fee 100644 --- a/src/View/Results/PlotAC.py +++ b/src/View/Results/PlotAC.py @@ -286,6 +286,10 @@ class PlotAC(PamhyrPlot): initial_z = np.asarray(profile.geometry.z(), dtype=float) result = self.results[self._current_res_id] + bed = self._get_bed_from_sections(profile, result, initial_z) + if bed is not None: + return bed + try: sediment_timestamps = sorted(result.get("sediment_timestamps")) except (KeyError, TypeError): @@ -324,6 +328,38 @@ class PlotAC(PamhyrPlot): return initial_z + def _get_bed_from_sections(self, profile, result, initial_z): + """Bed of the solvers that give the bed level of each section at + each time ('zfd') and exact sections at some times + ('bed_sections', RubarBE): the last exact section before the + current time, or the study geometry, with its moving points + ('bed_points') shifted by the change of the bed level since then. + None for the other results""" + try: + sections = result.get("bed_sections") + moving = result.get("bed_points")[profile.global_index] + except KeyError: + return None + + current = profile.get_ts_key(self._current_timestamp, "zfd") + if current is None: + return None + + index = profile.global_index + saved = [ts for ts in sorted(sections) + if ts <= self._current_timestamp and index in sections[ts]] + if saved: + base_timestamp = saved[-1] + base = np.asarray(sections[base_timestamp][index], dtype=float) + else: + base_timestamp = min(result.get("timestamps")) + base = initial_z + + level = profile.get_ts_key(base_timestamp, "zfd") + if level is None or len(base) != len(moving): + return None + return base + np.where(moving, current - level, 0.0) + def _get_wet_mask(self, profile, elevation): result = self.results[self._current_res_id]