optimization + debug

src/Model/Geometry/ProfileXYZ.py

@@ -20,6 +20,7 @@ import logging
 import numpy as np
 from typing import List
 from functools import reduce
+from dataclasses import dataclass
 
 from tools import timer
 from shapely import geometry
@@ -33,6 +34,13 @@ from Model.Geometry.Vector_1d import Vector1d
 logger = logging.getLogger()
 
 
+@dataclass
+class Tabulation:
+    z: np.array([])
+    A: np.array([])
+    L: np.array([])
+
+
 class ProfileXYZ(Profile, SQLSubModel):
     _sub_classes = [
         PointXYZ,
@@ -70,6 +78,14 @@ class ProfileXYZ(Profile, SQLSubModel):
             status=status,
         )
 
+        self.tab = Tabulation([], [], [])
+        self.tab_up_to_date = False
+        self.time_z = 0.0
+        self.time_A = 0.0
+        self.time_l = 0.0
+        self._station = []
+        self.station_up_to_date = False
+
     @classmethod
     def _db_create(cls, execute):
         execute("""
@@ -305,6 +321,8 @@ class ProfileXYZ(Profile, SQLSubModel):
             pt = PointXYZ(*point, profile=self, status=self._status)
             self.points.append(pt)
         self._status.modified()
+        self.tab_up_to_date = False
+        self.station_up_to_date = False
 
     def get_point_i(self, index: int) -> PointXYZ:
         """Get point at index.
@@ -366,6 +384,8 @@ class ProfileXYZ(Profile, SQLSubModel):
         point_xyz = PointXYZ(0., 0., 0., profile=self, status=self._status)
         self.points.append(point_xyz)
         self._status.modified()
+        self.tab_up_to_date = False
+        self.station_up_to_date = False
 
     def insert(self, index: int):
         """Insert a new point at index.
@@ -379,6 +399,8 @@ class ProfileXYZ(Profile, SQLSubModel):
         point = PointXYZ(0., 0., 0., profile=self, status=self._status)
         self.points.insert(index, point)
         self._status.modified()
+        self.tab_up_to_date = False
+        self.station_up_to_date = False
         return point
 
     def filter_isnan(self, lst):
@@ -410,7 +432,7 @@ class ProfileXYZ(Profile, SQLSubModel):
 
         return abs(rg.dist(rd))
 
-    def wet_width(self, z):
+    def compute_wet_width(self, z):
         start, end = self.get_all_water_limits_ac(z)
 
         if len(start) == 0:
@@ -421,6 +443,17 @@ class ProfileXYZ(Profile, SQLSubModel):
             length += abs(s - e)
         return length
 
+    def wet_width(self, z):
+        if self.tab_up_to_date:
+            if z > self.tab.z[-1]:
+                return self.tab.L[-1]
+            elif z < self.tab.z[0]:
+                return 0.0
+            else:
+                return np.interp(z, self.tab.z, self.tab.L)
+        else:
+            return self.compute_wet_width(z)
+
     def wet_perimeter(self, z):
         lines = self.wet_lines(z)
 
@@ -432,17 +465,46 @@ class ProfileXYZ(Profile, SQLSubModel):
             length += line.length
         return length
 
-    def wet_area(self, z):
+    def compute_wet_area(self, z):
-        lines = self.wet_lines(z)
-
-        if lines is None:
-            return 0
-
         area = 0.0
-        for line in lines:
+        if len(self.tab.L) > 0:
-            if len(line.coords) > 2:
+            if z < self.tab.z[0]:
-                poly = geometry.Polygon(line)
+                return 0.0
-                area += poly.area
+            i = np.searchsorted([z], self.tab.z, side='right')[0]
+            for j in range(i-1):
+                area += (self.tab.L[j] + self.tab.L[j+1]) * (
+                         self.tab.z[j+1] - self.tab.z[j]) / 2.0
+            area += (self.tab.L[i-1] + self.wet_width(z))
+        else:
+            lines = self.wet_lines(z)
+
+            if lines is None:
+                return 0.0
+            for line in lines:
+                if len(line.coords) > 2:
+                    poly = geometry.Polygon(line)
+                    area += poly.area
+        return area
+
+    def wet_area(self, z):
+        if self.tab_up_to_date:
+            if z > self.tab.z[-1]:
+                return self.tab.A[-1] + self.tab.L[-1] * (z - self.tab.z[-1])
+            elif z < self.tab.z[0]:
+                return 0.0
+            else:
+                return np.interp(z, self.tab.z, self.tab.A)
+        else:
+            lines = self.wet_lines(z)
+
+            if lines is None:
+                return 0.0
+
+            area = 0.0
+            for line in lines:
+                if len(line.coords) > 2:
+                    poly = geometry.Polygon(line)
+                    area += poly.area
         return area
 
     def wet_radius(self, z):
@@ -563,7 +625,7 @@ class ProfileXYZ(Profile, SQLSubModel):
             return None
 
         zz = list(map(lambda p: p.z, points))
-        station = self._get_station(points)
+        station = self.get_station()
 
         start = []
         if points[0].z <= z:
@@ -651,6 +713,19 @@ class ProfileXYZ(Profile, SQLSubModel):
 
         return pt_left, pt_right
 
+    def compute_tabulation(self):
+        sorted_points = sorted(self._points, key=lambda p: p.z)
+        self.tab.z = np.array([p.z for p in sorted_points], np.float64)
+        self.tab.L = np.zeros(len(self.tab.z), np.float64)
+        self.tab.A = np.zeros(len(self.tab.z), np.float64)
+        for i in range(1, len(self.tab.z)):
+            self.tab.L[i] = self.compute_wet_width(self.tab.z[i])
+            dx = (self.tab.L[i] + self.tab.L[i-1])/2
+            dz = self.tab.z[i] - self.tab.z[i-1]
+            self.tab.A[i] = self.tab.A[i-1] + dz * dx
+
+        self.tab_up_to_date = True
+
     def get_station(self):
         """Projection of the points of the profile on a plane.
 
@@ -663,7 +738,12 @@ class ProfileXYZ(Profile, SQLSubModel):
         if self.nb_points < 2:
             return [0.0]
         else:
-            return self._get_station(self.points)
+            if self.station_up_to_date:
+                return self._station
+            else:
+                self._station = self._get_station(self.points)
+                self.station_up_to_date = True
+                return self._station
 
     @timer
     def _get_station(self, points):

src/Solver/Mage.py

@@ -1001,7 +1001,6 @@ class Mage8(Mage):
 
                 # Add profile id correspondance to reach
                 key = (i1, i2)
-                print("mage keys reachs: ", key)
                 iprofiles[key] = r
 
                 # Profile ID offset
@@ -1034,8 +1033,6 @@ class Mage8(Mage):
             ]
             def ip_to_ri(r, i): return i - reach_offset[r]
 
