10b: Particle tracking with MODPATH
In this exercise, we will use MODPATH to simulate advective transport with the Freyberg flow model.
MODPATH is a particle tracking software that simulates the advective transport of particles using the output from a MODFLOW groundwater flow model. The MODPATH user guide, documentation, and executable can be found here. The MODPATH executable is downloaded with the stack of modflow related programs when get_modflow() was run in setup for this course. Students should already have this program on their machines.
For this exercise, we will use a quadtree version of the Freyberg model.
[1]:
from IPython.display import clear_output, display
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
from flopy.utils.gridintersect import GridIntersect
import flopy
import pandas as pd
The location of the contamination patch and the nodes that the define bounding cells of the patch are calculated below.
[2]:
# patch upper left and lower right
xmin, xmax = 250. * 1, 250. * 3
ymin, ymax = (40 - 14) * 250., (40 - 11) * 250.
csx, csy = [xmin, xmin, xmax, xmax, xmin], [ymin, ymax, ymax, ymin, ymin]
polygon = [list(zip(csx, csy))]
(xmin, ymax), (xmax, ymin)
[2]:
((250.0, 7250.0), (750.0, 6500.0))
Define the workspace and model names.
[3]:
load_ws = Path('data/quadtree')
ws = Path("temp/ex10a")
name = "project"
name_mp = f"{name}_mp"
exe_name = 'mf6'
Load the MODFLOW 6 Model
Load a simulation object using flopy.mf6.MFSimulation().load().
[4]:
%%capture
sim = flopy.mf6.MFSimulation.load(sim_name=name, exe_name=exe_name,
sim_ws=load_ws)
Load the groundwater flow model
[5]:
gwf = sim.get_model(name)
gwf.modelgrid
[5]:
xll:0.0; yll:0.0; rotation:0.0; units:meters; lenuni:2
Change the workspace
[6]:
sim.set_sim_path(ws)
Write the model files
[7]:
%%capture
sim.write_simulation()
Run the simulation.
[8]:
sim.run_simulation()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[8], line 1
----> 1 sim.run_simulation()
File ~/miniforge3/envs/pyclass/lib/python3.12/site-packages/flopy/mf6/mfsimbase.py:1820, in MFSimulationBase.run_simulation(self, silent, pause, report, processors, normal_msg, use_async, cargs, custom_print)
1818 else:
1819 silent = True
-> 1820 return run_model(
1821 self.exe_name,
1822 None,
1823 self.simulation_data.mfpath.get_sim_path(),
1824 silent=silent,
1825 pause=pause,
1826 report=report,
1827 processors=processors,
1828 normal_msg=normal_msg,
1829 use_async=use_async,
1830 cargs=cargs,
1831 custom_print=custom_print,
1832 )
File ~/miniforge3/envs/pyclass/lib/python3.12/site-packages/flopy/mbase.py:1692, in run_model(exe_name, namefile, model_ws, silent, pause, report, processors, normal_msg, use_async, cargs, custom_print)
1690 if exe_name is None:
1691 raise ValueError("An executable name or path must be provided")
-> 1692 exe_path = resolve_exe(exe_name)
1693 if not silent:
1694 print(
1695 "FloPy is using the following executable to run the model: "
1696 + flopy_io.relpath_safe(exe_path, model_ws)
1697 )
File ~/miniforge3/envs/pyclass/lib/python3.12/site-packages/flopy/mbase.py:110, in resolve_exe(exe_name, forgive)
104 warn(
105 f"The program {exe_name} does not exist or is not executable.",
106 category=UserWarning,
107 )
108 return None
--> 110 raise FileNotFoundError(
111 f"The program {exe_name} does not exist or is not executable."
112 )
FileNotFoundError: The program mf6 does not exist or is not executable.
Create and Run the MODPATH model
Lets plot the model grid and the location of the contamination patch.
[9]:
fig, ax = plt.subplots(figsize=(5, 9))
mm = flopy.plot.PlotMapView(gwf, layer=0, ax=ax)
mm.plot_bc('SFR', color="b", plotAll=True)
mm.plot_bc('WEL', plotAll=True)
mm.plot_inactive(alpha=0.75)
mm.plot_grid(lw=0.25, color='grey')
ax.fill(csx, csy, color='#e534eb');
Find the node numbers of the contamination patch using the groundwater model grid object and the GridIntersect utility. The GridIntersect utility has an intersect method that can be used to identified cellids from points, lines, and polygons.
[10]:
gx = GridIntersect(gwf.modelgrid)
results = gx.intersect(polygon, 'Polygon')
nodes = results.cellids
nodes
[10]:
array([1352, 1353, 1459, 1460, 1581, 1582], dtype=object)
Create the MODPATH particle data
[11]:
# create
sd = flopy.modpath.CellDataType()
p = flopy.modpath.NodeParticleData(subdivisiondata=[sd],
nodes=list(nodes))
# create forward particle group
fpth = name_mp + '.sloc'
pg = flopy.modpath.ParticleGroupNodeTemplate(particlegroupname='contaminant',
particledata=p,
filename=fpth)
[ ]:
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[25], line 1
----> 1 sd.data
AttributeError: 'CellDataType' object has no attribute 'data'
Create the MODPATH 7 files
[15]:
# create modpath files
mp = flopy.modpath.Modpath7(modelname=name_mp, flowmodel=gwf,
exe_name='mp7', model_ws=ws)
mpbas = flopy.modpath.Modpath7Bas(mp, porosity=0.1)
mpsim = flopy.modpath.Modpath7Sim(mp, simulationtype='pathline',
trackingdirection='forward',
weaksinkoption='pass_through',
weaksourceoption='pass_through',
referencetime=0.,
stoptimeoption='extend',
particlegroups=pg)
help(flopy.modpath.Modpath7Sim)
Help on class Modpath7Sim in module flopy.modpath.mp7sim:
class Modpath7Sim(flopy.pakbase.Package)
| Modpath7Sim(model, mpnamefilename=None, listingfilename=None, endpointfilename=None, pathlinefilename=None, timeseriesfilename=None, tracefilename=None, simulationtype='pathline', trackingdirection='forward', weaksinkoption='stop_at', weaksourceoption='stop_at', budgetoutputoption='no', traceparticledata=None, budgetcellnumbers=None, referencetime=None, stoptimeoption='extend', stoptime=None, timepointdata=None, zonedataoption='off', stopzone=None, zones=0, retardationfactoroption='off', retardation=1.0, particlegroups=None, extension='mpsim')
|
| MODPATH Simulation File Package Class.
|
| Parameters
| ----------
| model : model object
| The model object (of type :class:`flopy.modpath.Modpath7`) to
| which this package will be added.
| mpnamefilename : str
| Filename of the MODPATH 7 name file. If mpnamefilename is not
| defined it will be generated from the model name
| (default is None).
| listingfilename : str
| Filename of the MODPATH 7 listing file. If listingfilename is not
| defined it will be generated from the model name
| (default is None).
| endpointfilename : str
| Filename of the MODPATH 7 endpoint file. If endpointfilename is
| not defined it will be generated from the model name
| (default is None).
| pathlinefilename : str
| Filename of the MODPATH 7 pathline file. If pathlinefilename is
| not defined it will be generated from the model name
| (default is None).
| timeseriesfilename : str
| Filename of the MODPATH 7 timeseries file. If timeseriesfilename
| is not defined it will be generated from the model name
| (default is None).
| tracefilename : str
| Filename of the MODPATH 7 tracefile file. If tracefilename is not
| defined it will be generated from the model name
| (default is None).
| simulationtype : str
| MODPATH 7 simulation type. Valid simulation types are 'endpoint',
| 'pathline', 'timeseries', or 'combined' (default is 'pathline').
| trackingdirection : str
| MODPATH 7 tracking direction. Valid tracking directions are
| 'forward' or 'backward' (default os 'forward').
| weaksinkoption : str
| MODPATH 7 weak sink option. Valid weak sink options are
| 'pass_through' or 'stop_at' (default value is 'stop_at').
| weaksourceoption : str
| MODPATH 7 weak source option. Valid weak source options are
| 'pass_through' or 'stop_at' (default value is 'stop_at').
| budgetoutputoption : str
| MODPATH 7 budget output option. Valid budget output options are
| 'no' - individual cell water balance errors are not computed
| and budget record headers are not printed, 'summary' - a summary
| of individual cell water balance errors for each time step is
| printed in the listing file without record headers, or
| 'record_summary' - a summary of individual cell water balance
| errors for each time step is printed in the listing file with
| record headers (default is 'summary').
| traceparticledata : list or tuple
| List or tuple with two ints that define the particle group and
| particle id (zero-based) of the specified particle that is
| followed in detail. If traceparticledata is None, trace mode is
| off (default is None).
| budgetcellnumbers : int, list of ints, tuple of ints, or np.ndarray
| Cell numbers (zero-based) for which detailed water budgets are
| computed. If budgetcellnumbers is None, detailed water budgets are
| not calculated (default is None).
| referencetime : float, list, or tuple
| Specified reference time if a float or a list/tuple with a single
| float value is provided (reference time option 1). Otherwise a
| list or tuple with a zero-based stress period (int) and time
| step (int) and a float defining the relative time position in the
| time step is provided (reference time option 2). If referencetime
| is None, reference time is set to 0 (default is None).
| stoptimeoption : str
| String indicating how a particle tracking simulation is
| terminated based on time. If stop time option is 'total', particles
| will be stopped at the end of the final time step if 'forward'
| tracking is simulated or at the beginning of the first time step
| if backward tracking. If stop time option is 'extend', initial or
| final steady-state time steps will be extended and all particles
| will be tracked until they reach a termination location. If stop
| time option is 'specified', particles will be tracked until they
| reach a termination location or the specified stop time is reached
| (default is 'extend').
| stoptime : float
| User-specified value of tracking time at which to stop a particle
| tracking simulation. Stop time is only used if the stop time option
| is 'specified'. If stoptime is None and the stop time option is
| 'specified' particles will be terminated at the end of the last
| time step if 'forward' tracking or the beginning of the first time
| step if 'backward' tracking (default is None).
| timepointdata : list or tuple
| List or tuple with 2 items that is only used if simulationtype is
| 'timeseries' or 'combined'. If the second item is a float then the
| timepoint data corresponds to time point option 1 and the first
| entry is the number of time points (timepointcount) and the second
| entry is the time point interval. If the second item is a list,
| tuple, or np.ndarray then the timepoint data corresponds to time
| point option 2 and the number of time points entries
| (timepointcount) in the second item and the second item is an
| list, tuple, or array of user-defined time points. If Timepointdata
| is None, time point option 1 is specified and the total simulation
| time is split into 100 intervals (default is None).
| zonedataoption : str
| If zonedataoption is 'off', zone array data are not read and a zone
| value of 1 is applied to all cells. If zonedataoption is 'on',
| zone array data are read (default is 'off').
| stopzone : int
| A zero-based specified integer zone value that indicates an
| automatic stopping location for particles and is only used if
| zonedataoption is 'on'. A value of -1 indicates no automatic stop
| zone is used. Stopzone values less than -1 are not allowed. If
| stopzone is None, stopzone is set to -1 (default is None).
| zones : float or array of floats (nlay, nrow, ncol)
| Array of zero-based positive integer zones that are only used if
| zonedataoption is 'on' (default is 0).
| retardationfactoroption : str
| If retardationfactoroption is 'off', retardation array data are not
| read and a retardation factor of 1 is applied to all cells. If
| retardationfactoroption is 'on', retardation factor array data are
| read (default is 'off').
| retardation : float or array of floats (nlay, nrow, ncol)
| Array of retardation factors that are only used if
| retardationfactoroption is 'on' (default is 1).
| particlegroups : ParticleGroup or list of ParticleGroups
| ParticleGroup or list of ParticlesGroups that contain data for
| individual particle groups. If None is specified, a
| particle in the center of node 0 will be created (default is None).
| extension : string
| Filename extension (default is 'mpsim')
|
| Examples
| --------
|
| >>> import flopy
| >>> m = flopy.modflow.Modflow.load('mf2005.nam')
| >>> mp = flopy.modpath.Modpath7('mf2005_mp', flowmodel=m)
| >>> mpsim = flopy.modpath.Modpath7Sim(mp)
|
| Method resolution order:
| Modpath7Sim
| flopy.pakbase.Package
| flopy.pakbase.PackageInterface
| builtins.object
|
| Methods defined here:
|
| __init__(self, model, mpnamefilename=None, listingfilename=None, endpointfilename=None, pathlinefilename=None, timeseriesfilename=None, tracefilename=None, simulationtype='pathline', trackingdirection='forward', weaksinkoption='stop_at', weaksourceoption='stop_at', budgetoutputoption='no', traceparticledata=None, budgetcellnumbers=None, referencetime=None, stoptimeoption='extend', stoptime=None, timepointdata=None, zonedataoption='off', stopzone=None, zones=0, retardationfactoroption='off', retardation=1.0, particlegroups=None, extension='mpsim')
| Initialize self. See help(type(self)) for accurate signature.
|
| write_file(self, check=False)
| Write the package file
|
| Parameters
| ----------
| check : boolean
| Check package data for common errors. (default False)
|
| Returns
| -------
| None
|
| ----------------------------------------------------------------------
| Methods inherited from flopy.pakbase.Package:
|
| __getitem__(self, item)
|
| __repr__(self)
| Return repr(self).
|
| __setattr__(self, key, value)
| Implement setattr(self, name, value).
|
| __setitem__(self, key, value)
|
| export(self, f, **kwargs)
| Method to export a package to netcdf or shapefile based on the
| extension of the file name (.shp for shapefile, .nc for netcdf)
|
| Parameters
| ----------
| f : str
| filename
| kwargs : keyword arguments
| modelgrid : flopy.discretization.Grid instance
| user supplied modelgrid which can be used for exporting
| in lieu of the modelgrid associated with the model object
|
| Returns
| -------
| None or Netcdf object
|
| level1_arraylist(self, idx, v, name, txt)
|
| plot(self, **kwargs)
| Plot 2-D, 3-D, transient 2-D, and stress period list (MfList)
| package input data
|
| Parameters
| ----------
| **kwargs : dict
| filename_base : str
| Base file name that will be used to automatically generate file
| names for output image files. Plots will be exported as image
| files if file_name_base is not None. (default is None)
| file_extension : str
| Valid matplotlib.pyplot file extension for savefig(). Only used
| if filename_base is not None. (default is 'png')
| mflay : int
| MODFLOW zero-based layer number to return. If None, then all
| all layers will be included. (default is None)
| kper : int
| MODFLOW zero-based stress period number to return. (default is
| zero)
| key : str
| MfList dictionary key. (default is None)
|
| Returns
| -------
| axes : list
| Empty list is returned if filename_base is not None. Otherwise
| a list of matplotlib.pyplot.axis are returned.
|
| See Also
| --------
|
| Notes
| -----
|
| Examples
| --------
| >>> import flopy
| >>> ml = flopy.modflow.Modflow.load('test.nam')
| >>> ml.dis.plot()
|
| set_cbc_output_file(self, ipakcb, model, fname)
|
| to_shapefile(self, *args, **kwargs)
| Raises AttributeError, use :meth:`export`.
|
| webdoc(self)
| Open the web documentation.
|
| ----------------------------------------------------------------------
| Static methods inherited from flopy.pakbase.Package:
|
| add_to_dtype(dtype, field_names, field_types)
| Add one or more fields to a structured array data type
|
| Parameters
| ----------
| dtype : numpy.dtype
| Input structured array datatype to add to.
| field_names : str or list
| One or more field names.
| field_types : numpy.dtype or list
| One or more data types. If one data type is supplied, it is
| repeated for each field name.
|
| load(f: Union[str, bytes, os.PathLike], model, pak_type, ext_unit_dict=None, **kwargs)
| Default load method for standard boundary packages.
|
| ----------------------------------------------------------------------
| Readonly properties inherited from flopy.pakbase.Package:
|
| data_list
|
| package_type
|
| plottable
|
| ----------------------------------------------------------------------
| Data descriptors inherited from flopy.pakbase.Package:
|
| name
|
| parent
|
| ----------------------------------------------------------------------
| Methods inherited from flopy.pakbase.PackageInterface:
|
| check(self, f=None, verbose=True, level=1, checktype=None)
| Check package data for common errors.
|
| Parameters
| ----------
| f : str or file handle
| String defining file name or file handle for summary file
| of check method output. If a string is passed a file handle
| is created. If f is None, check method does not write
| results to a summary file. (default is None)
| verbose : bool
| Boolean flag used to determine if check method results are
| written to the screen
| level : int
| Check method analysis level. If level=0, summary checks are
| performed. If level=1, full checks are performed.
| checktype : check
| Checker type to be used. By default class check is used from
| check.py.
|
| Returns
| -------
| None
|
| Examples
| --------
|
| >>> import flopy
| >>> m = flopy.modflow.Modflow.load('model.nam')
| >>> m.dis.check()
|
| ----------------------------------------------------------------------
| Readonly properties inherited from flopy.pakbase.PackageInterface:
|
| has_stress_period_data
|
| ----------------------------------------------------------------------
| Data descriptors inherited from flopy.pakbase.PackageInterface:
|
| __dict__
| dictionary for instance variables
|
| __weakref__
| list of weak references to the object
/Users/aleaf/miniforge3/envs/pyclass/lib/python3.12/site-packages/flopy/mbase.py:104: UserWarning: The program mp7 does not exist or is not executable.
warn(
Write MODPATH 7 files and run the model
[16]:
# write modpath datasets
mp.write_input()
# run modpath
mp.run_model()
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[16], line 5
2 mp.write_input()
4 # run modpath
----> 5 mp.run_model()
File ~/miniforge3/envs/pyclass/lib/python3.12/site-packages/flopy/mbase.py:1382, in BaseModel.run_model(self, silent, pause, report, normal_msg)
1352 def run_model(
1353 self,
1354 silent=False,
(...) 1357 normal_msg="normal termination",
1358 ) -> tuple[bool, list[str]]:
1359 """
1360 This method will run the model using subprocess.Popen.
1361
(...) 1379
1380 """
-> 1382 return run_model(
1383 self.exe_name,
1384 self.namefile,
1385 model_ws=self.model_ws,
1386 silent=silent,
1387 pause=pause,
1388 report=report,
1389 normal_msg=normal_msg,
1390 )
File ~/miniforge3/envs/pyclass/lib/python3.12/site-packages/flopy/mbase.py:1691, in run_model(exe_name, namefile, model_ws, silent, pause, report, processors, normal_msg, use_async, cargs, custom_print)
1689 # make sure executable exists
1690 if exe_name is None:
-> 1691 raise ValueError("An executable name or path must be provided")
1692 exe_path = resolve_exe(exe_name)
1693 if not silent:
ValueError: An executable name or path must be provided
Post-Process the MODFLOW and MODPATH Results
Load MODFLOW and MODPATH results from the heads and pathline files
Load the MODFLOW heads
[17]:
hobj = gwf.output.head()
[18]:
hds = hobj.get_data()
Load the pathline file
[21]:
ppth =ws / f"{name_mp}.mppth"
p = flopy.utils.PathlineFile(ppth)
p0 = p.get_alldata()
[22]:
p0[0:2]
[22]:
[array([(0, 0, 0, 0, 0. , 291.66666, 7041.6665, 25.157747 , 0, 1352, 0.16666667, 0.16666667, 1.6666667e-01, 1, 1),
(0, 0, 0, 0, 1066.425, 292.2685 , 7053.4917, 24.259151 , 0, 1352, 0.16907398, 0.21396767, 0.0000000e+00, 1, 1),
(0, 0, 0, 0, 9177.14 , 296.90045, 7187.993 , 18.867577 , 1, 6460, 0.18760175, 0.7519722 , 0.0000000e+00, 1, 1),
(0, 0, 0, 0, 11623.706, 298.3129 , 7250. , 17.850622 , 2, 11568, 0.19325161, 1. , 8.1138092e-01, 1, 1),
(0, 0, 0, 0, 17445.506, 327.38663, 7500. , 17.049046 , 2, 11413, 0.3095466 , 1. , 4.6866903e-01, 1, 1),
(0, 0, 0, 0, 21365.98 , 351.95578, 7750. , 14.893128 , 2, 11258, 0.40782306, 1. , 3.3380401e-01, 1, 1),
(0, 0, 0, 0, 24204.107, 384.25403, 8000. , 14.742931 , 2, 11103, 0.53701615, 1. , 2.5803149e-01, 1, 1),
(0, 0, 0, 0, 26379.746, 437.96545, 8250. , 15.481973 , 2, 11002, 0.7518618 , 1. , 2.0690233e-01, 1, 1),
(0, 0, 0, 0, 27073.143, 500. , 8325.333 , 14.168972 , 2, 10883, 1. , 0.3013334 , 1.9355534e-01, 1, 1),
(0, 0, 0, 0, 28262.705, 692.84955, 8500. , 15.88128 , 2, 10884, 0.7713981 , 1. , 1.6877499e-01, 1, 1),
(0, 0, 0, 0, 28570.799, 750. , 8540.788 , 17.17807 , 2, 10735, 1. , 0.16315117, 1.6178566e-01, 1, 1),
(0, 0, 0, 0, 30038.275, 1000. , 8630.954 , 14.299529 , 2, 10736, 1. , 0.5238181 , 1.3712475e-01, 1, 1),
(0, 0, 0, 0, 31326.36 , 1250. , 8669.666 , 14.059713 , 2, 10737, 1. , 0.6786623 , 1.1771935e-01, 1, 1),
(0, 0, 0, 0, 32350.568, 1500. , 8683.666 , 15.130941 , 2, 10738, 1. , 0.73466474, 1.0161406e-01, 1, 1),
(0, 0, 0, 0, 33184.07 , 1750. , 8649.0205, 15.428246 , 2, 10739, 1. , 0.5960806 , 8.8422537e-02, 1, 1),
(0, 0, 0, 0, 34089.45 , 2000. , 8599.658 , 13.063103 , 2, 10740, 1. , 0.39863175, 7.7607140e-02, 1, 1),
(0, 0, 0, 0, 34740.848, 2250. , 8531.887 , 14.424319 , 2, 10741, 1. , 0.12754703, 6.8339482e-02, 1, 1),
(0, 0, 0, 0, 34987.312, 2322.559 , 8500. , 10.878742 , 2, 10742, 0.290236 , 0. , 6.5897785e-02, 1, 1),
(0, 0, 0, 0, 35640.72 , 2500. , 8422.05 , 9.153153 , 2, 10891, 1. , 0.68819803, 6.0578957e-02, 1, 1),
(0, 0, 0, 0, 36104.54 , 2589.1863 , 8375. , 3.7617726, 2, 10892, 0.35674474, 0.5 , 5.8074813e-02, 1, 1),
(0, 0, 0, 0, 36279.742, 2625. , 8357.09 , 3.7571237, 2, 10892, 0.5 , 0.42836112, 5.7156086e-02, 1, 1),
(0, 0, 0, 0, 36847.414, 2750. , 8296.136 , 3.74256 , 2, 10892, 1. , 0.1845447 , 5.4277979e-02, 1, 1),
(0, 0, 0, 0, 37268.14 , 2849.683 , 8250. , 3.6527355, 2, 10895, 0.79746544, 0. , 5.2135389e-02, 1, 1),
(0, 0, 0, 0, 37389.89 , 2875. , 8235.971 , 3.6129007, 2, 11012, 1. , 0.8877678 , 5.1527575e-02, 1, 1),
(0, 0, 0, 0, 37665.7 , 2937.5 , 8205.581 , 3.7412417, 2, 11013, 0.5 , 0.64464676, 5.0115835e-02, 1, 1),
(0, 0, 0, 0, 37816.004, 2972.0994 , 8187.5 , 3.7377973, 2, 11013, 0.77679414, 0.5 , 4.9362838e-02, 1, 1),
(0, 0, 0, 0, 37918.566, 3000. , 8175.3057, 3.7354765, 2, 11013, 1. , 0.4024466 , 4.8855513e-02, 1, 1),
(0, 0, 0, 0, 38132.367, 3062.5 , 8156.827 , 3.9864547, 2, 11018, 1. , 0.5092376 , 4.7761355e-02, 1, 1),
(0, 0, 0, 0, 38138.92 , 3064.4675 , 8156.25 , 4.1631 , 2, 11019, 0.0314803 , 0.5 , 4.7726989e-02, 1, 1),
(0, 0, 0, 0, 38233.55 , 3093.75 , 8148.09 , 4.161025 , 2, 11019, 0.5 , 0.36943898, 4.7233250e-02, 1, 1),
(0, 0, 0, 0, 38327.824, 3125. , 8139.738 , 4.15898 , 2, 11019, 1. , 0.23580766, 4.6746470e-02, 1, 1),
(0, 0, 0, 0, 38372.035, 3140.625 , 8136.8047, 4.203094 , 2, 11039, 0.5 , 0.37775394, 4.6436224e-02, 1, 1),
(0, 0, 0, 0, 38415.336, 3156.25 , 8133.7593, 4.2018466, 2, 11039, 1. , 0.28029665, 4.6134390e-02, 1, 1),
(0, 0, 0, 0, 38457.65 , 3171.875 , 8131.8296, 4.0527744, 2, 11042, 1. , 0.43707943, 4.2602737e-02, 1, 1),
(0, 0, 0, 0, 38564.637, 3187.5 , 8128.147 , 3.803116 , 2, 11043, 1. , 0.2013962 , 1.2858289e-03, 1, 1),
(0, 0, 0, 0, 38660.48 , 3191.2632 , 8125. , 3.708662 , 2, 11049, 0.12042002, 0. , 1.1944468e-03, 1, 1),
(0, 0, 0, 0, 38862.535, 3203.125 , 8118.368 , 3.697433 , 2, 11065, 0.5 , 0.7877901 , 1.0352990e-03, 1, 1),
(0, 0, 0, 0, 39095.668, 3218.75 , 8110.9272, 3.6967647, 2, 11065, 1. , 0.5496728 , 8.7783765e-04, 1, 1),
(0, 0, 0, 0, 39488.062, 3250. , 8098.543 , 3.5428562, 2, 11066, 1. , 0.15338022, 8.3847297e-04, 1, 1),
(0, 0, 0, 0, 39600.6 , 3259.7537 , 8093.75 , 3.3440962, 2, 11088, 0.1560589 , 0.5 , 8.2857930e-04, 1, 1),
(0, 0, 0, 0, 39877.094, 3281.25 , 8082.4824, 3.343992 , 2, 11088, 0.5 , 0.31971985, 8.0476422e-04, 1, 1),
(0, 0, 0, 0, 40233.7 , 3312.5 , 8070.1533, 3.3438623, 2, 11088, 1. , 0.12245228, 7.7505718e-04, 1, 1),
(0, 0, 0, 0, 40460.684, 3333.9407 , 8062.5 , 3.0362926, 2, 11089, 0.34305042, 0. , 7.5718301e-04, 1, 1),
(0, 0, 0, 0, 40903.15 , 3375. , 8047.6304, 3.1637266, 2, 11091, 1. , 0.7620828 , 7.2312116e-04, 1, 1),
(0, 0, 0, 0, 41541.184, 3437.5 , 8013.1743, 2.95539 , 2, 11092, 0.5 , 0.10539599, 6.7761948e-04, 1, 1),
(0, 0, 0, 0, 41800.92 , 3465.193 , 8000. , 2.9553099, 2, 11092, 0.7215439 , 0. , 6.5992685e-04, 1, 1),
(0, 0, 0, 0, 41982.45 , 3500. , 7990.8164, 3.210435 , 2, 11191, 1. , 0.9265324 , 6.4754800e-04, 1, 1),
(0, 0, 0, 0, 42502.145, 3625. , 7968.0464, 3.2902412, 2, 11206, 1. , 0.74437314, 6.1345613e-04, 1, 1),
(0, 0, 0, 0, 42670.46 , 3680.2456 , 7937.5 , 3.2615755, 2, 11207, 0.4419644 , 0.5 , 6.0290669e-04, 1, 1),
(0, 0, 0, 0, 42686.645, 3687.5 , 7934.9727, 3.261571 , 2, 11207, 0.5 , 0.47978276, 6.0190214e-04, 1, 1),
(0, 0, 0, 0, 42795.89 , 3750. , 7915.1133, 3.261541 , 2, 11207, 1. , 0.32090622, 5.9516344e-04, 1, 1),
(0, 0, 0, 0, 42827.008, 3773.8174 , 7906.25 , 3.237704 , 2, 11212, 0.38107696, 0.5 , 5.9033406e-04, 1, 1),
(0, 0, 0, 0, 42834.42 , 3781.25 , 7904.3936, 3.2376988, 2, 11212, 0.5 , 0.47029385, 5.8918889e-04, 1, 1),
(0, 0, 0, 0, 42858.8 , 3812.5 , 7896.461 , 3.237682 , 2, 11212, 1. , 0.34337416, 5.8543991e-04, 1, 1),
(0, 0, 0, 0, 42867.297, 3828.125 , 7891.5454, 3.2256145, 2, 11215, 0.5 , 0.5294564 , 5.5021694e-04, 1, 1),
(0, 0, 0, 0, 42868.305, 3830.3154 , 7890.625 , 3.2255962, 2, 11215, 0.5700975 , 0.5 , 5.4618990e-04, 1, 1),
(0, 0, 0, 0, 42872.234, 3843.75 , 7887.5947, 3.225526 , 2, 11215, 1. , 0.40302902, 5.3074013e-04, 1, 1),
(0, 0, 0, 0, 42874.312, 3859.375 , 7885.291 , 3.2016647, 2, 11219, 0. , 0.65862405, 3.6462079e-04, 1, 1),
(0, 0, 0, 0, 42874.312, 3859.375 , 7885.291 , 3.2188168, 2, 11218, 1. , 0.65862405, 3.6462079e-04, 1, 1)],
dtype=[('particleid', '<i4'), ('particlegroup', '<i4'), ('sequencenumber', '<i4'), ('particleidloc', '<i4'), ('time', '<f4'), ('x', '<f4'), ('y', '<f4'), ('z', '<f4'), ('k', '<i4'), ('node', '<i4'), ('xloc', '<f4'), ('yloc', '<f4'), ('zloc', '<f4'), ('stressperiod', '<i4'), ('timestep', '<i4')]),
array([(1, 0, 1, 1, 0. , 375. , 7041.6665, 25.157747 , 0, 1352, 0.5 , 0.16666667, 0.16666667, 1, 1),
(1, 0, 1, 1, 1066.425, 375.7526 , 7053.4917, 24.259151 , 0, 1352, 0.50301033, 0.21396767, 0. , 1, 1),
(1, 0, 1, 1, 9177.14 , 381.53244, 7187.993 , 18.867577 , 1, 6460, 0.5261298 , 0.7519722 , 0. , 1, 1),
(1, 0, 1, 1, 11623.706, 383.2931 , 7250. , 17.850622 , 2, 11568, 0.5331724 , 1. , 0.8113809 , 1, 1),
(1, 0, 1, 1, 17445.506, 406.1642 , 7500. , 17.049046 , 2, 11413, 0.62465686, 1. , 0.46866903, 1, 1),
(1, 0, 1, 1, 21365.98 , 429.19812, 7750. , 14.893128 , 2, 11258, 0.71679246, 1. , 0.333804 , 1, 1),
(1, 0, 1, 1, 24204.107, 465.9783 , 8000. , 14.742931 , 2, 11103, 0.86391324, 1. , 0.2580315 , 1, 1),
(1, 0, 1, 1, 25402.314, 500. , 8134.31 , 15.579909 , 2, 11002, 1. , 0.5372406 , 0.22848248, 1, 1),
(1, 0, 1, 1, 26073.605, 518.39496, 8250. , 15.031157 , 2, 11003, 0.07357987, 1. , 0.21383499, 1, 1),
(1, 0, 1, 1, 27358.307, 750. , 8465.95 , 15.943902 , 2, 10884, 1. , 0.8638008 , 0.18442652, 1, 1),
(1, 0, 1, 1, 28116.188, 894.3324 , 8500. , 12.073312 , 2, 10885, 0.5773297 , 1. , 0.17199726, 1, 1),
(1, 0, 1, 1, 28714.22 , 1000. , 8535.634 , 14.396254 , 2, 10736, 1. , 0.14253446, 0.1607872 , 1, 1),
(1, 0, 1, 1, 30002.307, 1250. , 8562.656 , 14.138697 , 2, 10737, 1. , 0.2506236 , 0.13803317, 1, 1),
(1, 0, 1, 1, 31026.514, 1500. , 8564.1875, 15.187186 , 2, 10738, 1. , 0.2567515 , 0.11914873, 1, 1),
(1, 0, 1, 1, 31860.018, 1750. , 8521.76 , 15.470393 , 2, 10739, 1. , 0.08703894, 0.10368085, 1, 1),
(1, 0, 1, 1, 32197.521, 1841.701 , 8500. , 13.130724 , 2, 10740, 0.36680412, 0. , 0.09875897, 1, 1),
(1, 0, 1, 1, 32723.895, 2000. , 8464.922 , 13.120698 , 2, 10889, 1. , 0.8596867 , 0.09138791, 1, 1),
(1, 0, 1, 1, 33319.28 , 2250. , 8386.608 , 14.194465 , 2, 10890, 1. , 0.5464344 , 0.08128989, 1, 1),
(1, 0, 1, 1, 34236.707, 2500. , 8257.554 , 9.194842 , 2, 10891, 1. , 0.0302132 , 0.07223024, 1, 1),
(1, 0, 1, 1, 34313.992, 2514.7893 , 8250. , 3.8308406, 2, 10892, 0.05915729, 0. , 0.07172394, 1, 1),
(1, 0, 1, 1, 34992.023, 2625. , 8182.07 , 3.7971873, 2, 11011, 0.5 , 0.72828007, 0.06734534, 1, 1),
(1, 0, 1, 1, 35506.47 , 2710.23 , 8125. , 3.781605 , 2, 11011, 0.84092045, 0.5 , 0.06420238, 1, 1),
(1, 0, 1, 1, 35697.39 , 2750. , 8103.5493, 3.776009 , 2, 11011, 1. , 0.41419792, 0.06307365, 1, 1),
(1, 0, 1, 1, 36268.164, 2875. , 8047.6475, 3.7971103, 2, 11014, 1. , 0.38118196, 0.05963819, 1, 1),
(1, 0, 1, 1, 36544.617, 2937.5 , 8022.5825, 3.747776 , 2, 11015, 0.5 , 0.1806616 , 0.05799387, 1, 1),
(1, 0, 1, 1, 36797.656, 2999.5725 , 8000. , 3.7411008, 2, 11015, 0.9965789 , 0. , 0.0565286 , 1, 1),
(1, 0, 1, 1, 36799.88 , 3000. , 7999.8086, 3.719146 , 2, 11114, 1. , 0.99846715, 0.05651587, 1, 1),
(1, 0, 1, 1, 37093.63 , 3062.5 , 7979.1167, 3.706241 , 2, 11117, 1. , 0.6658651 , 0.05481534, 1, 1),
(1, 0, 1, 1, 37227.023, 3093.75 , 7972.9053, 3.6043687, 2, 11118, 0.5 , 0.5664867 , 0.05405088, 1, 1),
(1, 0, 1, 1, 37310.426, 3113.6023 , 7968.75 , 3.6022983, 2, 11118, 0.81763643, 0.5 , 0.05357835, 1, 1),
(1, 0, 1, 1, 37355.496, 3125. , 7966.5186, 3.6011868, 2, 11118, 1. , 0.4642958 , 0.05332471, 1, 1),
(1, 0, 1, 1, 37416.727, 3140.625 , 7964.9746, 3.5058205, 2, 11126, 0.5 , 0.8791891 , 0.05288708, 1, 1),
(1, 0, 1, 1, 37478.727, 3156.25 , 7963.331 , 3.5038984, 2, 11126, 1. , 0.8265902 , 0.05244757, 1, 1),
(1, 0, 1, 1, 37537.65 , 3171.875 , 7961.9463, 3.5626056, 2, 11127, 1. , 0.5645512 , 0.0483855 , 1, 1),
(1, 0, 1, 1, 37646.598, 3187.5 , 7959.6675, 3.3983364, 2, 11128, 1. , 0.41871217, 0.0042506 , 1, 1),
(1, 0, 1, 1, 37893.9 , 3203.125 , 7954.166 , 3.4010422, 2, 11133, 0.5 , 0.5333076 , 0.00369481, 1, 1),
(1, 0, 1, 1, 37946.367, 3206.6306 , 7953.125 , 3.400573 , 2, 11133, 0.6121834 , 0.5 , 0.00358658, 1, 1),
(1, 0, 1, 1, 38120.16 , 3218.75 , 7949.732 , 3.3991148, 2, 11133, 1. , 0.39142844, 0.00325021, 1, 1),
(1, 0, 1, 1, 38526.08 , 3250. , 7943.0986, 3.4536526, 2, 11134, 1. , 0.17915015, 0.00310237, 1, 1),
(1, 0, 1, 1, 38884.27 , 3281.25 , 7937.916 , 3.6084507, 2, 11187, 0.5 , 0.00665233, 0.00298445, 1, 1),
(1, 0, 1, 1, 38915.113, 3284.1855 , 7937.5 , 3.6084082, 2, 11187, 0.5469703 , 0. , 0.00297451, 1, 1),
(1, 0, 1, 1, 39204.26 , 3312.5 , 7933.8066, 3.3178751, 2, 11189, 1. , 0.9409047 , 0.00288494, 1, 1),
(1, 0, 1, 1, 39746.887, 3375. , 7927.907 , 3.42754 , 2, 11190, 1. , 0.84651566, 0.00272271, 1, 1),
(1, 0, 1, 1, 40144.082, 3437.5 , 7925.4526, 3.2191124, 2, 11191, 0.5 , 0.403621 , 0.00261221, 1, 1),
(1, 0, 1, 1, 40462.62 , 3500. , 7913.7847, 3.2187355, 2, 11191, 1. , 0.3102781 , 0.00252683, 1, 1),
(1, 0, 1, 1, 40982.32 , 3625. , 7912.5156, 3.2981243, 2, 11206, 1. , 0.30012494, 0.0023938 , 1, 1),
(1, 0, 1, 1, 41153.348, 3687.5 , 7899.7046, 3.26941 , 2, 11207, 0.5 , 0.19763616, 0.00235198, 1, 1),
(1, 0, 1, 1, 41262.594, 3750. , 7891.862 , 3.2692919, 2, 11207, 1. , 0.13489589, 0.00232565, 1, 1),
(1, 0, 1, 1, 41298.973, 3781.25 , 7888.7163, 3.245443 , 2, 11212, 0.5 , 0.21946454, 0.0023036 , 1, 1),
(1, 0, 1, 1, 41323.348, 3812.5 , 7885.8555, 3.2453768, 2, 11212, 1. , 0.17368884, 0.00228894, 1, 1),
(1, 0, 1, 1, 41331.047, 3828.125 , 7884.7173, 3.232934 , 2, 11215, 0.5 , 0.31095538, 0.00216386, 1, 1),
(1, 0, 1, 1, 41335.785, 3843.75 , 7883.209 , 3.2326002, 2, 11215, 1. , 0.2626827 , 0.00209027, 1, 1),
(1, 0, 1, 1, 41337.863, 3859.375 , 7883.1943, 3.206547 , 2, 11219, 0. , 0.52444786, 0.00143603, 1, 1),
(1, 0, 1, 1, 41337.863, 3859.375 , 7883.1943, 3.223687 , 2, 11218, 1. , 0.52444786, 0.00143603, 1, 1)],
dtype=[('particleid', '<i4'), ('particlegroup', '<i4'), ('sequencenumber', '<i4'), ('particleidloc', '<i4'), ('time', '<f4'), ('x', '<f4'), ('y', '<f4'), ('z', '<f4'), ('k', '<i4'), ('node', '<i4'), ('xloc', '<f4'), ('yloc', '<f4'), ('zloc', '<f4'), ('stressperiod', '<i4'), ('timestep', '<i4')])]
Plot the heads and pathlines
[23]:
fig, ax = plt.subplots(figsize=(5, 9))
mm = flopy.plot.PlotMapView(model=gwf, layer=0, ax=ax)
mm.plot_array(hds, masked_values=[1e30])
mm.plot_bc('SFR', color='b', plotAll=True)
mm.plot_bc('WEL', plotAll=True)
mm.plot_ibound()
mm.plot_pathline(p0, layer='all', color='blue', lw=1)
mm.plot_grid(lw=0.2, color="0.5")
ax = plt.gca()
ax.fill(csx, csy, color='#e534eb', zorder=100, alpha=.75);
[19]:
fig, ax = plt.subplots(figsize=(5, 9))
mm = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=0)
mm.plot_array(hds, masked_values=[1e30])
mm.plot_bc('SFR', color='b', plotAll=True)
mm.plot_bc('WEL', plotAll=True)
mm.plot_ibound()
mm.plot_grid(lw=0.2, color="0.5")
p = flopy.utils.PathlineFile(ppth)
p0 = p.get_alldata()
df = pd.DataFrame.from_records(p0[0])
for ix, recarray in enumerate(p0):
dft = pd.DataFrame.from_records(recarray)
df = pd.concat((df, dft), ignore_index=True)
vmin, vmax = df.time.min(), df.time.max()
times = list(range(0, 74001, 1000))
for ix in range(1, len(times)):
tmp = df[(df.time >= times[ix - 1]) & (df.time < times[ix])]
s = ax.scatter(tmp.x.values, tmp.y.values, c=tmp.time.values, vmin=vmin, vmax=vmax, cmap="magma")
ax.set_title(f"{times[ix - 1]} - {times[ix]} days")
if ix == 1:
cbar = fig.colorbar(s, shrink=0.7)
cbar.set_label('particle travel time from source', rotation=270, labelpad=14)
display(fig)
clear_output(wait=True)
plt.pause(0.1)
Export the pathlines to a shapefile
[ ]:
spth = Path(ws / 'pathline.shp')
p.write_shapefile(p0, mg=gwf.modelgrid, one_per_particle=False, shpname=spth)
[('particleid', '<i4'), ('particlegroup', '<i4'), ('sequencenumber', '<i4'), ('particleidloc', '<i4'), ('time', '<f4'), ('x', '<f4'), ('y', '<f4'), ('z', '<f4'), ('k', '<i4'), ('node', '<i4'), ('xloc', '<f4'), ('yloc', '<f4'), ('zloc', '<f4'), ('stressperiod', '<i4'), ('timestep', '<i4')]
C:\Users\jlarsen\Documents\anaconda3\envs\pyclass3\Lib\site-packages\flopy\export\shapefile_utils.py:496: UserWarning: Truncating shapefile fieldname particlegroup to partiroup_
warn(f"Truncating shapefile fieldname {s} to {name}")
C:\Users\jlarsen\Documents\anaconda3\envs\pyclass3\Lib\site-packages\flopy\export\shapefile_utils.py:496: UserWarning: Truncating shapefile fieldname sequencenumber to sequember_
warn(f"Truncating shapefile fieldname {s} to {name}")
C:\Users\jlarsen\Documents\anaconda3\envs\pyclass3\Lib\site-packages\flopy\export\shapefile_utils.py:496: UserWarning: Truncating shapefile fieldname particleidloc to partidloc_
warn(f"Truncating shapefile fieldname {s} to {name}")
C:\Users\jlarsen\Documents\anaconda3\envs\pyclass3\Lib\site-packages\flopy\export\shapefile_utils.py:496: UserWarning: Truncating shapefile fieldname stressperiod to stresriod_
warn(f"Truncating shapefile fieldname {s} to {name}")
No CRS information for writing a .prj file.
Supply an valid coordinate system reference to the attached modelgrid object or .export() method.
Load the shapefile into geopandas
[21]:
rshp = gpd.read_file(spth)
Show the dataframe
[22]:
rshp
[22]:
| particleid | partiroup_ | sequember_ | partidloc_ | time | x | y | z | k | node | xloc | yloc | zloc | stresriod_ | timestep | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1 | 1 | 1066.444824 | 292.266235 | 7053.407715 | 24.259151 | 1 | 1353 | 0.169065 | 0.213630 | 0.000000 | 1 | 1 | LINESTRING (291.667 7041.667, 292.266 7053.408) |
| 1 | 1 | 1 | 1 | 1 | 9176.973633 | 296.880524 | 7186.960938 | 18.867577 | 2 | 6461 | 0.187522 | 0.747843 | 0.000000 | 1 | 1 | LINESTRING (292.266 7053.408, 296.881 7186.961) |
| 2 | 1 | 1 | 1 | 1 | 11675.894531 | 298.317780 | 7250.000000 | 17.831028 | 3 | 11569 | 0.193271 | 1.000000 | 0.807747 | 1 | 1 | LINESTRING (296.881 7186.961, 298.318 7250) |
| 3 | 1 | 1 | 1 | 1 | 17508.531250 | 327.402008 | 7500.000000 | 17.036358 | 3 | 11414 | 0.309608 | 1.000000 | 0.466072 | 1 | 1 | LINESTRING (298.318 7250, 327.402 7500) |
| 4 | 1 | 1 | 1 | 1 | 21433.228516 | 351.983856 | 7750.000000 | 14.882585 | 3 | 11259 | 0.407935 | 1.000000 | 0.331823 | 1 | 1 | LINESTRING (327.402 7500, 351.984 7750) |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 7326 | 162 | 1 | 162 | 162 | 49935.636719 | 3312.500000 | 5125.286133 | 1.553995 | 3 | 12772 | 1.000000 | 0.004577 | 0.009084 | 1 | 1 | LINESTRING (3305.245 5125, 3312.5 5125.286) |
| 7327 | 162 | 1 | 162 | 162 | 49954.375000 | 3328.125000 | 5128.049805 | 1.580439 | 3 | 12775 | 0.500000 | 0.097600 | 0.008830 | 1 | 1 | LINESTRING (3312.5 5125.286, 3328.125 5128.05) |
| 7328 | 162 | 1 | 162 | 162 | 49966.113281 | 3343.750000 | 5128.823730 | 1.579678 | 3 | 12775 | 1.000000 | 0.122366 | 0.008674 | 1 | 1 | LINESTRING (3328.125 5128.05, 3343.75 5128.824) |
| 7329 | 162 | 1 | 162 | 162 | 49971.105469 | 3359.375000 | 5131.083496 | 1.578258 | 3 | 12778 | 1.000000 | 0.389351 | 0.005555 | 1 | 1 | LINESTRING (3343.75 5128.824, 3359.375 5131.083) |
| 7330 | 162 | 1 | 162 | 162 | 49971.105469 | 3359.375000 | 5131.083496 | 1.587265 | 3 | 12779 | 0.000000 | 0.389351 | 0.005555 | 1 | 1 | LINESTRING (3359.375 5131.083, 3359.375 5131.083) |
7331 rows × 16 columns
Plot the geopandas dataframe on the modelgrid
[23]:
fig, ax = plt.subplots(figsize=(6, 10))
pmv = flopy.plot.PlotMapView(gwf, ax=ax)
pmv.plot_grid(lw=0.5)
pmv.plot_shapes(rshp.geometry.to_list(), edgecolor="grey") #, edgecolor="b")
pmv.plot_ibound();
Export the data to VTK (Visualization ToolKit) format for a 3d representation
FloPy’s Vtk module allows us to create Visualization ToolKit (VTK) files that can be opened and explored with pyVISTAs or external software such as ParaView.
[24]:
from flopy.export.vtk import Vtk
vtk = Vtk(model=gwf, binary=False, vertical_exageration=50, smooth=False)
vtk.add_model(gwf)
vtk.add_pathline_points(p0)
Write VTK file for use in ParaView or other VTK software
[ ]:
pth = Path("temp")
pth.mkdir(exist_ok=True)
vtk.write(pth / "freyberg")
Alternatively we can visualize the VTK representation using pyvista. To activate these cells, set plotvtk=True
[26]:
plotvtk = True
if plotvtk:
import pyvista as pv
Convert the VTK to pyvista meshes and rotate the meshes to match the orientation of our matplotlib plots
[27]:
if plotvtk:
grid, pathlines = vtk.to_pyvista()
[28]:
if plotvtk:
axes = pv.Axes(show_actor=True, actor_scale=2.0, line_width=5)
grid.rotate_z(160, point=axes.origin, inplace=True)
pathlines.rotate_z(160, point=axes.origin, inplace=True)
Select particle release locations and build a dictionary of particle tracks (pathlines). This will be used below for particle labelling, as well as for animation.
Note: while below we construct pathlines manually from data read from the exported VTK files, pathlines may also be read directly from the MODPATH 7 pathline output file (provided the simulation was run in pathline or combined mode, as this one was).
[29]:
if plotvtk:
tracks = {}
particle_ids = set()
release_locs = list()
for i, t in enumerate(pathlines["time"]):
pid = str(round(float(pathlines["particleid"][i])))
loc = pathlines.points[i]
if pid not in tracks:
tracks[pid] = []
particle_ids.add(pid)
release_locs.append(loc)
# store the particle location in the corresponding track
tracks[pid].append((loc, t))
release_locs = np.array(release_locs)
tracks = {k: np.array(v, dtype=object) for k, v in tracks.items()}
max_track_len = max([len(v) for v in tracks.values()])
View the grid and pathlines with PyVista, with particle tracks/locations colored by time. Also add particle ID labels to a few particles’ release locations.
[30]:
if plotvtk:
pv.set_plot_theme("document")
pv.set_jupyter_backend("static")
# create the plot and add the grid and pathline meshes
p = pv.Plotter()
p.add_mesh(grid, opacity=0.05)
p.add_mesh(pathlines, scalars="time")
# add a particle ID label to each 20th particle's starting point
label_coords = []
start_labels = []
for pid, track in tracks.items():
if int(pid) % 20 == 0:
label_coords.append(track[0][0])
start_labels.append(f"Particle {pid}")
p.add_point_labels(
label_coords,
start_labels,
font_size=10,
point_size=15,
point_color="black",
)
# zoom in and show the plot
p.camera.zoom(2.4)
p.show()
Create an animated GIF of the particles traveling along their pathlines, with particles colored by time.
[31]:
if plotvtk:
# create plotter
p = pv.Plotter(notebook=False, off_screen=True)
# open GIF file
gif_path = pth / f"freyberg_tracks.gif"
p.open_gif(str(gif_path))
# create mesh from release locations
spls = pv.PolyData(release_locs)
spls.point_data["time"] = np.zeros(len(spls.points))
# add the underlying grid mesh and particle data, then zoom in
p.add_mesh(grid, opacity=0.05)
p.add_mesh(spls, clim=[0, 1.23e09])
p.camera.zoom(2.4)
# cycle through time steps and update particle location
for i in range(1, max_track_len):
pts = []
times = []
segments = []
for pid in particle_ids:
track = tracks[pid]
npts = len(track)
# use last locn if particle has already terminated
loc, t = track[i] if i < npts else track[npts - 1]
pts.append(loc)
times.append(t)
if i < npts:
segments.append(track[i - 1][0])
segments.append(loc)
spls.points = np.vstack(pts)
spls.point_data["time"] = np.array(times)
p.add_lines(np.array(segments), width=1, color="black")
p.write_frame() # write frame to file
# close the plotter and the GIF file
p.close()
[32]:
if plotvtk:
from IPython.core.display import Image
display(Image(data=open(gif_path, "rb").read(), format="gif"))
<IPython.core.display.Image object>
[ ]: