update swi with new tva subpackage for reading time variable saltwate…

…r heads

autotest/test_gwf_swi01.py

@@ -173,20 +173,21 @@ def make_head_plot(sim, idx, title):
     plt.close("all")
 
 
+def plot_output(idx, test):
+    title = {
+        0: "Case 1a -- Confined",
+        1: "Case 1b -- Unconfined",
+        2: "Case 1c -- Unconfined Newton",
+    }
+    sim = test.sims[0]
+    make_cross_section_plot(sim, idx, title[idx])
+    make_head_plot(sim, idx, title[idx])
+
+
 def check_output(idx, test):
     # get the flopy sim object
     sim = test.sims[0]
 
-    makeplot = True
-    if makeplot:
-        title = {
-            0: "Case 1a -- Confined",
-            1: "Case 1b -- Unconfined",
-            2: "Case 1c -- Unconfined Newton",
-        }
-        make_cross_section_plot(sim, idx, title[idx])
-        make_head_plot(sim, idx, title[idx])
-
     xanalytical, hanalytical, hmodel = get_heads_for_comparison(sim, idx)
     diff = hmodel - hanalytical
     print(f"\nSimulation idx {idx}")
@@ -209,12 +210,13 @@ def check_output(idx, test):
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets):
+def test_mf6model(idx, name, function_tmpdir, targets, plot):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         build=lambda t: build_models(idx, t),
         check=lambda t: check_output(idx, t),
+        plot=lambda t: plot_output(idx, t) if plot else None,
         targets=targets,
     )
     test.run()

autotest/test_gwf_swi02.py

@@ -111,10 +111,11 @@ def build_models(idx, test):
     return sim, None
 
 
-def make_plot(sim, idx):
+def plot_output(idx, test):
     import matplotlib.pyplot as plt
 
-    ws = pl.Path(sim.sim_path)
+    ws = test.workspace
+    sim = test.sims[0]
     gwf = sim.gwf[0]
     x = gwf.modelgrid.xcellcenters.flatten()
     fpth = pl.Path(ws) / f"{gwf.name}.zta"
@@ -138,10 +139,6 @@ def check_output(idx, test):
     # get the flopy sim object
     sim = test.sims[0]
 
-    makeplot = True
-    if makeplot:
-        make_plot(sim, idx)
-
     ws = pl.Path(sim.sim_path)
     gwf = sim.gwf[0]
     x = gwf.modelgrid.xcellcenters.flatten()
@@ -158,12 +155,13 @@ def check_output(idx, test):
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets):
+def test_mf6model(idx, name, function_tmpdir, targets, plot):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         build=lambda t: build_models(idx, t),
         check=lambda t: check_output(idx, t),
+        plot=lambda t: plot_output(idx, t) if plot else None,
         targets=targets,
     )
     test.run()

autotest/test_gwf_swi03.py

@@ -131,10 +131,11 @@ def build_models(idx, test):
     return sim, None
 
 
-def make_plot(sim, idx):
+def plot_output(idx, test):
     import matplotlib.pyplot as plt
 
-    ws = pl.Path(sim.sim_path)
+    ws = test.workspace
+    sim = test.sims[0]
     gwf = sim.gwf[0]
     x = gwf.modelgrid.xcellcenters.flatten()
     fpth = pl.Path(ws) / f"{gwf.name}.zta"
@@ -158,10 +159,6 @@ def check_output(idx, test):
     # get the flopy sim object
     sim = test.sims[0]
 
-    makeplot = True
-    if makeplot:
-        make_plot(sim, idx)
-
     ws = pl.Path(sim.sim_path)
     gwf = sim.gwf[0]
     x = gwf.modelgrid.xcellcenters.flatten()
@@ -178,12 +175,14 @@ def check_output(idx, test):
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets):
+def test_mf6model(idx, name, function_tmpdir, targets, plot):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         build=lambda t: build_models(idx, t),
         check=lambda t: check_output(idx, t),
+        plot=lambda t: plot_output(idx, t) if plot else None,
         targets=targets,
     )
     test.run()
+

autotest/test_gwf_swi04.py

@@ -115,12 +115,12 @@ def build_models(idx, test):
     return sim, None
 
 
-def make_cross_section_plot(test, sim, idx, title):
+def plot_output(idx, test):
     import matplotlib.pyplot as plt
 
-    ws = pl.Path(sim.sim_path)
+    ws = test.workspace
+    sim = test.sims[0]
     gwf = sim.gwf[0]
-    ws = sim.sim_path
     head = gwf.output.head().get_data().flatten()
     zeta = gwf.swi.output.zeta().get_data().flatten()
     pxs = flopy.plot.PlotCrossSection(gwf, line={"row": 0})
@@ -135,6 +135,7 @@ def make_cross_section_plot(test, sim, idx, title):
     pxs.plot_fill_between(zeta, head=head, colors=colors, ax=ax, edgecolors="none")
 
     pxs.plot_grid()
+    title = f"saltwater head = {saltwater_head[idx]}"
     ax.set_title(title)
     ax.set_ylim(botm_aquifer, top_island)
     plt.savefig(ws / "zeta.png")
@@ -145,19 +146,15 @@ def check_output(idx, test):
     # get the flopy sim object
     sim = test.sims[0]
 
-    makeplot = False
-    if makeplot:
-        title = f"saltwater head = {saltwater_head[idx]}"
-        make_cross_section_plot(test, sim, idx, title)
-
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets):
+def test_mf6model(idx, name, function_tmpdir, targets, plot):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         build=lambda t: build_models(idx, t),
         check=lambda t: check_output(idx, t),
+        plot=lambda t: plot_output(idx, t) if plot else None,
         targets=targets,
     )
     test.run()

