Feat snf

autotest/test_snf_disl.py

@@ -0,0 +1,186 @@
+"""
+
+Simple 4 segment network with 5 vertices
+
+
+zero-based diagram below
+v1
+o
+ |
+  |
+   |  s0
+    |
+     |
+      |
+o------o------o------o
+v0     v2     v3     v4
+    s1     s2     s3
+
+ia  ja
+0   0 2
+2   1 2
+4   2 0 1 3
+8   3 2
+10
+
+"""
+
+import os
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+from simulation import TestSimulation
+
+ex = ["snf-disl01",]
+
+
+def build_model(idx, dir):
+
+    sim_ws = dir
+    name = ex[idx]
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, version="mf6", exe_name="mf6", sim_ws=sim_ws,
+        memory_print_option='all',
+    )
+
+    tdis = flopy.mf6.ModflowTdis(sim)
+    ems = flopy.mf6.ModflowEms(sim)
+    snf = flopy.mf6.ModflowSnf(sim, modelname=name, save_flows=True)
+
+    vertices = [
+        [0, 0., 0., 0.], 
+        [1, 0., 1., 0.],
+        [2, 1., 0., 0.],
+        [3, 2., 0., 0.],
+        [4, 3., 0., 0.],
+    ]
+    # icell1d fdc ncvert icvert
+    cell2d = [
+        [0, 0.5, 2, 1, 2],
+        [1, 0.5, 2, 0, 2],
+        [2, 0.5, 2, 2, 3],
+        [3, 0.5, 2, 3, 4],
+    ]
+
+    nodes = len(cell2d)
+    nvert = len(vertices)
+
+    disl = flopy.mf6.ModflowSnfdisl(
+        snf, 
+        nodes=nodes, 
+        nvert=nvert,
+        segment_length=1000.0,
+        tosegment=[2, 2, 3, -1],   # -1 gives 0 in one-based, which means outflow cell
+        idomain=1, 
+        vertices=vertices, 
+        cell2d=cell2d,
+    )
+
+    # note: for specifying zero-based reach number, put reach number in tuple
+    fname = f"{name}.mmr.obs.csv"
+    mmr_obs = {
+        fname: [
+            ("OUTFLOW", "EXT-OUTFLOW", (nodes - 1,)),
+        ],
+        "digits": 10,
+    }
+
+    mmr = flopy.mf6.ModflowSnfmmr(
+        snf, 
+        observations=mmr_obs,
+        print_flows=True,
+        save_flows=True,
+        iseg_order=list(range(nodes)),
+        qoutflow0=0.0,
+        k_coef=0.001, 
+        x_coef=0.2
+    )
+
+    # output control
+    oc = flopy.mf6.ModflowSnfoc(
+        snf,
+        budget_filerecord=f"{name}.bud",
+        saverecord=[("BUDGET", "ALL"), ],
+        printrecord=[("BUDGET", "ALL"), ],
+    )
+
+    flw_spd = {0: [[0, 1000.0], [1, 500.]]}
+    flw = flopy.mf6.ModflowSnfflw(
+        snf,
+        print_input=True,
+        print_flows=True,
+        stress_period_data=flw_spd,
+    )
+
+    return sim, None
+
+
+def eval_model(sim):
+    print("evaluating model...")
+
+    # read the observation output
+    name = ex[sim.idxsim]
+    fpth = os.path.join(sim.simpath, f"{name}.mmr.obs.csv")
+    obsvals = np.genfromtxt(fpth, names=True, delimiter=",")
+
+    # read the binary grid file
+    fpth = os.path.join(sim.simpath, f"{name}.disl.grb")
+    grb = flopy.mf6.utils.MfGrdFile(fpth)
+    ia = grb.ia
+    ja = grb.ja
+    assert ia.shape[0] == grb.nodes + 1, "ia in grb file is not correct size"
+
+    # read the budget file
+    fpth = os.path.join(sim.simpath, f"{name}.bud")
+    budobj = flopy.utils.binaryfile.CellBudgetFile(fpth)
+    flowja = budobj.get_data(text="FLOW-JA-FACE")
+    qstorage = budobj.get_data(text="STORAGE")
+    qflw = budobj.get_data(text="FLW")
+    qextoutflow = budobj.get_data(text="EXT-OUTFLOW")
+    qresidual = np.zeros(grb.nodes)
+
+    # check budget terms
+    for itime in range(len(flowja)):
+        print (f"evaluating timestep {itime}")
+
+        fja = flowja[itime].flatten()
+        for n in range(grb.nodes):
+            ipos = ia[n]
+            qresidual[n] = fja[ipos]
+        assert np.allclose(qresidual, 0.), "residual in flowja diagonal is not zero"
+
+        for n in range(grb.nodes):
+            qs = qstorage[itime].flatten()[n]
+            if n + 1 in qflw[itime]["node"]:
+                idx, = np.where(qflw[itime]["node"] == n + 1)
+                idx = idx[0]
+                qf = qflw[itime].flatten()["q"][idx]
+            else:
+                qf = 0.
+            qe = qextoutflow[itime].flatten()[n]
+            qdiag = fja[ia[n]]
+            print(f"{n=} {qs=} {qf=} {qe=} {qdiag=}")
+            for ipos in range(ia[n] + 1, ia[n + 1]):
+                j = ja[ipos]
+                q = fja[ipos]
+                print(f"  {ipos=} {j=} {q=}")        
+
+    return
+
+
+@pytest.mark.parametrize(
+    "idx, name",
+    list(enumerate(ex)),
+)
+def test_mf6model(idx, name, function_tmpdir, targets):
+    ws = str(function_tmpdir)
+    test = TestFramework()
+    test.build(build_model, idx, ws)
+    test.run(
+        TestSimulation(
+            name=name, exe_dict=targets, exfunc=eval_model, idxsim=idx
+        ),
+        ws,
+    )

