Skip to content

Commit

Permalink
ci(plotting): use plot hook in prt tests
Browse files Browse the repository at this point in the history
  • Loading branch information
wpbonelli committed Dec 6, 2024
1 parent 27b7820 commit c441801
Show file tree
Hide file tree
Showing 14 changed files with 1,128 additions and 803 deletions.
163 changes: 100 additions & 63 deletions autotest/test_prt_budget.py
Original file line number Diff line number Diff line change
Expand Up @@ -195,16 +195,10 @@ def check_output(idx, test):
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()
assert (prt_ws / prt_track_csv_file).is_file()
assert (prt_ws / prp_track_file).is_file()
assert (prt_ws / prp_track_csv_file).is_file()

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

# load mf6 pathline results
mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)

# load mp7 pathline results
plf = PathlineFile(mp7_ws / mp7_pathline_file)
Expand All @@ -216,8 +210,21 @@ def check_output(idx, test):
mp7_pls["node"] = mp7_pls["node"] + 1
mp7_pls["k"] = mp7_pls["k"] + 1

# load mf6 pathline results
mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
# extract head, budget, and specific discharge results from GWF model
hds = HeadFile(gwf_ws / gwf_head_file).get_data()
bud = gwf.output.budget()
spdis = bud.get_data(text="DATA-SPDIS")[0]
qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)

assert (gwf_ws / gwf_budget_file).is_file()
assert (gwf_ws / gwf_head_file).is_file()
assert (prt_ws / prt_track_file).is_file()
assert (prt_ws / prt_track_csv_file).is_file()
assert (prt_ws / prp_track_file).is_file()
assert (prt_ws / prp_track_csv_file).is_file()

# check mp7 output files exist
assert (mp7_ws / mp7_pathline_file).is_file()

# make sure pathline df has "name" (boundname) column and default values
assert "name" in mf6_pls
Expand Down Expand Up @@ -260,57 +267,6 @@ def check_output(idx, test):
track_csv=track_csv,
)

# extract head, budget, and specific discharge results from GWF model
hds = HeadFile(gwf_ws / gwf_head_file).get_data()
bud = gwf.output.budget()
spdis = bud.get_data(text="DATA-SPDIS")[0]
qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)

# setup plot
plot_results = False
if plot_results:
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))
for a in ax:
a.set_aspect("equal")

# plot mf6 pathlines in map view
pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0])
pmv.plot_grid()
pmv.plot_array(hds[0], alpha=0.1)
pmv.plot_vector(qx, qy, normalize=True, color="white")
mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"])
for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
pl.plot(
title="MF6 pathlines",
kind="line",
x="x",
y="y",
ax=ax[0],
legend=False,
color=cm.plasma(ipl / len(mf6_plines)),
)

# plot mp7 pathlines in map view
pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[1])
pmv.plot_grid()
pmv.plot_array(hds[0], alpha=0.1)
pmv.plot_vector(qx, qy, normalize=True, color="white")
mp7_plines = mp7_pls.groupby(["particleid"])
for ipl, (pid, pl) in enumerate(mp7_plines):
pl.plot(
title="MP7 pathlines",
kind="line",
x="x",
y="y",
ax=ax[1],
legend=False,
color=cm.plasma(ipl / len(mp7_plines)),
)

# view/save plot
plt.show()
plt.savefig(gwf_ws / f"test_{simname}.png")

# convert mf6 pathlines to mp7 format
mf6_pls = to_mp7_pathlines(mf6_pls)

Expand All @@ -335,13 +291,94 @@ def check_output(idx, test):
assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)


def plot_output(idx, test):
name = test.name
gwf_ws = test.workspace
prt_ws = test.workspace / "prt"
mp7_ws = test.workspace / "mp7"
gwf_name = get_model_name(name, "gwf")
prt_name = get_model_name(name, "prt")
mp7_name = get_model_name(name, "mp7")
gwf_sim = test.sims[0]
gwf = gwf_sim.get_model(gwf_name)
mg = gwf.modelgrid

# check mf6 output files exist
gwf_head_file = f"{gwf_name}.hds"
prt_track_csv_file = f"{prt_name}.trk.csv"
mp7_pathline_file = f"{mp7_name}.mppth"

# load mf6 pathline results
mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)

# load mp7 pathline results
plf = PathlineFile(mp7_ws / mp7_pathline_file)
mp7_pls = pd.DataFrame(
plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
)
# convert zero-based to one-based indexing in mp7 results
mp7_pls["particlegroup"] = mp7_pls["particlegroup"] + 1
mp7_pls["node"] = mp7_pls["node"] + 1
mp7_pls["k"] = mp7_pls["k"] + 1

# extract head, budget, and specific discharge results from GWF model
hds = HeadFile(gwf_ws / gwf_head_file).get_data()
bud = gwf.output.budget()
spdis = bud.get_data(text="DATA-SPDIS")[0]
qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)

# set up plot
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))
for a in ax:
a.set_aspect("equal")

# plot mf6 pathlines in map view
pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0])
pmv.plot_grid()
pmv.plot_array(hds[0], alpha=0.1)
pmv.plot_vector(qx, qy, normalize=True, color="white")
mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"])
for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
pl.plot(
title="MF6 pathlines",
kind="line",
x="x",
y="y",
ax=ax[0],
legend=False,
color=cm.plasma(ipl / len(mf6_plines)),
)

# plot mp7 pathlines in map view
pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[1])
pmv.plot_grid()
pmv.plot_array(hds[0], alpha=0.1)
pmv.plot_vector(qx, qy, normalize=True, color="white")
mp7_plines = mp7_pls.groupby(["particleid"])
for ipl, (pid, pl) in enumerate(mp7_plines):
pl.plot(
title="MP7 pathlines",
kind="line",
x="x",
y="y",
ax=ax[1],
legend=False,
color=cm.plasma(ipl / len(mp7_plines)),
)

# view/save plot
plt.show()
plt.savefig(prt_ws / f"test_{simname}.png")


@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,
compare=None,
)
Expand Down
Loading

0 comments on commit c441801

Please sign in to comment.