autotest/test_gwf_swi05.py

@@ -104,12 +104,12 @@ def build_models(idx, test):
     return sim, None
 
 
-def make_cross_section_plot(test, sim, idx):
+def plot_output(idx, test):
     import matplotlib.pyplot as plt
 
-    ws = pl.Path(sim.sim_path)
+    ws = test.workspace
+    sim = test.sims[0]
     gwf = sim.gwf[0]
-    ws = sim.sim_path
     head = gwf.output.head().get_data().flatten()
     zeta = gwf.swi.output.zeta().get_data().flatten()
     volume_fresh = (head - zeta) * delr
@@ -132,18 +132,15 @@ def check_output(idx, test):
     # get the flopy sim object
     sim = test.sims[0]
 
-    makeplot = False
-    if makeplot:
-        make_cross_section_plot(test, sim, idx)
-
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets):
+def test_mf6model(idx, name, function_tmpdir, targets, plot):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         build=lambda t: build_models(idx, t),
         check=lambda t: check_output(idx, t),
+        plot=lambda t: plot_output(idx, t) if plot else None,
         targets=targets,
     )
     test.run()

autotest/test_gwf_swi06.py

@@ -0,0 +1,163 @@
+"""
+Test of the seawater intrusion (SWI) package for a sloping offshore ghb
+boundary.
+
+"""
+
+import pathlib as pl
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+
+cases = [
+    "swi06a",
+]
+
+Lx = 10000. # x position of right edge of model domain
+delr = 100.
+delc = 1.
+nlay = 1
+nrow = 1
+ncol = int(Lx / delr)
+top_island = 50.
+botm_aquifer = -50.
+recharge = 0.0001
+hydraulic_conductivity = 100.0
+ghb_cond_fact = 1.0
+perlen = 10 * 365.
+nper = 10
+perioddata = nper * [(perlen, 1, 1.0)]
+sealevel_start = 0.
+sealevel_stop = 40.
+
+
+def build_models(idx, test):
+    ws = test.workspace
+    name = "mymodel"
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name,
+        sim_ws=ws,
+        exe_name="mf6",
+        memory_print_option="all",
+        print_input=True,
+    )
+    tdis = flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=perioddata)
+    ims = flopy.mf6.ModflowIms(
+        sim,
+        print_option="summary",
+        no_ptcrecord=True,
+        outer_maximum=500,
+        linear_acceleration="bicgstab",
+    )
+    gwf = flopy.mf6.ModflowGwf(
+        sim,
+        modelname=name,
+        save_flows=True,
+        newtonoptions=True,
+    )
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top_island,
+        botm=botm_aquifer,
+    )
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=0)
+    npf = flopy.mf6.ModflowGwfnpf(
+        gwf,
+        save_specific_discharge=True,
+        # alternative_cell_averaging="amt-hmk",  # "harmonic",
+        icelltype=1,
+        k=hydraulic_conductivity,
+    )
+    zeta_file = name + ".zta"
+    swi = flopy.mf6.ModflowGwfswi(
+        gwf, 
+        zeta_filerecord=zeta_file,
+        tva6_filerecord=f"{name}.swi.tva",
+        zetastrt=-10,
+    )
+
+    # initialize the tva subpackage with the saltwater head
+    sl = np.linspace(sealevel_start, sealevel_stop, nper)
+    aux = {}
+    for kper in range(nper):
+        # this list has 2 entries, one for hsalt and one for myauxvar
+        aux[kper] = [sl[kper], np.ones((nlay, nrow, ncol))]
+    swi.tva.initialize(auxiliary=["hsalt", "myauxvar"], aux=aux)
+
+    cghb = hydraulic_conductivity * ghb_cond_fact * delr * delc / 1
+    ghbspd = {}
+    for kper in range(nper):
+        ghb_list = []
+        for j in [0, ncol - 1]:
+            freshwater_head = sl[kper]
+            ghb_list.append([0, 0, j, freshwater_head, cghb])
+        ghbspd[kper] = ghb_list
+
+    ghb = flopy.mf6.ModflowGwfghb(
+        gwf,
+        stress_period_data=ghbspd,
+    )
+    rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=recharge)
+    budget_file = name + ".bud"
+    head_file = name + ".hds"
+    oc = flopy.mf6.ModflowGwfoc(
+        gwf,
+        budget_filerecord=budget_file,
+        head_filerecord=head_file,
+        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+        printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+
+    return sim, None
+
+
+def plot_output(idx, test):
+    import matplotlib.pyplot as plt
+
+    sim = test.sims[0]
+    ws = test.workspace
+    gwf = sim.gwf[0]
+    head = gwf.output.head().get_data().flatten()
+    zeta = gwf.swi.output.zeta().get_data().flatten()
+    volume_fresh = (head - zeta) * delr
+    volume_fresh = volume_fresh.sum()
+    pxs = flopy.plot.PlotCrossSection(gwf, line={"row": 0})
+    ax = pxs.ax
+
+    zobj = gwf.swi.output.zeta()
+    times = zobj.times
+    for t in times:
+        zeta = zobj.get_data(totim=t)
+        pxs.plot_surface(zeta)
+
+    pxs.plot_grid()
+    title = f"Rising Sea Level"
+    ax.set_title(title)
+    ax.set_ylim(botm_aquifer, top_island)
+    plt.savefig(ws / "zeta.png")
+    plt.close("all")
+
+
+def check_output(idx, test):
+    # get the flopy sim object
+    sim = test.sims[0]
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets, plot):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        plot=lambda t: plot_output(idx, t) if plot else None,
+        targets=targets,
+    )
+    test.run()