-            print("mage iprofiles: ", iprofiles)
-
             ts = set()
             end = False
             while not end:
@@ -1071,7 +1068,23 @@ class Mage8(Mage):
 
             logger.debug(reachs[0].profiles[0]._data)
             results.set("timestamps", ts)
-            print("set timestamps mage: ", ts)
+            ts_list = sorted(ts)
+            logger.info(f"compute tab...")
+            for r in reachs:
+                for p in r.profiles:
+                    if not p.geometry.tab_up_to_date:
+                        p.geometry.compute_tabulation()
+
+            logger.info(f"compute velocily...")
+
+            for r in reachs:
+                for t in ts_list:
+                    for i, p in enumerate(r.profiles):
+                        v = p.geometry.speed(
+                                p.get_ts_key(t, "Q"),
+                                p.get_ts_key(t, "Z")
+                        )
+                        r.set(i, t, "V", v)
             logger.info(f"read_bin: ... end with {len(ts)} timestamp read")
 
     @timer

src/View/About/Window.py

@@ -55,7 +55,7 @@ class AboutWindow(PamhyrDialog):
 
             label = self.get_label_text("label_version")
             label = label.replace("@version", version)
-            label = label.replace("@codename", "(Tahiti)")
+            label = label.replace("@codename", "(Adis-TS)")
             self.set_label_text("label_version", label)
 
         # Authors

