feat(gwf): add ghbg package (#2247)

* add ghba package
* add ghba dfn
* cleanup unused variables
* nodelist contains nodesuser nodes
* add structured netcdf grid write support
* add netcdf read and write support
* rebuild makefiles
* rebuild makefiles
* store dynamic netcdf input as timeseries
* add tests
* mark gridded array input packages as dev feature
* netcdf timeseries export array uses dis length units
* review updates
* add developmode marker to ghba tests
* apply maxbound to nodelist
* consolidate ghb package code
* update makefile
* rebuild makefiles
* rebuild makefiles
* restore makedefaults
* add mesh global attribute
* modflow_input attr should contain tag name
---------
Co-authored-by: mjreno <mjreno@IGSAAA071L01144.gs.doi.net>

Author: mjreno

autotest/test_gwf_disv_uzf.py

@@ -14,7 +14,7 @@
 import numpy as np
 import pytest
 from flopy.utils.gridutil import get_disv_kwargs
-from framework import TestFramework
+from framework import DNODATA, TestFramework
 
 cases = ["disv_with_uzf"]
 nlay = 5
@@ -109,20 +109,21 @@
     uzf_spd.update({t: spd})
 
 
-# Work up the GHB boundary
+# Work up the GHB / GHBG boundary
 ghb_ids = [(ncol - 1) + i * ncol for i in range(nrow)]
 ghb_spd = []
+abhead = np.full((nlay, ncpl), DNODATA, dtype=float)
+acond = np.full((nlay, ncpl), DNODATA, dtype=float)
 cond = 1e4
 for k in np.arange(3, 5, 1):
     for i in ghb_ids:
         ghb_spd.append([(k, i), 14.0, cond])
+        abhead[k, i] = 14.0
+        acond[k, i] = cond
 
 
-def build_models(idx, test):
-    name = cases[idx]
-
+def get_model(ws, name, array_input=False):
     # build MODFLOW 6 files
-    ws = test.workspace
     sim = flopy.mf6.MFSimulation(
         sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws
     )
@@ -171,7 +172,12 @@ def build_models(idx, test):
     sto = flopy.mf6.ModflowGwfsto(gwf, iconvert=1, ss=1e-5, sy=0.2, transient=True)
 
     # general-head boundary
-    ghb = flopy.mf6.ModflowGwfghb(gwf, print_flows=True, stress_period_data=ghb_spd)
+    if array_input:
+        ghb = flopy.mf6.ModflowGwfghbg(
+            gwf, print_flows=True, maxbound=20, bhead=abhead, cond=acond
+        )
+    else:
+        ghb = flopy.mf6.ModflowGwfghb(gwf, print_flows=True, stress_period_data=ghb_spd)
 
     # unsaturated-zone flow
     etobs = []
@@ -220,18 +226,18 @@ def build_models(idx, test):
     obs_dict = {f"{name}.obs.csv": obs_lst}
     obs = flopy.mf6.ModflowUtlobs(gwf, pname="head_obs", digits=20, continuous=obs_dict)
 
-    return sim, None
+    return sim
 
 
-def check_output(idx, test):
+def check_output(ws, name):
     # Next, get the binary printed heads
-    fpth = os.path.join(test.workspace, test.name + ".hds")
+    fpth = os.path.join(ws, name + ".hds")
     hobj = flopy.utils.HeadFile(fpth, precision="double")
     hds = hobj.get_alldata()
     hds = hds.reshape((np.sum(nstp), 5, 10, 10))
 
     # Get the MF6 cell-by-cell fluxes
-    bpth = os.path.join(test.workspace, test.name + ".cbc")
+    bpth = os.path.join(ws, name + ".cbc")
     bobj = flopy.utils.CellBudgetFile(bpth, precision="double")
     bobj.get_unique_record_names()
     # '          STO-SS'
@@ -249,7 +255,7 @@ def check_output(idx, test):
     gwet = gwetv.reshape((np.sum(nstp), 5, 10, 10))
 
     # Also retrieve the binary UZET output
-    uzpth = os.path.join(test.workspace, test.name + ".uzf.bud")
+    uzpth = os.path.join(ws, name + ".uzf.bud")
     uzobj = flopy.utils.CellBudgetFile(uzpth, precision="double")
     uzobj.get_unique_record_names()
     #  b'    FLOW-JA-FACE',
@@ -353,14 +359,41 @@ def check_output(idx, test):
     print("Finished running checks")
 
 
+def build_models(idx, test):
+    # build MODFLOW 6 files
+    ws = test.workspace
+    name = cases[idx]
+    sim = get_model(ws, name)
+
+    # build comparison array_input model
+    ws = os.path.join(test.workspace, "mf6")
+    mc = get_model(ws, name, array_input=True)
+
+    return sim, mc
+
+
+def check_outputs(idx, test):
+    name = cases[idx]
+
+    # check output MODFLOW 6 files
+    ws = test.workspace
+    check_output(ws, name)
+
+    # check output comparison array_input model
+    ws = os.path.join(test.workspace, "mf6")
+    check_output(ws, name)
+
+
 @pytest.mark.slow
+@pytest.mark.developmode
 @pytest.mark.parametrize("idx, name", enumerate(cases))
 def test_mf6model(idx, name, function_tmpdir, targets):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         targets=targets,
         build=lambda t: build_models(idx, t),
-        check=lambda t: check_output(idx, t),
+        check=lambda t: check_outputs(idx, t),
+        compare="mf6",
     )
     test.run()

