RubarBE: show bed evolution in the results window (cross section and longitudinal profile), from the bed level in profil and the sections of geomac-N

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=""):

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"

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]