diff --git a/src/Model/Geometry/ProfileXYZ.py b/src/Model/Geometry/ProfileXYZ.py index bb3faa0b7eba7908f72b5eb4d004e94ec24cbab0..916d579f0283bcaf53535fb8726c44661d03cc9d 100644 --- a/src/Model/Geometry/ProfileXYZ.py +++ b/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 compute_wet_area(self, z): + area = 0.0 + if len(self.tab.L) > 0: + if z < self.tab.z[0]: + return 0.0 + 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): - lines = self.wet_lines(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 + 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 + 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): diff --git a/src/Solver/Mage.py b/src/Solver/Mage.py index 62a906b7c57e873303f10b3299e152c239224927..414eb649deb7b956bf733d4bb760e38d6aadc4b3 100644 --- a/src/Solver/Mage.py +++ b/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 diff --git a/src/View/About/Window.py b/src/View/About/Window.py index a3605ea859e4115343fbd048060c347100e76f66..a110c2ea88062ccbd48950a7b70818e67c5a556d 100644 --- a/src/View/About/Window.py +++ b/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 diff --git a/src/View/DIFAdisTS/Table.py b/src/View/DIFAdisTS/Table.py index 27cfad4ff1a7a623b5a9d2c07f99d7855b808623..5c24b02702b197faf1aa298923aa43964e3c9d08 100644 --- a/src/View/DIFAdisTS/Table.py +++ b/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, diff --git a/src/View/DIFAdisTS/Window.py b/src/View/DIFAdisTS/Window.py index 4100125b19b54d287057bde96bf4036d861c2c91..ed16ae3b9b5b6864f84a3913b011438b4fef6928 100644 --- a/src/View/DIFAdisTS/Window.py +++ b/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) diff --git a/src/View/MainWindow.py b/src/View/MainWindow.py index 0fd12fcf6421c29baa9a10b972a8e5fadbddda8f..3fe55953d93cf0c33fa709e15cf3439b9c6295d4 100644 --- a/src/View/MainWindow.py +++ b/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): diff --git a/src/View/Results/CustomPlot/Plot.py b/src/View/Results/CustomPlot/Plot.py index 0d77df14f94d7c7e4924706557ec08f44a69be02..446c661a442c3cc4e67bb51e3e6f979db2dfa9bb 100644 --- a/src/View/Results/CustomPlot/Plot.py +++ b/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( - lambda p: list(map( - lambda q, z: - p.geometry.speed(q, z), - p.get_key("Q"), p.get_key("Z") - )), reach.profiles + velocities = list( + map( + lambda p: p.get_key("V"), + reach.profiles ) ) @@ -440,9 +436,7 @@ class CustomPlot(PamhyrPlot): fr = list( map( lambda p: - p.geometry.speed( - p.get_ts_key(self._timestamp, "Q"), - p.get_ts_key(self._timestamp, "Z")) / + p.get_ts_key(self._timestamp, "V") / 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( - 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) + self.lines["velocity"][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, "Q"), - p.get_ts_key(self._timestamp, "Z")) / + p.get_ts_key(self._timestamp, "V") / 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: - profile.geometry.speed(q, z) / + map(lambda z, v: + 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: - profile.geometry.speed(q, z) / + map(lambda z, v: + 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) diff --git a/src/View/Results/Table.py b/src/View/Results/Table.py index 66c1419c21bcfca9ce49f33225bfa46c22ee5641..229178b3da3246cc1d80c6ce8f260ef58b4620a7 100644 --- a/src/View/Results/Table.py +++ b/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)) diff --git a/src/View/Results/Window.py b/src/View/Results/Window.py index 13f15138fc6650e66f892f166807a9cd2a0c20fd..5f52f061878dc87174621b8b11c1afe4f5e96976 100644 --- a/src/View/Results/Window.py +++ 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( - p.get_ts_key(timestamp, "Q"), - p.get_ts_key(timestamp, "Z")), + lambda p: p.get_ts_key(timestamp, "V"), reach.profiles ) ) if envelop: velocities = list(map( - lambda p: list(map( - lambda q, z: - p.geometry.speed(q, z), - p.get_key("Q"), p.get_key("Z") - )), reach.profiles + lambda p: p.get_key("V"), + 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, "Q"), - p.get_ts_key(timestamp, "Z")) / + p.get_ts_key(timestamp, "V") / 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( - map( - lambda q, z: profile.geometry.speed(q, z), - q, z - ) - ) + my_dict[dict_y["velocity"]] = v 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: - profile.geometry.speed(q, z) / - sqrt(9.81 * ( + map(lambda z, v: + v / 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( diff --git a/src/View/Results/WindowAdisTS.py b/src/View/Results/WindowAdisTS.py index c7d9ab8bfadf668e94ed3bac56f632ec51f2ed7a..9d9988268c238e5ae874fb0dac12d8976ff2f367 100644 --- a/src/View/Results/WindowAdisTS.py +++ b/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} | " +