diff --git a/autotest/test_prt_budget.py b/autotest/test_prt_budget.py index 52fbe339401..acf9cdb5d96 100644 --- a/autotest/test_prt_budget.py +++ b/autotest/test_prt_budget.py @@ -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) @@ -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 @@ -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) @@ -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, ) diff --git a/autotest/test_prt_disv1.py b/autotest/test_prt_disv1.py index 98d8bbf65d8..d9d62e7793a 100644 --- a/autotest/test_prt_disv1.py +++ b/autotest/test_prt_disv1.py @@ -302,56 +302,6 @@ def build_models(idx, test): return gwf_sim, prt_sim, mp7_sim -def plot_output(name, gwf, head, spdis, mf6_pls, mp7_pls, fpath): - # setup 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=gwf.modelgrid, ax=ax[0]) - pmv.plot_grid() - pmv.plot_array(head[0], alpha=0.2) - pmv.plot_vector(*spdis, normalize=True, color="white") - - # plot labeled nodes and vertices - plot_nodes_and_vertices(gwf, gwf.modelgrid, None, gwf.modelgrid.ncpl, ax[0]) - mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"]) - for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): - pl.plot( - title="MF6 pathlines", - linestyle="None" if "trst" in name else "--", - marker="o", - markersize=2, - 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=gwf.modelgrid, ax=ax[1]) - pmv.plot_grid() - pmv.plot_array(head[0], alpha=0.2) - pmv.plot_vector(*spdis, 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(fpath) - - def compare_output(name, mf6_pls, mp7_pls, mp7_eps): mf6_eps = mf6_pls[mf6_pls.ireason == 3] # get prt endpoints mp7_eps = to_mp7_pathlines(mp7_eps) # convert mp7 pathlines to mp7 format @@ -405,6 +355,13 @@ def check_output(idx, test): prt = prt_sim.get_model(prt_name) mg = gwf.modelgrid + 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" + mp7_pathline_file = f"{mp7_name}.mppth" + mp7_endpoint_file = f"{mp7_name}.mpend" + # if invalid release points, check for error message if "bprp" in name: buff = test.buffs[1] @@ -412,10 +369,6 @@ def check_output(idx, test): return # check mf6 output files exist - 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" assert (gwf_ws / gwf_budget_file).is_file() assert (gwf_ws / gwf_head_file).is_file() assert (prt_ws / prt_track_file).is_file() @@ -427,8 +380,6 @@ def check_output(idx, test): assert (prt_ws / prp_track_csv_file).is_file() # check mp7 output files exist - mp7_pathline_file = f"{mp7_name}.mppth" - mp7_endpoint_file = f"{mp7_name}.mpend" assert (mp7_ws / mp7_pathline_file).is_file() assert (mp7_ws / mp7_endpoint_file).is_file() @@ -486,25 +437,114 @@ def check_output(idx, test): spdis = bud.get_data(text="DATA-SPDIS")[0] qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) - # plot results if enabled - plot = False - if "bprp" not in name and plot: - plot_output( - name, gwf, head, (qx, qy), mf6_pls, mp7_pls, gwf_ws / f"test_{simname}.png" - ) - # compare prt and mp7 results compare_output(name, mf6_pls, mp7_pls, mp7_eps) +def plot_output(idx, test): + name = test.name + gwf_ws = test.workspace + prt_ws = test.workspace / "prt" + mp7_ws = test.workspace / "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 + + gwf_head_file = f"{gwf_name}.hds" + prt_track_csv_file = f"{prt_name}.trk.csv" + mp7_pathline_file = f"{mp7_name}.mppth" + mp7_endpoint_file = f"{mp7_name}.mpend" + + # extract head, budget, and specific discharge results from GWF model + head = 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) + + # 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 + + # load mp7 endpoint results + epf = EndpointFile(mp7_ws / mp7_endpoint_file) + mp7_eps = pd.DataFrame(epf.get_destination_endpoint_data(range(mg.nnodes))) + # convert zero-based to one-based indexing in mp7 results + mp7_eps["particlegroup"] = mp7_eps["particlegroup"] + 1 + mp7_eps["node"] = mp7_eps["node"] + 1 + mp7_eps["k"] = mp7_eps["k"] + 1 + + # setup 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=gwf.modelgrid, ax=ax[0]) + pmv.plot_grid() + pmv.plot_array(head[0], alpha=0.2) + pmv.plot_vector(qx, qy, normalize=True, color="white") + + # plot labeled nodes and vertices + plot_nodes_and_vertices(gwf, gwf.modelgrid, None, gwf.modelgrid.ncpl, ax[0]) + mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"]) + for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): + pl.plot( + title="MF6 pathlines", + linestyle="None" if "trst" in name else "--", + marker="o", + markersize=2, + 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=gwf.modelgrid, ax=ax[1]) + pmv.plot_grid() + pmv.plot_array(head[0], alpha=0.2) + 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"{name}.png") + + @pytest.mark.developmode @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 and "bprp" not in name) else None, targets=targets, compare=None, xfail=[False, "bprp" in name, False], diff --git a/autotest/test_prt_drape.py b/autotest/test_prt_drape.py index 7a7cf453f1b..703f4ef0d13 100644 --- a/autotest/test_prt_drape.py +++ b/autotest/test_prt_drape.py @@ -244,6 +244,16 @@ 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" + + # 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) + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + assert (gwf_ws / gwf_budget_file).is_file() assert (gwf_ws / gwf_head_file).is_file() assert (prt_ws / prt_track_file).is_file() @@ -251,16 +261,10 @@ def check_output(idx, test): assert (prt_ws / prp_track_file).is_file() assert (prt_ws / prp_track_csv_file).is_file() - # load mf6 pathline results - mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) - # make sure all mf6 pathline data have correct model and PRP index (1) assert all_equal(mf6_pls["imdl"], 1) assert all_equal(mf6_pls["iprp"], 1) - # check budget data were written to mf6 prt list file - # check_budget_data(ws / f"{name}_prt.lst", perlen, nper) - # check mf6 prt particle track data were written to binary/CSV files # and that different formats are equal for track_csv in [ @@ -273,39 +277,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=1, figsize=(10, 10)) - ax.set_aspect("equal") - - # plot mf6 pathlines in map view - pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) - 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=f"MF6 pathlines{' (drape)' if drape else ''}", - kind="line", - x="x", - y="y", - ax=ax, - legend=False, - color=cm.plasma(ipl / len(mf6_plines)), - ) - - # view/save plot - plt.show() - plt.savefig(gwf_ws / f"test_{simname}.png") - if drape: assert mf6_pls.shape[0] == 36 else: @@ -315,13 +286,64 @@ def check_output(idx, test): assert mf6_pls.istatus.eq(8).all() +def plot_output(idx, test): + name = test.name + gwf_ws = test.workspace + prt_ws = test.workspace / "prt" + gwf_name = get_model_name(name, "gwf") + prt_name = get_model_name(name, "prt") + gwf_sim = test.sims[0] + gwf = gwf_sim.get_model(gwf_name) + mg = gwf.modelgrid + drape = "drp" in name + + # check mf6 output files exist + gwf_head_file = f"{gwf_name}.hds" + prt_track_csv_file = f"{prt_name}.trk.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) + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + + # set up plot + fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 10)) + ax.set_aspect("equal") + + # plot mf6 pathlines in map view + pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) + 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=f"MF6 pathlines{' (drape)' if drape else ''}", + kind="line", + x="x", + y="y", + ax=ax, + legend=False, + color=cm.plasma(ipl / len(mf6_plines)), + ) + + # view/save plot + plt.show() + plt.savefig(prt_ws / f"{name}.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, ) diff --git a/autotest/test_prt_exg.py b/autotest/test_prt_exg.py index 88a71d5cd78..8c7176fb803 100644 --- a/autotest/test_prt_exg.py +++ b/autotest/test_prt_exg.py @@ -194,6 +194,13 @@ def check_output(idx, test): mp7_pathline_file = f"{mp7_name}.mppth" assert (mp7_ws / mp7_pathline_file).is_file() + # load head, budget, and specific discharge results from GWF model + gwf = sim.get_model(gwf_name) + 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) + # load mp7 pathline results plf = PathlineFile(mp7_ws / mp7_pathline_file) mp7_pls = pd.DataFrame( @@ -234,58 +241,6 @@ def check_output(idx, test): track_csv=gwf_ws / prt_track_csv_file, ) - # extract head, budget, and specific discharge results from GWF model - gwf = sim.get_model(gwf_name) - 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_{name}.png") - # convert mf6 pathlines to mp7 format mf6_pls = to_mp7_pathlines(mf6_pls) @@ -310,13 +265,103 @@ 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 = Path(test.workspace) + mp7_ws = gwf_ws / "mp7" + + # model names + gwf_name = get_model_name(idx, "gwf") + prt_name = get_model_name(idx, "prt") + mp7_name = get_model_name(idx, "mp7") + + # extract model objects + sim = test.sims[0] + gwf = sim.get_model(gwf_name) + + # extract model grid + mg = gwf.modelgrid + + gwf_budget_file = f"{gwf_name}.bud" + gwf_head_file = f"{gwf_name}.hds" + prt_track_file = f"{prt_name}.trk" + prt_track_csv_file = f"{prt_name}.trk.csv" + mp7_pathline_file = f"{mp7_name}.mppth" + + # load head, budget, and specific discharge results from GWF model + gwf = sim.get_model(gwf_name) + 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) + + # 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 + mp7_pls["particlegroup"] = mp7_pls["particlegroup"] + 1 + mp7_pls["node"] = mp7_pls["node"] + 1 + mp7_pls["k"] = mp7_pls["k"] + 1 + + # load mf6 pathline results + mf6_pls = pd.read_csv(gwf_ws / prt_track_csv_file).replace( + r"^\s*$", np.nan, regex=True + ) + + # setup 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(gwf_ws / f"{name}.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, ) diff --git a/autotest/test_prt_fmi.py b/autotest/test_prt_fmi.py index c92dacd383a..8ebe9f2b80d 100644 --- a/autotest/test_prt_fmi.py +++ b/autotest/test_prt_fmi.py @@ -227,16 +227,13 @@ 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() + + # 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) # load mp7 pathline results plf = PathlineFile(mp7_ws / mp7_pathline_file) @@ -251,10 +248,22 @@ def check_output(idx, test): # load mf6 pathline results mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + prt_budget_file = prt_ws / f"{prt_name}.bud" + prt_bud = flopy.utils.CellBudgetFile(prt_budget_file, precision="double") + prt_bud_data = prt_bud.get_data(kstpkper=(0, 0)) + # make sure pathline df has "name" (boundname) column and default values assert "name" in mf6_pls assert has_default_boundnames(mf6_pls) + 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() + assert (mp7_ws / mp7_pathline_file).is_file() + # make sure all mf6 pathline data have correct model and PRP index (1) assert all_equal(mf6_pls["imdl"], 1) assert all_equal(mf6_pls["iprp"], 1) @@ -267,9 +276,7 @@ def check_output(idx, test): ) # check cell-by-cell particle mass budget file - prt_budget_file = prt_ws / f"{prt_name}.bud" - prt_bud = flopy.utils.CellBudgetFile(prt_budget_file, precision="double") - prt_bud_data = prt_bud.get_data(kstpkper=(0, 0)) + assert len(prt_bud_data) == 2 assert prt_bud_data[0].shape == (1, 1, 460) assert prt_bud_data[1].shape == (9,) @@ -286,57 +293,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_data = False - if plot_data: - 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") - if "noext" in name: # maximum tracking time should be simulation stop time assert mf6_pls.t.max() == FlopyReadmeCase.nper * FlopyReadmeCase.perlen @@ -365,13 +321,97 @@ 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" + + # 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) + + # 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 + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + + prt_budget_file = prt_ws / f"{prt_name}.bud" + prt_bud = flopy.utils.CellBudgetFile(prt_budget_file, precision="double") + + # 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"{name}.png") + + @pytest.mark.parametrize("idx, name", enumerate(cases)) -def test_mf6model(idx, name, function_tmpdir, targets, benchmark): +def test_mf6model(idx, name, function_tmpdir, targets, benchmark, 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, xfail=[False, "bprp" in name, False], diff --git a/autotest/test_prt_quad_refinement.py b/autotest/test_prt_quad_refinement.py index 8813875c9cb..1e372ce4be4 100644 --- a/autotest/test_prt_quad_refinement.py +++ b/autotest/test_prt_quad_refinement.py @@ -181,24 +181,30 @@ def check_output(idx, test, snapshot): name = test.name gwf_ws = Path(test.workspace) mp7_ws = gwf_ws / "mp7" - - # model names gwf_name = get_model_name(idx, "gwf") prt_name = get_model_name(idx, "prt") - - # extract model objects sim = test.sims[0] gwf = sim.get_model(gwf_name) prt = sim.get_model(prt_name) - - # extract model grid mg = gwf.modelgrid - - # check mf6 output files exist gwf_budget_file = f"{gwf_name}.bud" gwf_head_file = f"{gwf_name}.hds" prt_track_file = f"{prt_name}.trk" prt_track_csv_file = f"{prt_name}.trk.csv" + + # load mf6 pathline results + mf6_pls = pd.read_csv(gwf_ws / prt_track_csv_file).replace( + r"^\s*$", np.nan, regex=True + ) + + # extract head, budget, and specific discharge results from GWF model + gwf = sim.get_model(gwf_name) + 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) + + # check mf6 output files exist assert (gwf_ws / gwf_budget_file).is_file() assert (gwf_ws / gwf_head_file).is_file() assert (gwf_ws / prt_track_file).is_file() @@ -218,15 +224,27 @@ def check_output(idx, test, snapshot): track_csv=gwf_ws / prt_track_csv_file, ) + # extract endpoints and compare to snapshot + mf6_eps = mf6_pls[mf6_pls.ireason == 3] + assert snapshot == mf6_eps.round(2).to_records(index=False) + + +def plot_output(idx, test): + name = test.name + gwf_ws = Path(test.workspace) + gwf_name = get_model_name(idx, "gwf") + prt_name = get_model_name(idx, "prt") + sim = test.sims[0] + gwf = sim.get_model(gwf_name) + mg = gwf.modelgrid + gwf_head_file = f"{gwf_name}.hds" + prt_track_csv_file = f"{prt_name}.trk.csv" + # load mf6 pathline results mf6_pls = pd.read_csv(gwf_ws / prt_track_csv_file).replace( r"^\s*$", np.nan, regex=True ) - # extract endpoints and compare to snapshot - mf6_eps = mf6_pls[mf6_pls.ireason == 3] - assert snapshot == mf6_eps.round(2).to_records(index=False) - # extract head, budget, and specific discharge results from GWF model gwf = sim.get_model(gwf_name) hds = HeadFile(gwf_ws / gwf_head_file).get_data() @@ -234,44 +252,44 @@ def check_output(idx, test, snapshot): 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=1, figsize=(10, 10)) - ax.set_aspect("equal") - - # plot mf6 pathlines in map view - pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) - 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, - legend=False, - color=cm.plasma(ipl / len(mf6_plines)), - ) - - # plot nodes - xc, yc = (mg.get_xcellcenters_for_layer(0), mg.get_ycellcenters_for_layer(0)) - for i in range(mg.ncpl): - x, y = xc[i], yc[i] - ax.annotate(str(i + 1), (x, y), color="grey", alpha=0.5) - - # view/save plot - plt.show() - plt.savefig(gwf_ws / f"test_{name}.png") + # set up plot + fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 10)) + ax.set_aspect("equal") + + # plot mf6 pathlines in map view + pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) + 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, + legend=False, + color=cm.plasma(ipl / len(mf6_plines)), + ) + + # plot nodes + xc, yc = (mg.get_xcellcenters_for_layer(0), mg.get_ycellcenters_for_layer(0)) + for i in range(mg.ncpl): + x, y = xc[i], yc[i] + ax.annotate(str(i + 1), (x, y), color="grey", alpha=0.5) + + # view/save plot + plt.show() + plt.savefig(gwf_ws / f"test_{name}.png") @pytest.mark.parametrize("idx, name", enumerate(cases)) @pytest.mark.parametrize("levels", [1, 2]) @pytest.mark.parametrize("method", ["pollock", "ternary"]) -def test_mf6model(idx, name, function_tmpdir, targets, levels, method, array_snapshot): +def test_mf6model( + idx, name, function_tmpdir, targets, levels, method, array_snapshot, plot +): test = TestFramework( name=name, workspace=function_tmpdir, @@ -284,6 +302,7 @@ def test_mf6model(idx, name, function_tmpdir, targets, levels, method, array_sna smoothing_level_horizontal=levels, ), check=lambda t: check_output(idx, t, array_snapshot), + plot=lambda t: plot_output(idx, t) if plot else None, targets=targets, compare=None, ) diff --git a/autotest/test_prt_release_timing.py b/autotest/test_prt_release_timing.py index 6fb4a49f528..e343f0e1bfd 100644 --- a/autotest/test_prt_release_timing.py +++ b/autotest/test_prt_release_timing.py @@ -262,56 +262,11 @@ def build_models(test): return gwf_sim, prt_sim, mp7_sim -def plot_output(grid, head, spdis, mf6_pls, mp7_pls, fpath): - 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=grid, ax=ax[0]) - pmv.plot_grid() - pmv.plot_array(head[0], alpha=0.1) - pmv.plot_vector(*spdis[:2], 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=grid, ax=ax[1]) - pmv.plot_grid() - pmv.plot_array(head[0], alpha=0.1) - pmv.plot_vector(*spdis[:2], 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(fpath) - - def check_output(test, snapshot): from flopy.plot.plotutil import to_mp7_pathlines name = test.name ws = test.workspace - prt_ws = test.workspace / "prt" mp7_ws = test.workspace / "mp7" gwf_name = get_model_name(name, "gwf") @@ -321,22 +276,43 @@ def check_output(test, snapshot): gwf = gwf_sim.get_model(gwf_name) mg = gwf.modelgrid - # check mf6 output files exist gwf_budget_file = f"{gwf_name}.bud" 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" + mp7_pathline_file = f"{mp7_name}.mppth" + + # load head, budget and intercell flows from gwf model + head = gwf.output.head().get_data() + bud = gwf.output.budget() + spdis = bud.get_data(text="DATA-SPDIS")[0] + qx, qy, _ = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) + + # 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 + + # apply reference time to mp7 results (mp7 reports relative times) + mp7_pls["time"] = mp7_pls["time"] + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + + # check output files exist assert (ws / gwf_budget_file).is_file() assert (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() # check list file for logged release configuration @@ -351,22 +327,6 @@ def check_output(test, snapshot): li = lines.index("PARTICLE RELEASE FOR PRP 1") assert "RELEASE SCHEDULE:" in lines[li + 1] - # 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 - - # apply reference time to mp7 results (mp7 reports relative times) - mp7_pls["time"] = mp7_pls["time"] - - # load mf6 pathline results - mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) - # make sure pathline df has "name" (boundname) column and empty values assert "name" in mf6_pls assert (mf6_pls["name"] == "").all() @@ -394,24 +354,6 @@ def check_output(test, snapshot): track_csv=track_csv, ) - # load head, budget and intercell flows from gwf model - head = gwf.output.head().get_data() - bud = gwf.output.budget() - spdis = bud.get_data(text="DATA-SPDIS")[0] - qx, qy, _ = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) - - # plot results if enabled - plot = False - if plot: - plot_output( - grid=gwf.modelgrid, - head=head, - spdis=(qx, qy), - mf6_pls=mf6_pls, - mp7_pls=mp7_pls, - fpath=ws / f"test_{simname}.png", - ) - # compare pathlines with snapshot assert snapshot == mf6_pls.drop("name", axis=1).round(3).to_records(index=False) @@ -457,14 +399,95 @@ def check_output(test, snapshot): assert np.allclose(mf6_pls, mp7_pls, atol=1e-3) +def plot_output(test): + name = test.name + 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 + + prt_track_csv_file = f"{prt_name}.trk.csv" + mp7_pathline_file = f"{mp7_name}.mppth" + + # load head, budget and intercell flows from gwf model + head = gwf.output.head().get_data() + bud = gwf.output.budget() + spdis = bud.get_data(text="DATA-SPDIS")[0] + qx, qy, _ = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) + + # 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 + + # apply reference time to mp7 results (mp7 reports relative times) + mp7_pls["time"] = mp7_pls["time"] + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + + # 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(head[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(head[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"{name}.png") + + @requires_pkg("syrupy") @pytest.mark.parametrize("name", cases) -def test_mf6model(name, function_tmpdir, targets, array_snapshot): +def test_mf6model(name, function_tmpdir, targets, array_snapshot, plot): test = TestFramework( name=name, workspace=function_tmpdir, build=lambda t: build_models(t), check=lambda t: check_output(t, array_snapshot), + plot=lambda t: plot_output(t) if plot else None, targets=targets, compare=None, # expect case using FRACTION to fail diff --git a/autotest/test_prt_stop_zones.py b/autotest/test_prt_stop_zones.py index 158e2625d16..0bc8d0d84e6 100644 --- a/autotest/test_prt_stop_zones.py +++ b/autotest/test_prt_stop_zones.py @@ -229,19 +229,20 @@ def check_output(idx, test): gwf = gwf_sim.get_model(gwf_name) mg = gwf.modelgrid - # check mf6 output files exist gwf_budget_file = f"{gwf_name}.bud" gwf_head_file = f"{gwf_name}.hds" prt_track_file = f"{prt_name}.trk" prt_track_csv_file = f"{prt_name}.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() - - # check mp7 output files exist mp7_pathline_file = f"{mp7_name}.mppth" - assert (mp7_ws / mp7_pathline_file).is_file() + + # get head, budget, and spdis 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) + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file) # load mp7 pathline results plf = PathlineFile(mp7_ws / mp7_pathline_file) @@ -253,8 +254,12 @@ 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) + # check output files exist + 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 (mp7_ws / mp7_pathline_file).is_file() # check budget data were written to mf6 prt list file check_budget_data( @@ -270,101 +275,6 @@ def check_output(idx, test): track_csv=prt_ws / prt_track_csv_file, ) - # get head, budget, and spdis 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 map view 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)), - ) - - def sort_square_verts(verts): - """Sort 4 or more points on a square in clockwise order, - starting with the top-left point - """ - - # sort by y coordinate - verts.sort(key=lambda v: v[1], reverse=True) - - # separate top and bottom rows - y0 = verts[0][1] - t = [v for v in verts if v[1] == y0] - b = verts[len(t) :] - - # sort top and bottom rows by x coordinate - t.sort(key=lambda v: v[0]) - b.sort(key=lambda v: v[0]) - - # return vertices in clockwise order - return t + list(reversed(b)) - - def plot_stop_zone(nn, ax): - ifaces = [] - iverts = mg.iverts[nn] - - # sort vertices of well cell in clockwise order - verts = [tuple(mg.verts[v]) for v in iverts] - sorted_verts = sort_square_verts(list(set(verts.copy()))) - for i in range(len(sorted_verts) - 1): - if i == 0: - p0 = sorted_verts[-1] - p1 = sorted_verts[i] - ifaces.append([p0, p1]) - p0 = sorted_verts[i] - p1 = sorted_verts[(i + 1)] - ifaces.append([p0, p1]) - - lc = LineCollection(ifaces, color="red", lw=4) - ax.add_collection(lc) - - # plot stop zones - for iz in stopzone_cells: - for a in ax: - plot_stop_zone(mg.get_node([iz])[0], a) - - # view/save plot - plt.show() - plt.savefig(gwf_ws / f"test_{name}_map.png") - # check that cell numbers are correct for i, row in list(mf6_pls.iterrows()): # todo debug final cell number disagreement @@ -405,13 +315,137 @@ def plot_stop_zone(nn, ax): 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 + + gwf_head_file = f"{gwf_name}.hds" + prt_track_csv_file = f"{prt_name}.trk.csv" + mp7_pathline_file = f"{mp7_name}.mppth" + + # get head, budget, and spdis 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) + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file) + + # 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 + mp7_pls["particlegroup"] = mp7_pls["particlegroup"] + 1 + mp7_pls["node"] = mp7_pls["node"] + 1 + mp7_pls["k"] = mp7_pls["k"] + 1 + + # 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)), + ) + + def sort_square_verts(verts): + """Sort 4 or more points on a square in clockwise order, + starting with the top-left point + """ + + # sort by y coordinate + verts.sort(key=lambda v: v[1], reverse=True) + + # separate top and bottom rows + y0 = verts[0][1] + t = [v for v in verts if v[1] == y0] + b = verts[len(t) :] + + # sort top and bottom rows by x coordinate + t.sort(key=lambda v: v[0]) + b.sort(key=lambda v: v[0]) + + # return vertices in clockwise order + return t + list(reversed(b)) + + def plot_stop_zone(nn, ax): + ifaces = [] + iverts = mg.iverts[nn] + + # sort vertices of well cell in clockwise order + verts = [tuple(mg.verts[v]) for v in iverts] + sorted_verts = sort_square_verts(list(set(verts.copy()))) + for i in range(len(sorted_verts) - 1): + if i == 0: + p0 = sorted_verts[-1] + p1 = sorted_verts[i] + ifaces.append([p0, p1]) + p0 = sorted_verts[i] + p1 = sorted_verts[(i + 1)] + ifaces.append([p0, p1]) + + lc = LineCollection(ifaces, color="red", lw=4) + ax.add_collection(lc) + + # plot stop zones + for iz in stopzone_cells: + for a in ax: + plot_stop_zone(mg.get_node([iz])[0], a) + + # view/save plot + plt.show() + plt.savefig(prt_ws / f"{name}.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, ) diff --git a/autotest/test_prt_ternary_methods.py b/autotest/test_prt_ternary_methods.py index ef635917031..5d7ce756527 100644 --- a/autotest/test_prt_ternary_methods.py +++ b/autotest/test_prt_ternary_methods.py @@ -132,7 +132,27 @@ def build_models(idx, test, exit_solve_tolerance=1e-7): return gwf_sim, prt_sim -def plot_output(name, gwf, head, spdis, pls, fpath): +def plot_output(idx, test): + name = test.name + prt_ws = test.workspace / "prt" + gwf_name = get_model_name(name, "gwf") + prt_name = get_model_name(name, "prt") + gwf_sim = test.sims[0] + gwf = gwf_sim.get_model(gwf_name) + + # get gwf output + gwf = gwf_sim.get_model() + head = gwf.output.head().get_data() + bdobj = gwf.output.budget() + spdis = bdobj.get_data(text="DATA-SPDIS")[0] + qx, qy, _ = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) + + # get prt output + prt_name = get_model_name(name, "prt") + prt_track_csv_file = f"{prt_name}.prp.trk.csv" + pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + endpts = pls.sort_values("t").groupby(["imdl", "iprp", "irpt", "trelease"]).tail(1) + # plot in 2d with mpl fig = plt.figure(figsize=(16, 10)) ax = plt.subplot(1, 1, 1, aspect="equal") @@ -161,7 +181,7 @@ def plot_output(name, gwf, head, spdis, pls, fpath): ) ) ax.legend(handles=handles, loc="lower right") - pmv.plot_vector(*spdis, normalize=True, alpha=0.25) + pmv.plot_vector(qx, qy, normalize=True, alpha=0.25) if "wel" in name: pmv.plot_bc(ftype="WEL") mf6_plines = pls.groupby(["iprp", "irpt", "trelease"]) @@ -199,7 +219,7 @@ def plot_output(name, gwf, head, spdis, pls, fpath): ax.annotate(str(i + 1), (x, y), color="grey", alpha=0.5) plt.show() - plt.savefig(fpath) + plt.savefig(prt_ws / f"{name}.png") def check_output(idx, test, snapshot): @@ -228,22 +248,18 @@ def check_output(idx, test, snapshot): assert endpts.shape == (2, 16) assert set(endpts.icell) == {111, 112} - # plot results if enabled - plot = False - if plot: - plot_output(name, gwf, head, (qx, qy), pls, fpath=prt_ws / f"{name}.png") - # check pathlines against snapshot assert snapshot == pls.drop("name", axis=1).round(3).to_records(index=False) @pytest.mark.parametrize("idx, name", enumerate(cases)) -def test_mf6model(idx, name, function_tmpdir, targets, benchmark, array_snapshot): +def test_mf6model(idx, name, function_tmpdir, targets, benchmark, array_snapshot, plot): test = TestFramework( name=name, workspace=function_tmpdir, build=lambda t: build_models(idx, t), check=lambda t: check_output(idx, t, array_snapshot), + plot=lambda t: plot_output(idx, t) if plot else None, targets=targets, compare=None, ) diff --git a/autotest/test_prt_track_events.py b/autotest/test_prt_track_events.py index 4de9bd5f8f2..c83dc83b175 100644 --- a/autotest/test_prt_track_events.py +++ b/autotest/test_prt_track_events.py @@ -321,19 +321,20 @@ def check_output(idx, test): gwf = gwf_sim.get_model(gwf_name) mg = gwf.modelgrid - # check mf6 output files exist gwf_budget_file = f"{gwf_name}.bud" gwf_head_file = f"{gwf_name}.hds" prt_track_file = f"{prt_name}.trk" prt_track_csv_file = f"{prt_name}.trk.csv" mp7_pathline_file = f"{mp7_name}.mppth" - 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() - # check mp7 output files exist - assert (mp7_ws / mp7_pathline_file).is_file() + # get head, budget, and spdis 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) + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file) # load mp7 pathline results plf = PathlineFile(mp7_ws / mp7_pathline_file) @@ -345,8 +346,12 @@ 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) + # check output files exist + 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 (mp7_ws / mp7_pathline_file).is_file() # check pathlines total size expected_len = 0 @@ -396,64 +401,6 @@ def all_equal(col, val): track_csv=prt_ws / prt_track_csv_file, ) - # check that particle names are particle indices - # assert len(mf6_pldata) == len( - # mf6_pldata[mf6_pldata["irpt"].astype(str).eq(mf6_pldata["name"])] - # ) - - # get head, budget, and spdis 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", - marker="o", - markersize=2, - linestyle="None", - 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") - # check that cell numbers are correct for i, row in list(mf6_pls.iterrows()): # todo debug final cell number disagreement @@ -497,13 +444,97 @@ def all_equal(col, val): 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 + + gwf_budget_file = f"{gwf_name}.bud" + gwf_head_file = f"{gwf_name}.hds" + prt_track_file = f"{prt_name}.trk" + prt_track_csv_file = f"{prt_name}.trk.csv" + mp7_pathline_file = f"{mp7_name}.mppth" + + # get head, budget, and spdis 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) + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file) + + # 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 + mp7_pls["particlegroup"] = mp7_pls["particlegroup"] + 1 + mp7_pls["node"] = mp7_pls["node"] + 1 + mp7_pls["k"] = mp7_pls["k"] + 1 + + # set up plots + 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", + marker="o", + markersize=2, + linestyle="None", + 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"{name}.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, ) diff --git a/autotest/test_prt_triangle.py b/autotest/test_prt_triangle.py index 9186f3a0537..2f7bdcc3e22 100644 --- a/autotest/test_prt_triangle.py +++ b/autotest/test_prt_triangle.py @@ -225,13 +225,35 @@ def build_models(idx, test): return gwf_sim, prt_sim -def plot_output(name, grid, head, spdis, pls): +def plot_output(idx, test): + name = test.name + prt_ws = test.workspace / "prt" + gwf_name = get_model_name(name, "gwf") + prt_name = get_model_name(name, "prt") + gwf_sim = test.sims[0] + gwf = gwf_sim.get_model(gwf_name) + grid = gwf.modelgrid + + # get gwf output + gwf = gwf_sim.get_model() + head = gwf.output.head().get_data() + bdobj = gwf.output.budget() + spdis = bdobj.get_data(text="DATA-SPDIS")[0] + qx, qy, _ = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) + + # get prt output + prt_name = get_model_name(name, "prt") + prt_track_csv_file = f"{prt_name}.trk.csv" + pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + endpts = pls[pls.ireason == 3] # termination + + # set up plot fig = plt.figure(figsize=(10, 10)) ax = plt.subplot(1, 1, 1, aspect="equal") - pmv = flopy.plot.PlotMapView(model=grid, ax=ax) + pmv = flopy.plot.PlotMapView(model=gwf, ax=ax) pmv.plot_grid() pmv.plot_array(head, cmap="Blues", alpha=0.25) - pmv.plot_vector(*spdis, normalize=True, alpha=0.25) + pmv.plot_vector(qx, qy, normalize=True, alpha=0.25) mf6_plines = pls.groupby(["iprp", "irpt", "trelease"]) for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): pl.plot( @@ -248,7 +270,9 @@ def plot_output(name, grid, head, spdis, pls): x, y = xc[i], yc[i] ax.plot(x, y, "o", color="grey", alpha=0.25, ms=2) ax.annotate(str(i + 1), (x, y), color="grey", alpha=0.5) + plt.show() + plt.savefig(prt_ws / f"{name}.png") def check_output(idx, test, snapshot): @@ -275,10 +299,6 @@ def check_output(idx, test, snapshot): # check termination points against snapshot assert snapshot == endpts.drop("name", axis=1).round(3).to_records(index=False) - plot_debug = False - if plot_debug: - plot_output(name, gwf.modelgrid, head, (qx, qy), pls) - if "r2l" in name: assert pls.shape == (164, 16) assert (pls.z == 0.5).all() @@ -309,12 +329,13 @@ def check_output(idx, test, snapshot): @pytest.mark.parametrize("idx, name", enumerate(cases)) -def test_mf6model(idx, name, function_tmpdir, targets, array_snapshot): +def test_mf6model(idx, name, function_tmpdir, targets, array_snapshot, plot): test = TestFramework( name=name, workspace=function_tmpdir, build=lambda t: build_models(idx, t), check=lambda t: check_output(idx, t, array_snapshot), + plot=lambda t: plot_output(idx, t) if plot else None, targets=targets, compare=None, ) diff --git a/autotest/test_prt_voronoi1.py b/autotest/test_prt_voronoi1.py index 2bdb18cfb23..16e0ebf1d65 100644 --- a/autotest/test_prt_voronoi1.py +++ b/autotest/test_prt_voronoi1.py @@ -24,7 +24,6 @@ from flopy.utils.voronoi import VoronoiGrid from framework import TestFramework from modflow_devtools.markers import requires_pkg -from modflow_devtools.misc import is_in_ci from prt_test_utils import get_model_name from shapely.geometry import LineString, Point @@ -245,7 +244,24 @@ def build_models(idx, test): return gwf_sim, prt_sim -def plot_output(name, gwf, head, spdis, pls, fpath): +def plot_output(idx, test): + name = test.name + prt_ws = test.workspace / "prt" + prt_name = get_model_name(name, "prt") + gwfsim = test.sims[0] + + # get gwf output + gwf = gwfsim.get_model() + head = gwf.output.head().get_data() + bdobj = gwf.output.budget() + spdis = bdobj.get_data(text="DATA-SPDIS")[0] + qx, qy, _ = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) + + # get prt output + prt_track_csv_file = f"{prt_name}.prp.trk.csv" + pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + endpts = pls[pls.ireason == 3] # termination + # plot in 2d with mpl fig = plt.figure(figsize=(16, 10)) ax = plt.subplot(1, 1, 1, aspect="equal") @@ -274,7 +290,7 @@ def plot_output(name, gwf, head, spdis, pls, fpath): ) ) ax.legend(handles=handles, loc="lower right") - pmv.plot_vector(*spdis, normalize=True, alpha=0.25) + pmv.plot_vector(qx, qy, normalize=True, alpha=0.25) if "wel" in name: pmv.plot_bc(ftype="WEL") mf6_plines = pls.groupby(["iprp", "irpt", "trelease"]) @@ -312,40 +328,7 @@ def plot_output(name, gwf, head, spdis, pls, fpath): ax.annotate(str(i + 1), (x, y), color="grey", alpha=0.5) plt.show() - plt.savefig(fpath) - - # plot in 3d with pyvista (via vtk) - import pyvista as pv - from flopy.export.vtk import Vtk - from flopy.plot.plotutil import to_mp7_pathlines - - def get_meshes(model, pathlines): - vtk = Vtk(model=model, binary=False, smooth=False) - vtk.add_model(model) - vtk.add_pathline_points(to_mp7_pathlines(pathlines.to_records(index=False))) - grid_mesh, path_mesh = vtk.to_pyvista() - grid_mesh.rotate_x(-100, point=axes.origin, inplace=True) - grid_mesh.rotate_z(90, point=axes.origin, inplace=True) - grid_mesh.rotate_y(120, point=axes.origin, inplace=True) - path_mesh.rotate_x(-100, point=axes.origin, inplace=True) - path_mesh.rotate_z(90, point=axes.origin, inplace=True) - path_mesh.rotate_y(120, point=axes.origin, inplace=True) - return grid_mesh, path_mesh - - def callback(mesh, value): - sub = pls[pls.t <= value] - gm, pm = get_meshes(gwf, sub) - mesh.shallow_copy(pm) - - pv.set_plot_theme("document") - axes = pv.Axes(show_actor=True, actor_scale=2.0, line_width=5) - p = pv.Plotter(notebook=False) - grid_mesh, path_mesh = get_meshes(gwf, pls) - p.add_mesh(grid_mesh, scalars=head[0], cmap="Blues", opacity=0.5) - p.add_mesh(path_mesh, label="Time", style="points", color="black") - p.camera.zoom(1) - p.add_slider_widget(lambda v: callback(path_mesh, v), [0, 30202]) - p.show() + plt.savefig(prt_ws / f"{name}.png") def check_output(idx, test, snapshot): @@ -369,24 +352,20 @@ def check_output(idx, test, snapshot): # compare pathlines with snapshot. particles shouldn't # have moved vertically. round for cross-platform error. # skip macos-14 in CI because grid is slightly different - if not (is_in_ci() and system() == "Darwin" and processor() == "arm"): + if not (system() == "Darwin" and processor() == "arm"): assert snapshot == endpts.drop("name", axis=1).round(1).to_records(index=False) - # plot results if enabled - plot = False - if plot: - plot_output(name, gwf, head, (qx, qy), pls, fpath=prt_ws / f"{name}.png") - @requires_pkg("syrupy") @pytest.mark.slow @pytest.mark.parametrize("idx, name", enumerate(cases)) -def test_mf6model(idx, name, function_tmpdir, targets, benchmark, array_snapshot): +def test_mf6model(idx, name, function_tmpdir, targets, benchmark, array_snapshot, plot): test = TestFramework( name=name, workspace=function_tmpdir, build=lambda t: build_models(idx, t), check=lambda t: check_output(idx, t, array_snapshot), + plot=lambda t: plot_output(idx, t) if plot else None, targets=targets, compare=None, ) diff --git a/autotest/test_prt_voronoi2.py b/autotest/test_prt_voronoi2.py index 30cc846b926..5ec99a808d3 100644 --- a/autotest/test_prt_voronoi2.py +++ b/autotest/test_prt_voronoi2.py @@ -227,96 +227,89 @@ def check_output(idx, test): prt_track_csv_file = f"{prt_name}.prp.trk.csv" pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) - plot_2d = False - if plot_2d: - # plot in 2d with mpl - fig = plt.figure(figsize=(16, 10)) - ax = plt.subplot(1, 1, 1, aspect="equal") - pmv = flopy.plot.PlotMapView(model=gwf, ax=ax) - pmv.plot_grid(alpha=0.25) - pmv.plot_ibound(alpha=0.5) - # headmesh = pmv.plot_array(head, alpha=0.25) - # headctr = pmv.contour_array(head, levels=np.linspace(0, 1, 9), colors="black") - # plt.clabel(headctr) - # plt.colorbar(headmesh, shrink=0.25, ax=ax, label="Head", location="right") - concmesh = pmv.plot_array(conc, cmap="jet") - concctr = pmv.contour_array(conc, levels=(0.0001, 0.001, 0.01, 0.1), colors="y") - plt.clabel(concctr) - plt.colorbar( - concmesh, shrink=0.25, ax=ax, label="Concentration", location="right" - ) + # pathlines should be west -> east and no elevation change + assert len(pls.irpt.unique()) == 3 + assert np.allclose(pls.y, 967.556473) + assert pls.x.min() < 250 + assert pls.x.max() > 1950 + assert np.allclose(pls[pls.irpt == 1].z, 0.166666) + assert np.allclose(pls[pls.irpt == 2].z, 0.5) + assert np.allclose(pls[pls.irpt == 3].z, 0.833333) + + +def plot_output(idx, test): + name = test.name + prt_ws = test.workspace / "prt" + prt_name = get_model_name(name, "prt") + gwfsim, gwtsim, prtsim = test.sims + + # get gwf output + gwf = gwfsim.get_model() + head = gwf.output.head().get_data() + bdobj = gwf.output.budget() + spdis = bdobj.get_data(text="DATA-SPDIS")[0] + qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) + + # get gwt output + gwt = gwtsim.get_model() + conc = gwt.output.concentration().get_data() - handles = [ - mpl.lines.Line2D( - [0], - [0], - marker=">", - linestyle="", - label="Specific discharge", - color="grey", - markerfacecolor="gray", - ), - ] - ax.legend(handles=handles, loc="lower right") - pmv.plot_vector(qx, qy, normalize=True, alpha=0.25) - mf6_plines = pls.groupby(["iprp", "irpt", "trelease"]) - for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): - title = "DISV voronoi grid particle tracks" - pl.plot( - title=title, - kind="line", - x="x", - y="y", - ax=ax, - legend=False, - color="black", - ) - plt.show() - plt.savefig(prt_ws / f"{name}.png") - - plot_3d = False - if plot_3d: - # plot in 3d with pyvista (via vtk) - import pyvista as pv - from flopy.export.vtk import Vtk - from flopy.plot.plotutil import to_mp7_pathlines - - def get_meshes(model, pathlines): - vtk = Vtk(model=model, binary=False, smooth=False) - vtk.add_model(model) - vtk.add_pathline_points(to_mp7_pathlines(pathlines.to_records(index=False))) - grid_mesh, path_mesh = vtk.to_pyvista() - grid_mesh.rotate_x(-100, point=axes.origin, inplace=True) - grid_mesh.rotate_z(90, point=axes.origin, inplace=True) - grid_mesh.rotate_y(120, point=axes.origin, inplace=True) - path_mesh.rotate_x(-100, point=axes.origin, inplace=True) - path_mesh.rotate_z(90, point=axes.origin, inplace=True) - path_mesh.rotate_y(120, point=axes.origin, inplace=True) - return grid_mesh, path_mesh - - def callback(mesh, value): - sub = pls[pls.t <= value] - gm, pm = get_meshes(gwf, sub) - mesh.shallow_copy(pm) - - pv.set_plot_theme("document") - axes = pv.Axes(show_actor=True, actor_scale=2.0, line_width=5) - p = pv.Plotter(notebook=False) - grid_mesh, path_mesh = get_meshes(gwf, pls) - p.add_mesh(grid_mesh, scalars=head[0], cmap="Blues", opacity=0.5) - p.add_mesh(path_mesh, label="Time", style="points", color="black") - p.camera.zoom(1) - p.add_slider_widget(lambda v: callback(path_mesh, v), [0, 30202]) - p.show() + # get prt output + prt_track_csv_file = f"{prt_name}.prp.trk.csv" + pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) + + # plot in 2d with mpl + fig = plt.figure(figsize=(16, 10)) + ax = plt.subplot(1, 1, 1, aspect="equal") + pmv = flopy.plot.PlotMapView(model=gwf, ax=ax) + pmv.plot_grid(alpha=0.25) + pmv.plot_ibound(alpha=0.5) + # headmesh = pmv.plot_array(head, alpha=0.25) + # headctr = pmv.contour_array(head, levels=np.linspace(0, 1, 9), colors="black") + # plt.clabel(headctr) + # plt.colorbar(headmesh, shrink=0.25, ax=ax, label="Head", location="right") + concmesh = pmv.plot_array(conc, cmap="jet") + concctr = pmv.contour_array(conc, levels=(0.0001, 0.001, 0.01, 0.1), colors="y") + plt.clabel(concctr) + plt.colorbar(concmesh, shrink=0.25, ax=ax, label="Concentration", location="right") + + handles = [ + mpl.lines.Line2D( + [0], + [0], + marker=">", + linestyle="", + label="Specific discharge", + color="grey", + markerfacecolor="gray", + ), + ] + ax.legend(handles=handles, loc="lower right") + pmv.plot_vector(qx, qy, normalize=True, alpha=0.25) + mf6_plines = pls.groupby(["iprp", "irpt", "trelease"]) + for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): + title = "DISV voronoi grid particle tracks" + pl.plot( + title=title, + kind="line", + x="x", + y="y", + ax=ax, + legend=False, + color="black", + ) + plt.show() + plt.savefig(prt_ws / f"{name}.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, ) diff --git a/autotest/test_prt_weak_sinks.py b/autotest/test_prt_weak_sinks.py index 545b82b9750..78b08ac5223 100644 --- a/autotest/test_prt_weak_sinks.py +++ b/autotest/test_prt_weak_sinks.py @@ -220,19 +220,20 @@ def check_output(idx, test): gwf = gwf_sim.get_model(gwf_name) mg = gwf.modelgrid - # check mf6 output files exist gwf_budget_file = f"{gwf_name}.bud" gwf_head_file = f"{gwf_name}.hds" prt_track_file = f"{prt_name}.trk" prt_track_csv_file = f"{prt_name}.trk.csv" mp7_pathline_file = f"{mp7_name}.mppth" - 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() - # check mp7 output files exist - assert (mp7_ws / mp7_pathline_file).is_file() + # 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) + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file) # load mp7 pathline results plf = PathlineFile(mp7_ws / mp7_pathline_file) @@ -244,8 +245,12 @@ 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) + # check output files exist + 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 (mp7_ws / mp7_pathline_file).is_file() # if STOP_AT_WEAK_SINK disabled, # check for an extra datum when particle exited weak sink @@ -278,70 +283,6 @@ def all_equal(col, val): track_csv=prt_ws / prt_track_csv_file, ) - # 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(model=gwf, ax=ax[0]) - pmv.plot_grid() - pmv.plot_array(hds[0], alpha=0.1) - pmv.plot_vector(qx, qy, normalize=True, color="white") - pmv.plot_bc("WEL") - 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(model=gwf, ax=ax[1]) - pmv.plot_grid() - pmv.plot_array(hds[0], alpha=0.1) - pmv.plot_vector(qx, qy, normalize=True, color="white") - pmv.plot_bc("WEL") - 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)), - ) - - # plot cell centers - # xc, yc = mg.get_xcellcenters_for_layer(0), mg.get_ycellcenters_for_layer(0) - # xc = xc.flatten() - # yc = yc.flatten() - # for i in range(mg.ncpl): - # x, y = xc[i], yc[i] - # nn = mg.get_node(mg.intersect(x, y, 0))[0] - # for a in ax: - # a.plot(x, y, "ro") - # a.annotate(str(nn + 1), (x, y), color="r") - - # 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) @@ -370,13 +311,97 @@ def all_equal(col, val): 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 + + gwf_budget_file = f"{gwf_name}.bud" + gwf_head_file = f"{gwf_name}.hds" + prt_track_file = f"{prt_name}.trk" + prt_track_csv_file = f"{prt_name}.trk.csv" + mp7_pathline_file = f"{mp7_name}.mppth" + + # 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) + + # load mf6 pathline results + mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file) + + # 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 + + # 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(model=gwf, ax=ax[0]) + pmv.plot_grid() + pmv.plot_array(hds[0], alpha=0.1) + pmv.plot_vector(qx, qy, normalize=True, color="white") + pmv.plot_bc("WEL") + 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(model=gwf, ax=ax[1]) + pmv.plot_grid() + pmv.plot_array(hds[0], alpha=0.1) + pmv.plot_vector(qx, qy, normalize=True, color="white") + pmv.plot_bc("WEL") + 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"{name}.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, )