src/View/DIFAdisTS/Table.py

@@ -187,7 +187,6 @@ class DIFTableModel(PamhyrTableModel):
                     )
                 )
             elif self._headers[column] == "reach":
-                print(self._river.edge(value).id)
                 self._undo.push(
                     SetCommandSpec(
                         self._lst,

src/View/DIFAdisTS/Window.py

@@ -287,9 +287,7 @@ class DIFAdisTSWindow(PamhyrWindow):
             self._table_spec.add(rows[0])
 
     def delete(self):
-        print("del")
         rows = self.index_selected_rows()
         if len(rows) == 0:
-            print("len 0")
             return
         self._table_spec.delete(rows)

src/View/MainWindow.py

@@ -186,10 +186,10 @@ class ApplicationWindow(QMainWindow, ListedSubWindow, WindowToolKit):
         title = "(dbg) " if self.conf.debug else ""
 
         if self._study is not None:
-            title += f"Pamhyr2 (Tahiti 🌴) - {self._study.name}"
+            title += f"Pamhyr2 - {self._study.name}"
             self.setWindowTitle(title)
         else:
-            title += "Pamhyr2 (Tahiti 🌴)"
+            title += "Pamhyr2"
             self.setWindowTitle(title)
 
     def setup_tab(self):

src/View/Results/CustomPlot/Plot.py

@@ -187,6 +187,12 @@ class CustomPlot(PamhyrPlot):
                 reach.profiles
             )
         )
+        v = list(
+            map(
+                lambda p: p.get_ts_key(self._timestamp, "V"),
+                reach.profiles
+            )
+        )
 
         shift = 0
         compt = 0
@@ -327,14 +333,6 @@ class CustomPlot(PamhyrPlot):
         if "velocity" in self._y:
 
             ax = self._axes[unit["velocity"]]
