debug compare results with bedload

master
Theophile Terraz 2026-04-03 13:51:15 +02:00
parent 9670bde56f
commit 0886956ec8
1 changed files with 46 additions and 84 deletions

View File

@ -88,86 +88,27 @@ class CustomPlot(PamhyrPlot):
self.lines = {} self.lines = {}
def draw_bottom_with_bedload(self, reach): def draw_bottom_with_bedload(self, reach):
self._bedrock = self.sl_compute_bedrock(reach)
rk = reach.geometry.get_rk() 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( z = list(
map( map(
lambda z, sl: reduce( lambda p: p.get_ts_key(
lambda z, h: z + h[0], self._current_timestamp, "zfd"
sl, z
), ),
z_br, # Bedrock elevation reach.profiles
sl # Current sediment layers
) )
) )
return z 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].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): def get_ts_zmin(self, profile, res_id):
results = self.data[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) reach = results.river.reach(self._current_reach)
berdrock = self.sl_compute_bedrock(reach) # berdrock = self.sl_compute_bedrock(reach)
sl = reach.profile(profile).get_key("sl") # sl = reach.profile(profile).get_key("sl")
zfd = reach.profile(profile).get_key("zfd")
ts_z_bedrock = [berdrock[profile]]*nt return zfd
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): def _draw_rk(self):
results = self.data[self._current_res_id] results = self.data[self._current_res_id]
@ -176,7 +117,7 @@ class CustomPlot(PamhyrPlot):
reach1 = self.data[0].river.reach(self._current_reach) reach1 = self.data[0].river.reach(self._current_reach)
reach2 = self.data[1].river.reach(self._current_reach) reach2 = self.data[1].river.reach(self._current_reach)
rk = reach.geometry.get_rk() rk = reach.geometry.get_rk()
if reach.has_sediment(): if reach.has_bedload():
z_min = self.draw_bottom_with_bedload(reach) z_min = self.draw_bottom_with_bedload(reach)
else: else:
z_min = reach.geometry.get_z_min() z_min = reach.geometry.get_z_min()
@ -216,14 +157,32 @@ class CustomPlot(PamhyrPlot):
ax = self._axes[unit["bed_elevation"]] ax = self._axes[unit["bed_elevation"]]
if self._current_res_id < 2: if self._current_res_id < 2:
if reach.has_bedload():
dz = self.draw_bottom_with_bedload(reach)
else:
dz = z_min
line = ax.plot( line = ax.plot(
rk, z_min, rk, dz,
color='grey', lw=1., color='grey', lw=1.,
) )
else: else:
if reach.has_sediment(): if reach.has_bedload():
z_min1 = self.draw_bottom_with_bedload(reach1) z_min1 = list(
z_min2 = self.draw_bottom_with_bedload(reach2) 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: else:
z_min1 = reach1.geometry.get_z_min() z_min1 = reach1.geometry.get_z_min()
z_min2 = reach2.geometry.get_z_min() z_min2 = reach2.geometry.get_z_min()
@ -241,7 +200,7 @@ class CustomPlot(PamhyrPlot):
self.lines["bed_elevation"] = line self.lines["bed_elevation"] = line
if (self._envelop and if (self._envelop and
reach.has_sediment() and reach.has_bedload() and
self._current_res_id < 2): self._current_res_id < 2):
ax = self._axes[unit["bed_elevation_envelop"]] ax = self._axes[unit["bed_elevation_envelop"]]
@ -618,7 +577,7 @@ class CustomPlot(PamhyrPlot):
rk = reach.geometry.get_rk() rk = reach.geometry.get_rk()
z_min = reach.geometry.get_z_min() z_min = reach.geometry.get_z_min()
if reach.has_sediment(): if reach.has_bedload():
z_min = self.draw_bottom_with_bedload(reach) z_min = self.draw_bottom_with_bedload(reach)
else: else:
z_min = reach.geometry.get_z_min() z_min = reach.geometry.get_z_min()
@ -643,9 +602,12 @@ class CustomPlot(PamhyrPlot):
) )
if "bed_elevation" in self._y: if "bed_elevation" in self._y:
if self._current_res_id < 2: if self._current_res_id < 2:
if reach.has_bedload():
dz = self.draw_bottom_with_bedload(reach)
else:
dz = z_min dz = z_min
else: else:
if reach.has_sediment(): if reach.has_bedload():
z_min1 = self.draw_bottom_with_bedload(reach1) z_min1 = self.draw_bottom_with_bedload(reach1)
z_min2 = self.draw_bottom_with_bedload(reach2) z_min2 = self.draw_bottom_with_bedload(reach2)
else: else:
@ -904,7 +866,7 @@ class CustomPlot(PamhyrPlot):
v = profile.get_key("V") v = profile.get_key("V")
z_min = profile.geometry.z_min() z_min = profile.geometry.z_min()
if self._current_res_id < 2: if self._current_res_id < 2:
if reach.has_sediment(): if reach.has_bedload():
ts_z_min = self.get_ts_zmin( ts_z_min = self.get_ts_zmin(
self._current_profile_id, self._current_res_id) self._current_profile_id, self._current_res_id)
else: else:
@ -922,9 +884,9 @@ class CustomPlot(PamhyrPlot):
ax = self._axes[unit["bed_elevation"]] ax = self._axes[unit["bed_elevation"]]
if self._current_res_id == 2: if self._current_res_id == 2:
if reach.has_sediment(): if reach.has_bedload():
ts_z_min1 = self.get_ts_zmin(self._current_profile_id1, 0) ts_z_min1 = self.get_ts_zmin(self._current_profile_id, 0)
ts_z_min2 = self.get_ts_zmin(self._current_profile_id2, 1) ts_z_min2 = self.get_ts_zmin(self._current_profile_id, 1)
ts_z_min = list( ts_z_min = list(
map( map(
lambda x, y: x - y, lambda x, y: x - y,
@ -1140,7 +1102,7 @@ class CustomPlot(PamhyrPlot):
z = profile.get_key("Z") z = profile.get_key("Z")
v = profile.get_key("V") v = profile.get_key("V")
if self._current_res_id < 2: if self._current_res_id < 2:
if reach.has_sediment(): if reach.has_bedload():
ts_z_min = self.get_ts_zmin( ts_z_min = self.get_ts_zmin(
self._current_profile_id, self._current_res_id) self._current_profile_id, self._current_res_id)
else: else:
@ -1153,9 +1115,9 @@ class CustomPlot(PamhyrPlot):
) )
if "bed_elevation" in self._y: if "bed_elevation" in self._y:
if self._current_res_id == 2: if self._current_res_id == 2:
if reach.has_sediment(): if reach.has_bedload():
ts_z_min1 = self.get_ts_zmin(self._current_profile_id1, 0) ts_z_min1 = self.get_ts_zmin(self._current_profile_id, 0)
ts_z_min2 = self.get_ts_zmin(self._current_profile_id2, 1) ts_z_min2 = self.get_ts_zmin(self._current_profile_id, 1)
ts_z_min = list( ts_z_min = list(
map( map(
lambda x, y: x - y, lambda x, y: x - y,