Skip to content

Commit

Permalink
test: move prt autotests to folder, deduplicate model setup
Browse files Browse the repository at this point in the history
  • Loading branch information
wpbonelli committed Oct 6, 2023
1 parent d556e20 commit 2681687
Show file tree
Hide file tree
Showing 9 changed files with 336 additions and 661 deletions.
110 changes: 101 additions & 9 deletions autotest/prt_test_utils.py → autotest/prt/prt_test_utils.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import os
from typing import Union
from types import SimpleNamespace
from typing import Tuple, Union

import flopy
import numpy as np
Expand Down Expand Up @@ -39,7 +40,7 @@ def check_track_data(
# check each column separately to avoid:
# TypeError: The DType <class 'numpy._FloatAbstractDType'> could not be promoted by <class 'numpy.dtype[void]'>
for k in data_bin.dtype.names:
if k == 'name':
if k == "name":
continue
assert np.allclose(data_bin[k], data_csv[k], equal_nan=True)

Expand Down Expand Up @@ -93,18 +94,23 @@ def get_output_event(case_name):
if "wksk" in case_name
else "TERMINATE"
if "terminate" in case_name
else "ALL" # default
else "ALL" # default
)


def get_ireason_code(output_event):
return (
0 if output_event == "RELEASE"
else 1 if output_event == "TRANSIT"
else 2 if output_event == "TIMESTEP"
else 3 if output_event == "TERMINATE"
else 4 if output_event == "WEAKSINK"
else -1 # default
0
if output_event == "RELEASE"
else 1
if output_event == "TRANSIT"
else 2
if output_event == "TIMESTEP"
else 3
if output_event == "TERMINATE"
else 4
if output_event == "WEAKSINK"
else -1 # default
)


Expand Down Expand Up @@ -140,3 +146,89 @@ def get_partdata(grid, rpts):
timeoffset=0,
drape=0,
)


def build_gwf_sim(
name, ws, mf6
) -> Tuple[flopy.mf6.MFSimulation, SimpleNamespace]:
ctx = SimpleNamespace(
nlay=1,
nrow=10,
ncol=10,
top=1.0,
botm=[0.0],
nper=1,
perlen=1.0,
nstp=1,
tsmult=1.0,
porosity=0.1,
)

# create simulation
sim = flopy.mf6.MFSimulation(
sim_name=name,
exe_name=mf6,
version="mf6",
sim_ws=ws,
)

# create tdis package
flopy.mf6.modflow.mftdis.ModflowTdis(
sim,
pname="tdis",
time_units="DAYS",
nper=ctx.nper,
perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
)

# create gwf model
gwfname = f"{name}_gwf"
gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True)

# create gwf discretization
flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
gwf,
pname="dis",
nlay=ctx.nlay,
nrow=ctx.nrow,
ncol=ctx.ncol,
)

# create gwf initial conditions package
flopy.mf6.modflow.mfgwfic.ModflowGwfic(gwf, pname="ic")

# create gwf node property flow package
flopy.mf6.modflow.mfgwfnpf.ModflowGwfnpf(
gwf,
pname="npf",
save_saturation=True,
save_specific_discharge=True,
)

# create gwf chd package
spd = {
0: [[(0, 0, 0), 1.0, 1.0], [(0, 9, 9), 0.0, 0.0]],
1: [[(0, 0, 0), 0.0, 0.0], [(0, 9, 9), 1.0, 2.0]],
}
chd = flopy.mf6.ModflowGwfchd(
gwf,
pname="CHD-1",
stress_period_data=spd,
auxiliary=["concentration"],
)

# create gwf output control package
# output file names
gwf_budget_file = f"{gwfname}.bud"
gwf_head_file = f"{gwfname}.hds"
oc = flopy.mf6.ModflowGwfoc(
gwf,
budget_filerecord=gwf_budget_file,
head_filerecord=gwf_head_file,
saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
)

# create iterative model solution for gwf model
ims = flopy.mf6.ModflowIms(sim)

return sim, ctx
62 changes: 21 additions & 41 deletions autotest/test_prt_exg01.py → autotest/prt/test_prt_exg01.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
Test GWF and PRT models in the same simulation
with an exchange.
The grid is a 10x10 square with a single layer.
The same flow system shown on the FloPy readme.
The grid is a 10x10 square with a single layer,
the same flow system shown on the FloPy readme.
Particles are released from the top left cell.
Results are compared against a MODPATH 7 model.
Expand All @@ -21,39 +21,21 @@
import numpy as np
import pandas as pd
import pytest
from flopy.plot.plotutil import to_mp7_pathlines
from flopy.utils import PathlineFile
from flopy.utils.binaryfile import HeadFile
from flopy.plot.plotutil import to_mp7_pathlines
from prt_test_utils import (all_equal, build_gwf_sim, check_budget_data,
check_track_data, has_default_boundnames)