autotest/test_snf_muskingum.py

@@ -0,0 +1,187 @@
+"""
+
+http://uon.sdsu.edu/enghydro/engineering_hydrology_09.php#muskingum
+Ponce -- Muskingum routing example
+single reach
+"""
+
+import os
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+from simulation import TestSimulation
+
+
+# data from ponce Table 9-1
+# time(d) coi2(m4/s) c1i1(m3/s) c2i1(m3/s) outflow(m3/s)
+ponce_data = """0  352.0  0  0  0  352.0
+1  587.0  76.6  107.1  199.0  382.7
+2  1353.0  176.5  178.6  216.3  571.4
+3  2725.0  355.4  411.8  323.0  1090.2
+4  4408.5  575.0  829.4  616.2  2020.6
+5  5987.0  780.9  1341.7  1142.1  3264.7
+6  6704.0  874.4  1822.1  1845.3  4541.8
+7  6951.0  906.7  2040.3  2567.1  5514.1
+8  6839.0  892.0  2115.5  3116.7  6124.2
+9  6207.0  809.6  2081.5  3461.5  6352.6
+10  5346.0  697.3  1889.1  3590.6  6177.0
+11  4560.0  594.8  1627.0  3491.4  5713.2
+12  3861.5  503.7  1387.8  3229.2  5120.7
+13  3007.0  392.2  1175.2  2894.3  4461.7
+14  2357.5  307.5  915.2  2521.8  3744.5
+15  1779.0  232.0  717.5  2116.5  3066.0
+16  1405.0  183.3  541.4  1733.0  2457.7
+17  1123.0  146.5  427.6  1389.1  1963.2
+18  952.5  124.2  341.8  1109.6  1575.6
+19  730.0  95.2  289.9  890.6  1275.7
+20  605.0  78.9  222.2  721.0  1022.1
+21  514.0  67.1  184.1  577.7  828.9
+22  422.0  55.1  156.4  468.5  680.0
+23  352.0  45.9  128.4  384.4  558.7
+24  352.0  45.9  107.1  315.8  468.8
+25  352.0  45.9  107.1  265.0  418.0"""
+
+
+def get_ponce_data():
+    qinflow = []
+    qoutflow = []
+    time_days = []
+    for line in ponce_data.split("\n"):
+        itime, inflow, c0, c1, c2, outflow = line.strip().split("  ")
+        time_days.append(float(itime))
+        qinflow.append(float(inflow))
+        qoutflow.append(float(outflow))
+    return time_days, qinflow, qoutflow
+
+
+ex = ["ponce01",]
+
+
+def build_model(idx, dir):
+
+    sim_ws = dir
+    name = ex[idx]
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, version="mf6", exe_name="mf6", sim_ws=sim_ws,
+        memory_print_option='all',
+    )
+    nper = 25
+    tdis_rc = [(1., 1, 1.0) for ispd in range(nper)]
+    tdis = flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=tdis_rc)
+    ems = flopy.mf6.ModflowEms(sim)
+    snf = flopy.mf6.ModflowSnf(sim, modelname=name, save_flows=True)
+
+    vertices = None
+    cell2d = None
+    nvert = None
+
+    nodes = 1
+    channel_length = 500. # miles
+    channel_length = channel_length * 5280. / 3.2808 # meters
+    segment_length = channel_length * 5280. / nodes
+
+    k_coef = 2. # days
+    x_coef = 0.1 # dimensionless
+
+    time_days, qinflow, qoutflow = get_ponce_data()
+
+    disl = flopy.mf6.ModflowSnfdisl(
+        snf, 
+        nodes=nodes, 
+        nvert=nvert,
+        segment_length=segment_length,
+        tosegment=[(-1,)],   # (-1,) gives 0 in one-based, which means outflow cell
+        idomain=1,
+        vertices=vertices, 
+        cell2d=cell2d,
+    )
+
+
+    # note: for specifying reach number, use fortran indexing!
+    fname = f"{name}.mmr.obs.csv"
+    mmr_obs = {
+        fname: [
+            ("OUTFLOW1", "EXT-OUTFLOW", 1),
+        ],
+        "digits": 10,
+    }
+
+    mmr = flopy.mf6.ModflowSnfmmr(
+        snf, 
+        print_flows=True,
+        observations=mmr_obs,
+        iseg_order=list(range(nodes)),
+        qoutflow0=qinflow[0],
+        k_coef=k_coef, 
+        x_coef=x_coef
+    )
+
+    # output control
+    oc = flopy.mf6.ModflowSnfoc(
+        snf,
+        budget_filerecord=f"{name}.bud",
+        saverecord=[("BUDGET", "ALL"), ],
+        printrecord=[("BUDGET", "ALL"), ],
+    )
+
+    inflow = qinflow[1:]
+    flw_spd = {ispd: [[0, inflow[ispd]]] for ispd in range(nper)}
+    flw = flopy.mf6.ModflowSnfflw(
+        snf,
+        print_input=True,
+        print_flows=True,
+        stress_period_data=flw_spd,
+    )
+
+    return sim, None
+
+
+def eval_model(sim):
+    print("evaluating model...")
+
+    # get back the ponce data for comparison
+    time_days, qinflow, qoutflow = get_ponce_data()
+
+    # read the observation output
+    name = ex[sim.idxsim]
+    fpth = os.path.join(sim.simpath, f"{name}.mmr.obs.csv")
+    obsvals = np.genfromtxt(fpth, names=True, delimiter=",")
+
+    # compare output with known result
+    time_days = time_days[1:]
+    qoutflow = -np.array(qoutflow[1:])
+    atol = 0.06
+    success = np.allclose(qoutflow, obsvals["OUTFLOW1"], atol=atol)
+    if not success:
+        for i, t in enumerate(time_days):
+            qa = qoutflow[i]
+            qs = obsvals["OUTFLOW1"][i]
+            d = qa - qs
+            if i == 0:
+                maxdiff = d
+            else:
+                if abs(d) > maxdiff:
+                    maxdiff = abs(d)
+            print(t, qa, qs, d)
+        print(f"maximum difference is {maxdiff}")
+    assert success
+
+    return
+
+
+@pytest.mark.parametrize(
+    "idx, name",
+    list(enumerate(ex)),
+)
+def test_mf6model(idx, name, function_tmpdir, targets):
+    ws = str(function_tmpdir)
+    test = TestFramework()
+    test.build(build_model, idx, ws)
+    test.run(
+        TestSimulation(
+            name=name, exe_dict=targets, exfunc=eval_model, idxsim=idx
+        ),
+        ws,
+    )