Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

update weigths between zmin and zmax in uv points #200

Merged
merged 7 commits into from
Sep 4, 2024
Merged

update weigths between zmin and zmax in uv points #200

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
2d32528
update weigths between zmin and zmax in uv points
maartenvanormondt May 15, 2024
33be56e
run pre-commit
roeldegoede Aug 26, 2024
6c28129
Merge branch 'main' of https://github.com/Deltares/hydromt_sfincs int…
roeldegoede Aug 26, 2024
0c2eb7d
added q_table_option for backward compatability and changed max_gradi…
roeldegoede Aug 28, 2024
4200bcd
add valueerror for max_gradient when using old q_table_option
roeldegoede Aug 29, 2024
cbc5923
Update sfincs.py
Leynse Sep 4, 2024
9c422ab
Update sfincs.py
Leynse Sep 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 16 additions & 1 deletion hydromt_sfincs/sfincs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -619,8 +619,10 @@ def setup_subgrid(
nbins: int = None,
nr_subgrid_pixels: int = 20,
nrmax: int = 2000, # blocksize
max_gradient: float = 5.0,
max_gradient: float = 99999.0,
z_minimum: float = -99999.0,
huthresh: float = 0.01,
q_table_option: int = 2,
manning_land: float = 0.04,
manning_sea: float = 0.02,
rgh_lev_land: float = 0.0,
Expand Down Expand Up @@ -710,6 +712,12 @@ def setup_subgrid(
to prevent numerical stability problems, by default 5.0
z_minimum : float, optional
Minimum depth in the subgrid tables, by default -99999.0
huthresh : float, optional
Threshold depth in SFINCS model, by default 0.01 m
q_table_option : int, optional
Option for the computation of the representative roughness and conveyance depth at u/v points, by default 2.
1: "old" weighting method, compliant with SFINCS < v2.1.1, taking the avarage of the adjacent cells
2: "improved" weighting method, recommended for SFINCS >= v2.1.1, that takes into account the wet fractions of the adjacent cells
manning_land, manning_sea : float, optional
Constant manning roughness values for land and sea, by default 0.04 and 0.02 s.m-1/3
Note that these values are only used when no Manning's n datasets are provided,
Expand Down Expand Up @@ -753,6 +761,11 @@ def setup_subgrid(
)
nlevels = nbins

if q_table_option == 1 and max_gradient > 20.0:
raise ValueError(
"For the old q_table_option, a max_gradient of 5.0 is recommended to improve numerical stability"
)

if self.grid_type == "regular":
self.reggrid.subgrid.build(
da_mask=self.mask,
Expand All @@ -765,6 +778,8 @@ def setup_subgrid(
nrmax=nrmax,
max_gradient=max_gradient,
z_minimum=z_minimum,
huthresh=huthresh,
q_table_option=q_table_option,
manning_land=manning_land,
manning_sea=manning_sea,
rgh_lev_land=rgh_lev_land,
Expand Down
105 changes: 68 additions & 37 deletions hydromt_sfincs/subgrid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -371,9 +371,10 @@ def build(
nlevels: int = 10,
nr_subgrid_pixels: int = 20,
nrmax: int = 2000,
max_gradient: float = 5.0,
max_gradient: float = 99999.0,
z_minimum: float = -99999.0,
huthresh: float = 0.01,
q_table_option: int = 2,
manning_land: float = 0.04,
manning_sea: float = 0.02,
rgh_lev_land: float = 0.0,
Expand Down Expand Up @@ -428,6 +429,10 @@ def build(
Minimum depth in the subgrid tables, by default -99999.0
huthresh : float, optional
Threshold depth in SFINCS model, by default 0.01 m
q_table_option : int, optional
Option for the computation of the representative roughness and conveyance depth at u/v points, by default 2.
1: "old" weighting method, compliant with SFINCS < v2.1.0, taking the avarage of the adjecent cells
2: "improved" weighting method, recommended for SFINCS >= v2.1.0, that takes into account the wet fractions of the adjacent cells
manning_land, manning_sea : float, optional
Constant manning roughness values for land and sea,
by default 0.04 and 0.02 s.m-1/3
Expand Down Expand Up @@ -716,6 +721,7 @@ def build(
yg,
max_gradient,
huthresh,
q_table_option,
da_mask.raster.crs.is_geographic,
)

Expand Down Expand Up @@ -797,6 +803,7 @@ def process_tile_regular(
yg,
max_gradient,
huthresh,
q_table_option,
is_geographic=False,
):
"""calculate subgrid properties for a single tile"""
Expand Down Expand Up @@ -864,7 +871,7 @@ def process_tile_regular(
manning = manning_grid[nn : nn + refi, mm : mm + refi]
manning = np.transpose(manning)
zmin, zmax, havg, nrep, pwet, ffit, navg, zz = subgrid_q_table(
zgu.flatten(), manning.flatten(), nlevels, huthresh
zgu.flatten(), manning.flatten(), nlevels, huthresh, q_table_option
)
u_zmin[n, m] = zmin
u_zmax[n, m] = zmax
Expand All @@ -880,7 +887,7 @@ def process_tile_regular(
zgu = zg[nn : nn + refi, mm : mm + refi]
manning = manning_grid[nn : nn + refi, mm : mm + refi]
zmin, zmax, havg, nrep, pwet, ffit, navg, zz = subgrid_q_table(
zgu.flatten(), manning.flatten(), nlevels, huthresh
zgu.flatten(), manning.flatten(), nlevels, huthresh, q_table_option
)
v_zmin[n, m] = zmin
v_zmax[n, m] = zmax
Expand Down Expand Up @@ -952,7 +959,7 @@ def subgrid_v_table(
zvolmin: float
minimum elevation value to use for volume calculation (typically -20 m)
max_gradient: float
maximum gradient to use for volume calculation (typically 0.1)
maximum gradient to use for volume calculation

Return
------
Expand Down Expand Up @@ -1000,7 +1007,11 @@ def subgrid_v_table(

@njit
def subgrid_q_table(
elevation: np.ndarray, manning: np.ndarray, nlevels: int, huthresh: float
elevation: np.ndarray,
manning: np.ndarray,
nlevels: int,
huthresh: float,
option: int,
):
"""
map vector of elevation values into a hypsometric hydraulic radius - depth relationship for one u/v point
Expand All @@ -1010,6 +1021,8 @@ def subgrid_q_table(
manning : np.ndarray (nr of pixels in one cell) containing subgrid manning roughness values for one grid cell [s m^(-1/3)]
nlevels : int, number of vertical levels [-]
huthresh : float, threshold depth [m]
option : int, option to use "old" or "new" method for computing conveyance depth at u/v points

Returns
-------
zmin : float, minimum elevation [m]
Expand Down Expand Up @@ -1053,46 +1066,62 @@ def subgrid_q_table(
# Determine level size (metres)
dlevel = (zmax - zmin) / (nlevels - 1)

# Grid mean roughness
navg = np.mean(manning)
# Option can be either 1 ("old, compliant with SFINCS < v2.1.0") or 2 ("new", recommended SFINCS >= v2.1.0)
option = option

# Loop through levels
for ilevel in range(nlevels):
# Top of level
zlevel = zmin + ilevel * dlevel
zz[ilevel] = zlevel

# ibelow = np.where(elevation<=zlevel) # index of pixels below level level
h = np.maximum(zlevel - elevation, 0.0) # water depth in each pixel
iwet = np.where(zlevel - elevation > -1.0e-6)[0] # indices of wet pixels
hmean = np.mean(h)
havg[ilevel] = hmean # conveyance depth
pwet[ilevel] = len(iwet) / n # wet fraction
for ibin in range(nlevels):
# Top of bin
zbin = zmin + ibin * dlevel
zz[ibin] = zbin

h = np.maximum(zbin - elevation, 0.0) # water depth in each pixel

pwet[ibin] = (zbin - elevation > -1.0e-6).sum() / n

# Side A
h_a = np.maximum(
zlevel - dd_a, 0.0
zbin - dd_a, 0.0
) # Depth of all pixels (but set min pixel height to zbot). Can be negative, but not zero (because zmin = zbot + huthresh, so there must be pixels below zb).
q_a = h_a ** (5.0 / 3.0) / manning_a # Determine 'flux' for each pixel
q_a = np.mean(q_a) # Wet-average flux through all the pixels
q_a = np.mean(q_a) # Grid-average flux through all the pixels
h_a = np.mean(h_a) # Grid-average depth through all the pixels

# Side B
h_b = np.maximum(
zlevel - dd_b, 0.0
zbin - dd_b, 0.0
) # Depth of all pixels (but set min pixel height to zbot). Can be negative, but not zero (because zmin = zbot + huthresh, so there must be pixels below zb).
q_b = h_b ** (5.0 / 3.0) / manning_b # Determine 'flux' for each pixel
q_b = np.mean(q_b) # Wet-average flux through all the pixels

q_ab = np.minimum(q_a, q_b)

q_all = h ** (5.0 / 3.0) / manning # Determine 'flux' for each pixel
q_all = np.mean(q_all) # Wet-average flux through all the pixels

# Weighted average of q_ab and q_all
w = (ilevel) / (nlevels - 1)
q = (1.0 - w) * q_ab + w * q_all

nrep[ilevel] = hmean ** (5.0 / 3.0) / q # Representative n for qmean and hmean
q_b = np.mean(q_b) # Grid-average flux through all the pixels
h_b = np.mean(h_b) # Grid-average depth through all the pixels

# Compute q and h
q_all = np.mean(
h ** (5.0 / 3.0) / manning
) # Determine grid average 'flux' for each pixel
h_all = np.mean(h) # grid averaged depth of A and B combined
q_min = np.minimum(q_a, q_b)
h_min = np.minimum(h_a, h_b)

if option == 1:
# Use old 1 option (weighted average of q_ab and q_all) option (min at bottom bin, mean at top bin)
w = (ibin) / (
nlevels - 1
) # Weight (increase from 0 to 1 from bottom to top bin)
q = (1.0 - w) * q_min + w * q_all # Weighted average of q_min and q_all
hmean = h_all

elif option == 2:
# Use newer 2 option (minimum of q_a an q_b, minimum of h_a and h_b increasing to h_all, using pwet for weighting) option
pwet_a = (zbin - dd_a > -1.0e-6).sum() / (n / 2)
pwet_b = (zbin - dd_b > -1.0e-6).sum() / (n / 2)
# Weight increases linearly from 0 to 1 from bottom to top bin use percentage wet in sides A and B
w = 2 * np.minimum(pwet_a, pwet_b) / (pwet_a + pwet_b)
q = (1.0 - w) * q_min + w * q_all # Weighted average of q_min and q_all
hmean = (1.0 - w) * h_min + w * h_all # Weighted average of h_min and h_all

havg[ibin] = hmean # conveyance depth
nrep[ibin] = hmean ** (5.0 / 3.0) / q # Representative n for qmean and hmean

nrep_top = nrep[-1]
havg_top = havg[-1]
Expand All @@ -1101,14 +1130,16 @@ def subgrid_q_table(

# Determine nfit at zfit
zfit = zmax + zmax - zmin
h = np.maximum(zfit - elevation, 0.0) # water depth in each pixel
hfit = (
havg_top + zmax - zmin
) # mean water depth in cell as computed in SFINCS (assuming linear relation between water level and water depth above zmax)
q = h ** (5.0 / 3.0) / manning # unit discharge in each pixel
qmean = np.mean(q) # combined unit discharge for cell

nfit = hfit ** (5.0 / 3.0) / qmean
# Compute q and navg
h = np.maximum(zfit - elevation, 0.0) # water depth in each pixel
q = np.mean(h ** (5.0 / 3.0) / manning) # combined unit discharge for cell
navg = np.mean(manning)

nfit = hfit ** (5.0 / 3.0) / q

# Actually apply fit on gn2 (this is what is used in sfincs)
gnavg2 = 9.81 * navg**2
Expand Down
Loading