diff --git a/autotest/prt/test_prt_fmi08.py b/autotest/prt/test_prt_fmi08.py index 574707b5417..5857d2d778f 100644 --- a/autotest/prt/test_prt_fmi08.py +++ b/autotest/prt/test_prt_fmi08.py @@ -52,9 +52,13 @@ nodes = ncol * nrow porosity = 0.1 rpts = [ + [500, 100, 0.5], + [500, 350, 0.5], [500, 450, 0.5], [500, 500, 0.5], [500, 550, 0.5], + [500, 650, 0.5], + [500, 800, 0.5], ] @@ -279,55 +283,61 @@ def test_mf6model(name, function_tmpdir, targets): prt_track_csv_file = f"{prtname}.prp.trk.csv" pls = pd.read_csv(ws / prt_track_csv_file, na_filter=False) - # plot in 2d - fig = plt.figure(figsize=(10, 10)) - ax = plt.subplot(1, 1, 1, aspect="equal") - pmv = flopy.plot.PlotMapView(model=gwf, ax=ax) - pmv.plot_grid() - pmv.plot_array(head, cmap="jet", 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"]) - for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): - pl.plot( - title=f"MF6 pathlines ({name})", - kind="line", - x="x", - y="y", - ax=ax, - legend=False, - color=cm.plasma(ipl / len(mf6_plines)), - ) - # plt.show() + plot_debug = True + if plot_debug: + # plot in 2d with mpl + fig = plt.figure(figsize=(10, 10)) + ax = plt.subplot(1, 1, 1, aspect="equal") + pmv = flopy.plot.PlotMapView(model=gwf, ax=ax) + pmv.plot_grid() + pmv.plot_array(head, cmap="Blues", 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"]) + for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): + pl.plot( + title=f"MF6 pathlines ({name})", + kind="line", + x="x", + y="y", + ax=ax, + legend=False, + color="black", + ) + plt.show() - if "wel" in name: # plot in 3d with pyvista (via vtk) from flopy.export.vtk import Vtk from flopy.plot.plotutil import to_mp7_pathlines - - vtk = Vtk(model=gwf, binary=False, smooth=False) - vtk.add_model(gwf) - vtk.add_pathline_points(to_mp7_pathlines(pls.to_records(index=False))) - grid_mesh, path_mesh = vtk.to_pyvista() import pyvista as pv + 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) - 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) - p = pv.Plotter() - p.add_mesh(grid_mesh, scalars=head[0], cmap="Blues") - p.add_mesh( - path_mesh, - scalars="time", - cmap="Reds", - label="Time", - style="points", - ) - p.camera.zoom(2) - # p.show() + 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()