modflowpy / flopy

A Python package to create, run, and post-process MODFLOW-based models.
https://flopy.readthedocs.io
Other
506 stars 306 forks source link

bug: reverse Headfile fails for structured grids #2246

Closed cneyens closed 2 months ago

cneyens commented 2 months ago

Describe the bug We're doing some backward particle tracking using PRT in MODFLOW 6 on a structured grid. To reverse the flowfield, the reverse() methods for the head- and budgetfile are used. This returns an error when running MODFLOW 6. Looking in more detail, it seems that when writing the reversed head file to disk, only the first row is written for each layer. Seems the culprit is the second 0-index subsetting when getting the data to write to disk:

https://github.com/modflowpy/flopy/blob/2d372c653ad17aea0f93deb4a5d1e6b5d385c5d6/flopy/utils/binaryfile.py#L737

The budgetfile on the other hand is reversed as expected. For unstructured grids (DISV), everything works fine.

To Reproduce Reproducible example:

Fold code ```python import flopy import tempfile import os ws = tempfile.TemporaryDirectory().name name = 'mymodel' sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws, exe_name='mf6') tdis = flopy.mf6.ModflowTdis(sim) ims = flopy.mf6.ModflowIms(sim) gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) dis = flopy.mf6.ModflowGwfdis(gwf, nrow=10, ncol=10) ic = flopy.mf6.ModflowGwfic(gwf) npf = flopy.mf6.ModflowGwfnpf(gwf, save_specific_discharge=True) chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=[[(0, 0, 0), 1.], [(0, 9, 9), 0.]]) budget_file = name + '.bud' head_file = name + '.hds' oc = flopy.mf6.ModflowGwfoc(gwf, budget_filerecord=budget_file, head_filerecord=head_file, saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')]) sim.write_simulation(silent=True) success, buff = sim.run_simulation(silent=True) assert success, "Failed to converge" # connect to head file heads = flopy.utils.HeadFile(os.path.join(ws, head_file), tdis=tdis) # this creates a reversed head file, but only saves the the first row in each layer which will cause PRT to error out head_file_bckwrd = os.path.join(ws, head_file + '_backward') heads.reverse(head_file_bckwrd) # read to check dimensions heads_reverse = flopy.utils.HeadFile(head_file_bckwrd, tdis=tdis) # this throws a reshaping error heads_reverse.get_data() ```

Expected behavior All cells are included in the reversed head file on disk, for all types of grids.