autotest/test_gwf_sfr_inactive02.py

@@ -1,20 +1,19 @@
 # Test evap in SFR reaches (no interaction with gwf)
 
+import os
 
 import flopy
 import numpy as np
 import pytest
-from framework import TestFramework
+from framework import DNODATA, TestFramework
 
 HDRY, HNOFLO = -1e30, 1e30
 
 cases = ["sfr-inactive02"]
 
 
-def build_models(idx, test):
+def get_model(ws, name, array_input=False):
     # Base simulation and model name and workspace
-    ws = test.workspace
-    name = cases[idx]
 
     length_units = "m"
     time_units = "sec"
@@ -66,7 +65,15 @@ def build_models(idx, test):
         icelltype=1,  # >0 means saturated thickness varies with computed head
     )
     flopy.mf6.ModflowGwfic(gwf, strt=1.0)
-    flopy.mf6.ModflowGwfghb(gwf, stress_period_data=[((0, 0, 0), 1.0, 1e6)])
+    if array_input:
+        # if False:
+        bhead = np.full(nlay * nrow * ncol, DNODATA, dtype=float)
+        cond = np.full(nlay * nrow * ncol, DNODATA, dtype=float)
+        bhead[0] = 1.0
+        cond[0] = 1e6
+        flopy.mf6.ModflowGwfghbg(gwf, maxbound=1, bhead=bhead, cond=cond)
+    else:
+        flopy.mf6.ModflowGwfghb(gwf, stress_period_data=[((0, 0, 0), 1.0, 1e6)])
 
     # sfr data
     nreaches = 4
@@ -116,7 +123,7 @@ def build_models(idx, test):
         print_stage=True,
         print_flows=True,
         print_input=True,
-        stage_filerecord=f"{name}.sfr.hds",
+        stage_filerecord=f"{name}.sfr.stg",
         budget_filerecord=f"{name}.sfr.cbc",
         length_conversion=1.0,
         time_conversion=1.0,
@@ -150,11 +157,11 @@ def build_models(idx, test):
         saverecord=[("head", "all"), ("budget", "all")],
     )
 
-    return sim, None
+    return sim
 
 
-def check_output(idx, test):
-    sim = flopy.mf6.MFSimulation.load(sim_ws=test.workspace)
+def check_output(ws, name):
+    sim = flopy.mf6.MFSimulation.load(sim_ws=ws)
     gwf = sim.get_model()
     sfr = gwf.get_package("SFR-1")
     stage = sfr.output.stage().get_alldata().squeeze()
@@ -214,13 +221,40 @@ def check_output(idx, test):
     )
 
 
+def build_models(idx, test):
+    # build MODFLOW 6 files
+    ws = test.workspace
+    name = cases[idx]
+    sim = get_model(ws, name)
+
+    # build comparison array_input model
+    ws = os.path.join(test.workspace, "mf6")
+    mc = get_model(ws, name, array_input=True)
+
+    return sim, mc
+
+
+def check_outputs(idx, test):
+    name = cases[idx]
+
+    # check output MODFLOW 6 files
+    ws = test.workspace
+    check_output(ws, name)
+
+    # check output comparison array_input model
+    ws = os.path.join(test.workspace, "mf6")
+    check_output(ws, name)
+
+
+@pytest.mark.developmode
 @pytest.mark.parametrize("idx, name", enumerate(cases))
 def test_mf6model(idx, name, function_tmpdir, targets):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         targets=targets,
         build=lambda t: build_models(idx, t),
