From 55672ea79fb24c700eedb2655c79c31faa5abaac Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Thu, 22 Dec 2022 15:19:19 -0500 Subject: [PATCH] feat(utils): move cross section utils from modflow-devtools --- flopy/utils/crosssection.py | 276 ++++++++++++++++++++++++++++++++++++ 1 file changed, 276 insertions(+) create mode 100644 flopy/utils/crosssection.py diff --git a/flopy/utils/crosssection.py b/flopy/utils/crosssection.py new file mode 100644 index 0000000000..725c1e250d --- /dev/null +++ b/flopy/utils/crosssection.py @@ -0,0 +1,276 @@ +import numpy as np + +# power for Manning's hydraulic radius term +_mpow = 2.0 / 3.0 + + +def calculate_rectchan_mannings_discharge( + conversion_factor, roughness, slope, width, depth +): + """ + Calculate Manning's discharge for a rectangular channel. + + """ + area = width * depth + return conversion_factor * area * depth**_mpow * slope**0.5 / roughness + + +# n-point cross-section functions +def get_wetted_station( + x0, + x1, + h0, + h1, + depth, +): + """Get the wetted length in the x-direction""" + # -- calculate the minimum and maximum depth + hmin = min(h0, h1) + hmax = max(h0, h1) + + # -- if depth is less than or equal to the minimum value the + # station length (xlen) is zero + if depth <= hmin: + x1 = x0 + # -- if depth is between hmin and hmax, station length is less + # than h1 - h0 + elif depth < hmax: + xlen = x1 - x0 + dlen = h1 - h0 + if abs(dlen) > 0.0: + slope = xlen / dlen + else: + slope = 0.0 + if h0 > h1: + dx = (depth - h1) * slope + xt = x1 + dx + xt0 = xt + xt1 = x1 + else: + dx = (depth - h0) * slope + xt = x0 + dx + xt0 = x0 + xt1 = xt + x0 = xt0 + x1 = xt1 + return x0, x1 + + +def get_wetted_perimeter( + x0, + x1, + h0, + h1, + depth, +): + # -- calculate the minimum and maximum depth + hmin = min(h0, h1) + hmax = max(h0, h1) + + # -- calculate the wetted perimeter for the segment + xlen = x1 - x0 + if xlen > 0.0: + if depth > hmax: + dlen = hmax - hmin + else: + dlen = depth - hmin + else: + if depth > hmin: + dlen = min(depth, hmax) - hmin + else: + dlen = 0.0 + return np.sqrt(xlen**2.0 + dlen**2.0) + + +def get_wetted_area(x0, x1, h0, h1, depth): + # -- calculate the minimum and maximum depth + hmin = min(h0, h1) + hmax = max(h0, h1) + + # -- calculate the wetted area for the segment + xlen = x1 - x0 + area = 0.0 + if xlen > 0.0: + # -- add the area above hmax + if depth > hmax: + area = xlen * (depth - hmax) + # -- add the area below zmax + if hmax != hmin and depth > hmin: + area += 0.5 * (depth - hmin) + return area + + +def wetted_area( + x, + h, + depth, + verbose=False, +): + area = 0.0 + if x.shape[0] == 1: + area = x[0] * depth + else: + for idx in range(0, x.shape[0] - 1): + x0, x1 = x[idx], x[idx + 1] + h0, h1 = h[idx], h[idx + 1] + + # get station data + x0, x1 = get_wetted_station(x0, x1, h0, h1, depth) + + # get wetted area + a = get_wetted_area(x0, x1, h0, h1, depth) + area += a + + # write to screen + if verbose: + print( + f"{idx}->{idx + 1} ({x0},{x1}) - " + f"perimeter={x1 - x0} - area={a}" + ) + + return area + + +def wetted_perimeter( + x, + h, + depth, + verbose=False, +): + perimeter = 0.0 + if x.shape[0] == 1: + perimeter = x[0] + else: + for idx in range(0, x.shape[0] - 1): + x0, x1 = x[idx], x[idx + 1] + h0, h1 = h[idx], h[idx + 1] + + # get station data + x0, x1 = get_wetted_station(x0, x1, h0, h1, depth) + + # get wetted perimeter + perimeter += get_wetted_perimeter(x0, x1, h0, h1, depth) + + # write to screen + if verbose: + print(f"{idx}->{idx + 1} ({x0},{x1}) - perimeter={x1 - x0}") + + return perimeter + + +def manningsq( + x, + h, + depth, + roughness=0.01, + slope=0.001, + conv=1.0, +): + if isinstance(roughness, float): + roughness = np.ones(x.shape, dtype=float) * roughness + if x.shape[0] > 1: + q = 0.0 + for i0 in range(x.shape[0] - 1): + i1 = i0 + 1 + perimeter = get_wetted_perimeter(x[i0], x[i1], h[i0], h[i1], depth) + area = get_wetted_area(x[i0], x[i1], h[i0], h[i1], depth) + if perimeter > 0.0: + radius = area / perimeter + q += ( + conv * area * radius ** _mpow * slope ** 0.5 / roughness[i0] + ) + else: + perimeter = wetted_perimeter(x, h, depth) + area = wetted_area(x, h, depth) + radius = 0.0 + if perimeter > 0.0: + radius = area / perimeter + q = conv * area * radius ** _mpow * slope ** 0.5 / roughness[0] + return q + + +def get_depths( + flows, + x, + h, + roughness=0.01, + slope=0.001, + conv=1.0, + dd=1e-4, + verbose=False, +): + if isinstance(flows, float): + flows = np.array([flows], dtype=float) + if isinstance(roughness, float): + roughness = np.ones(x.shape, dtype=float) * roughness + depths = np.zeros(flows.shape, dtype=float) + for idx, q in enumerate(flows): + depths[idx] = qtodepth( + x, + h, + q, + roughness=roughness, + slope=slope, + conv=conv, + dd=dd, + verbose=False, + ) + + return depths + + +def qtodepth( + x, + h, + q, + roughness=0.01, + slope=0.001, + conv=1.0, + dd=1e-4, + verbose=False, +): + h0 = 0.0 + q0 = manningsq( + x, + h, + h0, + roughness=roughness, + slope=slope, + conv=conv, + ) + r = q0 - q + + iter = 0 + if verbose: + print(f"iteration {iter:>2d} - residual={r}") + while abs(r) > 1e-12: + q1 = manningsq( + x, + h, + h0 + dd, + roughness=roughness, + slope=slope, + conv=conv, + ) + dq = q1 - q0 + if dq != 0.0: + derv = dd / (q1 - q0) + else: + derv = 0.0 + h0 -= derv * r + q0 = manningsq( + x, + h, + h0, + roughness=roughness, + slope=slope, + conv=conv, + ) + r = q0 - q + + iter += 1 + if verbose: + print(f"iteration {iter:>2d} - residual={r}") + if iter > 100: + break + return h0