-            v = list(
-                map(
-                    lambda p: p.geometry.speed(
-                        p.get_ts_key(self._timestamp, "Q"),
-                        p.get_ts_key(self._timestamp, "Z")),
-                    reach.profiles
-                )
-            )
 
             line = ax.plot(
                 rk, v, lw=1.,
@@ -344,12 +342,10 @@ class CustomPlot(PamhyrPlot):
 
             if self._envelop:
 
-                velocities = list(map(
+                velocities = list(
-                    lambda p: list(map(
+                    map(
-                        lambda q, z:
+                        lambda p: p.get_key("V"),
-                            p.geometry.speed(q, z),
+                        reach.profiles
-                            p.get_key("Q"), p.get_key("Z")
-                        )), reach.profiles
                     )
                 )
 
@@ -440,9 +436,7 @@ class CustomPlot(PamhyrPlot):
             fr = list(
                 map(
                     lambda p:
-                        p.geometry.speed(
+                        p.get_ts_key(self._timestamp, "V") /
-                            p.get_ts_key(self._timestamp, "Q"),
-                            p.get_ts_key(self._timestamp, "Z")) /
                         sqrt(9.81 * (
                             p.geometry.wet_area(
                                 p.get_ts_key(self._timestamp, "Z")) /
@@ -506,6 +500,12 @@ class CustomPlot(PamhyrPlot):
                 reach.profiles
             )
         )
+        v = list(
+            map(
+                lambda p: p.get_ts_key(self._timestamp, "V"),
+                reach.profiles
+            )
+        )
         if "bed_elevation" in self._y:
             self.lines["bed_elevation"][0].set_ydata(z_min)
 
@@ -524,15 +524,7 @@ class CustomPlot(PamhyrPlot):
             self.lines["discharge"][0].set_ydata(q)
 
         if "velocity" in self._y:
-            v = list(
+            self.lines["velocity"][0].set_ydata(v)
-                map(
-                    lambda p: p.geometry.speed(
-                        p.get_ts_key(self._timestamp, "Q"),
-                        p.get_ts_key(self._timestamp, "Z")),
-                    reach.profiles
-                )
-            )
-            self.lines["discharge"][0].set_ydata(v)
 
         if "depth" in self._y:
             d = list(
@@ -558,9 +550,7 @@ class CustomPlot(PamhyrPlot):
             fr = list(
                 map(
                     lambda p:
-                        p.geometry.speed(
+                        p.get_ts_key(self._timestamp, "V") /
-                            p.get_ts_key(self._timestamp, "Q"),
-                            p.get_ts_key(self._timestamp, "Z")) /
                         sqrt(9.81 * (
                             p.geometry.wet_area(
                                 p.get_ts_key(self._timestamp, "Z")) /
@@ -642,6 +632,7 @@ class CustomPlot(PamhyrPlot):
 
         q = profile.get_key("Q")
         z = profile.get_key("Z")
+        v = profile.get_key("V")
         z_min = profile.geometry.z_min()
         if reach.has_sediment():
             ts_z_min = self.get_ts_zmin(self._profile)
@@ -693,12 +684,6 @@ class CustomPlot(PamhyrPlot):
         if "velocity" in self._y:
 
             ax = self._axes[unit["velocity"]]
-            v = list(
-                map(
-                    lambda q, z: profile.geometry.speed(q, z),
-                    q, z
-                )
-            )
 
             line = ax.plot(
                 ts, v, lw=1.,
@@ -736,12 +721,12 @@ class CustomPlot(PamhyrPlot):
 
             ax = self._axes[unit["froude"]]
             d = list(
-                map(lambda z, q:
+                map(lambda z, v:
-                    profile.geometry.speed(q, z) /
+                    v /
                     sqrt(9.81 * (
                         profile.geometry.wet_area(z) /
                         profile.geometry.wet_width(z))
-                    ), z, q)
+                    ), z, v)
                 )
 
             line = ax.plot(
@@ -782,6 +767,7 @@ class CustomPlot(PamhyrPlot):
 
         q = profile.get_key("Q")
         z = profile.get_key("Z")
+        v = profile.get_key("V")
         if reach.has_sediment():
             ts_z_min = self.get_ts_zmin(self._profile)
         else:
@@ -810,12 +796,6 @@ class CustomPlot(PamhyrPlot):
             self.lines["discharge"][0].set_ydata(q)
 
         if "velocity" in self._y:
-            v = list(
-                map(
-                    lambda q, z: profile.geometry.speed(q, z),
-                    q, z
-                )
-            )
             self.lines["velocity"][0].set_ydata(v)
 
         if "depth" in self._y:
@@ -832,12 +812,12 @@ class CustomPlot(PamhyrPlot):
 
         if "froude" in self._y:
             d = list(
-                map(lambda z, q:
+                map(lambda z, v:
-                    profile.geometry.speed(q, z) /
+                    v /
                     sqrt(9.81 * (
                         profile.geometry.wet_area(z) /
                         profile.geometry.wet_width(z))
-                    ), z, q)
+                    ), z, v)
                 )
             self.lines["froude"][0].set_ydata(d)
 

src/View/Results/Table.py

@@ -47,12 +47,8 @@ _translate = QCoreApplication.translate
 class TableModel(PamhyrTableModel):
     def _setup_lst(self):
         _river = self._data.river
-        print("data: ", self._data)
-        print("river: ", _river)
-        print("reaches: ", _river.reachs)
         if self._opt_data == "reach":
             self._lst = _river.reachs
-            print("optreach: ", self._lst)
         elif self._opt_data == "profile":
             self._lst = _river.reach(0).profiles
         elif self._opt_data == "raw_data":
@@ -95,7 +91,7 @@ class TableModel(PamhyrTableModel):
             elif self._headers[column] == "velocity":
                 q = self._lst[row].get_ts_key(self._timestamp, "Q")
                 z = self._lst[row].get_ts_key(self._timestamp, "Z")
-                v = self._lst[row].geometry.speed(q, z)
+                v = self._lst[row].get_ts_key(self._timestamp, "V")
                 return f"{v:.4f}"
             elif self._headers[column] == "width":
                 z = self._lst[row].get_ts_key(self._timestamp, "Z")
@@ -124,7 +120,7 @@ class TableModel(PamhyrTableModel):
             elif self._headers[column] == "froude":
                 q = self._lst[row].get_ts_key(self._timestamp, "Q")
                 z = self._lst[row].get_ts_key(self._timestamp, "Z")
-                v = self._lst[row].geometry.speed(q, z)
+                v = self._lst[row].get_ts_key(self._timestamp, "V")
                 a = self._lst[row].geometry.wet_area(z)
                 b = self._lst[row].geometry.wet_width(z)
                 froude = v / sqrt(9.81 * (a / b))

src/View/Results/Window.py

@@ -730,21 +730,16 @@ class ResultsWindow(PamhyrWindow):
         if "velocity" in y:
             my_dict[dict_y["velocity"]] = list(
                 map(
-                    lambda p: p.geometry.speed(
+                    lambda p: p.get_ts_key(timestamp, "V"),
-                        p.get_ts_key(timestamp, "Q"),
-                        p.get_ts_key(timestamp, "Z")),
                     reach.profiles
                 )
             )
             if envelop:
                 velocities = list(map(
-                    lambda p: list(map(
+                            lambda p: p.get_key("V"),
-                        lambda q, z:
+                            reach.profiles
-                            p.geometry.speed(q, z),
+                        )
-                            p.get_key("Q"), p.get_key("Z")
-                        )), reach.profiles
                     )
-                )
                 my_dict[dict_y["min_velocity"]] = [min(v) for v in velocities]
                 my_dict[dict_y["max_velocity"]] = [max(v) for v in velocities]
 
@@ -784,9 +779,7 @@ class ResultsWindow(PamhyrWindow):
             my_dict[dict_y["froude"]] = list(
                 map(
                     lambda p:
-                        p.geometry.speed(
+                        p.get_ts_key(timestamp, "V") /
-                            p.get_ts_key(timestamp, "Q"),
-                            p.get_ts_key(timestamp, "Z")) /
                         sqrt(9.81 * (
                             p.geometry.wet_area(
                                 p.get_ts_key(timestamp, "Z")) /
@@ -818,6 +811,7 @@ class ResultsWindow(PamhyrWindow):
         my_dict[dict_x["time"]] = ts
         z = profile.get_key("Z")
         q = profile.get_key("Q")
+        v = profile.get_key("V")
         if "bed_elevation" in y:
             my_dict[dict_y["bed_elevation"]] = [
                 profile.geometry.z_min()] * len(ts)
@@ -826,12 +820,7 @@ class ResultsWindow(PamhyrWindow):
         if "water_elevation" in y:
             my_dict[dict_y["water_elevation"]] = z
         if "velocity" in y:
-            my_dict[dict_y["velocity"]] = list(
+            my_dict[dict_y["velocity"]] = v
-                map(
-                    lambda q, z: profile.geometry.speed(q, z),
-                    q, z
-                )
-            )
         if "depth" in y:
             my_dict[dict_y["depth"]] = list(
                 map(lambda z: profile.geometry.max_water_depth(z), z)
@@ -842,12 +831,11 @@ class ResultsWindow(PamhyrWindow):
             )
         if "froude" in y:
             my_dict[dict_y["froude"]] = list(
-                map(lambda z, q:
+                map(lambda z, v:
-                    profile.geometry.speed(q, z) /
+                    v / sqrt(9.81 * (
-                    sqrt(9.81 * (
                         profile.geometry.wet_area(z) /
                         profile.geometry.wet_width(z))
-                    ), z, q)
+                    ), z, v)
                 )
         if "wet_area" in y:
             my_dict[dict_y["wet_area"]] = list(

src/View/Results/WindowAdisTS.py

@@ -376,7 +376,9 @@ class ResultsWindowAdisTS(PamhyrWindow):
         table = self.find(QTableView, f"tableView_pollutants")
         indexes = table.selectedIndexes()
         if len(indexes) != 0:
-            self.pollutant_label = [self._results.pollutants_list[i.row()+1] for i in indexes]
+            self.pollutant_label = [
+                self._results.pollutants_list[i.row()+1] for i in indexes
+            ]
 
         return (f"Reach: {reach.name} | " +
                 f"Profile: {pname} | " +