src/Model/Results/Results.py

@@ -219,11 +219,6 @@ class Results(SQLSubModel):
             qlog=None,
         )
 
-    def bufferize(self, key):
-        if self.is_valid:
-            for reach in self._river._reachs:
-                reach.bufferize(self._meta_data["timestamps"], key)
-
     def timestamps_to_struct(self):
         ts = self._meta_data["timestamps"]
         sf = ">" + ''.join(itertools.repeat("d", len(ts)))

src/Model/Results/River/River.py

@@ -17,7 +17,6 @@
 import struct
 import logging
 import itertools
-import numpy as np
 
 from tools import flatten
 from functools import reduce
@@ -77,12 +76,9 @@ class Profile(SQLSubModel):
                 return self._data[timestamp][key]
         return None
 
-    def has_sediment_layers(self):
+    def has_sediment(self):
         return any(map(lambda ts: "sl" in self._data[ts], self._data))
 
-    def has_bedload(self):
-        return any(map(lambda ts: "zfd" in self._data[ts], self._data))
-
     @classmethod
     def _db_create(cls, execute, ext=""):
         execute(f"""
@@ -279,7 +275,6 @@ class Reach(SQLSubModel):
                 lambda p: p.name[0:8] != 'interpol', self._profiles
             )
         )
-        self._buffers = {}
 
     def __len__(self):
         return len(self._profiles)
@@ -300,10 +295,6 @@ class Reach(SQLSubModel):
     def profile_mask(self):
         return self._profile_mask
 
-    @property
-    def buffers(self):
-        return self._buffers
-
     def profile(self, id):
         return self._profiles[id]
 
@@ -311,16 +302,7 @@ class Reach(SQLSubModel):
         self._profiles[profile_id].set(timestamp, key, data)
 
     def has_sediment(self):
-        return any(map(lambda profile: profile.has_sediment_layers(),
+        return any(map(lambda profile: profile.has_sediment(), self._profiles))
-                       self._profiles))
-
-    def has_bedload(self):
-        return any(map(lambda profile: profile.has_bedload(), self._profiles))
-
-    def bufferize(self, timestamps, key):
-        self._buffers[key] = np.zeros((len(timestamps), len(self)))
-        for i, p in enumerate(self._profiles):
-            self._buffers[key][:, i] = p.get_key(key)
 
     @classmethod
     def _db_create(cls, execute, ext=""):

src/Solver/Mage.py

@@ -19,7 +19,6 @@
 import os
 import logging
 import numpy as np
-from functools import reduce
 
 from itertools import chain
 from tools import timer, trace, logger_exception
@@ -1197,10 +1196,6 @@ class Mage8(Mage):
                         r.set(i, t, "V", v)
             logger.info(f"read_bin: ... end with {len(ts)} timestamp read")
 
-        results.bufferize("Z")
-        results.bufferize("Q")
-        results.bufferize("V")
-
     @timer
     def read_gra(self, study, repertory, results, qlog=None, name="0"):
         if not study.river.has_sediment():
@@ -1367,40 +1362,6 @@ class Mage8(Mage):
                 )
                 end = newline().size <= 0
 
-            ts_list = sorted(ts)
-
-            logger.info(f"compute river bed elevation...")
-
-            for r in reachs:
-                z_min = reach.geometry.get_z_min()
-                sls = list(map(
-                    lambda p: p.get_ts_key(ts_list[0], "sl")[0],
-                    r.profiles
-                    ))
-                z_br = list(map(
-                            lambda z, sl: reduce(
-                                lambda z, h: z - h[0],
-                                sl, z
-                            ),
-                            z_min,          # Original geometry
-                            sls             # Original sediment layers
-                            ))
-                for t in ts_list:
-                    sls = list(map(
-                        lambda p: p.get_ts_key(t, "sl")[0],
-                        r.profiles
-                        ))
-                    zfd = list(map(
-                                lambda z, sl: reduce(
-                                    lambda z, h: z + h[0],
-                                    sl, z
-                                ),
-                                z_br,       # bedrock
-                                sls         # Original sediment layers
-                            ))
-                    for i, p in enumerate(r.profiles):
-                        r.set(i, t, "zfd", zfd[i])
-
             results.set("sediment_timestamps", ts)
             logger.info(f"read_gra: ... end with {len(ts)} timestamp read")
 

src/View/MainWindow.py

@@ -1925,6 +1925,7 @@ class ApplicationWindow(QMainWindow, ListedSubWindow, WindowToolKit):
 
     @timer
     def _diff_results(self, solver1, solver2, solver3, result1, result2):
+        result3 = Results(study=self._study, solver=solver3)
         ts = sorted(
             list(
                 result1.get("timestamps").intersection(
@@ -1933,73 +1934,37 @@ class ApplicationWindow(QMainWindow, ListedSubWindow, WindowToolKit):
             )
         )
 
-        result3 = Results(study=self._study, solver=solver3)
-        result4 = Results(study=result1.study, solver=solver1)
-        result4._river = result1._river
-        result5 = Results(study=result2.study, solver=solver2)
-        result5._river = result2._river
-
         result3.set("nb_reach", result1.get("nb_reach"))
-        result4.set("nb_reach", result1.get("nb_reach"))
-        result5.set("nb_reach", result1.get("nb_reach"))
         result3.set("nb_profile", result1.get("nb_profile"))
-        result4.set("nb_profile", result1.get("nb_profile"))
-        result5.set("nb_profile", result1.get("nb_profile"))
         result3.set("timestamps", ts)
-        result4.set("timestamps", ts)
-        result5.set("timestamps", ts)
 
         for i in range(int(result1.get("nb_reach"))):
             r = result3.river.add(i)
 
-        for timestamp in ts:
+        for timestamp in result3.get("timestamps"):
             for r in range(int(result1.get("nb_reach"))):
+                reach1 = result1.river.reach(r)
+                reach2 = result2.river.reach(r)
                 reach3 = result3.river.reach(r)
-                reach1 = result4.river.reach(r)
-                reach2 = result5.river.reach(r)
 
-                for p, (profile1, profile2, profile3) in enumerate(zip(
+                for profile1, profile2, profile3 in zip(
                         reach1.profiles,
                         reach2.profiles,
-                        reach3.profiles)):
+                        reach3.profiles):
                     for key in ["Z", "Q", "V"]:
                         d1 = profile1.get_ts_key(timestamp, key)
                         d2 = profile2.get_ts_key(timestamp, key)
-                        d = d1-d2
+                        d = d1 - d2
-                        reach3.set(p, timestamp, key, d)
-                        reach1.set(p, timestamp, key, d1)
-                        reach2.set(p, timestamp, key, d2)
-                    if reach1.has_bedload():
-                        d1 = profile1.get_ts_key(timestamp, 'zfd')
-                        reach1.set(p, timestamp, 'zfd', d1)
-                    if reach2.has_bedload():
-                        d2 = profile2.get_ts_key(timestamp, 'zfd')
-                        reach2.set(p, timestamp, 'zfd', d2)
-                        if reach1.has_bedload():
-                            d3 = d1-d2
-                            reach3.set(p, timestamp, 'zfd', d3)
 
-                    limits = reach3.profile(p).geometry.get_water_limits(
+                        profile3.set(timestamp, key, d)
-                            reach3.profile(p).get_ts_key(timestamp, "Z")
+
+                    limits = profile3.geometry\
+                                     .get_water_limits(
+                                         profile3.get_ts_key(timestamp, "Z")
                                      )
-                    reach3.set(
+                    profile3.set(timestamp, "water_limits", limits)
-                            p, timestamp,
-                            "water_limits",
-                            limits
-                        )
-                    limits = profile1.get_ts_key(timestamp, "water_limits")
-                    reach1.set(p, timestamp, "water_limits", limits)
-                    limits = profile2.get_ts_key(timestamp, "water_limits")
-                    reach2.set(p, timestamp, "water_limits", limits)
 
-        for res in (result3, result4, result5):
+        return [result1, result2, result3]
-            for r in range(int(res.get("nb_reach"))):
-                for key in ["Z", "Q", "V"]:
-                    res.river.reach(r).bufferize(ts, key)
-                if res.river.reach(r).has_bedload():
-                    res.river.reach(r).bufferize(ts, "zfd")
-
-        return [result4, result5, result3]
 
     def open_results_adists(self):
         if self._study is None:

src/View/Results/CustomPlot/Plot.py

@@ -88,24 +88,86 @@ class CustomPlot(PamhyrPlot):
         self.lines = {}
 
     def draw_bottom_with_bedload(self, reach):
+        self._bedrock = self.sl_compute_bedrock(reach)
 
         rk = reach.geometry.get_rk()
+        z = self.sl_compute_current_z(reach)
+
+        return z
+
+    def sl_compute_current_z(self, reach):
+        z_br = self._bedrock
+        sl = self.sl_compute_current_rk(reach)
+
         z = list(
             map(
-                lambda p: p.get_ts_key(
+                lambda z, sl: reduce(
-                    self._current_timestamp, "zfd"
+                    lambda z, h: z + h[0],
+                    sl, z
                 ),
-                reach.profiles
+                z_br,           # Bedrock elevation
+                sl              # Current sediment layers
             )
         )
 
         return z
 
+    def sl_compute_bedrock(self, reach):
+        z_min = reach.geometry.get_z_min()
+        sl = self.sl_compute_initial(reach)
+
+        z = list(
+            map(
+                lambda z, sl: reduce(
+                    lambda z, h: z - h[0],
+                    sl, z
+                ),
+                z_min,          # Original geometry
+                sl              # Original sediment layers
+            )
+        )
+
+        return z
+
+    def sl_compute_initial(self, reach):
+        """
+        Get SL list for profile p at initial time (initial data)
+        """
+        t0 = min(list(self.data[self._current_res_id[0]].get("timestamps")))
+        return map(
+            lambda p: p.get_ts_key(t0, "sl")[0],
+            reach.profiles
+        )
+
+    def sl_compute_current_rk(self, reach):
+        """
+        Get SL list for profile p at current time
+        """
+        return map(
+            lambda p: p.get_ts_key(self._current_timestamp, "sl")[0],
+            reach.profiles
+        )
+
     def get_ts_zmin(self, profile, res_id):
         results = self.data[res_id]
+        nt = len(list(results.get("timestamps")))
         reach = results.river.reach(self._current_reach)
-        zfd = reach.profile(profile).get_key("zfd")
+        berdrock = self.sl_compute_bedrock(reach)
-        return zfd
+        sl = reach.profile(profile).get_key("sl")
+
+        ts_z_bedrock = [berdrock[profile]]*nt
+
+        ts_z_min = list(
+            map(
+                lambda z, sl: reduce(
+                    lambda z, h: z + h,
+                    sl, z
+                ),
+                ts_z_bedrock,           # Bedrock elevations
+                asarray(sl)[:, 0, :, 0]    # Sediment layers
+            )
+        )
+        return ts_z_min
 
     def _draw_rk(self):
         results = self.data[self._current_res_id]
@@ -114,7 +176,7 @@ class CustomPlot(PamhyrPlot):
             reach1 = self.data[0].river.reach(self._current_reach)
             reach2 = self.data[1].river.reach(self._current_reach)
         rk = reach.geometry.get_rk()
-        if reach.has_bedload():
+        if reach.has_sediment():
             z_min = self.draw_bottom_with_bedload(reach)
         else:
             z_min = reach.geometry.get_z_min()
@@ -154,32 +216,14 @@ class CustomPlot(PamhyrPlot):
 
             ax = self._axes[unit["bed_elevation"]]
             if self._current_res_id < 2:
-                if reach.has_bedload():
-                    dz = self.draw_bottom_with_bedload(reach)
-                else:
-                    dz = z_min
                 line = ax.plot(
-                    rk, dz,
+                    rk, z_min,
                     color='grey', lw=1.,
                 )
             else:
-                if reach.has_bedload():
+                if reach.has_sediment():
-                    z_min1 = list(
+                    z_min1 = self.draw_bottom_with_bedload(reach1)
-                            map(
+                    z_min2 = self.draw_bottom_with_bedload(reach2)
-                                lambda p: p.get_ts_key(
-                                    self._current_timestamp, "zfd"
-                                ),
-                                reach1.profiles
-                            )
-                        )
-                    z_min2 = list(
-                            map(
-                                lambda p: p.get_ts_key(
-                                    self._current_timestamp, "zfd"
-                                ),
-                                reach2.profiles
-                            )
-                        )
                 else:
                     z_min1 = reach1.geometry.get_z_min()
                     z_min2 = reach2.geometry.get_z_min()
@@ -197,7 +241,7 @@ class CustomPlot(PamhyrPlot):
             self.lines["bed_elevation"] = line
 
             if (self._envelop and
-                    reach.has_bedload() and
+                    reach.has_sediment() and
                     self._current_res_id < 2):
 
                 ax = self._axes[unit["bed_elevation_envelop"]]
@@ -566,7 +610,7 @@ class CustomPlot(PamhyrPlot):
 
         rk = reach.geometry.get_rk()
         z_min = reach.geometry.get_z_min()
-        if reach.has_bedload():
+        if reach.has_sediment():
             z_min = self.draw_bottom_with_bedload(reach)
         else:
             z_min = reach.geometry.get_z_min()
@@ -591,12 +635,9 @@ class CustomPlot(PamhyrPlot):
         )
         if "bed_elevation" in self._y:
             if self._current_res_id < 2:
-                if reach.has_bedload():
-                    dz = self.draw_bottom_with_bedload(reach)
-                else:
                 dz = z_min
             else:
-                if reach.has_bedload():
+                if reach.has_sediment():
                     z_min1 = self.draw_bottom_with_bedload(reach1)
                     z_min2 = self.draw_bottom_with_bedload(reach2)
                 else:
@@ -855,7 +896,7 @@ class CustomPlot(PamhyrPlot):
         v = profile.get_key("V")
         z_min = profile.geometry.z_min()
         if self._current_res_id < 2:
-            if reach.has_bedload():
+            if reach.has_sediment():
                 ts_z_min = self.get_ts_zmin(
                     self._current_profile_id, self._current_res_id)
             else:
@@ -873,9 +914,9 @@ class CustomPlot(PamhyrPlot):
             ax = self._axes[unit["bed_elevation"]]
 
             if self._current_res_id == 2:
-                if reach.has_bedload():
+                if reach.has_sediment():
-                    ts_z_min1 = self.get_ts_zmin(self._current_profile_id, 0)
+                    ts_z_min1 = self.get_ts_zmin(self._current_profile_id1, 0)
-                    ts_z_min2 = self.get_ts_zmin(self._current_profile_id, 1)
+                    ts_z_min2 = self.get_ts_zmin(self._current_profile_id2, 1)
                     ts_z_min = list(
                         map(
                             lambda x, y: x - y,
@@ -1085,7 +1126,7 @@ class CustomPlot(PamhyrPlot):
         z = profile.get_key("Z")
         v = profile.get_key("V")
         if self._current_res_id < 2:
-            if reach.has_bedload():
+            if reach.has_sediment():
                 ts_z_min = self.get_ts_zmin(
                     self._current_profile_id, self._current_res_id)
             else:
@@ -1098,9 +1139,9 @@ class CustomPlot(PamhyrPlot):
                 )
         if "bed_elevation" in self._y:
             if self._current_res_id == 2:
-                if reach.has_bedload():
+                if reach.has_sediment():
-                    ts_z_min1 = self.get_ts_zmin(self._current_profile_id, 0)
+                    ts_z_min1 = self.get_ts_zmin(self._current_profile_id1, 0)
-                    ts_z_min2 = self.get_ts_zmin(self._current_profile_id, 1)
+                    ts_z_min2 = self.get_ts_zmin(self._current_profile_id2, 1)
                     ts_z_min = list(
                         map(
                             lambda x, y: x - y,

src/View/Results/PlotRKC.py

@@ -96,29 +96,24 @@ class PlotRKC(PamhyrPlot):
         self._init = True
 
     def draw_bottom(self, reach):
-        if reach.has_bedload():
+        if reach.has_sediment():
             self.draw_bottom_with_bedload(reach)
         else:
             self.draw_bottom_geometry(reach)
 
     def draw_bottom_with_bedload(self, reach):
+        self._bedrock = self.sl_compute_bedrock(reach)
+
         rk = reach.geometry.get_rk()
-        zfd = list(
+        z = self.sl_compute_current_z(reach)
-            map(
-                lambda p: p.get_ts_key(
-                    self._current_timestamp, "zfd"
-                ),
-                reach.profiles
-            )
-        )
 
         self.line_bottom, = self.canvas.axes.plot(
-            rk, zfd,
+            rk, z,
             linestyle="solid", lw=1.,
             color=self.color_plot_river_bottom,
         )
 
-        self._river_bottom = zfd
+        self._river_bottom = z
 
     def draw_profiles_hs(self, reach):
         results = self.results[self._current_res_id]
@@ -152,6 +147,58 @@ class PlotRKC(PamhyrPlot):
                 fontsize=9, color=self.color_plot_previous,
             )
 
+    def sl_compute_bedrock(self, reach):
+        z_min = reach.geometry.get_z_min()
+        sl = self.sl_compute_initial(reach)
+
+        z = list(
+            map(
+                lambda z, sl: reduce(
+                    lambda z, h: z - h[0],
+                    sl, z
+                ),
+                z_min,          # Original geometry
+                sl              # Original sediment layers
+            )
+        )
+
+        return z
+
+    def sl_compute_current_z(self, reach):
+        z_br = self._bedrock
+        sl = self.sl_compute_current(reach)
+
+        z = list(
+            map(
+                lambda z, sl: reduce(
+                    lambda z, h: z + h[0],
+                    sl, z
+                ),
+                z_br,           # Bedrock elevation
+                sl              # Current sediment layers
+            )
+        )
+
+        return z
+
+    def sl_compute_initial(self, reach):
+        """
+        Get SL list for profile p at initial time (initial data)
+        """
+        return map(
+            lambda p: p.get_ts_key(min(self._timestamps), "sl")[0],
+            reach.profiles
+        )
+
+    def sl_compute_current(self, reach):
+        """
+        Get SL list for profile p at current time
+        """
+        return map(
+            lambda p: p.get_ts_key(self._current_timestamp, "sl")[0],
+            reach.profiles
+        )
+
     def draw_bottom_geometry(self, reach):
         rk = reach.geometry.get_rk()
         z_min = reach.geometry.get_z_min()
@@ -290,7 +337,7 @@ class PlotRKC(PamhyrPlot):
 
         results = self.results[self._current_res_id]
         reach = results.river.reach(self._current_reach_id)
-        if reach.has_bedload():
+        if reach.has_sediment():
             self.update_bottom_with_bedload()
 
         self.update_water_elevation()
@@ -350,23 +397,14 @@ class PlotRKC(PamhyrPlot):
         results = self.results[self._current_res_id]
         reach = results.river.reach(self._current_reach_id)
         rk = reach.geometry.get_rk()
-        # z = self.sl_compute_current_z(reach)
+        z = self.sl_compute_current_z(reach)
-
-        zfd = list(
-            map(
-                lambda p: p.get_ts_key(
-                    self._current_timestamp, "zfd"
-                ),
-                reach.profiles
-            )
-        )
 
         self.line_bottom.remove()
 
         self.line_bottom, = self.canvas.axes.plot(
-            rk, zfd,
+            rk, z,
             linestyle="solid", lw=1.,
             color=self.color_plot_river_bottom,
         )
 
-        self._river_bottom = zfd
+        self._river_bottom = z

src/View/Results/Window.py

@@ -558,7 +558,7 @@ class ResultsWindow(PamhyrWindow):
         table = self.find(QTableView, f"tableView_solver")
         indexes = table.selectedIndexes()
         if len(indexes) == 0:
-            return []
+            return
 
         return [i.row() for i in indexes]
 
@@ -846,17 +846,17 @@ class ResultsWindow(PamhyrWindow):
         my_dict = {}
         my_dict[dict_x["rk"]] = reach.geometry.get_rk()
         if "bed_elevation" in y:
-            if reach.has_bedload():
+            if self._current_results != 2:
+                zmin = reach.geometry.get_z_min()
+            else:
+                z_min1 = reach1.geometry.get_z_min()
+                z_min2 = reach2.geometry.get_z_min()
                 zmin = list(
                     map(
-                            lambda p: p.get_ts_key(
+                        lambda x, y: x - y,
-                                timestamp, "zfd"
+                        z_min2, z_min2
-                            ),
-                            reach.profiles
                         )
                     )
-            else:
-                zmin = reach.geometry.get_z_min()
             my_dict[dict_y["bed_elevation"]] = zmin
             # if envelop and reach.has_sediment():
         if "discharge" in y:
@@ -1091,16 +1091,29 @@ class ResultsWindow(PamhyrWindow):
         if self._current_results == 2:
             reach1 = self._results[0].river.reach(self._reach)
             reach2 = self._results[1].river.reach(self._reach)
-            reach3 = self._results[2].river.reach(self._reach)
             profile1 = reach1.profile(self._profile)
             profile2 = reach2.profile(self._profile)
-            profile3 = reach3.profile(self._profile)
+
+            q1 = profile1.get_key("Q")
+            z1 = profile1.get_key("Z")
+            v1 = profile1.get_key("V")
+
+            q2 = profile2.get_key("Q")
+            z2 = profile2.get_key("Z")
+            v2 = profile2.get_key("V")
 
         if "bed_elevation" in y:
-            if reach.has_bedload():
+            if self._current_results != 2:
-                z_min = profile.get_key("zfd")
-            else:
                 z_min = [profile.geometry.z_min()] * len(self._timestamps)
+            else:
+                z_min1 = profile1.geometry.z_min()
+                z_min2 = profile2.geometry.z_min()
+                z_min = list(
+                    map(
+                        lambda ts:  z_min1 - z_min2,
+                        self._timestamps
+                    )
+                )
             my_dict[dict_y["bed_elevation"]] = z_min
         if "discharge" in y:
             my_dict[dict_y["discharge"]] = q