debug compare results with bedload

dev_dylan
Theophile Terraz 2026-04-03 13:51:15 +02:00
parent 5b497d61f0
commit 2bcfa4d121
1 changed files with 46 additions and 84 deletions

View File

@ -88,86 +88,27 @@ 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")))
# 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
# berdrock = self.sl_compute_bedrock(reach)
# sl = reach.profile(profile).get_key("sl")
zfd = reach.profile(profile).get_key("zfd")
return zfd
def _draw_rk(self):
results = self.data[self._current_res_id]
@ -176,7 +117,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 +157,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 +200,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 +569,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 +594,12 @@ 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_sediment():
if reach.has_bedload():
z_min1 = self.draw_bottom_with_bedload(reach1)
z_min2 = self.draw_bottom_with_bedload(reach2)
else:
@ -896,7 +858,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 +876,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 +1088,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 +1101,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,