From c7af787110eb1d984fa660c49014c999d87b0774 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Wed, 1 May 2024 08:56:12 -0400 Subject: [PATCH] refactor(mp7particledata): match mp7 order in to_coords()/to_prp() (#2172) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * return the same particle order from to_coords() and to_prp() as generated by MP7 — preserves particle indexing for MF6 PRT, making it easier to compare results * refactoring/cleanup, factor out helper methods get_cell_release_points(), get_face_release_points() * defer materialization — return iterators not lists so the consumer can control memory use, previously done for to_coords() and to_prp() but not get_release_points(), might be relevant for very large release configurations * test-drive a pytest plugin for previously hardcoded snapshot tests, introduce some fixtures in autotest/conftest.py for nicer array snapshots than the default format, maybe these could live in devtools if they see wider use? --- DEVELOPER.md | 11 + ...test_lrcparticledata_to_prp_1_per_face.npy | Bin 0 -> 464 bytes ...edata_to_prp_divisions_defaults[False].npy | Bin 0 -> 3152 bytes ...ledata_to_prp_divisions_defaults[True].npy | Bin 0 -> 3152 bytes autotest/conftest.py | 105 +++++- autotest/test_particledata.py | 189 ++++++---- etc/environment.yml | 1 + flopy/modpath/mp7particledata.py | 350 ++++++++++-------- pyproject.toml | 1 + 9 files changed, 430 insertions(+), 227 deletions(-) create mode 100644 autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_1_per_face.npy create mode 100644 autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_divisions_defaults[False].npy create mode 100644 autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_divisions_defaults[True].npy diff --git a/DEVELOPER.md b/DEVELOPER.md index 8dc8dee0a5..ac775a7f01 100644 --- a/DEVELOPER.md +++ b/DEVELOPER.md @@ -28,6 +28,7 @@ This document describes how to set up a FloPy development environment, run the e - [Performance testing](#performance-testing) - [Benchmarking](#benchmarking) - [Profiling](#profiling) + - [Snapshot testing](#snapshot-testing) - [Branching model](#branching-model) - [Deprecation policy](#deprecation-policy) - [Miscellaneous](#miscellaneous) @@ -345,6 +346,16 @@ Profiling is [distinct](https://stackoverflow.com/a/39381805/6514033) from bench By default, `pytest-benchmark` will only print profiling results to `stdout`. If the `--benchmark-autosave` flag is provided, performance profile data will be included in the JSON files written to the `.benchmarks` save directory as described in the benchmarking section above. +### Snapshot testing + +Snapshot testing is a form of regression testing in which a "snapshot" of the results of some computation is verified and captured by the developer to be compared against when tests are subsequently run. This is accomplished with [`syrupy`](https://github.com/tophat/syrupy), which provides a `snapshot` fixture overriding the equality operator to allow comparison with e.g. `snapshot == result`. A few custom fixtures for snapshots of NumPy arrays are provided in `autotest/conftest.py`: + +- `array_snapshot`: saves an array in a binary file for compact storage, can be inspected programmatically with `np.load()` +- `text_array_snapshot`: flattens an array and stores it in a text file, compromise between readability and disk usage +- `readable_array_snapshot`: stores an array in a text file in its original shape, easy to inspect but largest on disk + +By default, tests run in comparison mode. This means a newly written test using any of the snapshot fixtures will fail until a snapshot is created. Snapshots can be created/updated by running pytest with the `--snapshot-update` flag. + ## Branching model This project follows the [git flow](https://nvie.com/posts/a-successful-git-branching-model/): development occurs on the `develop` branch, while `master` is reserved for the state of the latest release. Development PRs are typically squashed to `develop`, to avoid merge commits. At release time, release branches are merged to `master`, and then `master` is merged back into `develop`. diff --git a/autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_1_per_face.npy b/autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_1_per_face.npy new file mode 100644 index 0000000000000000000000000000000000000000..16a9342cb209db5d4fc7c842a42260ec6ef9d96d GIT binary patch literal 464 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$dItu2RItsN4WCJb+tl)z^L>QfZ02P;WK$V9Xg07wcst#QrOdOZKA5eQ> l@-X{2pyt8E(d~ojgNeiR2@t6d=1vi)dKgW}KA1ig2LPXIIxPSI literal 0 HcmV?d00001 diff --git a/autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_divisions_defaults[False].npy b/autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_divisions_defaults[False].npy new file mode 100644 index 0000000000000000000000000000000000000000..f9177363ee799e3dd665d7e50113f74b00dd6feb GIT binary patch literal 3152 zcmbW&y=oL;6vpw_POK~hEkuSTWTGsRXugaw6F(B8px8vPlErLDAx3t^MgnsW-bK8Q z+(&Aw#cCVD&7Ao=oR{-hI7MD|a(?qX%)2xD{q_Eo^!|8sFuXXwJNP&r&d!gI2A_ws&kHNd{bKEQzc?xW_qdY&c}%Jwrx(Ys zzw3T;T;Db}sbbFgsc@&So9q7;?y_|Bv@1?4-}9<;%(HrTO*(qU?uK;q#QVIP z(vcIpLOSMIeO^yGda_#nt#*D+I-BeE-Me||=ojZLNayFL?e~$Mccimk%(LS9($Uj< z++FGD={;^yI&xySBpvgtb{+EWNk>ms`@Zzr`TNqP=IHj_yJhL<7w4@=$2_?{Z&fh39rLU{?_uV6Pj~Nmn?K9mpPSOrGj@-p^ZoL?E$O(Q zmV2J{tfaGfp7Y#oOGnQzuhx#9&CloQ&#dQTnMcpi)!NZBbm@uAqbK`k`R#k6*3Rbm zojuFmyH8~vJ>xv>=xMp7|5N6*mJ+R-y~X-DSKlYKDH{nXCprd8`??HOnzw$io Y=*d2~`;uPDJj?0!-MLyjdggh505T0s;Q#;t literal 0 HcmV?d00001 diff --git a/autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_divisions_defaults[True].npy b/autotest/__snapshots__/test_particledata/test_lrcparticledata_to_prp_divisions_defaults[True].npy new file mode 100644 index 0000000000000000000000000000000000000000..d52b8f7a9f58e1706e5f8c0d3163db61957a047f GIT binary patch literal 3152 zcmbW%y-pNi6vpw_&R9`USXfL9VUTQ$D9EP_pn#f4Y%sAR8(2bNgzOp`Lh{}NcY)Wz zeQZrdsJ1Z?_RQbO*__9MX-;rf4|lK(d*oO-#cobbnZS+&CRxo`4_F?ym;vKB>nfgseWJHU$_6} z?(n**jNR08^bo!u2$$AI1?PPmp8rR<$I{W$&KXp`=a_WtGrb#^j-Ii5DjhxXKJS@y z$LpZA1xw%0s&?+v|_K=PT)Kzxkj0yfx|Q8Rli}=*j)}dFwKdo}tUy(KB>uL*~)b_p9F%Svz`$ zF1?p|^bB3rj-HmwdyM%(=Fu~BSvz`$E^W#@dh#5c&)V7EpbA~ul6mwDUDl4CJO}T$ zWgb04m$jp3=+Z}-M^DS``}>o%qi5*SCz(gj&}HrD$$eznIiF>o{oXux=dyP6WS{vL DQ(TyL literal 0 HcmV?d00001 diff --git a/autotest/conftest.py b/autotest/conftest.py index 0db0d3ec08..f1106b3f64 100644 --- a/autotest/conftest.py +++ b/autotest/conftest.py @@ -1,10 +1,12 @@ import re from importlib import metadata +from io import BytesIO, StringIO from pathlib import Path from platform import system -from typing import List +from typing import List, Optional import matplotlib.pyplot as plt +import numpy as np import pytest from modflow_devtools.misc import is_in_ci @@ -158,3 +160,104 @@ def pytest_report_header(config): if not_found: lines.append("optional packages not found: " + ", ".join(not_found)) return "\n".join(lines) + + +# snapshot classes and fixtures + +from syrupy.extensions.single_file import ( + SingleFileSnapshotExtension, + WriteMode, +) +from syrupy.types import ( + PropertyFilter, + PropertyMatcher, + SerializableData, + SerializedData, +) + + +def _serialize_bytes(data): + buffer = BytesIO() + np.save(buffer, data) + return buffer.getvalue() + + +class BinaryArrayExtension(SingleFileSnapshotExtension): + """ + Binary snapshot of a NumPy array. Can be read back into NumPy with + .load(), preserving dtype and shape. This is the recommended array + snapshot approach if human-readability is not a necessity, as disk + space is minimized. + """ + + _write_mode = WriteMode.BINARY + _file_extension = "npy" + + def serialize( + self, + data, + *, + exclude=None, + include=None, + matcher=None, + ): + return _serialize_bytes(data) + + +class TextArrayExtension(SingleFileSnapshotExtension): + """ + Text snapshot of a NumPy array. Flattens the array before writing. + Can be read back into NumPy with .loadtxt() assuming you know the + shape of the expected data and subsequently reshape it if needed. + """ + + _write_mode = WriteMode.TEXT + _file_extension = "txt" + + def serialize( + self, + data: "SerializableData", + *, + exclude: Optional["PropertyFilter"] = None, + include: Optional["PropertyFilter"] = None, + matcher: Optional["PropertyMatcher"] = None, + ) -> "SerializedData": + buffer = StringIO() + np.savetxt(buffer, data.ravel()) + return buffer.getvalue() + + +class ReadableArrayExtension(SingleFileSnapshotExtension): + """ + Human-readable snapshot of a NumPy array. Preserves array shape + at the expense of possible loss of precision (default 8 places) + and more difficulty loading into NumPy than TextArrayExtension. + """ + + _write_mode = WriteMode.TEXT + _file_extension = "txt" + + def serialize( + self, + data: "SerializableData", + *, + exclude: Optional["PropertyFilter"] = None, + include: Optional["PropertyFilter"] = None, + matcher: Optional["PropertyMatcher"] = None, + ) -> "SerializedData": + return np.array2string(data, threshold=np.inf) + + +@pytest.fixture +def array_snapshot(snapshot): + return snapshot.use_extension(BinaryArrayExtension) + + +@pytest.fixture +def text_array_snapshot(snapshot): + return snapshot.use_extension(TextArrayExtension) + + +@pytest.fixture +def readable_array_snapshot(snapshot): + return snapshot.use_extension(ReadableArrayExtension) diff --git a/autotest/test_particledata.py b/autotest/test_particledata.py index 991518dd8b..ed1f167ab8 100644 --- a/autotest/test_particledata.py +++ b/autotest/test_particledata.py @@ -1,11 +1,12 @@ from functools import reduce -from itertools import chain, repeat +from itertools import chain +from pprint import pformat import matplotlib.pyplot as plt import numpy as np import pandas as pd import pytest -from modflow_devtools.markers import requires_pkg +from modflow_devtools.markers import requires_exe, requires_pkg import flopy from autotest.test_grid_cases import GridCases @@ -22,9 +23,11 @@ Modpath7Sim, NodeParticleData, ParticleData, + ParticleGroupLRCTemplate, ParticleGroupNodeTemplate, ) -from flopy.utils.modpathfile import PathlineFile +from flopy.modpath.mp7particlegroup import ParticleGroup +from flopy.utils.modpathfile import EndpointFile, PathlineFile # utilities @@ -44,7 +47,7 @@ def flatten(a): ] -# test constructors +# test initializers structured_dtype = np.dtype( @@ -388,7 +391,7 @@ def test_particledata_to_prp_disv_9(localz): @pytest.mark.parametrize( "localz", [False, True] ) # whether to return local z coords -def test_lrcparticledata_to_prp_divisions_defaults(localz): +def test_lrcparticledata_to_prp_divisions_defaults(localz, array_snapshot): sd_data = CellDataType() regions = [[0, 0, 1, 0, 1, 1]] part_data = LRCParticleData( @@ -396,66 +399,6 @@ def test_lrcparticledata_to_prp_divisions_defaults(localz): ) grid = GridCases().structured_small() rpts_prt = flatten(list(part_data.to_prp(grid, localz=localz))) - rpts_exp = [ - [0, 0, 0, 1, 1.166666, 1.166666, 5.833333], - [1, 0, 0, 1, 1.166666, 1.166666, 7.5], - [2, 0, 0, 1, 1.166666, 1.166666, 9.166666], - [3, 0, 0, 1, 1.1666666, 1.5, 5.833333], - [4, 0, 0, 1, 1.1666666, 1.5, 7.5], - [5, 0, 0, 1, 1.1666666, 1.5, 9.166666], - [6, 0, 0, 1, 1.166666, 1.833333, 5.833333], - [7, 0, 0, 1, 1.166666, 1.833333, 7.5], - [8, 0, 0, 1, 1.166666, 1.833333, 9.166666], - [9, 0, 0, 1, 1.5, 1.166666, 5.833333], - [10, 0, 0, 1, 1.5, 1.166666, 7.5], - [11, 0, 0, 1, 1.5, 1.166666, 9.166666], - [12, 0, 0, 1, 1.5, 1.5, 5.833333], - [13, 0, 0, 1, 1.5, 1.5, 7.5], - [14, 0, 0, 1, 1.5, 1.5, 9.166666], - [15, 0, 0, 1, 1.5, 1.833333, 5.833333], - [16, 0, 0, 1, 1.5, 1.833333, 7.5], - [17, 0, 0, 1, 1.5, 1.833333, 9.166666], - [18, 0, 0, 1, 1.833333, 1.166666, 5.833333], - [19, 0, 0, 1, 1.833333, 1.166666, 7.5], - [20, 0, 0, 1, 1.833333, 1.166666, 9.166666], - [21, 0, 0, 1, 1.833333, 1.5, 5.833333], - [22, 0, 0, 1, 1.833333, 1.5, 7.5], - [23, 0, 0, 1, 1.833333, 1.5, 9.166666], - [24, 0, 0, 1, 1.833333, 1.833333, 5.833333], - [25, 0, 0, 1, 1.833333, 1.833333, 7.5], - [26, 0, 0, 1, 1.833333, 1.833333, 9.166666], - [27, 0, 1, 1, 1.166666, 0.166666, 5.833333], - [28, 0, 1, 1, 1.166666, 0.166666, 7.5], - [29, 0, 1, 1, 1.166666, 0.166666, 9.166666], - [30, 0, 1, 1, 1.166666, 0.5, 5.833333], - [31, 0, 1, 1, 1.166666, 0.5, 7.5], - [32, 0, 1, 1, 1.166666, 0.5, 9.166666], - [33, 0, 1, 1, 1.166666, 0.833333, 5.833333], - [34, 0, 1, 1, 1.166666, 0.833333, 7.5], - [35, 0, 1, 1, 1.166666, 0.833333, 9.166666], - [36, 0, 1, 1, 1.5, 0.166666, 5.833333], - [37, 0, 1, 1, 1.5, 0.166666, 7.5], - [38, 0, 1, 1, 1.5, 0.166666, 9.166666], - [39, 0, 1, 1, 1.5, 0.5, 5.833333], - [40, 0, 1, 1, 1.5, 0.5, 7.5], - [41, 0, 1, 1, 1.5, 0.5, 9.166666], - [42, 0, 1, 1, 1.5, 0.833333, 5.833333], - [43, 0, 1, 1, 1.5, 0.833333, 7.5], - [44, 0, 1, 1, 1.5, 0.833333, 9.166666], - [45, 0, 1, 1, 1.833333, 0.166666, 5.833333], - [46, 0, 1, 1, 1.833333, 0.166666, 7.5], - [47, 0, 1, 1, 1.833333, 0.166666, 9.166666], - [48, 0, 1, 1, 1.833333, 0.5, 5.833333], - [49, 0, 1, 1, 1.833333, 0.5, 7.5], - [50, 0, 1, 1, 1.833333, 0.5, 9.166666], - [51, 0, 1, 1, 1.833333, 0.833333, 5.833333], - [52, 0, 1, 1, 1.833333, 0.833333, 7.5], - [53, 0, 1, 1, 1.833333, 0.833333, 9.166666], - ] - if localz: - locz = [0.166666, 0.5, 0.833333] * 18 - rpts_exp = [(*rpt[0:-1], z) for rpt, z in zip(rpts_exp, locz)] - num_cells = reduce( sum, [ @@ -473,7 +416,7 @@ def test_lrcparticledata_to_prp_divisions_defaults(localz): * sd_data.layercelldivisions ) assert act_len == exp_len - assert np.allclose(rpts_prt, rpts_exp) + assert rpts_prt == array_snapshot def test_lrcparticledata_to_prp_divisions_custom(): @@ -581,7 +524,7 @@ def test_lrcparticledata_to_prp_top_bottom(): assert rpts_prt[1][6] == grid.top_botm[0, 1, 1] -def test_lrcparticledata_to_prp_1_per_face(): +def test_lrcparticledata_to_prp_1_per_face(array_snapshot): sddata = FaceDataType( horizontaldivisions1=1, verticaldivisions1=1, @@ -600,15 +543,6 @@ def test_lrcparticledata_to_prp_1_per_face(): data = LRCParticleData(subdivisiondata=[sddata], lrcregions=[lrcregions]) grid = GridCases().structured_small() rpts_prt = flatten(list(data.to_prp(grid))) - rpts_exp = [ - # irpt, k, i, j, x, y, z - [0, 0, 1, 1, 1.0, 0.5, 7.5], - [1, 0, 1, 1, 2.0, 0.5, 7.5], - [2, 0, 1, 1, 1.5, 0.0, 7.5], - [3, 0, 1, 1, 1.5, 1.0, 7.5], - [4, 0, 1, 1, 1.5, 0.5, 5.0], - [5, 0, 1, 1, 1.5, 0.5, 10.0], - ] num_cells = len( [ (lrc[3] - lrc[0]) * (lrc[4] - lrc[1]) * (lrc[5] - lrc[2]) @@ -616,7 +550,7 @@ def test_lrcparticledata_to_prp_1_per_face(): ] ) assert len(rpts_prt) == num_cells * 6 # 1 particle on each face - assert np.allclose(rpts_prt, rpts_exp) + assert rpts_prt == array_snapshot @pytest.mark.parametrize( @@ -810,7 +744,7 @@ def test_nodeparticledata_prp_disv_big(function_tmpdir): print(rpts_prt) -# test write +# test write() def test_lrcparticledata_write(function_tmpdir): @@ -841,3 +775,102 @@ def test_lrcparticledata_write(function_tmpdir): # check lines written lines = open(p).readlines() assert lines == ["1 1\n", "2 1 0\n", " 5 5 1\n", "1 3 3 3 4 4 \n"] + + +# To make it easier to compare MODFLOW 6 PRT and MODPATH 7 results, +# we want to_coords() and to_prp() to return release configurations +# in the same order as is generated by MODPATH 7. That is, an input +# file for MF6 PRT's PRP package should list particle release points +# in the same order that MODPATH 7 assigns particle IDs upon release +# from any input style. The tests below set up bare-bones MP7 models +# in endpoint mode and compare against their release points. + + +@pytest.fixture +def mf6_sim(function_tmpdir): + name = "tiny-gwf" + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=function_tmpdir, 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=1, ncol=1) + 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]]) + 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")], + ) + return sim + + +def get_mp7_sim(mf6_sim, groups): + mp = Modpath7( + modelname=f"{mf6_sim.name}_mp7", + flowmodel=mf6_sim.get_model(), + exe_name="mp7", + model_ws=mf6_sim.sim_path, + ) + mpbas = Modpath7Bas(mp) + mpsim = Modpath7Sim( + mp, + simulationtype="endpoint", + trackingdirection="forward", + particlegroups=groups, + ) + return mp + + +@requires_exe("mp7") +def test_lrcparticledata_celldatatype_to_coords_order(mf6_sim): + mf6_sim.write_simulation() + success, buff = mf6_sim.run_simulation() + assert success, pformat(buff) + + pdata = flopy.modpath.LRCParticleData() + pg = ParticleGroupLRCTemplate(particlegroupname="PG1", particledata=pdata) + mp7_sim = get_mp7_sim(mf6_sim, [pg]) + mp7_sim.write_input() + success, buff = mp7_sim.run_model() + assert success, pformat(buff) + + gwf = mf6_sim.get_model() + grid = gwf.modelgrid + ep_file = EndpointFile(mf6_sim.sim_path / f"{mp7_sim.name}.mpend") + expected = ep_file.get_destination_endpoint_data(range(grid.nnodes))[ + ["x0", "y0", "z0"] + ].tolist() + actual = list(pdata.to_coords(grid)) + assert len(expected) == len(actual) == 27 + assert np.allclose(expected, actual) + + +def test_lrcparticledata_facedatatype_to_coords_order(mf6_sim): + mf6_sim.write_simulation() + success, buff = mf6_sim.run_simulation() + assert success, pformat(buff) + + pdata = flopy.modpath.LRCParticleData( + subdivisiondata=[FaceDataType()], + ) + pg = ParticleGroupLRCTemplate(particlegroupname="PG1", particledata=pdata) + mp7_sim = get_mp7_sim(mf6_sim, [pg]) + mp7_sim.write_input() + success, buff = mp7_sim.run_model() + assert success, pformat(buff) + + gwf = mf6_sim.get_model() + grid = gwf.modelgrid + ep_file = EndpointFile(mf6_sim.sim_path / f"{mp7_sim.name}.mpend") + expected = ep_file.get_destination_endpoint_data(range(grid.nnodes))[ + ["x0", "y0", "z0"] + ].tolist() + actual = list(pdata.to_coords(grid)) + assert len(expected) == len(actual) == 54 + assert np.allclose(expected, actual) diff --git a/etc/environment.yml b/etc/environment.yml index d2e71205e9..f29f1ca406 100644 --- a/etc/environment.yml +++ b/etc/environment.yml @@ -28,6 +28,7 @@ dependencies: - pytest-cov - pytest-dotenv - pytest-xdist + - syrupy - virtualenv # optional diff --git a/flopy/modpath/mp7particledata.py b/flopy/modpath/mp7particledata.py index af42b19906..f905219398 100644 --- a/flopy/modpath/mp7particledata.py +++ b/flopy/modpath/mp7particledata.py @@ -5,6 +5,7 @@ """ +from collections import namedtuple from itertools import product from typing import Iterator, Tuple @@ -15,6 +16,17 @@ from ..utils.recarray_utils import create_empty_recarray +def reversed_product(*iterables, repeat=1): + """ + Like `itertools.product()`, but left-most elements advance first. + + Adapted from https://stackoverflow.com/a/32998481/6514033. + """ + + for t in product(*reversed(iterables), repeat=repeat): + yield tuple(reversed(t)) + + class ParticleData: """ Class to create the most basic particle data type (starting location @@ -328,7 +340,7 @@ def __init__( def write(self, f=None): """ - Write the particle data to disk. + Write the particle data template to a file. Parameters ---------- @@ -695,8 +707,6 @@ def write(self, f=None): ) f.write(line) - return - class CellDataType: """ @@ -753,15 +763,13 @@ def __init__( def write(self, f=None): """ + Write the cell data template to a file. Parameters ---------- f : fileobject Fileobject that is open with write access - Returns - ------- - """ # validate that a valid file object was passed if not hasattr(f, "write"): @@ -781,17 +789,23 @@ def write(self, f=None): f.write(line) -def get_release_points( - subdivisiondata, grid, k=None, i=None, j=None, nn=None, localz=False -): - if nn is None and (k is None or i is None or j is None): - raise ValueError( - "A cell (node) must be specified by indices (for structured grids) or node number (for vertex/unstructured)" - ) +Extent = namedtuple( + "Extent", + [ + "minx", + "maxx", + "miny", + "maxy", + "minz", + "maxz", + "xspan", + "yspan", + "zspan", + ], +) - rpts = [] - template = [k, i, j] if nn is None else [nn] +def get_extent(grid, k=None, i=None, j=None, nn=None, localz=False) -> Extent: # get cell coords and span in each dimension if not (k is None or i is None or j is None): verts = grid.get_cell_vertices(i, j) @@ -824,140 +838,187 @@ def get_release_points( xspan = maxx - minx yspan = maxy - miny zspan = maxz - minz + return Extent(minx, maxx, miny, maxy, minz, maxz, xspan, yspan, zspan) - if isinstance(subdivisiondata, FaceDataType): - # x1 (west) - if ( - subdivisiondata.verticaldivisions1 > 0 - and subdivisiondata.horizontaldivisions1 > 0 - ): - yincr = yspan / subdivisiondata.horizontaldivisions1 - ylocs = [ - (miny + (yincr * 0.5) + (yincr * d)) - for d in range(subdivisiondata.horizontaldivisions1) - ] - zincr = zspan / subdivisiondata.verticaldivisions1 - zlocs = [ - (minz + (zincr * 0.5) + (zincr * d)) - for d in range(subdivisiondata.verticaldivisions1) - ] - prod = list(product(*[ylocs, zlocs])) - rpts = rpts + [template + [minx, p[0], p[1]] for p in prod] - # x2 (east) - if ( - subdivisiondata.verticaldivisions2 > 0 - and subdivisiondata.horizontaldivisions2 > 0 - ): - yincr = yspan / subdivisiondata.horizontaldivisions2 - ylocs = [ - (miny + (yincr * 0.5) + (yincr * d)) - for d in range(subdivisiondata.horizontaldivisions2) - ] - zincr = zspan / subdivisiondata.verticaldivisions2 - zlocs = [ - (minz + (zincr * 0.5) + (zincr * d)) - for d in range(subdivisiondata.verticaldivisions2) - ] - prod = list(product(*[ylocs, zlocs])) - rpts = rpts + [template + [maxx, p[0], p[1]] for p in prod] +def get_face_release_points( + subdivisiondata, cellid, extent +) -> Iterator[tuple]: + """ + Get release points for MODPATH 7 input style 2, template + subdivision style 1, i.e. face (2D) subdivision, for the + given cell with the given extent. + """ - # y1 (south) - if ( - subdivisiondata.verticaldivisions3 > 0 - and subdivisiondata.horizontaldivisions3 > 0 - ): - xincr = xspan / subdivisiondata.horizontaldivisions3 - xlocs = [ - (minx + (xincr * 0.5) + (xincr * rd)) - for rd in range(subdivisiondata.horizontaldivisions3) - ] - zincr = zspan / subdivisiondata.verticaldivisions3 - zlocs = [ - (minz + (zincr * 0.5) + (zincr * d)) - for d in range(subdivisiondata.verticaldivisions3) - ] - prod = list(product(*[xlocs, zlocs])) - rpts = rpts + [template + [p[0], miny, p[1]] for p in prod] + # Product incrementing left elements first, to + # match the release point ordering used by MP7 + product = reversed_product - # y2 (north) - if ( - subdivisiondata.verticaldivisions4 > 0 - and subdivisiondata.horizontaldivisions4 > 0 - ): - xincr = xspan / subdivisiondata.horizontaldivisions4 - xlocs = [ - (minx + (xincr * 0.5) + (xincr * rd)) - for rd in range(subdivisiondata.horizontaldivisions4) - ] - zincr = zspan / subdivisiondata.verticaldivisions4 - zlocs = [ - (minz + (zincr * 0.5) + (zincr * d)) - for d in range(subdivisiondata.verticaldivisions4) - ] - prod = list(product(*[xlocs, zlocs])) - rpts = rpts + [template + [p[0], maxy, p[1]] for p in prod] + # x1 (west) + if ( + subdivisiondata.verticaldivisions1 > 0 + and subdivisiondata.horizontaldivisions1 > 0 + ): + yincr = extent.yspan / subdivisiondata.horizontaldivisions1 + ylocs = [ + (extent.miny + (yincr * 0.5) + (yincr * d)) + for d in range(subdivisiondata.horizontaldivisions1) + ] + zincr = extent.zspan / subdivisiondata.verticaldivisions1 + zlocs = [ + (extent.minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions1) + ] + for p in product(*[ylocs, zlocs]): + yield cellid + [extent.minx, p[0], p[1]] - # z1 (bottom) - if ( - subdivisiondata.rowdivisions5 > 0 - and subdivisiondata.columndivisions5 > 0 - ): - xincr = xspan / subdivisiondata.columndivisions5 - xlocs = [ - (minx + (xincr * 0.5) + (xincr * rd)) - for rd in range(subdivisiondata.columndivisions5) - ] - yincr = yspan / subdivisiondata.rowdivisions5 - ylocs = [ - (miny + (yincr * 0.5) + (yincr * rd)) - for rd in range(subdivisiondata.rowdivisions5) - ] - prod = list(product(*[xlocs, ylocs])) - rpts = rpts + [template + [p[0], p[1], minz] for p in prod] + # x2 (east) + if ( + subdivisiondata.verticaldivisions2 > 0 + and subdivisiondata.horizontaldivisions2 > 0 + ): + yincr = extent.yspan / subdivisiondata.horizontaldivisions2 + ylocs = [ + (extent.miny + (yincr * 0.5) + (yincr * d)) + for d in range(subdivisiondata.horizontaldivisions2) + ] + zincr = extent.zspan / subdivisiondata.verticaldivisions2 + zlocs = [ + (extent.minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions2) + ] + for p in product(*[ylocs, zlocs]): + yield cellid + [extent.maxx, p[0], p[1]] - # z2 (top) - if ( - subdivisiondata.rowdivisions6 > 0 - and subdivisiondata.columndivisions6 > 0 - ): - xincr = xspan / subdivisiondata.columndivisions6 - xlocs = [ - (minx + (xincr * 0.5) + (xincr * rd)) - for rd in range(subdivisiondata.columndivisions6) - ] - yincr = yspan / subdivisiondata.rowdivisions6 - ylocs = [ - (miny + (yincr * 0.5) + (yincr * rd)) - for rd in range(subdivisiondata.rowdivisions6) - ] - prod = list(product(*[xlocs, ylocs])) - rpts = rpts + [template + [p[0], p[1], maxz] for p in prod] - elif isinstance(subdivisiondata, CellDataType): - xincr = xspan / subdivisiondata.columncelldivisions + # y1 (south) + if ( + subdivisiondata.verticaldivisions3 > 0 + and subdivisiondata.horizontaldivisions3 > 0 + ): + xincr = extent.xspan / subdivisiondata.horizontaldivisions3 xlocs = [ - (minx + (xincr * 0.5) + (xincr * rd)) - for rd in range(subdivisiondata.columncelldivisions) + (extent.minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.horizontaldivisions3) ] - yincr = yspan / subdivisiondata.rowcelldivisions - ylocs = [ - (miny + (yincr * 0.5) + (yincr * d)) - for d in range(subdivisiondata.rowcelldivisions) + zincr = extent.zspan / subdivisiondata.verticaldivisions3 + zlocs = [ + (extent.minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions3) + ] + for p in product(*[xlocs, zlocs]): + yield cellid + [p[0], extent.miny, p[1]] + + # y2 (north) + if ( + subdivisiondata.verticaldivisions4 > 0 + and subdivisiondata.horizontaldivisions4 > 0 + ): + xincr = extent.xspan / subdivisiondata.horizontaldivisions4 + xlocs = [ + (extent.minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.horizontaldivisions4) ] - zincr = zspan / subdivisiondata.layercelldivisions + zincr = extent.zspan / subdivisiondata.verticaldivisions4 zlocs = [ - (minz + (zincr * 0.5) + (zincr * d)) - for d in range(subdivisiondata.layercelldivisions) + (extent.minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions4) + ] + for p in product(*[xlocs, zlocs]): + yield cellid + [p[0], extent.maxy, p[1]] + + # z1 (bottom) + if ( + subdivisiondata.rowdivisions5 > 0 + and subdivisiondata.columndivisions5 > 0 + ): + xincr = extent.xspan / subdivisiondata.columndivisions5 + xlocs = [ + (extent.minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.columndivisions5) + ] + yincr = extent.yspan / subdivisiondata.rowdivisions5 + ylocs = [ + (extent.miny + (yincr * 0.5) + (yincr * rd)) + for rd in range(subdivisiondata.rowdivisions5) ] - prod = list(product(*[xlocs, ylocs, zlocs])) - rpts = rpts + [template + [p[0], p[1], p[2]] for p in prod] + for p in product(*[xlocs, ylocs]): + yield cellid + [p[0], p[1], extent.minz] + + # z2 (top) + if ( + subdivisiondata.rowdivisions6 > 0 + and subdivisiondata.columndivisions6 > 0 + ): + xincr = extent.xspan / subdivisiondata.columndivisions6 + xlocs = [ + (extent.minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.columndivisions6) + ] + yincr = extent.yspan / subdivisiondata.rowdivisions6 + ylocs = [ + (extent.miny + (yincr * 0.5) + (yincr * rd)) + for rd in range(subdivisiondata.rowdivisions6) + ] + for p in product(*[xlocs, ylocs]): + yield cellid + [p[0], p[1], extent.maxz] + + +def get_cell_release_points( + subdivisiondata, cellid, extent +) -> Iterator[tuple]: + """ + Get release points for MODPATH 7 input style 2, template + subdivision type 2, i.e. cell (3D) subdivision, for the + given cell with the given extent. + """ + + # Product incrementing left elements first, to + # match the release point ordering used by MP7 + product = reversed_product + + xincr = extent.xspan / subdivisiondata.columncelldivisions + xlocs = [ + (extent.minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.columncelldivisions) + ] + yincr = extent.yspan / subdivisiondata.rowcelldivisions + ylocs = [ + (extent.miny + (yincr * 0.5) + (yincr * d)) + for d in range(subdivisiondata.rowcelldivisions) + ] + zincr = extent.zspan / subdivisiondata.layercelldivisions + zlocs = [ + (extent.minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.layercelldivisions) + ] + for p in product(*[xlocs, ylocs, zlocs]): + yield cellid + [p[0], p[1], p[2]] + + +def get_release_points( + subdivisiondata, grid, k=None, i=None, j=None, nn=None, localz=False +) -> Iterator[tuple]: + """ + Get MODPATH 7 release point tuples for the given cell. + """ + + if nn is None and (k is None or i is None or j is None): + raise ValueError( + "A cell (node) must be specified by indices (for structured grids) or node number (for vertex/unstructured)" + ) + + cellid = [k, i, j] if nn is None else [nn] + extent = get_extent(grid, k, i, j, nn, localz) + + if isinstance(subdivisiondata, FaceDataType): + return get_face_release_points(subdivisiondata, cellid, extent) + elif isinstance(subdivisiondata, CellDataType): + return get_cell_release_points(subdivisiondata, cellid, extent) else: raise ValueError( f"Unsupported subdivision data type: {type(subdivisiondata)}" ) - return rpts - class LRCParticleData: """ @@ -993,7 +1054,7 @@ def __init__(self, subdivisiondata=None, lrcregions=None): subdivisiondata = CellDataType() if lrcregions is None: - lrcregions = [[0, 0, 0, 0, 0, 0]] + lrcregions = [[[0, 0, 0, 0, 0, 0]]] if isinstance(subdivisiondata, (CellDataType, FaceDataType)): subdivisiondata = [subdivisiondata] @@ -1049,7 +1110,6 @@ def __init__(self, subdivisiondata=None, lrcregions=None): "{} columns".format(self.name, shapel[1]) ) - # totalcellregioncount = 0 for lrcregion in lrcregions: totalcellregioncount += lrcregion.shape[0] @@ -1062,15 +1122,13 @@ def __init__(self, subdivisiondata=None, lrcregions=None): def write(self, f=None): """ + Write the layer-row-column particle data template to a file. Parameters ---------- f : fileobject Fileobject that is open with write access - Returns - ------- - """ # validate that a valid file object was passed if not hasattr(f, "write"): @@ -1294,15 +1352,13 @@ def __init__(self, subdivisiondata=None, nodes=None): def write(self, f=None): """ + Write the node particle data template to a file. Parameters ---------- f : fileobject Fileobject that is open with write access - Returns - ------- - """ # validate that a valid file object was passed if not hasattr(f, "write"): @@ -1354,10 +1410,9 @@ def to_coords(self, grid, localz=False) -> Iterator[tuple]: for sd in self.subdivisiondata: for nd in self.nodedata: - rpts = get_release_points( + for rpt in get_release_points( sd, grid, nn=int(nd[0]), localz=localz - ) - for rpt in rpts: + ): yield (*rpt[1:4],) def to_prp(self, grid, localz=False) -> Iterator[tuple]: @@ -1383,10 +1438,9 @@ def to_prp(self, grid, localz=False) -> Iterator[tuple]: for sd in self.subdivisiondata: for nd in self.nodedata: - rpts = get_release_points( - sd, grid, nn=int(nd[0]), localz=localz - ) - for irpt, rpt in enumerate(rpts): + for irpt, rpt in enumerate( + get_release_points(sd, grid, nn=int(nd[0]), localz=localz) + ): row = [irpt] if grid.grid_type == "structured": k, i, j = grid.get_lrc([rpt[0]])[0] diff --git a/pyproject.toml b/pyproject.toml index 3b0f58ba42..3d9e0607ec 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,6 +55,7 @@ test = [ "pytest-dotenv", "pytest-xdist", "pyzmq >=25.1.2", + "syrupy", "virtualenv" ] optional = [