from framework import TestFramework
from prt_test_utils import (
all_equal,
check_budget_data,
check_track_data,
has_default_boundnames,
)
from simulation import TestSimulation


# simulation/model names
simname = "prtexg1"
simname = "prtexg01"

# test cases
ex = [simname, f"{simname}bnms"]


# model info
nlay = 1
nrow = 10
ncol = 10
top = 1.0
botm = [0.0]
nper = 1
perlen = 1.0
nstp = 1
tsmult = 1.0
porosity = 0.1

# release points
releasepts = [
# index, k, i, j, x, y, z
Expand All @@ -75,11 +57,9 @@ def get_model_name(idx, mdl):


def build_sim(idx, ws, mf6):
from test_prt_fmi01 import build_gwf_sim

# create simulation
name = ex[idx]
sim = build_gwf_sim(name, ws, mf6)
sim, ctx = build_gwf_sim(name, ws, mf6)

# create prt model
prtname = get_model_name(idx, "prt")
Expand All @@ -89,13 +69,13 @@ def build_sim(idx, ws, mf6):
flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
prt,
pname="dis",
nlay=nlay,
nrow=nrow,
ncol=ncol,
nlay=ctx.nlay,
nrow=ctx.nrow,
ncol=ctx.ncol,
)

# create mip package
flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)

# create prp package
rpts = (
Expand Down Expand Up @@ -124,7 +104,7 @@ def build_sim(idx, ws, mf6):
)

# create a flow model interface
# todo Fienen's report (crash when FMI created but not needed)
# todo Mike Fienen's report (crash when FMI created but not needed)
# flopy.mf6.ModflowPrtfmi(
# prt,
# packagedata=[
Expand All @@ -151,10 +131,10 @@ def build_sim(idx, ws, mf6):
)
sim.register_solution_package(ems, [prt.name])

return sim
return sim, ctx


def build_mp7_sim(idx, ws, mp7, gwf):
def build_mp7_sim(ctx, idx, ws, mp7, gwf):
partdata = flopy.modpath.ParticleData(
partlocs=[p[0] for p in releasepts_mp7],
localx=[p[1] for p in releasepts_mp7],
Expand All @@ -177,7 +157,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
)
mpbas = flopy.modpath.Modpath7Bas(
mp,
porosity=porosity,
porosity=ctx.porosity,
)
mpsim = flopy.modpath.Modpath7Sim(
mp,
Expand All @@ -191,12 +171,12 @@ def build_mp7_sim(idx, ws, mp7, gwf):
return mp


def eval_results(sim):
def eval_results(ctx, sim):
print(f"Evaluating results for sim {sim.name}")
simpath = Path(sim.simpath)

# check budget data
check_budget_data(simpath / f"{sim.name}_prt.lst", perlen, nper)
check_budget_data(simpath / f"{sim.name}_prt.lst", ctx.perlen, ctx.nper)

# check particle track data
prt_track_file = simpath / f"{sim.name}_prt.trk"
Expand All @@ -215,15 +195,15 @@ def eval_results(sim):
@pytest.mark.parametrize("idx, name", enumerate(ex))
def test_mf6model(idx, name, function_tmpdir, targets):
ws = function_tmpdir
sim = build_sim(idx, str(ws), targets.mf6)
sim, ctx = build_sim(idx, str(ws), targets.mf6)
sim.write_simulation()

test = TestFramework()
test.run(
TestSimulation(
name=name,
exe_dict=targets,
exfunc=eval_results,
exfunc=lambda s: eval_results(ctx, s),
idxsim=0,
make_comparison=False,
),
Expand All @@ -243,7 +223,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
mg = gwf.modelgrid

# build mp7 model
mp7sim = build_mp7_sim(idx, ws, targets.mp7, gwf)
mp7sim = build_mp7_sim(ctx, idx, ws, targets.mp7, gwf)

# run mp7 model
mp7sim.write_input()
Expand Down Expand Up @@ -293,7 +273,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
assert pd.isna(mf6_pls["name"]).all()

# check budget data were written to mf6 prt list file
check_budget_data(ws / f"{name}_prt.lst", perlen, nper)
check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)

# check mf6 prt particle track data were written to binary/CSV files
check_track_data(
Expand Down
Loading

0 comments on commit 2681687

Please sign in to comment.