-        check=lambda t: check_output(idx, t),
+        check=lambda t: check_outputs(idx, t),
+        compare="mf6",
     )
     test.run()

autotest/test_gwf_uzf01.py

@@ -8,17 +8,15 @@
 import flopy
 import numpy as np
 import pytest
-from framework import TestFramework
+from framework import DNODATA, TestFramework
 
 cases = ["gwf_uzf01a"]
 nlay, nrow, ncol = 100, 1, 1
 
 crs = "EPSG:26916"
 
 
-def build_models(idx, test):
-    name = cases[idx]
-
+def get_model(ws, name, array_input=False):
     perlen = [500.0]
     nper = len(perlen)
     nstp = [10]
@@ -39,7 +37,6 @@ def build_models(idx, test):
         tdis_rc.append((perlen[i], nstp[i], tsmult[i]))
 
     # build MODFLOW 6 files
-    ws = test.workspace
     sim = flopy.mf6.MFSimulation(
         sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws
     )
@@ -104,16 +101,40 @@ def build_models(idx, test):
         transient={0: True},
     )
 
-    # ghb
-    ghbspdict = {
-        0: [[(nlay - 1, 0, 0), 1.5, 1.0]],
-    }
-    ghb = flopy.mf6.ModflowGwfghb(
-        gwf,
+    # ghb / ghbg
+    if array_input:
+        ghb_obs = {f"{name}.ghb.obs.csv": [("100_1_1", "GHB", (99, 0, 0))]}
+        bhead = np.full(nlay * nrow * ncol, DNODATA, dtype=float)
+        cond = np.full(nlay * nrow * ncol, DNODATA, dtype=float)
+        bhead[nlay - 1] = 1.5
+        cond[nlay - 1] = 1.0
+        ghb = flopy.mf6.ModflowGwfghbg(
+            gwf,
+            print_input=True,
+            print_flows=True,
+            maxbound=1,
+            bhead=bhead,
+            cond=cond,
+            save_flows=False,
+        )
+    else:
+        ghb_obs = {f"{name}.ghb.obs.csv": [("100_1_1", "GHB", (99, 0, 0))]}
+        ghbspdict = {
+            0: [[(nlay - 1, 0, 0), 1.5, 1.0]],
+        }
+        ghb = flopy.mf6.ModflowGwfghb(
+            gwf,
+            print_input=True,
+            print_flows=True,
+            stress_period_data=ghbspdict,
+            save_flows=False,
+        )
+
+    ghb.obs.initialize(
+        filename=f"{name}.ghb.obs",
+        digits=20,
         print_input=True,
-        print_flows=True,
-        stress_period_data=ghbspdict,
-        save_flows=False,
+        continuous=ghb_obs,
     )
 
     # note: for specifying lake number, use fortran indexing!
@@ -174,13 +195,10 @@ def build_models(idx, test):
     obs_dict = {f"{name}.obs.csv": obs_lst}
     obs = flopy.mf6.ModflowUtlobs(gwf, pname="head_obs", digits=20, continuous=obs_dict)
 
-    return sim, None
+    return sim
 
 
-def check_output(idx, test):
-    name = test.name
-    ws = test.workspace
-
+def check_output(ws, name):
     # check binary grid file
     fname = os.path.join(ws, name + ".dis.grb")
     grbobj = flopy.mf6.utils.MfGrdFile(fname)
@@ -228,13 +246,40 @@ def check_output(idx, test):
         )
 
 
+def build_models(idx, test):
+    # build MODFLOW 6 files
+    ws = test.workspace
+    name = cases[idx]
+    sim = get_model(ws, name)
+
+    # build comparison array_input model
+    ws = os.path.join(test.workspace, "mf6")
+    mc = get_model(ws, name, array_input=True)
+
+    return sim, mc
+
+
+def check_outputs(idx, test):
+    name = cases[idx]
+
+    # check output MODFLOW 6 files
+    ws = test.workspace
+    check_output(ws, name)
+
+    # check output comparison array_input model
+    ws = os.path.join(test.workspace, "mf6")
+    check_output(ws, name)
+
+
+@pytest.mark.developmode
 @pytest.mark.parametrize("idx, name", enumerate(cases))
 def test_mf6model(idx, name, function_tmpdir, targets):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         build=lambda t: build_models(idx, t),
-        check=lambda t: check_output(idx, t),
+        check=lambda t: check_outputs(idx, t),
         targets=targets,
+        compare="mf6",
     )
     test.run()