Skip to content

Commit

Permalink
run mp7 comparison models with test framework
Browse files Browse the repository at this point in the history
  • Loading branch information
wpbonelli committed Dec 31, 2023
1 parent 091eb8f commit 39fe521
Show file tree
Hide file tree
Showing 12 changed files with 360 additions and 449 deletions.
30 changes: 17 additions & 13 deletions autotest/framework.py
Original file line number Diff line number Diff line change
Expand Up @@ -376,10 +376,19 @@ def _compare_heads(self, cpth=None, extensions="hds") -> bool:
else:
files2.append(None)

if self.cmp_namefile is None:
# todo: clean up namfile path detection?
nf = next(iter(get_namefiles(cpth)), None)
cmp_namefile = (
None
if "mf6" in self.compare or "libmf6" in self.compare
else os.path.basename(nf)
if nf
else None
)
if cmp_namefile is None:
pth = None
else:
pth = os.path.join(cpth, self.cmp_namefile)
pth = os.path.join(cpth, cmp_namefile)

for i in range(len(files1)):
file1 = files1[i]
Expand Down Expand Up @@ -642,16 +651,6 @@ def run_sim_or_model(
if self.verbose:
print(f"Running {target} in {workspace}")

# needed in _compare_heads()... todo: inject explicitly?
nf = next(iter(get_namefiles(workspace)), None)
self.cmp_namefile = (
None
if "mf6" in target.name or "libmf6" in target.name
else os.path.basename(nf)
if nf
else None
)

# run the model
try:
# via MODFLOW API
Expand Down Expand Up @@ -694,8 +693,13 @@ def run_sim_or_model(
else:
# non-MF6 model
try:
nf_ext = ".mpsim" if "mp7" in target.name else ".nam"
namefile = next(iter(workspace.glob(f"*{nf_ext}")), None)
assert (
namefile
), f"Control file with extension {nf_ext} not found"
success, buff = flopy.run_model(
target, self.cmp_namefile, workspace, report=True
target, namefile, workspace, report=True
)
except Exception:
warn(f"{target} model failed:\n{format_exc()}")
Expand Down
106 changes: 48 additions & 58 deletions autotest/prt/test_prt_disv.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,14 +19,9 @@
from flopy.utils import PathlineFile
from flopy.utils.binaryfile import HeadFile
from flopy.utils.gridutil import get_disv_kwargs
from prt_test_utils import (
all_equal,
check_budget_data,
check_track_data,
get_partdata,
has_default_boundnames,
plot_nodes_and_vertices,
)
from prt_test_utils import (all_equal, check_budget_data, check_track_data,
get_partdata, has_default_boundnames,
plot_nodes_and_vertices)

from framework import TestFramework

Expand Down Expand Up @@ -72,9 +67,9 @@


def build_gwf_sim(idx, ws, mf6):
gwfname = f"{cases[idx]}_gwf"
gwf_name = f"{cases[idx]}_gwf"
sim = flopy.mf6.MFSimulation(
sim_name=gwfname, version="mf6", exe_name=mf6, sim_ws=ws
sim_name=gwf_name, version="mf6", exe_name=mf6, sim_ws=ws
)

# create tdis package
Expand All @@ -84,7 +79,7 @@ def build_gwf_sim(idx, ws, mf6):

# create gwf model
gwf = flopy.mf6.ModflowGwf(
sim, modelname=gwfname, newtonoptions="NEWTON", save_flows=True
sim, modelname=gwf_name, newtonoptions="NEWTON", save_flows=True
)

# create iterative model solution and register the gwf model with it
Expand Down Expand Up @@ -133,12 +128,12 @@ def build_gwf_sim(idx, ws, mf6):
# output control
oc = flopy.mf6.ModflowGwfoc(
gwf,
budget_filerecord="{}.cbc".format(gwfname),
head_filerecord="{}.hds".format(gwfname),
budget_filerecord="{}.cbc".format(gwf_name),
head_filerecord="{}.hds".format(gwf_name),
headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
filename="{}.oc".format(gwfname),
filename="{}.oc".format(gwf_name),
)

# Print human-readable heads
Expand All @@ -147,7 +142,7 @@ def build_gwf_sim(idx, ws, mf6):
for i in np.arange(40, 50, 1):
obs_lst.append(["obs_" + str(i + 1), "head", (k, i)])

obs_dict = {f"{gwfname}.obs.csv": obs_lst}
obs_dict = {f"{gwf_name}.obs.csv": obs_lst}
obs = flopy.mf6.ModflowUtlobs(
gwf, pname="head_obs", digits=20, continuous=obs_dict
)
Expand All @@ -171,8 +166,8 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
)

# create prt model
prtname = f"{cases[idx]}_prt"
prt = flopy.mf6.ModflowPrt(sim, modelname=prtname)
prt_name = f"{cases[idx]}_prt"
prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name)

# create prt discretization
disv = flopy.mf6.ModflowGwfdisv(prt, **disvkwargs)
Expand All @@ -188,12 +183,12 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
releasepts[0] = (0, (0, 1), 0.5, 9.5, 22.5)

# create prp package
prp_track_file = f"{prtname}.prp.trk"
prp_track_csv_file = f"{prtname}.prp.trk.csv"
prp_track_file = f"{prt_name}.prp.trk"
prp_track_csv_file = f"{prt_name}.prp.trk.csv"
flopy.mf6.ModflowPrtprp(
prt,
pname="prp1",
filename=f"{prtname}_1.prp",
filename=f"{prt_name}_1.prp",
nreleasepts=len(releasepts),
packagedata=releasepts,
perioddata={0: ["FIRST"]},
Expand All @@ -204,8 +199,8 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
)

