debug

src/Model/Results/Results.py

@@ -219,6 +219,11 @@ 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,6 +17,7 @@
 import struct
 import logging
 import itertools
+import numpy as np
 
 from tools import flatten
 from functools import reduce
@@ -76,9 +77,12 @@ class Profile(SQLSubModel):
                 return self._data[timestamp][key]
         return None
 
-    def has_sediment(self):
+    def has_sediment_layers(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"""
@@ -275,6 +279,7 @@ class Reach(SQLSubModel):
                 lambda p: p.name[0:8] != 'interpol', self._profiles
             )
         )
+        self._buffers = {}
 
     def __len__(self):
         return len(self._profiles)
@@ -295,6 +300,10 @@ 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]
 
@@ -302,7 +311,16 @@ class Reach(SQLSubModel):
         self._profiles[profile_id].set(timestamp, key, data)
 
     def has_sediment(self):
-        return any(map(lambda profile: profile.has_sediment(), self._profiles))
+        return any(map(lambda profile: profile.has_sediment_layers(),
+                       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,6 +19,7 @@
 import os
 import logging
 import numpy as np
+from functools import reduce
 
 from itertools import chain
 from tools import timer, trace, logger_exception
@@ -1196,6 +1197,10 @@ 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():
@@ -1362,6 +1367,40 @@ 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,7 +1925,6 @@ 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(
@@ -1934,37 +1933,73 @@ 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 result3.get("timestamps"):
+        for timestamp in ts:
             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 profile1, profile2, profile3 in zip(
+                for p, (profile1, profile2, profile3) in enumerate(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)
 
-                        profile3.set(timestamp, key, d)
+                    limits = reach3.profile(p).geometry.get_water_limits(
+                            reach3.profile(p).get_ts_key(timestamp, "Z")
+                    )
+                    reach3.set(
+                            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)
 
-                    limits = profile3.geometry\
-                                     .get_water_limits(
-                                         profile3.get_ts_key(timestamp, "Z")
-                                     )
-                    profile3.set(timestamp, "water_limits", limits)
+        for res in (result3, result4, result5):
+            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 [result1, result2, result3]
+        return [result4, result5, result3]
 
     def open_results_adists(self):
         if self._study is None:

src/View/Results/CustomPlot/Plot.py

@@ -88,86 +88,24 @@ 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 z, sl: reduce(
-                    lambda z, h: z + h[0],
-                    sl, z
+                lambda p: p.get_ts_key(
+                    self._current_timestamp, "zfd"
                 ),
-                z_br,           # Bedrock elevation
-                sl              # Current sediment layers
+                reach.profiles
             )
         )
 
         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)
-        berdrock = self.sl_compute_bedrock(reach)
-        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
+        zfd = reach.profile(profile).get_key("zfd")
+        return zfd
 
     def _draw_rk(self):
         results = self.data[self._current_res_id]
@@ -176,7 +114,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_sediment():
+        if reach.has_bedload():
             z_min = self.draw_bottom_with_bedload(reach)
         else:
             z_min = reach.geometry.get_z_min()
@@ -216,14 +154,32 @@ 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, z_min,
+                    rk, dz,
                     color='grey', lw=1.,
                 )
             else:
-                if reach.has_sediment():
-                    z_min1 = self.draw_bottom_with_bedload(reach1)
-                    z_min2 = self.draw_bottom_with_bedload(reach2)
+                if reach.has_bedload():
+                    z_min1 = list(
+                            map(
+                                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()
@@ -241,7 +197,7 @@ class CustomPlot(PamhyrPlot):
             self.lines["bed_elevation"] = line
 
             if (self._envelop and
-                    reach.has_sediment() and
+                    reach.has_bedload() and
                     self._current_res_id < 2):
 
                 ax = self._axes[unit["bed_elevation_envelop"]]
@@ -610,7 +566,7 @@ class CustomPlot(PamhyrPlot):
 
         rk = reach.geometry.get_rk()
         z_min = reach.geometry.get_z_min()
-        if reach.has_sediment():
+        if reach.has_bedload():
             z_min = self.draw_bottom_with_bedload(reach)
         else:
             z_min = reach.geometry.get_z_min()
@@ -635,9 +591,12 @@ class CustomPlot(PamhyrPlot):
         )
         if "bed_elevation" in self._y:
             if self._current_res_id < 2:
-                dz = z_min
+                if reach.has_bedload():
+                    dz = self.draw_bottom_with_bedload(reach)
+                else:
+                    dz = z_min
             else:
-                if reach.has_sediment():
+                if reach.has_bedload():
                     z_min1 = self.draw_bottom_with_bedload(reach1)
                     z_min2 = self.draw_bottom_with_bedload(reach2)
                 else:
@@ -896,7 +855,7 @@ class CustomPlot(PamhyrPlot):
         v = profile.get_key("V")
         z_min = profile.geometry.z_min()
         if self._current_res_id < 2:
-            if reach.has_sediment():
+            if reach.has_bedload():
                 ts_z_min = self.get_ts_zmin(
                     self._current_profile_id, self._current_res_id)
             else:
@@ -914,9 +873,9 @@ class CustomPlot(PamhyrPlot):
             ax = self._axes[unit["bed_elevation"]]
 
             if self._current_res_id == 2:
-                if reach.has_sediment():
-                    ts_z_min1 = self.get_ts_zmin(self._current_profile_id1, 0)
-                    ts_z_min2 = self.get_ts_zmin(self._current_profile_id2, 1)
+                if reach.has_bedload():
+                    ts_z_min1 = self.get_ts_zmin(self._current_profile_id, 0)
+                    ts_z_min2 = self.get_ts_zmin(self._current_profile_id, 1)
                     ts_z_min = list(
                         map(
                             lambda x, y: x - y,
@@ -1126,7 +1085,7 @@ class CustomPlot(PamhyrPlot):
         z = profile.get_key("Z")
         v = profile.get_key("V")
         if self._current_res_id < 2:
-            if reach.has_sediment():
+            if reach.has_bedload():
                 ts_z_min = self.get_ts_zmin(
                     self._current_profile_id, self._current_res_id)
             else:
@@ -1139,9 +1098,9 @@ class CustomPlot(PamhyrPlot):
                 )
         if "bed_elevation" in self._y:
             if self._current_res_id == 2:
-                if reach.has_sediment():
-                    ts_z_min1 = self.get_ts_zmin(self._current_profile_id1, 0)
-                    ts_z_min2 = self.get_ts_zmin(self._current_profile_id2, 1)
+                if reach.has_bedload():
+                    ts_z_min1 = self.get_ts_zmin(self._current_profile_id, 0)
+                    ts_z_min2 = self.get_ts_zmin(self._current_profile_id, 1)
                     ts_z_min = list(
                         map(
                             lambda x, y: x - y,

src/View/Results/PlotRKC.py

@@ -96,24 +96,29 @@ class PlotRKC(PamhyrPlot):
         self._init = True
 
     def draw_bottom(self, reach):
-        if reach.has_sediment():
+        if reach.has_bedload():
             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()
-        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, = self.canvas.axes.plot(
-            rk, z,
+            rk, zfd,
             linestyle="solid", lw=1.,
             color=self.color_plot_river_bottom,
         )
 
-        self._river_bottom = z
+        self._river_bottom = zfd
 
     def draw_profiles_hs(self, reach):
         results = self.results[self._current_res_id]
@@ -147,58 +152,6 @@ 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()
@@ -337,7 +290,7 @@ class PlotRKC(PamhyrPlot):
 
         results = self.results[self._current_res_id]
         reach = results.river.reach(self._current_reach_id)
-        if reach.has_sediment():
+        if reach.has_bedload():
             self.update_bottom_with_bedload()
 
         self.update_water_elevation()
@@ -397,14 +350,23 @@ 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, z,
+            rk, zfd,
             linestyle="solid", lw=1.,
             color=self.color_plot_river_bottom,
         )
 
-        self._river_bottom = z
+        self._river_bottom = zfd

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 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()
+            if reach.has_bedload():
                 zmin = list(
-                    map(
-                        lambda x, y: x - y,
-                        z_min2, z_min2
+                        map(
+                            lambda p: p.get_ts_key(
+                                timestamp, "zfd"
+                            ),
+                            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,29 +1091,16 @@ 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)
-
-            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")
+            profile3 = reach3.profile(self._profile)
 
         if "bed_elevation" in y:
-            if self._current_results != 2:
-                z_min = [profile.geometry.z_min()] * len(self._timestamps)
+            if reach.has_bedload():
+                z_min = profile.get_key("zfd")
             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
-                    )
-                )
+                z_min = [profile.geometry.z_min()] * len(self._timestamps)
             my_dict[dict_y["bed_elevation"]] = z_min
         if "discharge" in y:
             my_dict[dict_y["discharge"]] = q