better plots, cleanup

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()