diff --git a/autotest/conftest.py b/autotest/conftest.py index 5f0073b18c2..27ff70c453d 100644 --- a/autotest/conftest.py +++ b/autotest/conftest.py @@ -1,8 +1,11 @@ import platform +import sys from pathlib import Path +from os import PathLike +from typing import Dict import pytest -from modflow_devtools.executables import Executables, build_default_exe_dict +from modflow_devtools.executables import Executables, get_suffixes pytest_plugins = ["modflow_devtools.fixtures"] project_root_path = Path(__file__).resolve().parent.parent @@ -45,7 +48,42 @@ def libmf6_path(bin_path) -> Path: @pytest.fixture(scope="session") def targets(bin_path) -> Executables: - return Executables(**build_default_exe_dict(bin_path)) + exe_ext, lib_ext = get_suffixes(sys.platform) + dl_bin = bin_path / "downloaded" + rb_bin = bin_path / "rebuilt" + tgts = dict() + + # downloaded executables + tgts["mf2000"] = dl_bin / f"mf2000{exe_ext}" + tgts["mf2005"] = dl_bin / f"mf2005dbl{exe_ext}" + tgts["mfnwt"] = dl_bin / f"mfnwtdbl{exe_ext}" + tgts["mfusg"] = dl_bin / f"mfusgdbl{exe_ext}" + tgts["mflgr"] = dl_bin / f"mflgrdbl{exe_ext}" + tgts["mf2005s"] = dl_bin / f"mf2005{exe_ext}" + tgts["mt3dms"] = dl_bin / f"mt3dms{exe_ext}" + tgts["crt"] = dl_bin / f"crt{exe_ext}" + tgts["gridgen"] = dl_bin / f"gridgen{exe_ext}" + tgts["mp6"] = dl_bin / f"mp6{exe_ext}" + tgts["mp7"] = dl_bin / f"mp7{exe_ext}" + tgts["swtv4"] = dl_bin / f"swtv4{exe_ext}" + tgts["sutra"] = dl_bin / f"sutra{exe_ext}" + tgts["triangle"] = dl_bin / f"triangle{exe_ext}" + tgts["vs2dt"] = dl_bin / f"vs2dt{exe_ext}" + tgts["zonbudusg"] = dl_bin / f"zonbudusg{exe_ext}" + + # binaries rebuilt from last release + tgts["mf6_regression"] = rb_bin / f"mf6{exe_ext}" + tgts["libmf6_regression"] = rb_bin / f"libmf6{lib_ext}" + tgts["mf5to6_regression"] = rb_bin / f"mf5to6{exe_ext}" + tgts["zbud6_regression"] = rb_bin / f"zbud6{exe_ext}" + + # local development binaries + tgts["mf6"] = bin_path / f"mf6{exe_ext}" + tgts["libmf6"] = bin_path / f"libmf6{lib_ext}" + tgts["mf5to6"] = bin_path / f"mf5to6{exe_ext}" + tgts["zbud6"] = bin_path / f"zbud6{exe_ext}" + + return Executables(**tgts) @pytest.fixture diff --git a/autotest/test_gwf.py b/autotest/test_gwf.py deleted file mode 100644 index 59e0d18f6ac..00000000000 --- a/autotest/test_gwf.py +++ /dev/null @@ -1,27 +0,0 @@ -from modflow_devtools.executables import Executables -from pytest_cases import parametrize_with_cases -from simulation import TestSimulation -from test_gwf_maw04 import GwfMaw04Cases -from test_gwf_maw_cases import GwfMawCases - - -@parametrize_with_cases("case", cases=[GwfMawCases, GwfMaw04Cases]) -def test_gwf_models(case, targets: Executables): - data, sim, cmp, exfunc = case - sim.write_simulation() - if cmp: - cmp.write_simulation() - - test = TestSimulation( - name=data.name, - exe_dict=targets, - exfunc=exfunc, - idxsim=0, # TODO: remove parameter from TestSimulation - mf6_regression=True, - require_failure=data.xfail, - make_comparison=data.compare, - ) - - test.set_model(sim.simulation_data.mfpath.get_sim_path(), testModel=False) - test.run() - test.compare() diff --git a/autotest/test_gwf_maw01.py b/autotest/test_gwf_maw01.py new file mode 100644 index 00000000000..7db88abe711 --- /dev/null +++ b/autotest/test_gwf_maw01.py @@ -0,0 +1,213 @@ +import os +from types import SimpleNamespace as Case + +import flopy +import numpy as np +import pytest + +from framework import TestFramework +from simulation import TestSimulation + +budtol = 1e-2 +bud_lst = ["GWF_IN", "GWF_OUT", "RATE_IN", "RATE_OUT"] + +well1 = Case( + observations={"maw_obs.csv": [("mh1", "head", 1)]}, + packagedata=[[0, 0.1, 50.0, 100.0, "THIEM", 1]], + connectiondata=[[0, 0, (0, 0, 1), 100.0, 50.0, 1.0, 0.1]], + perioddata=[[0, "rate", 0.0]], +) + +ex = ["maw01", "maw01nwt", "maw01nwtur"] +krylov = ["CG", "BICGSTAB", "BICGSTAB"] +newton = [None, "NEWTON", "NEWTON UNDER_RELAXATION"] +nlay = 1 +nrow = 1 +ncol = 3 +nper = 3 +delr = 300 +delc = 300 +perlen = 3 * [1] +nstp = 3 * [1] +tsmult = 3 * [1] +well = well1 +strt = 100 +hk = 1 +nouter = 100 +ninner = 300 +hclose = 1e-9 +rclose = 1e-3 +relaxation_factor = 1 +compare = False + + +def build_model(idx, ws, mf6): + name = ex[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, + version="mf6", + exe_name=mf6, + sim_ws=ws, + ) + + # create tdis package + tdis_rc = [(perlen[i], nstp[i], tsmult[i]) for i in range(nper)] + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # create gwf model + gwf = flopy.mf6.MFModel( + sim, + model_type="gwf6", + modelname=name, + model_nam_file=f"{name}.nam", + ) + gwf.name_file.newtonoptions = newton[idx] + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration=krylov[idx], + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relaxation_factor, + ) + sim.register_ims_package(ims, [gwf.name]) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=100.0, + botm=0.0, + idomain=1, + filename=f"{name}.dis", + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{name}.ic") + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=1, + k=hk, + k33=hk, + filename=f"{name}.npf", + ) + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=1, + ss=0.0, + sy=0.1, + steady_state={0: True}, + # transient={1: False}, + filename=f"{name}.sto", + ) + + # chd files + chdlist0 = [] + chdlist0.append([(0, 0, 0), 100.0]) + chdlist0.append([(0, 0, 2), 100.0]) + + chdlist1 = [] + chdlist1.append([(0, 0, 0), 25.0]) + chdlist1.append([(0, 0, 2), 25.0]) + + chdspdict = {0: chdlist0, 1: chdlist1, 2: chdlist0} + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdspdict, + save_flows=False, + filename=f"{name}.chd", + ) + + # wel files + # wel = flopy.mf6.ModflowGwfwel(gwf, print_input=True, print_flows=True, + # maxbound=len(ws), + # periodrecarray=wd6, + # save_flows=False) + # MAW + maw = flopy.mf6.ModflowGwfmaw( + gwf, + filename=f"{name}.maw", + print_input=True, + print_head=True, + print_flows=True, + save_flows=True, + observations=well.observations, + packagedata=well.packagedata, + connectiondata=well.connectiondata, + perioddata=well.perioddata, + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{name}.oc", + ) + + return sim, None + + +def eval_results(sim): + print("evaluating MAW heads...") + + # MODFLOW 6 maw results + fpth = os.path.join(sim.simpath, "maw_obs.csv") + tc = np.genfromtxt(fpth, names=True, delimiter=",") + + # create known results array + tc0 = np.array([100.0, 25.0, 100.0]) + + # calculate maximum absolute error + diff = tc["MH1"] - tc0 + diffmax = np.abs(diff).max() + dtol = 1e-9 + msg = f"maximum absolute maw head difference ({diffmax}) " + + if diffmax > dtol: + sim.success = False + msg += f"exceeds {dtol}" + assert diffmax < dtol, msg + else: + sim.success = True + print(" " + msg) + + +@pytest.mark.parametrize("idx, name", list(enumerate(ex))) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + sim, _ = build_model(idx, ws, targets.mf6) + sim.write_simulation() + sim.run_simulation() + test = TestFramework() + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + ), + ws, + ) diff --git a/autotest/test_gwf_maw02.py b/autotest/test_gwf_maw02.py new file mode 100644 index 00000000000..09eda8d61be --- /dev/null +++ b/autotest/test_gwf_maw02.py @@ -0,0 +1,312 @@ +import os +from types import SimpleNamespace as Case + +import flopy +import numpy as np +import pytest + +from framework import TestFramework +from simulation import TestSimulation + +budtol = 1e-2 +bud_lst = ["GWF_IN", "GWF_OUT", "RATE_IN", "RATE_OUT"] + +well2 = Case( + observations={"maw_obs.csv": [("mh1", "head", 1)]}, + packagedata=[ + [0, 0.1, 0.0, 100.0, "THIEM", 1], + [1, 0.1, 0.0, 100.0, "THIEM", 1], + ], + connectiondata=[ + [0, 0, (0, 0, 1), 100.0, 0.0, 1.0, 0.1], + [1, 0, (0, 0, 1), 100.0, 0.0, 1.0, 0.1], + ], + perioddata={ + 0: [ + [0, "rate", -20.0], + [0, "status", "inactive"], + [0, "rate_scaling", 1.0, 15.0], + [1, "rate", -30.0], + [1, "status", "inactive"], + [1, "rate_scaling", 5.0, 15.0], + ], + 1: [ + [0, "rate", -110.0], + [0, "status", "active"], + [1, "rate", -130.0], + [1, "status", "active"], + ], + 3: [[0, "status", "inactive"]], + 4: [[0, "status", "active"]], + }, +) + +ex = ["maw02"] +krylov = "CG" +nlay = 1 +nrow = 1 +ncol = 3 +nper = 5 +delr = 300 +delc = 300 +perlen = 5 * [1] +nstp = 5 * [1] +tsmult = 5 * [1] +well = well2 +strt = 100 +hk = 1 +nouter = 100 +ninner = 300 +hclose = 1e-9 +rclose = 1e-3 +relaxation_factor = 1 +compare = False + + +def build_model(idx, ws, mf6): + name = ex[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name=mf6, sim_ws=ws + ) + + # create tdis package + tdis_rc = [(perlen[i], nstp[i], tsmult[i]) for i in range(nper)] + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # create gwf model + gwf = flopy.mf6.MFModel( + sim, + model_type="gwf6", + modelname=name, + model_nam_file=f"{name}.nam", + ) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration=krylov, + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relaxation_factor, + ) + sim.register_ims_package(ims, [gwf.name]) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=100.0, + botm=0.0, + idomain=1, + filename=f"{name}.dis", + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{name}.ic") + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=1, + k=hk, + k33=hk, + filename=f"{name}.npf", + ) + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=1, + ss=0.0, + sy=0.1, + steady_state={0: True}, + # transient={1: False}, + filename=f"{name}.sto", + ) + + # chd files + chdlist0 = [] + chdlist0.append([(0, 0, 0), 100.0]) + chdlist0.append([(0, 0, 2), 100.0]) + + chdlist1 = [] + chdlist1.append([(0, 0, 0), 25.0]) + chdlist1.append([(0, 0, 2), 25.0]) + + chdspdict = {0: chdlist0, 1: chdlist1, 2: chdlist0} + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdspdict, + save_flows=False, + filename=f"{name}.chd", + ) + + # MAW + maw = flopy.mf6.ModflowGwfmaw( + gwf, + filename=f"{name}.maw", + budget_filerecord=f"{name}.maw.cbc", + print_input=True, + print_head=True, + print_flows=True, + save_flows=True, + observations=well.observations, + packagedata=well.packagedata, + connectiondata=well.connectiondata, + perioddata=well.perioddata, + pname="MAW-1", + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{name}.oc", + ) + + return sim, None + + +def eval_results(sim): + print("evaluating MAW budgets...") + + shape3d = (nlay, nrow, ncol) + size3d = nlay * nrow * ncol + + # get results from listing file + fpth = os.path.join(sim.simpath, f"{os.path.basename(sim.name)}.lst") + budl = flopy.utils.Mf6ListBudget( + fpth, budgetkey="MAW-1 BUDGET FOR ENTIRE MODEL AT END OF TIME STEP" + ) + names = list(bud_lst) + d0 = budl.get_budget(names=names)[0] + dtype = d0.dtype + nbud = d0.shape[0] + + # get results from cbc file + cbc_bud = ["GWF", "RATE"] + d = np.recarray(nbud, dtype=dtype) + for key in bud_lst: + d[key] = 0.0 + fpth = os.path.join(sim.simpath, f"{os.path.basename(sim.name)}.maw.cbc") + cobj = flopy.utils.CellBudgetFile(fpth, precision="double") + kk = cobj.get_kstpkper() + times = cobj.get_times() + cbc_vals = [] + for idx, (k, t) in enumerate(zip(kk, times)): + for text in cbc_bud: + qin = 0.0 + qout = 0.0 + v = cobj.get_data(kstpkper=k, text=text)[0] + if isinstance(v, np.recarray): + vt = np.zeros(size3d, dtype=float) + wq = [] + for jdx, node in enumerate(v["node"]): + vt[node - 1] += v["q"][jdx] + wq.append(v["q"][jdx]) + v = vt.reshape(shape3d) + if text == cbc_bud[-1]: + cbc_vals.append(wq) + for kk in range(v.shape[0]): + for ii in range(v.shape[1]): + for jj in range(v.shape[2]): + vv = v[kk, ii, jj] + if vv < 0.0: + qout -= vv + else: + qin += vv + d["totim"][idx] = t + d["time_step"][idx] = k[0] + d["stress_period"] = k[1] + key = f"{text}_IN" + d[key][idx] = qin + key = f"{text}_OUT" + d[key][idx] = qout + + maw_vals = [ + [0.000, 0.000], + [-106.11303563809453, -96.22598985147631], + [-110.000, -130.000], + [0.0, -130.000], + [-110.000, -130.000], + ] + + # evaluate if well rates in cbc file are equal to expected values + diffv = [] + for ovs, svs in zip(maw_vals, cbc_vals): + for ov, sv in zip(ovs, svs): + diffv.append(ov - sv) + diffv = np.abs(np.array(diffv)).max() + msg = f"\nmaximum absolute maw rate difference ({diffv})\n" + + # calculate difference between water budget items in the lst and cbc files + diff = np.zeros((nbud, len(bud_lst)), dtype=float) + for idx, key in enumerate(bud_lst): + diff[:, idx] = d0[key] - d[key] + diffmax = np.abs(diff).max() + msg += f"maximum absolute total-budget difference ({diffmax}) " + + # write summary + fpth = os.path.join( + sim.simpath, f"{os.path.basename(sim.name)}.bud.cmp.out" + ) + f = open(fpth, "w") + for i in range(diff.shape[0]): + if i == 0: + line = f"{'TIME':>10s}" + for idx, key in enumerate(bud_lst): + line += f"{key + '_LST':>25s}" + line += f"{key + '_CBC':>25s}" + line += f"{key + '_DIF':>25s}" + f.write(line + "\n") + line = f"{d['totim'][i]:10g}" + for idx, key in enumerate(bud_lst): + line += f"{d0[key][i]:25g}" + line += f"{d[key][i]:25g}" + line += f"{diff[i, idx]:25g}" + f.write(line + "\n") + f.close() + + if diffmax > budtol or diffv > budtol: + sim.success = False + msg += f"\n...exceeds {budtol}" + assert diffmax < budtol, msg + else: + sim.success = True + print(" " + msg) + + +@pytest.mark.parametrize("idx, name", list(enumerate(ex))) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + sim, _ = build_model(idx, ws, targets.mf6) + sim.write_simulation() + sim.run_simulation() + test = TestFramework() + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + ), + ws, + ) diff --git a/autotest/test_gwf_maw03.py b/autotest/test_gwf_maw03.py new file mode 100644 index 00000000000..4a4c6ae71a8 --- /dev/null +++ b/autotest/test_gwf_maw03.py @@ -0,0 +1,229 @@ +import os +from types import SimpleNamespace as Case + +import flopy +import numpy as np +import pytest + +from framework import TestFramework +from simulation import TestSimulation + +budtol = 1e-2 +bud_lst = ["GWF_IN", "GWF_OUT", "RATE_IN", "RATE_OUT"] + + +def well3(name): + perioddata = { + "maw03a": [ + (0, "rate", 2000.0), + ], + "maw03b": [(0, "rate", 2000.0), (0, "head_limit", 0.4)], + "maw03c": [(0, "rate", 2000.0), (0, "rate_scaling", 0.0, 1.0)], + } + wellbottom = -1000 + return Case( + observations={ + f"{name}.maw.obs.csv": [ + ("m1head", "head", (0,)), + ("m1rate", "rate", (0,)), + ] # is this index one-based? Not if in a tuple + }, + packagedata=[[0, 0.15, wellbottom, 0.0, "THIEM", 1]], + connectiondata=[[0, 0, (0, 50, 50), 0.0, wellbottom, 0.0, 0.0]], + perioddata=perioddata[name], + ) + + +ex = ["maw03a", "maw03b", "maw03c"] +krylov = "CG" +nlay = 1 +nrow = 101 +ncol = 101 +nper = 1 +delr = 142 +delc = 142 +perlen = [1000] +nstp = [50] +tsmult = [1.2] +strt = 0 +hk = 10 +nouter = 100 +ninner = 100 +hclose = 1e-6 +rclose = 1e-6 +relaxation_factor = 1 +compare = False + + +def build_model(idx, ws, mf6): + top = 0.0 + botm = [-1000.0] + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = ex[idx] + sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws, exe_name=mf6) + + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # create gwf model + gwf = flopy.mf6.MFModel( + sim, + model_type="gwf6", + modelname=name, + model_nam_file=f"{name}.nam", + ) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration=krylov, + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relaxation_factor, + ) + sim.register_ims_package(ims, [gwf.name]) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=1, + filename=f"{name}.dis", + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{name}.ic") + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=1, + k=hk, + k33=hk, + filename=f"{name}.npf", + ) + + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=0, + ss=1.0e-5, + sy=0.1, + steady_state={0: False}, + transient={0: True}, + filename=f"{name}.sto", + ) + + # MAW + well = well3(name) + maw = flopy.mf6.ModflowGwfmaw( + gwf, + filename=f"{name}.maw", + print_input=True, + print_head=True, + print_flows=True, + save_flows=True, + observations=well.observations, + packagedata=well.packagedata, + connectiondata=well.connectiondata, + perioddata=well.perioddata, + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + headprintrecord=[ + ("COLUMNS", ncol, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("HEAD", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{name}.oc", + ) + + # head observations + obs_data0 = [("head_well_cell", "HEAD", (0, 0, 0))] + obs_recarray = {f"{name}.obs.csv": obs_data0} + obs = flopy.mf6.ModflowUtlobs( + gwf, + pname="head_obs", + filename=f"{name}.obs", + digits=15, + print_input=True, + continuous=obs_recarray, + ) + + return sim, None + + +def eval_results(sim): + print("evaluating MAW heads...") + + # MODFLOW 6 maw results + test_name = sim.name + fpth = os.path.join(sim.simpath, f"{test_name}.maw.obs.csv") + tc = np.genfromtxt(fpth, names=True, delimiter=",") + + if test_name.endswith("a"): + # M1RATE should be 2000. + msg = "The injection rate should be 2000. for all times" + assert tc["M1RATE"].min() == tc["M1RATE"].max() == 2000, msg + + elif test_name.endswith("b"): + # M1RATE should have a minimum value less than 200 and + # M1HEAD should not exceed 0.400001 + msg = ( + "Injection rate should fall below 200 and the head should not" + "exceed 0.4" + ) + assert tc["M1RATE"].min() < 200.0, msg + assert tc["M1HEAD"].max() < 0.400001, msg + + elif test_name.endswith("c"): + # M1RATE should have a minimum value less than 800 + # M1HEAD should not exceed 1.0. + msg = ( + "Min injection rate should be less than 800 and well " + "head should not exceed 1.0" + ) + assert tc["M1RATE"].min() < 800.0 and tc["M1HEAD"].max() < 1.0, msg + + +@pytest.mark.parametrize("idx, name", list(enumerate(ex))) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + sim, _ = build_model(idx, ws, targets.mf6) + sim.write_simulation() + sim.run_simulation() + test = TestFramework() + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + ), + ws, + ) diff --git a/autotest/test_gwf_maw04.py b/autotest/test_gwf_maw04.py index e77656f0c68..462a3afa735 100644 --- a/autotest/test_gwf_maw04.py +++ b/autotest/test_gwf_maw04.py @@ -1,10 +1,12 @@ import os -from types import SimpleNamespace +from types import SimpleNamespace as Case import flopy import numpy as np -from modflow_devtools.case import Case -from pytest_cases import parametrize +import pytest + +from framework import TestFramework +from simulation import TestSimulation # temporal discretization nper = 2 @@ -92,7 +94,7 @@ def well4(label): [0, 1, (1, nhalf, nhalf), botm[0], botm[1], hks, sradius[label]], ] perioddata = {1: [[0, "RATE", wellq]]} - return SimpleNamespace( + return Case( print_input=True, no_well_storage=True, packagedata=packagedata, @@ -101,198 +103,188 @@ def well4(label): ) -case = Case( - name="maw_iss305", - nlay=nlay, - nrow=nrow, - ncol=ncol, - nper=nper, - delr=delr, - perlen=perlen, - nstp=nstp, - tsmult=tsmult, - steady=steady, - strt=0, - hk=10, - nouter=100, - ninner=100, - hclose=1e-9, - rclose=1e-6, - relax=1, - top=top, - botm=botm, - confined=confined, - ss=ss, - chd_spd=chd_spd, - chd5_spd=chd5_spd, - nhalf=nhalf, - radius=radius, - wellq=wellq, - compare=False, -) -cases = [case.copy_update(name=case.name + "a", well=well4("a"),)] + [ - case.copy_update(name=case.name + label, well=well4(label), xfail=True) - for label in [ - "b", - # "c", # todo: this one passes when it should fail - "d", - "e", - "f", - ] -] +# npf data +strt = 0 +hk = 10 +# solver +nouter = 100 +ninner = 100 +hclose = 1e-9 +rclose = 1e-6 +relax = 1 -class GwfMaw04Cases: - @parametrize(data=cases, ids=[c.name for c in cases]) - def case_4(self, function_tmpdir, targets, data): - name = data.name - ws = str(function_tmpdir) +# subproblems +subprobs = ["a", "b", "c", "d", "e", "f"] +ex = [f"maw_iss305{sp}" for sp in subprobs] +wells = [well4(sp) for sp in subprobs] +xfail = [False, True, True, True, True, True] - # build MODFLOW 6 files - sim = flopy.mf6.MFSimulation( - sim_name=name, version="mf6", exe_name=targets["mf6"], sim_ws=ws - ) - # create tdis package - tdis_rc = [] - for idx in range(data.nper): - tdis_rc.append( - (data.perlen[idx], data.nstp[idx], data.tsmult[idx]) - ) - tdis = flopy.mf6.ModflowTdis( - sim, time_units="DAYS", nper=data.nper, perioddata=tdis_rc - ) - # create iterative model solution - ims = flopy.mf6.ModflowIms( - sim, - inner_dvclose=data.hclose, - rcloserecord=data.rclose, - outer_dvclose=data.hclose, - ) +def build_model(idx, ws, mf6): + name = ex[idx] + well = wells[idx] + # build MODFLOW 6 files + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name=mf6, sim_ws=ws + ) + # create tdis package + tdis_rc = [] + for kper in range(nper): + tdis_rc.append((perlen[kper], nstp[kper], tsmult[kper])) + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) - # create gwf model - gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) + # create iterative model solution + ims = flopy.mf6.ModflowIms( + sim, + inner_dvclose=hclose, + rcloserecord=rclose, + outer_dvclose=hclose, + ) - # discretization - dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=data.nlay, - nrow=data.nrow, - ncol=data.ncol, - delr=data.delr, - delc=data.delr, - top=data.top, - botm=data.botm, - ) - # initial conditions - ic = flopy.mf6.ModflowGwfic(gwf, strt=data.strt) + # create gwf model + gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) - # node property flow - npf = flopy.mf6.ModflowGwfnpf( - gwf, save_flows=False, icelltype=data.confined, k=data.hk - ) - # storage - sto = flopy.mf6.ModflowGwfsto( - gwf, - save_flows=False, - iconvert=data.confined, - ss=data.ss, - steady_state={0: True}, - transient={1: True}, + # discretization + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delr, + top=top, + botm=botm, + ) + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt) + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, save_flows=False, icelltype=confined, k=hk + ) + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=False, + iconvert=confined, + ss=ss, + steady_state={0: True}, + transient={1: True}, + ) + # constant head + chd = flopy.mf6.ModflowGwfchd( + gwf, stress_period_data=chd_spd, save_flows=False + ) + # multi-aquifer well + maw = flopy.mf6.ModflowGwfmaw( + gwf, + print_input=well.print_input, + no_well_storage=well.no_well_storage, + packagedata=well.packagedata, + connectiondata=well.connectiondata, + perioddata=well.perioddata, + ) + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + # build MODFLOW-2005 files + if xfail[idx]: + mc = None + else: + cmppth = "mf2005" + ws = os.path.join(ws, cmppth) + mc = flopy.modflow.Modflow(name, model_ws=ws, version=cmppth) + dis = flopy.modflow.ModflowDis( + mc, + nlay=nlay, + nrow=nrow, + ncol=ncol, + nper=nper, + perlen=perlen, + nstp=nstp, + tsmult=tsmult, + steady=steady, + delr=delr, + delc=delr, + top=top, + botm=botm, ) - # constant head - chd = flopy.mf6.ModflowGwfchd( - gwf, stress_period_data=data.chd_spd, save_flows=False + bas = flopy.modflow.ModflowBas(mc, strt=strt) + lpf = flopy.modflow.ModflowLpf( + mc, + laytyp=confined, + hk=hk, + vka=hk, + ss=ss, + sy=0, ) - # multi-aquifer well - maw = flopy.mf6.ModflowGwfmaw( - gwf, - print_input=data.well.print_input, - no_well_storage=data.well.no_well_storage, - packagedata=data.well.packagedata, - connectiondata=data.well.connectiondata, - perioddata=data.well.perioddata, + chd = flopy.modflow.ModflowChd(mc, stress_period_data=chd5_spd) + # mnw2 + # empty mnw2 file to create recarrays + mnw2 = flopy.modflow.ModflowMnw2(mc) + node_data = mnw2.get_empty_node_data(2) + node_data["ztop"] = np.array([top, botm[0]]) + node_data["zbotm"] = np.array([botm[0], botm[1]]) + node_data["i"] = np.array([nhalf, nhalf]) + node_data["j"] = np.array([nhalf, nhalf]) + node_data["wellid"] = np.array(["well1", "well1"]) + node_data["losstype"] = np.array(["skin", "skin"]) + node_data["rw"] = np.array([radius, radius]) + node_data["rskin"] = np.array([sradius[name[-1]], sradius[name[-1]]]) + hks = hk * skin_mult[name[-1]] + node_data["kskin"] = np.array([hks, hks]) + dtype = [("wellid", np.unicode_, 20), ("qdes", " dtol: - config.success = False - msg += f"exceeds {dtol}" - assert diffmax < dtol, msg - else: - config.success = True - print(" " + msg) - - case2 = Case( - name="maw02", - krylov="CG", - nlay=1, - nrow=1, - ncol=3, - nper=5, - delr=300, - delc=300, - perlen=5 * [1], - nstp=5 * [1], - tsmult=5 * [1], - well=well2, - strt=100, - hk=1, - nouter=100, - ninner=300, - hclose=1e-9, - rclose=1e-3, - relaxation_factor=1, - compare=False, - ) - cases2 = [case2] - - @parametrize(data=cases2, ids=[c.name for c in cases2]) - def case_2(self, function_tmpdir, targets, data): - name = data.name - ws = str(function_tmpdir) - sim = flopy.mf6.MFSimulation( - sim_name=name, version="mf6", exe_name=targets["mf6"], sim_ws=ws - ) - - # create tdis package - tdis_rc = [ - (data.perlen[i], data.nstp[i], data.tsmult[i]) - for i in range(data.nper) - ] - tdis = flopy.mf6.ModflowTdis( - sim, time_units="DAYS", nper=data.nper, perioddata=tdis_rc - ) - - # create gwf model - gwf = flopy.mf6.MFModel( - sim, - model_type="gwf6", - modelname=name, - model_nam_file=f"{name}.nam", - ) - - # create iterative model solution and register the gwf model with it - ims = flopy.mf6.ModflowIms( - sim, - print_option="SUMMARY", - outer_dvclose=data.hclose, - outer_maximum=data.nouter, - under_relaxation="NONE", - inner_maximum=data.ninner, - inner_dvclose=data.hclose, - rcloserecord=data.rclose, - linear_acceleration=data.krylov, - scaling_method="NONE", - reordering_method="NONE", - relaxation_factor=data.relaxation_factor, - ) - sim.register_ims_package(ims, [gwf.name]) - - dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=data.nlay, - nrow=data.nrow, - ncol=data.ncol, - delr=data.delr, - delc=data.delc, - top=100.0, - botm=0.0, - idomain=1, - filename=f"{name}.dis", - ) - - # initial conditions - ic = flopy.mf6.ModflowGwfic(gwf, strt=data.strt, filename=f"{name}.ic") - - # node property flow - npf = flopy.mf6.ModflowGwfnpf( - gwf, - save_flows=True, - icelltype=1, - k=data.hk, - k33=data.hk, - filename=f"{name}.npf", - ) - # storage - sto = flopy.mf6.ModflowGwfsto( - gwf, - save_flows=True, - iconvert=1, - ss=0.0, - sy=0.1, - steady_state={0: True}, - # transient={1: False}, - filename=f"{name}.sto", - ) - - # chd files - chdlist0 = [] - chdlist0.append([(0, 0, 0), 100.0]) - chdlist0.append([(0, 0, 2), 100.0]) - - chdlist1 = [] - chdlist1.append([(0, 0, 0), 25.0]) - chdlist1.append([(0, 0, 2), 25.0]) - - chdspdict = {0: chdlist0, 1: chdlist1, 2: chdlist0} - chd = flopy.mf6.ModflowGwfchd( - gwf, - stress_period_data=chdspdict, - save_flows=False, - filename=f"{name}.chd", - ) - - # MAW - maw = flopy.mf6.ModflowGwfmaw( - gwf, - filename=f"{name}.maw", - budget_filerecord=f"{name}.maw.cbc", - print_input=True, - print_head=True, - print_flows=True, - save_flows=True, - observations=data.well.observations, - packagedata=data.well.packagedata, - connectiondata=data.well.connectiondata, - perioddata=data.well.perioddata, - pname="MAW-1", - ) - - # output control - oc = flopy.mf6.ModflowGwfoc( - gwf, - budget_filerecord=f"{name}.cbc", - head_filerecord=f"{name}.hds", - headprintrecord=[ - ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") - ], - saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - filename=f"{name}.oc", - ) - - return data, sim, None, self.eval_2 - - def eval_2(self, config, data): - print("evaluating MAW budgets...") - - shape3d = (data.nlay, data.nrow, data.ncol) - size3d = data.nlay * data.nrow * data.ncol - - # get results from listing file - fpth = os.path.join( - config.simpath, f"{os.path.basename(config.name)}.lst" - ) - budl = flopy.utils.Mf6ListBudget( - fpth, budgetkey="MAW-1 BUDGET FOR ENTIRE MODEL AT END OF TIME STEP" - ) - names = list(self.bud_lst) - d0 = budl.get_budget(names=names)[0] - dtype = d0.dtype - nbud = d0.shape[0] - - # get results from cbc file - cbc_bud = ["GWF", "RATE"] - d = np.recarray(nbud, dtype=dtype) - for key in self.bud_lst: - d[key] = 0.0 - fpth = os.path.join( - config.simpath, f"{os.path.basename(config.name)}.maw.cbc" - ) - cobj = flopy.utils.CellBudgetFile(fpth, precision="double") - kk = cobj.get_kstpkper() - times = cobj.get_times() - cbc_vals = [] - for idx, (k, t) in enumerate(zip(kk, times)): - for text in cbc_bud: - qin = 0.0 - qout = 0.0 - v = cobj.get_data(kstpkper=k, text=text)[0] - if isinstance(v, np.recarray): - vt = np.zeros(size3d, dtype=float) - wq = [] - for jdx, node in enumerate(v["node"]): - vt[node - 1] += v["q"][jdx] - wq.append(v["q"][jdx]) - v = vt.reshape(shape3d) - if text == cbc_bud[-1]: - cbc_vals.append(wq) - for kk in range(v.shape[0]): - for ii in range(v.shape[1]): - for jj in range(v.shape[2]): - vv = v[kk, ii, jj] - if vv < 0.0: - qout -= vv - else: - qin += vv - d["totim"][idx] = t - d["time_step"][idx] = k[0] - d["stress_period"] = k[1] - key = f"{text}_IN" - d[key][idx] = qin - key = f"{text}_OUT" - d[key][idx] = qout - - maw_vals = [ - [0.000, 0.000], - [-106.11303563809453, -96.22598985147631], - [-110.000, -130.000], - [0.0, -130.000], - [-110.000, -130.000], - ] - - # evaluate if well rates in cbc file are equal to expected values - diffv = [] - for ovs, svs in zip(maw_vals, cbc_vals): - for ov, sv in zip(ovs, svs): - diffv.append(ov - sv) - diffv = np.abs(np.array(diffv)).max() - msg = f"\nmaximum absolute maw rate difference ({diffv})\n" - - # calculate difference between water budget items in the lst and cbc files - diff = np.zeros((nbud, len(self.bud_lst)), dtype=float) - for idx, key in enumerate(self.bud_lst): - diff[:, idx] = d0[key] - d[key] - diffmax = np.abs(diff).max() - msg += f"maximum absolute total-budget difference ({diffmax}) " - - # write summary - fpth = os.path.join( - config.simpath, f"{os.path.basename(config.name)}.bud.cmp.out" - ) - f = open(fpth, "w") - for i in range(diff.shape[0]): - if i == 0: - line = f"{'TIME':>10s}" - for idx, key in enumerate(self.bud_lst): - line += f"{key + '_LST':>25s}" - line += f"{key + '_CBC':>25s}" - line += f"{key + '_DIF':>25s}" - f.write(line + "\n") - line = f"{d['totim'][i]:10g}" - for idx, key in enumerate(self.bud_lst): - line += f"{d0[key][i]:25g}" - line += f"{d[key][i]:25g}" - line += f"{diff[i, idx]:25g}" - f.write(line + "\n") - f.close() - - if diffmax > self.budtol or diffv > self.budtol: - config.success = False - msg += f"\n...exceeds {self.budtol}" - assert diffmax < self.budtol, msg - else: - config.success = True - print(" " + msg) - - case3 = Case( - name="maw03", - krylov="CG", - nlay=1, - nrow=101, - ncol=101, - nper=1, - delr=142, - delc=142, - perlen=[1000], - nstp=[50], - tsmult=[1.2], - strt=0, - hk=10, - nouter=100, - ninner=100, - hclose=1e-6, - rclose=1e-6, - relaxation_factor=1, - compare=False, - ) - cases3 = [ - case3.copy_update( - name="maw03a", - well=well3("maw03a"), - ), - case3.copy_update(name="maw03b", well=well3("maw03b")), - case3.copy_update(name="maw03c", well=well3("maw03c")), - ] - - @parametrize(data=cases3, ids=[c.name for c in cases3]) - def case_3(self, function_tmpdir, targets, data): - top = 0.0 - botm = [-1000.0] - - tdis_rc = [] - for i in range(data.nper): - tdis_rc.append((data.perlen[i], data.nstp[i], data.tsmult[i])) - - name = data.name - ws = str(function_tmpdir) - sim = flopy.mf6.MFSimulation( - sim_name=name, sim_ws=ws, exe_name=targets["mf6"] - ) - - # create tdis package - tdis = flopy.mf6.ModflowTdis( - sim, time_units="DAYS", nper=data.nper, perioddata=tdis_rc - ) - - # create gwf model - gwf = flopy.mf6.MFModel( - sim, - model_type="gwf6", - modelname=name, - model_nam_file=f"{name}.nam", - ) - - # create iterative model solution and register the gwf model with it - ims = flopy.mf6.ModflowIms( - sim, - print_option="SUMMARY", - outer_dvclose=data.hclose, - outer_maximum=data.nouter, - under_relaxation="NONE", - inner_maximum=data.ninner, - inner_dvclose=data.hclose, - rcloserecord=data.rclose, - linear_acceleration=data.krylov, - scaling_method="NONE", - reordering_method="NONE", - relaxation_factor=data.relaxation_factor, - ) - sim.register_ims_package(ims, [gwf.name]) - - dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=data.nlay, - nrow=data.nrow, - ncol=data.ncol, - delr=data.delr, - delc=data.delc, - top=top, - botm=botm, - idomain=1, - filename=f"{name}.dis", - ) - - # initial conditions - ic = flopy.mf6.ModflowGwfic(gwf, strt=data.strt, filename=f"{name}.ic") - - # node property flow - npf = flopy.mf6.ModflowGwfnpf( - gwf, - save_flows=True, - icelltype=1, - k=data.hk, - k33=data.hk, - filename=f"{name}.npf", - ) - - # storage - sto = flopy.mf6.ModflowGwfsto( - gwf, - save_flows=True, - iconvert=0, - ss=1.0e-5, - sy=0.1, - steady_state={0: False}, - transient={0: True}, - filename=f"{name}.sto", - ) - - # MAW - maw = flopy.mf6.ModflowGwfmaw( - gwf, - filename=f"{name}.maw", - print_input=True, - print_head=True, - print_flows=True, - save_flows=True, - observations=data.well.observations, - packagedata=data.well.packagedata, - connectiondata=data.well.connectiondata, - perioddata=data.well.perioddata, - ) - - # output control - oc = flopy.mf6.ModflowGwfoc( - gwf, - budget_filerecord=f"{name}.cbc", - head_filerecord=f"{name}.hds", - headprintrecord=[ - ("COLUMNS", data.ncol, "WIDTH", 15, "DIGITS", 6, "GENERAL") - ], - saverecord=[("HEAD", "ALL")], - printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - filename=f"{name}.oc", - ) - - # head observations - obs_data0 = [("head_well_cell", "HEAD", (0, 0, 0))] - obs_recarray = {f"{name}.obs.csv": obs_data0} - obs = flopy.mf6.ModflowUtlobs( - gwf, - pname="head_obs", - filename=f"{name}.obs", - digits=15, - print_input=True, - continuous=obs_recarray, - ) - - return data, sim, None, self.eval_3 - - def eval_3(self, config, data): - print("evaluating MAW heads...") - - # MODFLOW 6 maw results - test_name = config.name - case_name = data.name - fpth = os.path.join(config.simpath, f"{test_name}.maw.obs.csv") - tc = np.genfromtxt(fpth, names=True, delimiter=",") - - if case_name == "a": - - # M1RATE should be 2000. - msg = "The injection rate should be 2000. for all times" - assert tc["M1RATE"].min() == tc["M1RATE"].max() == 2000, msg - - elif case_name == "b": - - # M1RATE should have a minimum value less than 200 and - # M1HEAD should not exceed 0.400001 - msg = ( - "Injection rate should fall below 200 and the head should not" - "exceed 0.4" - ) - assert tc["M1RATE"].min() < 200.0, msg - assert tc["M1HEAD"].max() < 0.400001, msg - - elif case_name == "c": - - # M1RATE should have a minimum value less than 800 - # M1HEAD should not exceed 1.0. - msg = ( - "Min injection rate should be less than 800 and well " - "head should not exceed 1.0" - ) - assert tc["M1RATE"].min() < 800.0 and tc["M1HEAD"].max() < 1.0, msg - - # TODO - # case4 = Case( - # name='maw_iss305', - # krylov='CG', - # nlay=2, - # nrow=101, - # ncol=101, - # nper=2, - # delr=142, - # delc=142, - # perlen=[0.0, 365.0], - # nstp=[1, 25], - # tsmult=[1.0, 1.1], - # steady=[True, False], - # well=None, - # strt=0, - # hk=10, - # nouter=100, - # ninner=100, - # hclose=1e-9, - # rclose=1e-6, - # relax=1, - # require_failure=True - # ) - # cases4 = [ - # case4.copy_update({ - # 'name': "maw_iss305a", - # 'well': well4(case4, "CUMULATIVE"), - # 'require_failure': False - # }), - # case4.copy_update({ - # 'name': "maw_iss305b", - # 'well': well4(case4, "SKIN") - # }), - # case4.copy_update({ - # 'name': "maw_iss305c", - # 'well': well4(case4, "SKIN") - # }), - # case4.copy_update({ - # 'name': "maw_iss305d", - # 'well': well4(case4, "SKIN") - # }), - # case4.copy_update({ - # 'name': "maw_iss305e", - # 'well': well4(case4, "SPECIFIED") - # }), - # case4.copy_update({ - # 'name': "maw_iss305f", - # 'well': well4(case4, "CUMULATIVE") - # }) - # ] - - # @parametrize(data=cases4, ids=[c['name'] for c in cases4]) - # def case_4(self, function_tmpdir, targets, data): - # pass - - # def eval_4(self, sim, data): - # pass diff --git a/autotest/test_gwf_maw_obs.py b/autotest/test_gwf_maw_obs.py index 7455949fc37..050de4f96b9 100644 --- a/autotest/test_gwf_maw_obs.py +++ b/autotest/test_gwf_maw_obs.py @@ -14,7 +14,6 @@ def build_model(dir, exe): - nlay, nrow, ncol = 1, 1, 3 nper = 3 perlen = [1.0, 1.0, 1.0]