Repository navigation
feature: Modpath7 zones defined for DISV grid #2612
Description
I am trying to run particle tracking on modflow 6 model that has a unstructued grid (DISV package). When i run the modpath7 for it. It runs smoothly but when i want to define zones, it gives error and ask for a 3d array (layer, row, col). First, I thought that modpath7 is not recognizing that its an unstructured grid but then I checked the type of data allowed in defining zones in flopy.modpath.mp7sim module and i understood that there is no option for unstructured grid array, i.e 2d array of (layer, nodes) .
I would really like if someone can explain if there is something that i am understanding wrong, or there is another way to put zones and stopzones in modpath7 which i am unaware of.
Activity
wpbonelli commentedon Sep 27, 2025
I think you should be able to provide either a 2D array of shape
(nlay, ncpl)or a 3d array of shape(nlay, 1, ncpl). Do you see an error with either of these?The docstring does make it seem like it must be a structured grid cell index.
hafmood95 commentedon Sep 27, 2025
I provide 2d array (nlay, ncpl) and it gives error saying 3d array required. I will try doing (nlay, 1, ncpl) and can give an update if it works
hafmood95 commentedon Sep 28, 2025
Update: I have tried both. it does not work and gives the same error. I am attaching the code i am runing and the error. When i remove the stopzones, zones and zonedata option in the code the model run just fine.
npf = gwf.get_package("npf")
##get hydraulic conductivity array
k_horizontal = npf.k.array
##Store vertical hydraulic conductivity (K33) as a variable
k_vertical = npf.k33.array
#defining zones starting for my 12 layers [100,101,102,103,104,105,106,107,108,109,110 111]
layer_vals = np.arange(100, 100 + nlay).reshape(nlay, 1)
##broadcast to full array
zones = np.broadcast_to(layer_vals, (nlay, ncpl)) #resultant array (nlay,ncpl)
zones_ = np.expand_dims(zones, axis=1) #resultant array (nlay,1, ncpl)
aquifer_stopzones=[100,101,103,105,107,109,111]
pdata = ParticleData(partlocs=cellid_loc_iol,
particleids=particl_id_iol, structured=False,
localx=x_loc_iol, localy=y_loc_iol, localz=z_loc_iol, drape=1)
pg = ParticleGroup(particlegroupname="my_particles", particledata=pdata, filename="sgo.sloc")
mp7_path=r"C:\als\ALS_092125\mp7.exe"
mp = flopy.modpath.Modpath7(
modelname='mp',
flowmodel=gwf,
exe_name=mp7_path,
model_ws=sim_ws)
mpbas = flopy.modpath.Modpath7Bas(mp, porosity=0.3)
for ij in aquifer_stopzones:
mpsim = flopy.modpath.Modpath7Sim(
mp,endpointfilename='zone'+str(ij)+'.mpend',
simulationtype="endpoint",
trackingdirection="forward",
weaksinkoption="pass_through",
weaksourceoption="pass_through",
budgetoutputoption="summary",
zones=zones_,
referencetime=[0, 0, 0.0],
stoptimeoption="extend",
zonedataoption ='on',
stopzone=ij,
particlegroups=pg)
#write modpath datasets
mp.write_input()
ERROR*********
Traceback (most recent call last):
File "", line 41, in
particlegroups=pg)
File "C:\anaconda\envs\ropen\lib\site-packages\flopy\modpath\mp7sim.py", line 531, in init
locat=self.unit_number[0],
File "C:\anaconda\envs\ropen\lib\site-packages\flopy\utils\util_array.py", line 560, in init
"Util3d: expected 3 dimensions, found shape {0}".format(shape)
ValueError: Util3d: expected 3 dimensions, found shape (12, 218793)
zones_.shape
Out[15]: (12, 1, 218793)
wpbonelli commentedon Sep 28, 2025
It looks like you have a pretty old version of flopy. This issue seems to have been fixed in #1415 a few years ago. Does updating solve it?