# create output control package
prt_track_file = f"{prtname}.trk"
prt_track_csv_file = f"{prtname}.trk.csv"
prt_track_file = f"{prt_name}.trk"
prt_track_csv_file = f"{prt_name}.trk.csv"
flopy.mf6.ModflowPrtoc(
prt,
pname="oc",
Expand All @@ -214,9 +209,9 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
)

# create the flow model interface
gwfname = f"{cases[idx]}_gwf"
gwf_budget_file = gwf_ws / f"{gwfname}.cbc"
gwf_head_file = gwf_ws / f"{gwfname}.hds"
gwf_name = f"{cases[idx]}_gwf"
gwf_budget_file = gwf_ws / f"{gwf_name}.cbc"
gwf_head_file = gwf_ws / f"{gwf_name}.hds"
flopy.mf6.ModflowPrtfmi(
prt,
packagedata=[
Expand All @@ -229,7 +224,7 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
ems = flopy.mf6.ModflowEms(
sim,
pname="ems",
filename=f"{prtname}.ems",
filename=f"{prt_name}.ems",
)
sim.register_solution_package(ems, [prt.name])

Expand All @@ -238,17 +233,19 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):

def build_mp7_sim(idx, ws, mp7, gwf):
partdata = get_partdata(gwf.modelgrid, releasepts_mp7)
mp7name = f"{cases[idx]}_mp7"
mp7_name = f"{cases[idx]}_mp7"
pg = flopy.modpath.ParticleGroup(
particlegroupname="G1",
particledata=partdata,
filename=f"{mp7name}.sloc",
filename=f"{mp7_name}.sloc",
)
mp = flopy.modpath.Modpath7(
modelname=mp7name,
modelname=mp7_name,
flowmodel=gwf,
exe_name=mp7,
model_ws=ws,
headfilename=f"{gwf.name}.hds",
budgetfilename=f"{gwf.name}.cbc",
)
mpbas = flopy.modpath.Modpath7Bas(
mp,
Expand All @@ -267,50 +264,43 @@ def build_mp7_sim(idx, ws, mp7, gwf):


def build_models(idx, test):
gwfsim = build_gwf_sim(idx, test.workspace, test.targets.mf6)
prtsim = build_prt_sim(
gwf_sim = build_gwf_sim(idx, test.workspace, test.targets.mf6)
prt_sim = build_prt_sim(
idx, test.workspace, test.workspace / "prt", test.targets.mf6
)
return gwfsim, prtsim
mp7_sim = build_mp7_sim(
idx, test.workspace / "mp7", test.targets.mp7, gwf_sim.get_model()
)
return gwf_sim, prt_sim, mp7_sim


def check_output(idx, test):
name = test.name
gwf_ws = test.workspace
prt_ws = test.workspace / "prt"
mp7_ws = test.workspace / "mp7"
gwfname = f"{name}_gwf"
prtname = f"{name}_prt"
mp7name = f"{name}_mp7"
gwf_name = f"{name}_gwf"
prt_name = f"{name}_prt"
mp7_name = f"{name}_mp7"
gwf_sim = test.sims[0]
prt_sim = test.sims[1]
gwf = gwf_sim.get_model(gwf_name)
prt = prt_sim.get_model(prt_name)
mg = gwf.modelgrid

# if invalid release points, check for error message
if "bprp" in name:
buff = test.buffs[1]
assert any("Error: release point" in l for l in buff)
return

# extract mf6 simulations/models and grid
gwfsim = test.sims[0]
prtsim = test.sims[1]
gwf = gwfsim.get_model(gwfname)
prt = prtsim.get_model(prtname)
mg = gwf.modelgrid

# todo build mp7 model
mp7sim = build_mp7_sim(idx, mp7_ws, test.targets.mp7, gwf)

# todo run mp7 model
mp7sim.write_input()
success, buff = mp7sim.run_model(report=True)
assert success, pformat(buff)

# check mf6 output files exist
gwf_budget_file = f"{gwfname}.cbc"
gwf_head_file = f"{gwfname}.hds"
prt_track_file = f"{prtname}.trk"
prt_track_csv_file = f"{prtname}.trk.csv"
prp_track_file = f"{prtname}.prp.trk"
prp_track_csv_file = f"{prtname}.prp.trk.csv"
gwf_budget_file = f"{gwf_name}.cbc"
gwf_head_file = f"{gwf_name}.hds"
prt_track_file = f"{prt_name}.trk"
prt_track_csv_file = f"{prt_name}.trk.csv"
prp_track_file = f"{prt_name}.prp.trk"
prp_track_csv_file = f"{prt_name}.prp.trk.csv"
assert (gwf_ws / gwf_budget_file).is_file()
assert (gwf_ws / gwf_head_file).is_file()
assert (prt_ws / prt_track_file).is_file()
Expand All @@ -319,7 +309,7 @@ def check_output(idx, test):
assert (prt_ws / prp_track_csv_file).is_file()

# check mp7 output files exist
mp7_pathline_file = f"{mp7name}.mppth"
mp7_pathline_file = f"{mp7_name}.mppth"
assert (mp7_ws / mp7_pathline_file).is_file()

# load mp7 pathline results
Expand Down Expand Up @@ -446,6 +436,6 @@ def test_mf6model(idx, name, function_tmpdir, targets):
check=lambda t: check_output(idx, t),
targets=targets,
compare=None,
xfail=[False, "bprp" in name],
xfail=[False, "bprp" in name, False],
)
test.run()
Loading

0 comments on commit 39fe521

Please sign in to comment.