diff --git a/.github/workflows/ci_examples.yaml b/.github/workflows/ci_examples.yaml index 0feec3db..fa098b99 100644 --- a/.github/workflows/ci_examples.yaml +++ b/.github/workflows/ci_examples.yaml @@ -50,6 +50,23 @@ jobs: cache-environment: true cache-downloads: true + - name: Checkout MODFLOW 6 + uses: actions/checkout@v4 + with: + repository: MODFLOW-USGS/modflow6 + ref: develop + path: modflow6 + + - name: Update flopy MODFLOW 6 classes + working-directory: modflow6/autotest + run: | + python update_flopy.py + + - name: Install mf6 nightly build binaries + uses: modflowpy/install-modflow-action@v1 + with: + repo: modflow6-nightly-build + - name: Install error reporter run: | pip install pytest-github-actions-annotate-failures diff --git a/autotest/test_mmr_to_mf6_dfw.py b/autotest/test_mmr_to_mf6_dfw.py index b61711f6..7c019cbb 100644 --- a/autotest/test_mmr_to_mf6_dfw.py +++ b/autotest/test_mmr_to_mf6_dfw.py @@ -15,14 +15,14 @@ # Getting the gis files causes problems in parallel in CI. See regression test. # See below to check if these answers are still up-to-date with -# mf6/autotest/test_swf_dfw.py +# mf6/autotest/test_chf_dfw.py # Note the answers are from the binary FLW files in -# mf6/autotest/test_swf_dfw.py, if you switch to text files, the line should be +# mf6/autotest/test_chf_dfw.py, if you switch to text files, the line should be # flw_list = [ # (int(binary), 100), # ] # one-based cell numbers here if binary, zero-based if text -answers_swf_dfw = { +answers_chf_dfw = { "ia": np.array([0, 2, 5, 7]), "ja": np.array([0, 1, 1, 0, 2, 2, 1], dtype=np.int32), "stage": np.array([[[[1.00196123, 1.00003366, 1.0]]]]), @@ -52,9 +52,9 @@ @pytest.mark.skipif(mf6_bin_unavailable, reason="mf6 binary not available") @pytest.mark.domainless @pytest.mark.parametrize("binary_flw", [True, False]) -def test_mmr_to_mf6_swf_dfw(tmp_path, binary_flw): +def test_mmr_to_mf6_chf_dfw(tmp_path, binary_flw): # The point of this test is to reproduce the - # modflow6/autotest/test_swf_dfw.py + # modflow6/autotest/test_chf_dfw.py # Here we supply "seg_mid_elevation" in the parameter data and # "stress_period_data" in chd options. This bypasses the calculation of @@ -68,8 +68,8 @@ def test_mmr_to_mf6_swf_dfw(tmp_path, binary_flw): # stream segments, the other is about the units of volume/flow being in # cubicfeet. - name = "swf-dfw01" - output_dir = tmp_path / "test_swf_dfw01" + name = "chf-dfw01" + output_dir = tmp_path / "test_chf_dfw01" save_flows = True print_flows = True @@ -123,16 +123,16 @@ def test_mmr_to_mf6_swf_dfw(tmp_path, binary_flw): # vertices could also be supplied by shapefiles vertices = [] vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] - cell2d = [] + cell1d = [] for j in range(nreach): - cell2d.append([j, 0.5, 2, j, j + 1]) - nodes = len(cell2d) + cell1d.append([j, 0.5, 2, j, j + 1]) + nodes = len(cell1d) nvert = len(vertices) disv1d_options = { "nodes": nodes, "nvert": nvert, "vertices": vertices, - "cell2d": cell2d, + "cell1d": cell1d, } # dfw @@ -246,18 +246,22 @@ def test_mmr_to_mf6_swf_dfw(tmp_path, binary_flw): # Checks assert success - # one can verify the answers match the current mf6 results for test_swf_dfw - # by placing the path to its output here + # one can verify the answers match the current mf6 results for test_chf_dfw + # by placing the path to its output here. Run the following lines to + # generate the mf6 results: + # cd modflow6_for_pws_ci/autotest # or your mf6 repo location + # pytest -s -vv test_chf_dfw.py --keep=keepers # output_dir = pl.Path( - # "../../modflow6/autotest/keepers/test_mf6model[0-swf-dfw01]0" + # "../../modflow6_for_pws_ci/autotest/keepers/" + # "test_mf6model[0-chf-dfw01]0" # ) # check binary grid file grb = flopy.mf6.utils.MfGrdFile(output_dir / f"{name}.disv1d.grb") ia = grb.ia ja = grb.ja - assert (answers_swf_dfw["ia"] == ia).all() - assert (answers_swf_dfw["ja"] == ja).all() + assert (answers_chf_dfw["ia"] == ia).all() + assert (answers_chf_dfw["ja"] == ja).all() assert ia.shape[0] == grb.nodes + 1, "ia in grb file is not correct size" # check stage file @@ -265,7 +269,7 @@ def test_mmr_to_mf6_swf_dfw(tmp_path, binary_flw): output_dir / f"{name}.stage", precision="double", text="STAGE" ) stage = qobj.get_alldata() - assert ((answers_swf_dfw["stage"] - stage) < 1.0e-7).all() + assert ((answers_chf_dfw["stage"] - stage) < 1.0e-7).all() # read the budget file budobj = flopy.utils.binaryfile.CellBudgetFile(output_dir / f"{name}.bud") @@ -275,17 +279,17 @@ def test_mmr_to_mf6_swf_dfw(tmp_path, binary_flw): qchd = np.array(budobj.get_data(text="CHD")[0].tolist()[0]) qresidual = np.zeros(grb.nodes)[0] - assert (answers_swf_dfw["flowja"] - flowja < 1.0e-7).all() - assert (answers_swf_dfw["qstorage"] - qstorage < 1.0e-7).all() - assert (answers_swf_dfw["qflw"] - qflw < 1.0e-7).all() - assert (answers_swf_dfw["qchd"] - qchd < 1.0e-7).all() - assert (answers_swf_dfw["qresidual"] - qresidual < 1.0e-7).all() + assert (answers_chf_dfw["flowja"] - flowja < 1.0e-7).all() + assert (answers_chf_dfw["qstorage"] - qstorage < 1.0e-7).all() + assert (answers_chf_dfw["qflw"] - qflw < 1.0e-7).all() + assert (answers_chf_dfw["qchd"] - qchd < 1.0e-7).all() + assert (answers_chf_dfw["qresidual"] - qresidual < 1.0e-7).all() # < answers_regression_means = { - "stage_all": 1.03667372881148, - "flow_all": 44.685014111989425, + "stage_all": 1.0372047024253908, + "flow_all": 44.70210250910929, } @@ -468,4 +472,5 @@ def get_outflow(itime): for kk, vv in answers_regression_means.items(): abs_diff = abs(locals()[kk].mean() - vv) + assert abs_diff < 1e-5, f"results for {kk} are not close" diff --git a/autotest/test_obsin_flow_node.py b/autotest/test_obsin_flow_node.py index 984d439a..39a905f3 100644 --- a/autotest/test_obsin_flow_node.py +++ b/autotest/test_obsin_flow_node.py @@ -1,6 +1,6 @@ import numpy as np import pytest -from pyPRMS.DataFile import DataFile as PRMSStreamflowData +from pyPRMS import DataFile as PRMSStreamflowData from pywatershed import PRMSChannel from pywatershed.base.adapter import Adapter, AdapterNetcdf, adapter_factory diff --git a/examples/06_flow_graph_starfit.ipynb b/examples/06_flow_graph_starfit.ipynb index 8c4a4dd4..16b02643 100644 --- a/examples/06_flow_graph_starfit.ipynb +++ b/examples/06_flow_graph_starfit.ipynb @@ -48,7 +48,11 @@ "plot_height = 600\n", "plot_width = 1000\n", "\n", - "jupyter_black.load()" + "jupyter_black.load()\n", + "\n", + "# to remove:\n", + "import warnings\n", + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -292,26 +296,8 @@ " .sel(nhm_seg=44438)\n", " .rename(\"modeled\")\n", " .to_dataframe()[\"modeled\"]\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3b385029-6554-4f4d-9fd9-4c255d70e40d", - "metadata": {}, - "outputs": [], - "source": [ - "obs_all = pyPRMS.DataFile(domain_dir / \"sf_data\").data_by_variable(\"runoff\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "54058459-fed5-4cc3-bed2-8a3edfb22a13", - "metadata": {}, - "outputs": [], - "source": [ + ")\n", + "obs_all = pyPRMS.DataFile(domain_dir / \"sf_data\").data_by_variable(\"runoff\")\n", "wh_poi_obs = np.where(params.parameters[\"poi_gage_segment\"] == 184)\n", "gage_id = params.parameters[\"poi_gage_id\"][wh_poi_obs][0]\n", "obs = obs_all[f\"runoff_{gage_id}\"]\n", diff --git a/examples/mmr_to_mf6_dfw.ipynb b/examples/07_mmr_to_mf6_chf_dfw.ipynb similarity index 96% rename from examples/mmr_to_mf6_dfw.ipynb rename to examples/07_mmr_to_mf6_chf_dfw.ipynb index 61495c25..f1ad7bc9 100644 --- a/examples/mmr_to_mf6_dfw.ipynb +++ b/examples/07_mmr_to_mf6_chf_dfw.ipynb @@ -7,9 +7,9 @@ "source": [ "# MMR To MF6 DFW\n", "\n", - "This notebook performs MF6 diffussive wave routing (DFW) simulaton based on PRMS MMR (Muskingum-Mann Routing) inputs and boundary flows from the PRMS model. We'll compare and contrast the MF6 DFW and PRMS MMR simulations and delve into their differences at the largest gaged flows in the domain (that are not affected by tides). \n", + "This notebook performs 1-D diffussive wave (DFW) routing in MODFLOW 6 using the CHF (channel flow) model. The static and time-varying boundary conditions, come from PRMS. Specifically, PRMS's MMR (Muskingum-Mann Routing) parameters and its channel inflows calculated during a pywatershed run. We'll compare and contrast the MF6 CHF-DFW and PRMS MMR simulations and delve into their differences at the largest gaged flows in the domain (that are not affected by tides). \n", "\n", - "This notebook requries the develop branch of MF6 and a flopy upto date with this branch by running `python update_flopy.py` in modflow6/autotest. You will also need to add the MF6 executable to your path below. You will also need to have run \n", + "This notebook requries the develop branch of MF6 and a flopy upto date with this branch by running `python update_flopy.py` in modflow6/autotest. You will also need to add the MF6 executable to your path below. \n", "\n", "## User configuration" ] @@ -21,16 +21,20 @@ "metadata": {}, "outputs": [], "source": [ - "# Set YOUR path to MF6 here\n", + "# Set YOUR path to MF6 in this block\n", "import pathlib as pl\n", "\n", "mf6_bin = pl.Path(\"../../modflow6/bin/mf6\")\n", + "# double check\n", "msg = \"A build of mf6/develop branch required for DFW simulation\"\n", "assert mf6_bin.exists, msg\n", "\n", "# Rerun MF6 model below or use an existing run?\n", "rerun_mf6 = True\n", - "rerun_prms = True" + "rerun_prms = True\n", + "\n", + "# Perform the full, 2 year run period or just the first 45 days?\n", + "full_run_period = False" ] }, { @@ -96,7 +100,7 @@ "id": "e6a36353-3a61-4b57-8db5-e48293ffc0d1", "metadata": {}, "source": [ - "## Run PRMS\n", + "## Run PRMS NHM configuration using pywatershed\n", "Running PRMS gives the boundary conditions for the MF6 DFW run and it also produces its own streamflow simulation with its Muskingum-Mann routing method. " ] }, @@ -107,11 +111,9 @@ "metadata": {}, "outputs": [], "source": [ - "%time\n", - "\n", - "prms_run_dir = repo_root_dir / \"examples/mmr_to_mf6_dfw/prms_run\"\n", + "prms_run_dir = repo_root_dir / \"examples/07_mmr_to_mf6_chf_dfw/prms_run\"\n", "\n", - "if rerun_prms:\n", + "if rerun_prms and prms_run_dir.exists():\n", " shutil.rmtree(prms_run_dir)\n", "\n", "if not prms_run_dir.exists():\n", @@ -186,12 +188,10 @@ "control_file = domain_dir / \"nhm.control\"\n", "control = pws.Control.load_prms(control_file)\n", "ndays_run = control.n_times\n", - "# Could shorten the run duration\n", - "# Subtract one from ndays_run becase end day/time would be included in the PRMS run\n", - "# ndays_run = 45\n", - "# control.edit_end_time(\n", - "# control.start_time + ((ndays_run - 1) * control.time_step)\n", - "# )" + "\n", + "if not full_run_period:\n", + " ndays_run = 45\n", + " control.edit_n_time_steps(ndays_run)" ] }, { @@ -482,7 +482,10 @@ "metadata": {}, "outputs": [], "source": [ - "tt = 3000\n", + "if ndays_run < 3000:\n", + " tt = ndays_run\n", + "else:\n", + " tt = 3000\n", "zoom = 1.0\n", "figsize = (7 * zoom, 10 * zoom)\n", "dt = tdis_perlen / tdis_nstp\n", diff --git a/pywatershed/utils/mmr_to_mf6_dfw.py b/pywatershed/utils/mmr_to_mf6_dfw.py index 7f14318a..024d2722 100644 --- a/pywatershed/utils/mmr_to_mf6_dfw.py +++ b/pywatershed/utils/mmr_to_mf6_dfw.py @@ -67,7 +67,7 @@ class MmrToMf6Dfw: Units of meters from metadata (these are converted to units requested for modflow, which we currently hardcode to meters and seconds.) * seg_slope: - Passed directly to SWF DFW. Also used to calculate seg_mid_elevation if + Passed directly to CHF DFW. Also used to calculate seg_mid_elevation if not present in parameters. See its description below. Unitless slope. * hru_segment: Used to calculate seg_mid_elevation if not present in parameters. See @@ -81,7 +81,7 @@ class MmrToMf6Dfw: This is (apparently) the NHD bank-full width present in the PRMS parameter files but unused in PRMS. Units are in meters. * mann_n: - Mannings N coefficient for MMR. Passed directly to SWF DFW. Units in + Mannings N coefficient for MMR. Passed directly to CHF DFW. Units in seconds / meter ** (1/3) according to the metadata. The following optional parameters can be user-supplied but are NOT part of @@ -246,7 +246,7 @@ def __init__( self._handle_control_parameters() self._set_sim() self._set_tdis() - self._set_swf() + self._set_chf() self._set_disv1d() self._set_flw() self._set_ims() @@ -389,8 +389,8 @@ def _set_tdis(self): perioddata=tdis_rc, ) - def _set_swf(self): - self._swf = flopy.mf6.ModflowSwf( + def _set_chf(self): + self._chf = flopy.mf6.ModflowChf( self._sim, modelname=self._sim_name, save_flows=self._save_flows ) @@ -405,16 +405,6 @@ def _set_disv1d(self): if "idomain" not in self._disv1d_options.keys(): self._disv1d_options["idomain"] = 1 - if "length" not in self._disv1d_options.keys(): - # meters - # unit-ed quantities - segment_units = self._units(meta.parameters["seg_length"]["units"]) - self._segment_length = parameters["seg_length"] - self._segment_length = self._segment_length * segment_units - self._disv1d_options["length"] = ( - self._segment_length.to_base_units().magnitude - ) - if "width" not in self._disv1d_options.keys(): # meters, per metadata self._disv1d_options["width"] = parameters["seg_width"] @@ -432,13 +422,13 @@ def _set_disv1d(self): # < self._disv1d_options["bottom"] = self._seg_mid_elevation - # < vertices and cell2d + # < vertices and cell1d if self._segment_shp_file is not None: opts_set = [ "nodes", "nvert", "vertices", - "cell2d", + "cell1d", ] for oo in opts_set: self._warn_option_overwrite(oo, opt_dict_name, method_name) @@ -486,9 +476,9 @@ def _set_disv1d(self): for vv in range(len(inds)): vertices += [[inds[vv], geoms[vv][0], geoms[vv][1]]] - # cell2d: for each reach, if it has a "to", get the to's first ind + # cell1d: for each reach, if it has a "to", get the to's first ind # as the last vertex ind for the reach - cell2d = [] + cell1d = [] for rr in range(nreach): reach_id = parameters["nhm_seg"][rr] icvert = reach_vert_inds[reach_id] @@ -498,7 +488,7 @@ def _set_disv1d(self): icvert += [to_reach_start_ind] # < ncvert = len(icvert) - cell2d += [[rr, 0.5, ncvert, *icvert]] + cell1d += [[rr, 0.5, ncvert, *icvert]] # < Just a check check_connectivity = True @@ -507,8 +497,8 @@ def _set_disv1d(self): reach_id = parameters["nhm_seg"][rr] reach_ind = rr - reach_cell2d_start_ind = cell2d[reach_ind][3] - reach_cell2d_end_ind = cell2d[reach_ind][-1] + reach_cell1d_start_ind = cell1d[reach_ind][3] + reach_cell1d_end_ind = cell1d[reach_ind][-1] wh_geom = np.where( segment_gdf.nsegment_v.values == reach_id @@ -524,9 +514,9 @@ def _set_disv1d(self): to_reach_ind = np.where( parameters["nhm_seg"] == to_reach_id )[0][0] - to_reach_cell2d_start_ind = cell2d[to_reach_ind][3] + to_reach_cell1d_start_ind = cell1d[to_reach_ind][3] assert ( - reach_cell2d_end_ind == to_reach_cell2d_start_ind + reach_cell1d_end_ind == to_reach_cell1d_start_ind ) # wh_down_geom = np.where( @@ -548,9 +538,9 @@ def _set_disv1d(self): from_reach_ind = np.where( parameters["nhm_seg"] == from_reach_id )[0][0] - from_reach_cell2d_end_ind = cell2d[from_reach_ind][-1] + from_reach_cell1d_end_ind = cell1d[from_reach_ind][-1] assert ( - from_reach_cell2d_end_ind == reach_cell2d_start_ind + from_reach_cell1d_end_ind == reach_cell1d_start_ind ) # wh_up_geom = np.where( @@ -561,7 +551,7 @@ def _set_disv1d(self): # ).tolist() # assert reach_up_geom[-1] == reach_geom[0] - nnodes = len(cell2d) + nnodes = len(cell1d) assert nnodes == self._nsegment nvert = len(vertices) @@ -571,11 +561,11 @@ def _set_disv1d(self): self._disv1d_options["nodes"] = self._nsegment self._disv1d_options["nvert"] = nvert self._disv1d_options["vertices"] = vertices - self._disv1d_options["cell2d"] = cell2d + self._disv1d_options["cell1d"] = cell1d # < - self._disv1d = flopy.mf6.ModflowSwfdisv1D( - self._swf, **self._disv1d_options + self._disv1d = flopy.mf6.ModflowChfdisv1D( + self._chf, **self._disv1d_options ) def _set_flw(self, **kwargs): @@ -586,7 +576,7 @@ def _set_flw(self, **kwargs): # aggregate inflows over the contributing fluxes # For non-binary data, this method has to be called after - # flopy.mf6.Swfdisl is set on self._swf + # flopy.mf6.Chfdisl is set on self._chf parameters = self.parameters.parameters @@ -706,7 +696,7 @@ def read_inflow(vv, start_time, end_time): self.control.start_time + ispd * self.control.time_step ) bin_name = ( - f"swf_flw_bc/" + f"chf_flw_bc/" f"flw_{flow_name}_{i_time_str.replace(':', '_')}.bin" ) bin_name_pl = self._output_dir / bin_name @@ -727,8 +717,8 @@ def read_inflow(vv, start_time, end_time): if key not in self._flw_options.keys(): self._flw_options[key] = True - _ = flopy.mf6.ModflowSwfflw( - self._swf, + _ = flopy.mf6.ModflowChfflw( + self._chf, save_flows=self._save_flows, stress_period_data=flw_spd, maxbound=self._nsegment, @@ -752,13 +742,13 @@ def _set_dfw(self): # seconds / meter ** (1/3) self._dfw_options["manningsn"] = parameters["mann_n"] - _ = flopy.mf6.ModflowSwfdfw(self._swf, **self._dfw_options) + _ = flopy.mf6.ModflowChfdfw(self._chf, **self._dfw_options) def _set_sto(self): if "save_flows" not in self._sto_options: self._sto_options["save_flows"] = self._save_flows - _ = flopy.mf6.ModflowSwfsto( - self._swf, + _ = flopy.mf6.ModflowChfsto( + self._chf, **self._sto_options, ) @@ -771,20 +761,20 @@ def _set_ic(self): "strt": self._seg_mid_elevation + self._ic_options["strt"] } - self._ic = flopy.mf6.ModflowSwfic(self._swf, **ic_options) + self._ic = flopy.mf6.ModflowChfic(self._chf, **ic_options) def _set_cxs(self): if self._cxs_options is None: pass else: - self._cxs = flopy.mf6.ModflowSwfcxs(self._swf, **self._cxs_options) + self._cxs = flopy.mf6.ModflowChfcxs(self._chf, **self._cxs_options) def _set_oc(self): if self._oc_options is None: self._oc_options = {} - self._oc = flopy.mf6.ModflowSwfoc( - self._swf, + self._oc = flopy.mf6.ModflowChfoc( + self._chf, budget_filerecord=f"{self._sim_name}.bud", stage_filerecord=f"{self._sim_name}.stage", # qoutflow_filerecord=f"{self._sim_name}.qoutflow", @@ -814,7 +804,7 @@ def _set_chd(self): self._chd_options["maxbound"] = len( self._chd_options["stress_period_data"] ) - self._chd = flopy.mf6.ModflowSwfchd(self._swf, **self._chd_options) + self._chd = flopy.mf6.ModflowChfchd(self._chf, **self._chd_options) def _calculate_seg_mid_elevations(self, check=False): nseg = self._nsegment