feat(channel_utils): move open channel flow utils from modflow6

Author: wpbonelli

autotest/test_channel_utils.py

@@ -0,0 +1,275 @@
+import math
+
+import numpy as np
+import pytest
+
+from flopy.utils.channel_utils import (
+    get_discharge,
+    get_discharge_rect,
+    get_segment_wetted_area,
+    get_segment_wetted_perimeter,
+    get_segment_wetted_station,
+    get_wetted_area,
+    get_wetted_perimeter,
+)
+
+
+def test_get_segment_wetted_station_all_dry():
+    depth = 10
+    p0 = (0, 12)
+    p1 = (10, 11)
+    x0, x1 = get_segment_wetted_station(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    # zero-length segment at x0
+    assert x0 == x1 == p0[0]
+
+
+def test_get_segment_wetted_station_partial():
+    depth = 10
+
+    # left bank (sloping downwards to the right)
+    p0 = (0, 12)
+    p1 = (10, 8)
+    x0, x1 = get_segment_wetted_station(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    # left endpoint should be moved to the right
+    assert x0 != x1
+    assert x0 != p0[0]
+    assert x1 == p1[0]
+    assert x0 == (p1[0] - p0[0]) / 2
+
+    # right bank (sloping upwards to the right)
+    p0 = (0, 8)
+    p1 = (10, 12)
+    x0, x1 = get_segment_wetted_station(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    # right endpoint should be moved to the left
+    assert x0 != x1
+    assert x0 == p0[0]
+    assert x1 != p1[0]
+    assert x1 == (p1[0] - p0[0]) / 2
+
+
+def test_get_segment_wetted_station_submerged():
+    depth = 13
+    p0 = (0, 12)
+    p1 = (10, 11)
+    x0, x1 = get_segment_wetted_station(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    # entire segment should be returned
+    assert x0 == p0[0]
+    assert x1 == p1[0]
+
+
+def test_get_segment_wetted_perimeter_dry():
+    depth = 10
+    p0 = (0, 12)
+    p1 = (10, 11)
+    perim = get_segment_wetted_perimeter(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    assert perim == 0
+
+
+point_lists = [
+    # single segments
+    [(0, -1), (1, -1)],
+    [(0, 0), (1, 1)],
+    [(0, 1), (2, -1)],
+    [(0, -1), (2, 1)],
+    [(0, 1), (1, 0)],
+    # channels with multiple segments
+    [(0, -1), (1, -1), (2, -1)],  # flat
+    [(0, 0), (1, -1), (2, 0)],  # triangular
+    [(0, -1), (1, -2), (2, -2), (3, -1)],  # trapezoidal
+    [(0, -1), (1, -2), (2, -4), (3, -4), (4, -1)],  # complex
+]
+
+
+@pytest.mark.parametrize(
+    "depth, p0, p1",
+    [
+        (0, point_lists[0][0], point_lists[0][1]),
+        (0, point_lists[1][0], point_lists[1][1]),
+        (0, point_lists[2][0], point_lists[2][1]),
+        (3, point_lists[3][0], point_lists[3][1]),
+    ],
+)
+def test_get_segment_wetted_perimeter(depth, p0, p1):
+    wp = get_segment_wetted_perimeter(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+
+    xlen = abs(p1[0] - p0[0])
+    hlen = abs(p1[1] - p0[1])
+    hmax = max([p0[1], p1[1]])
+    hmin = min([p0[1], p1[1]])
+    seg_len = math.sqrt(hlen**2 + xlen**2)
+
+    # if segment is fully wetted, wetted perimeter is just its length
+    if depth >= hmax:
+        # expect perimeter 0 if the water surface is level with a flat segment
+        if depth == hmin == hmax:
+            assert wp == 0
+        else:
+            assert wp == seg_len
+
+    # if segment is partially submerged, wetted perimeter should be
+    # less than the length of the segment but greater than zero
+    elif depth > hmin:
+        assert wp > 0
+        assert wp < seg_len
+
+    # if segment is completely dry, wetted perimeter should be zero
+    else:
+        assert wp == 0
+
+
+@pytest.mark.parametrize(
+    "depth, p0, p1",
+    [
+        (0, point_lists[0][0], point_lists[0][1]),
+        (0, point_lists[1][0], point_lists[1][1]),
+        (0, point_lists[2][0], point_lists[2][1]),
+        (3, point_lists[3][0], point_lists[3][1]),
+    ],
+)
+def test_get_segment_wetted_area(depth, p0, p1):
+    wa = get_segment_wetted_area(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+
+    xlen = abs(p1[0] - p0[0])
+    hlen = abs(p1[1] - p0[1])
+    hmax = max([p0[1], p1[1]])
+    hmin = min([p0[1], p1[1]])
+    seg_len = math.sqrt(hlen**2 + xlen**2)
+    tri_area = 0 if hlen == 0 else (0.5 * hlen * xlen)
+
+    # if segment is submerged wetted area is that of quadrilateral
+    # formed by the segment, the water surface, and vertical sides
+    if depth > hmax:
+        rect_area = xlen * (depth - hmax)
+        assert wa == (rect_area + tri_area)
+
+    # if segment is partially submerged, wetted area should be
+    # less than the area of the triangle formed by the segment
+    # and the water surface, with one vertical side
+    elif depth > hmin:
+        assert wa > 0
+        assert wa < tri_area
+
+    # if segment is completely dry, wetted area should be zero
+    else:
+        assert wa == 0
+
+
+@pytest.mark.parametrize("points", [np.array(pts) for pts in point_lists])
+@pytest.mark.parametrize("depth", [0, 1, -1, -2])
+def test_get_wetted_perimeter(depth, points):
+    def total_perim(pts):
+        return sum(
+            [
+                math.sqrt(
+                    (pts[i][0] - pts[i - 1][0]) ** 2
+                    + (pts[i][1] - pts[i - 1][1]) ** 2
+                )
+                for i in range(1, len(pts))
+            ]
+        )
+
+    wp = get_wetted_perimeter(
+        x=points[:, 0], h=points[:, 1], depth=depth, verbose=True
+    )
+
+    hmax = max(points[:, 1])
+    hmin = min(points[:, 1])
+    total_perim = total_perim(points)
+
+    # if all segments are submerged, wetted perimeter is total perimeter
+    if depth >= hmax:
+        # expect perimeter 0 if the water surface is level with a flat channel
+        if all(p == depth for p in points[:, 1]):
+            assert wp == 0
+        else:
+            assert wp == total_perim
+
+    # if at least some segments are at least partially submerged...
+    elif depth > hmin:
+        assert wp > 0
+        assert wp < total_perim
+
+        def patch(x0, x1, h0, h1, depth):
+            # TODO: refactor get_segment_wetted_perimeter() to handle partial wetting
+            #  internally so separate get_segment_wetted_station() is no longer needed?
+
+            x0, x1 = get_segment_wetted_station(
+                x0=x0, x1=x1, h0=h0, h1=h1, depth=depth
+            )
+            return get_segment_wetted_perimeter(
+                x0=x0, x1=x1, h0=h0, h1=h1, depth=depth
+            )
+
+        assert np.isclose(
+            wp,
+            sum(
+                [
+                    patch(
+                        x0=points[i][0],
+                        x1=points[i + 1][0],
+                        h0=points[i][1],
+                        h1=points[i + 1][1],
+                        depth=depth,
+                    )
+                    for i in range(0, len(points) - 1)
+                ]
+            ),
+        )
+
+    # if all segments are dry, wetted perimeter is zero
+    else:
+        assert wp == 0
+
+
+@pytest.mark.skip(reason="todo")
+def test_get_wetted_area():
+    pass
+
+
+def test_get_discharge_rect():
+    width = 0.5
+    depth = 0.25
+    roughness = 0.022
+    slope = 0.005
+    conv = 1.0
+    expected = 0.1
+    actual = get_discharge_rect(width, depth, roughness, slope, conv)
+    assert np.isclose(expected, actual, rtol=1e-2)
+
+
+@pytest.mark.xfail(reason="debug")
+def test_get_discharge():
+    width = 0.5
+    depth = 0.25
+    roughness = 0.022
+    slope = 0.005
+    conv = 1.0
+    expected = 0.1
+    x = np.array([width])
+    h = np.array([depth])
+    actual = get_discharge(x, h, depth, roughness, slope, conv)
+    assert np.isclose(expected, actual, rtol=1e-2)
+
+
+@pytest.mark.skip(reason="todo")
+def test_get_depth():
+    pass
+
+
+@pytest.mark.skip(reason="todo")
+def test_get_depths():
+    pass