fix(BudgetTerm.f90): ensure that the size of an array is >0 before deallocating

autotest/test_gwf_mvr01.py

@@ -374,48 +374,29 @@ def check_output(idx, test):
     times = bobj.get_times()
     records = bobj.get_data(totim=times[-1])
     adt = [("node", "<i4"), ("node2", "<i4"), ("q", "<f8")]
-    assert len(records) == 25
-    assert records[0].shape == (0,)
+    assert len(records) == 5
 
-    assert records[1].shape == (2,)
+    assert records[0].shape == (2,)
     a = np.array([(1, 2, -0.0), (2, 2, -0.0)], dtype=adt)
-    assert np.array_equal(records[1], a)
+    assert np.array_equal(records[0], a)
 
-    assert records[2].shape == (2,)
+    assert records[1].shape == (2,)
     a = np.array([(1, 1, -0.0), (2, 1, -0.0)], dtype=adt)
-    assert np.array_equal(records[2], a)
+    assert np.array_equal(records[1], a)
 
-    assert records[3].shape == (2,)
+    assert records[2].shape == (2,)
     a = np.array([(1, 3, -0.00545875), (2, 3, -0.00468419)], dtype=adt)
-    assert np.allclose(records[3]["node"], a["node"])
-    assert np.allclose(records[3]["node2"], a["node2"])
-    assert np.allclose(records[3]["q"], a["q"], atol=0.001), "{}\n{}".format(
-        records[3]["q"], a["q"]
+    assert np.allclose(records[2]["node"], a["node"])
+    assert np.allclose(records[2]["node2"], a["node2"])
+    assert np.allclose(records[2]["q"], a["q"], atol=0.001), "{}\n{}".format(
+        records[2]["q"], a["q"]
     )
 
-    assert records[4].shape == (0,)
-    assert records[5].shape == (0,)
-    assert records[6].shape == (0,)
-    assert records[7].shape == (0,)
-    assert records[8].shape == (3,)
+    assert records[3].shape == (3,)
     a = np.array([(1, 1, -0.0), (1, 2, -0.0005), (1, 3, -0.0)], dtype=adt)
-    assert np.array_equal(records[8], a)
+    assert np.array_equal(records[3], a)
 
-    assert records[9].shape == (0,)
-    assert records[10].shape == (0,)
-    assert records[11].shape == (0,)
-    assert records[12].shape == (0,)
-    assert records[13].shape == (0,)
-    assert records[14].shape == (0,)
-    assert records[15].shape == (0,)
-    assert records[16].shape == (0,)
-    assert records[17].shape == (0,)
-    assert records[18].shape == (0,)
-    assert records[19].shape == (0,)
-    assert records[20].shape == (0,)
-    assert records[21].shape == (0,)
-    assert records[22].shape == (0,)
-    assert records[23].shape == (9,)
+    assert records[4].shape == (9,)
     a = np.array(
         [
             (1, 1, -1.0e-04),
@@ -430,14 +411,12 @@ def check_output(idx, test):
         ],
         dtype=adt,
     )
-    assert np.allclose(records[23]["node"], a["node"])
-    assert np.allclose(records[23]["node2"], a["node2"])
-    assert np.allclose(records[23]["q"], a["q"], atol=0.001), "{}\n{}".format(
-        records[23]["q"], a["q"]
+    assert np.allclose(records[4]["node"], a["node"])
+    assert np.allclose(records[4]["node2"], a["node2"])
+    assert np.allclose(records[4]["q"], a["q"], atol=0.001), "{}\n{}".format(
+        records[4]["q"], a["q"]
     )
 
-    assert records[24].shape == (0,)
-
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
 def test_mf6model(idx, name, function_tmpdir, targets):

autotest/test_gwf_sfrsft_gwdischrg.py

@@ -0,0 +1,256 @@
+# Test groundwater discharge to a stream and then go on to test
+# that a transport model with a single reach works.
+
+import pathlib as pl
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+
+cases = ["sfr-gwfout", "sfr-gwf-trnsprt"]
+
+length_units = "m"
+time_units = "sec"
+
+nrow = 1
+ncol = 1
+nlay = 1
+delr = delc = 1.0
+
+nouter, ninner = 100, 300
+dvclose, rclose, relax = 1e-10, 1e-10, 1.0
+
+
+def add_gwt_model(sim, gwtname):
+    gwt = flopy.mf6.ModflowGwt(sim, modelname=gwtname, model_nam_file=f"{gwtname}.nam")
+    gwt.name_file.save_flows = True
+
+    imsgwt = flopy.mf6.ModflowIms(
+        sim,
+        print_option="SUMMARY",
+        outer_dvclose=dvclose,
+        outer_maximum=nouter,
+        under_relaxation="NONE",
+        inner_maximum=ninner,
+        inner_dvclose=dvclose,
+        rcloserecord=rclose,
+        linear_acceleration="BICGSTAB",
+        scaling_method="NONE",
+        reordering_method="NONE",
+        relaxation_factor=relax,
+        filename=f"{gwtname}.ims",
+    )
+    sim.register_ims_package(imsgwt, [gwt.name])
+
+    flopy.mf6.ModflowGwtdis(
+        gwt,
+        length_units=length_units,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=0.0,
+        botm=-100.0,
+    )
+
+    flopy.mf6.ModflowGwtic(gwt, strt=1.0)
+    flopy.mf6.ModflowGwtmst(gwt, porosity=0.2, filename=f"{gwtname}.mst")
+
+    sourcerecarray = [("GHB", "AUX", "CONCENTRATION")]
+    flopy.mf6.ModflowGwtssm(
+        gwt,
+        sources=sourcerecarray,
+        pname="SSM",
+        filename=f"{gwtname}.ssm",
+    )
+
+    # Instantiate Streamflow Transport package
+    sft_packagedata = []
+    t = (0, 1.0)
+    sft_packagedata.append(t)
+
+    sft_perioddata = []
+    sft_perioddata.append((0, "INFLOW", 0.0))
+    flwpckname = "SFR"
+
+    flopy.mf6.modflow.ModflowGwtsft(
+        gwt,
+        boundnames=False,
+        save_flows=True,
+        print_input=True,
+        print_flows=True,
+        print_concentration=True,
+        concentration_filerecord=gwtname + ".sft.bin",
+        budget_filerecord=gwtname + ".sft.bud",
+        packagedata=sft_packagedata,
+        reachperioddata=sft_perioddata,
+        flow_package_name=flwpckname,
+        pname="SFT",
+        filename=f"{gwtname}.sft",
+    )
+
+    flopy.mf6.ModflowGwtoc(
+        gwt,
+        budget_filerecord=f"{gwtname}.cbc",
+        concentration_filerecord=f"{gwtname}.ucn",
+        concentrationprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
+        saverecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")],
+        printrecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")],
+    )
+
+    return sim
+
+
+def build_models(idx, test):
+    # Base simulation and model name and workspace
+    ws = test.workspace
+    name = cases[idx]
+
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6"
+    )
+    flopy.mf6.ModflowTdis(sim, time_units=time_units)
+    flopy.mf6.ModflowIms(
+        sim,
+        inner_dvclose=1e-5,
+        inner_hclose=1e-6,
+    )
+
+    gwf = flopy.mf6.ModflowGwf(
+        sim,
+        modelname=name,
+    )
+    flopy.mf6.ModflowGwfdis(
+        gwf,
+        length_units=length_units,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=0.0,
+        botm=-100.0,
+    )
+    flopy.mf6.ModflowGwfnpf(
+        gwf,
+        icelltype=1,  # >0 means saturated thickness varies with computed head
+    )
+    flopy.mf6.ModflowGwfic(gwf, strt=1.0)
+    flopy.mf6.ModflowGwfghb(
+        gwf,
+        auxiliary="CONCENTRATION",
+        stress_period_data=[((0, 0, 0), 1.0, 1e6, 1.0)],
+        pname="GHB",
+    )
+
+    # sfr data
+    # <ifno> <cellid> <rlen> <rwid> <rgrd> <rtp> <rbth> <rhk> <man> <ncon> <ustrf> <ndv>
+    package_data = [(0, (0, 0, 0), delr, 1.0, 1e-3, 0.0, 1.0, 1.0, 0.001, 0, 0.0, 0)]
+    connection_data = [(0)]
+
+    sfr_obs = {
+        f"{name}.sfr.obs.csv": [
+            ("gwf", "sfr", (0,)),
+            ("outflow", "ext-outflow", (0,)),
+            ("depth", "depth", (0,)),
+        ],
+        "filename": name + ".sfr.obs",
+    }
+
+    flopy.mf6.ModflowGwfsfr(
+        gwf,
+        save_flows=True,
+        print_stage=True,
+        print_flows=True,
+        print_input=True,
+        length_conversion=1.0,
+        time_conversion=1.0,
+        nreaches=1,
+        packagedata=package_data,
+        connectiondata=connection_data,
+        observations=sfr_obs,
+        pname="SFR",
+    )
+
+    if idx > 0:
+        gwtname = "gwt-sft"
+        sim = add_gwt_model(sim, gwtname)
+
+        # Add the flow-transport exchanges
+        flopy.mf6.ModflowGwfgwt(
+            sim,
+            exgtype="GWF6-GWT6",
+            exgmnamea=name,
+            exgmnameb=gwtname,
+            pname="GWFGWT1",
+            filename=f"{gwtname}.gwfgwt1",
+        )
+
+    return sim, None
+
+
+def check_output(idx, test):
+    answer = np.array(
+        [
+            1.0,
+            -0.92094535738673577,
+            -0.92094535738673577,
+            0.79053721667952215e-1,
+        ]
+    )
+    obs_pth = pl.Path(f"{test.workspace}/{cases[idx]}.sfr.obs.csv")
+    sim_data = flopy.utils.Mf6Obs(obs_pth).get_data()
+    data_names = sim_data.dtype.names
+    for ct, name in enumerate(data_names):
+        assert np.allclose(
+            sim_data[name][0], answer[ct]
+        ), f"simulated sfr {name} results do not match answer"
+
+    if idx > 0:
+        sft_ans = np.array([[[[1.0]]]])
+
+        gwtname = "gwt-sft"
+        sft_obs_fl = pl.Path(f"{test.workspace}/{gwtname}.sft.bin")
+        try:
+            # load simulated concentration in SFT
+            cobj = flopy.utils.HeadFile(
+                sft_obs_fl,
+                text="CONCENTRATION",  # precision="double"
+            )
+            sim_conc_sft = cobj.get_alldata()
+        except:
+            assert False, f'could not load concentration data from "{sft_obs_fl}"'
+
+        gwt_sim_conc = pl.Path(f"{test.workspace}/{gwtname}.ucn")
+        try:
+            # load simulated concentration of groundwater
+            cobj = flopy.utils.HeadFile(
+                gwt_sim_conc,
+                text="CONCENTRATION",  # precision="double"
+            )
+            conc_gw = cobj.get_alldata()
+        except:
+            assert False, f'could not load temperature data from "{sft_obs_fl}"'
+
+        msg0 = "The simulation is not matching the established answer"
+        msg1 = (
+            "Groundwater discharge is the only source of flow in the "
+            "channel. Thus, the gw and sft water should have the same "
+            "concentration, but don't"
+        )
+        assert np.allclose(sft_ans, sim_conc_sft), msg0
+        assert np.allclose(conc_gw, sim_conc_sft), msg1
+
+
+@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),
+    )
+    test.run()