From 623091b13a3f5052a6b31d0bc6b3209a7dcce0ad Mon Sep 17 00:00:00 2001 From: "Leaf, Andrew T" Date: Thu, 13 May 2021 12:57:39 -0500 Subject: [PATCH] refactor: harmonize OC setup between MF6 and MFNWT; add support for empty periods to turn off output writing --- docs/source/concepts/index.rst | 1 - docs/source/index.rst | 3 +- docs/source/input/index.rst | 7 + docs/source/input/oc.rst | 48 +++++++ mfsetup/mf6_defaults.yml | 2 +- mfsetup/mf6model.py | 23 +--- mfsetup/mfnwt_defaults.yml | 4 +- mfsetup/mfnwtmodel.py | 17 ++- mfsetup/oc.py | 120 ++++++++++++++++++ mfsetup/tdis.py | 2 +- mfsetup/tests/data/shellmound.yml | 2 +- mfsetup/tests/test_mf6_shellmound.py | 27 ++-- mfsetup/tests/test_mf6_shellmound_rot_grid.py | 2 +- mfsetup/tests/test_oc.py | 64 ++++++++++ mfsetup/tests/test_pfl_mfnwt_inset.py | 20 ++- mfsetup/tests/test_pleasant_mfnwt_inset.py | 6 + 16 files changed, 302 insertions(+), 46 deletions(-) create mode 100644 docs/source/input/index.rst create mode 100644 docs/source/input/oc.rst create mode 100644 mfsetup/oc.py create mode 100644 mfsetup/tests/test_oc.py diff --git a/docs/source/concepts/index.rst b/docs/source/concepts/index.rst index 7461b8b5..a69951a5 100644 --- a/docs/source/concepts/index.rst +++ b/docs/source/concepts/index.rst @@ -4,6 +4,5 @@ Modflow-setup concepts and methods .. toctree:: :maxdepth: 1 - :caption: Concepts and methods Specifying perimeter boundary conditions diff --git a/docs/source/index.rst b/docs/source/index.rst index c360deb0..a6df98ab 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -27,7 +27,8 @@ Modflow-setup is intended to facilitate automated setup of MODFLOW models, from Basic program structure and usage The configuration file - Concepts and Methods + Concepts and methods + Input instructions by package Troubleshooting diff --git a/docs/source/input/index.rst b/docs/source/input/index.rst new file mode 100644 index 00000000..dccb5339 --- /dev/null +++ b/docs/source/input/index.rst @@ -0,0 +1,7 @@ +=========================================================== +Input instructions by package +=========================================================== +.. toctree:: + :maxdepth: 1 + + Output Control diff --git a/docs/source/input/oc.rst b/docs/source/input/oc.rst new file mode 100644 index 00000000..075a3e6c --- /dev/null +++ b/docs/source/input/oc.rst @@ -0,0 +1,48 @@ +=========================================================== +MODFLOW Output Control +=========================================================== + +Stress period input format +-------------------------- +Regardless of the model version (MODFLOW-2005-style or MODFLOW 6), output control can be specified in a format similar to native MODFLOW 6 input: + +.. code-block:: yaml + + oc: + period_options: + 0: ['save head last', 'save budget last'] + 10: [] + 15: ['save head last', 'save budget last'] + +The above ``period_options:`` block would save the head and cell budget output on the last timestep of stress periods 0 through 9, and from 15 on, but turn off output saving for stress periods 10 through 14. This behavior is consistent with MODFLOW 6 but differs from MODFLOW-2005, where each stress periods must be explicitly included for output to be written. Other options besides ``'last'`` include ``all, first, frequency , and steps ``; see the MODFLOW 6 input instructions for more details. + +Output filenames and other arguments +------------------------------------ +For MODFLOW 6 models, the ``head_fileout_fmt`` and ``budget_fileout_fmt`` arguments can also be supplied to tell Flopy where to save the head and cell budget files, and how to name them. Modflow-setup fills any format specifiers (``'{}'``) with the model name, and passes the resulting strings to the ``head_filerecord`` and ``budget_filerecord`` arguments to the :py:class:`flopy.mf6.ModflowGwfoc ` constructor. + +.. code-block:: yaml + + oc: + head_fileout_fmt: '{}.hds' + budget_fileout_fmt: '{}.cbc' + period_options: + 0: ['save head last', 'save budget last'] + +Any other valid arguments to the :py:class:`flopy.mf6.ModflowGwfoc ` and :py:class:`flopy.modflow.ModflowOc ` constructors can be supplied as keys in the ``oc:`` dictionary block. For example: + + .. code-block:: yaml + + oc: + unitnumber: [14, 51, 52, 53, 0] + +would set the unit numbers for the head, drawdown, budget, and ibound output files. See the Flopy documentation for more details. Invalid arguments are filtered out prior to calling the constructor. + +Alternative stress period input formats +---------------------------------------- +As with other arguments, stress period input can also be directly specified in the Flopy input formats. For example, ``stress_period_data`` could be supplied for a MODFLOW-2005 model as it would be supplied to the :py:class:`flopy.modflow.ModflowOc ` constructor: + +.. code-block:: yaml + + oc: + stress_period_data: + (0, 1): ['save head', 'save budget'] diff --git a/mfsetup/mf6_defaults.yml b/mfsetup/mf6_defaults.yml index 087e1653..aa84f1d9 100644 --- a/mfsetup/mf6_defaults.yml +++ b/mfsetup/mf6_defaults.yml @@ -183,7 +183,7 @@ chd: oc: head_fileout_fmt: '{}.hds' budget_fileout_fmt: '{}.cbc' - # example of using MODFLOW 6 text input + # example of using MODFLOW 6-style text input period_options: {0: ['save head last', 'save budget last'] } diff --git a/mfsetup/mf6model.py b/mfsetup/mf6model.py index 511215ba..99ef9443 100644 --- a/mfsetup/mf6model.py +++ b/mfsetup/mf6model.py @@ -40,6 +40,7 @@ from mfsetup.mfmodel import MFsetupMixin from mfsetup.mover import get_mover_sfr_package_input from mfsetup.obs import setup_head_observations +from mfsetup.oc import parse_oc_period_input from mfsetup.tdis import add_date_comments_to_tdis from mfsetup.tmr import TmrNew from mfsetup.units import convert_time_units @@ -911,25 +912,9 @@ def setup_oc(self): kwargs = self.cfg[package] kwargs['budget_filerecord'] = self.cfg[package]['budget_fileout_fmt'].format(self.name) kwargs['head_filerecord'] = self.cfg[package]['head_fileout_fmt'].format(self.name) - # parse both flopy and mf6-style input into flopy input - for rec in ['printrecord', 'saverecord']: - if rec in kwargs: - data = kwargs[rec] - mf6_input = {} - for kper, words in data.items(): - mf6_input[kper] = [] - for var, instruction in words.items(): - mf6_input[kper].append((var, instruction)) - kwargs[rec] = mf6_input - elif 'period_options' in kwargs: - mf6_input = defaultdict(list) - for kper, options in kwargs['period_options'].items(): - for words in options: - type, var, instruction = words.split() - if type == rec.replace('record', ''): - mf6_input[kper].append((var, instruction)) - if len(mf6_input) > 0: - kwargs[rec] = mf6_input + + period_input = parse_oc_period_input(kwargs) + kwargs.update(period_input) kwargs = get_input_arguments(kwargs, mf6.ModflowGwfoc) oc = mf6.ModflowGwfoc(self, **kwargs) diff --git a/mfsetup/mfnwt_defaults.yml b/mfsetup/mfnwt_defaults.yml index dd371b7c..12cd7ccb 100644 --- a/mfsetup/mfnwt_defaults.yml +++ b/mfsetup/mfnwt_defaults.yml @@ -123,8 +123,8 @@ mnw: oc: head_fileout_fmt: '{}.hds' budget_fileout_fmt: '{}.cbc' - period_options: {0: ['save head', - 'save budget'] + period_options: {0: ['save head last', + 'save budget last'] } hyd: diff --git a/mfsetup/mfnwtmodel.py b/mfsetup/mfnwtmodel.py index 11ada547..b7cea264 100644 --- a/mfsetup/mfnwtmodel.py +++ b/mfsetup/mfnwtmodel.py @@ -38,6 +38,7 @@ ) from mfsetup.mfmodel import MFsetupMixin from mfsetup.obs import read_observation_data, setup_head_observations +from mfsetup.oc import parse_oc_period_input from mfsetup.tdis import get_parent_stress_periods, setup_perioddata_group from mfsetup.tmr import TmrNew from mfsetup.units import convert_length_units, itmuni_text, lenuni_text @@ -324,12 +325,18 @@ def setup_oc(self): package = 'oc' print('\nSetting up {} package...'.format(package.upper())) t0 = time.time() - stress_period_data = {} - for i, r in self.perioddata.iterrows(): - stress_period_data[(r.per, r.nstp -1)] = r.oc - + #stress_period_data = {} + #for i, r in self.perioddata.iterrows(): + # stress_period_data[(r.per, r.nstp -1)] = r.oc + + # use stress_period_data if supplied + # (instead of period_input defaults) + if 'stress_period_data' in self.cfg['oc']: + del self.cfg['oc']['period_options'] kwargs = self.cfg['oc'] - kwargs['stress_period_data'] = stress_period_data + period_input = parse_oc_period_input(kwargs, nstp=self.perioddata.nstp, + output_fmt='mfnwt') + kwargs.update(period_input) kwargs = get_input_arguments(kwargs, fm.ModflowOc) oc = fm.ModflowOc(model=self, **kwargs) print("finished in {:.2f}s\n".format(time.time() - t0)) diff --git a/mfsetup/oc.py b/mfsetup/oc.py new file mode 100644 index 00000000..b913f7a2 --- /dev/null +++ b/mfsetup/oc.py @@ -0,0 +1,120 @@ +"""Functions for handling MODFLOW output control +""" +from collections import defaultdict + + +def parse_oc_period_input(period_input, nstp=None, output_fmt='mf6'): + """Parse both flopy and mf6-style stress period output control input + into flopy input. + + Parameters + ---------- + period_input : dict + Dictionary of stress period input (see examples) + nstp : list-like + Number of timesteps in each stress period + output_fmt : str + 'mf6' for MODFLOW 6 style input (to :py:func:`flopy.mf6.ModflowGwfoc`), otherwise, + input for :py:func:`flopy.modflow.ModflowOc` is produced. + + Returns + ------- + flopy_input : dict + Input to the flopy output control package constructor. + + Examples + -------- + >>> period_input = {'saverecord': {0: {'head': 'last', 'budget': 'last'}} + {0: [('head', 'last'), ('budget', 'last')]} + """ + if nstp is not None: + nstp = list(nstp) + + flopy_input = {} + mf6_flopy_input = {} + for rec in ['printrecord', 'saverecord']: + if rec in period_input: + if output_fmt != 'mf6': + msg = ("MODFLOW 6 Flopy-style OC input (printrecord or " + "saverecord arguments) only supported for MODFLOW 6 models.") + raise NotImplementedError(msg) + data = period_input[rec] + mf6_record_input = {} + for kper, words in data.items(): + mf6_record_input[kper] = [] + for var, instruction in words.items(): + mf6_record_input[kper].append((var, instruction)) + mf6_flopy_input[rec] = mf6_record_input + elif 'period_options' in period_input: + mf6_record_input = defaultdict(list) + mf_record_input = defaultdict(list) + for kper, options in period_input['period_options'].items(): + # empty period for turning off output + if len(options) == 0: + mf6_record_input[kper] = [] + mf_record_input[(kper, 0)] = [] + else: + for words in options: + type, var, *instruction = words.split() + if type == rec.replace('record', ''): + if output_fmt == 'mf6': + mf6_record_input[kper].append((var, *instruction)) + else: + if nstp is None: + raise ValueError("MODFLOW 2005-style OC input requires " + "timestep information.") + # parse MF6-style instructions + kstp = 0 + nstep_idx = kper if kper < len(nstp) else -1 + instruction, *values = instruction + instruction = instruction.lower() + if instruction == 'all': + for kstp in range(nstp[nstep_idx]): + mf_record_input[(kper, kstp)].append(f"{type} {var}") + elif 'frequency' in instruction: + if len(values) == 0: + raise ValueError("mfsetup.oc.parse_oc: " + "'frequency' instruction needs a value") + freq = int(values[0]) + steps = list(range(nstp[nstep_idx]))[::freq] + for kstp in steps: + mf_record_input[(kper, kstp)].append(f"{type} {var}") + elif 'steps' in instruction: + if len(values) == 0: + raise ValueError("mfsetup.oc.parse_oc: " + "'steps' instruction needs one or more values") + for kstp in values: + mf_record_input[(kper, int(kstp))].append(f"{type} {var}") + elif instruction == 'first': + mf_record_input[(kper, 0)].append(f"{type} {var}") + elif instruction == 'last': + kstp = nstp[nstep_idx] - 1 + mf_record_input[(kper, int(kstp))].append(f"{type} {var}") + else: + raise ValueError("mfsetup.oc.parse_oc: instruction " + f"'{instruction}' not understood") + if len(mf6_record_input) > 0: + mf6_flopy_input[rec] = dict(mf6_record_input) + if len(mf_record_input) > 0: + mf_record_input = fill_oc_stress_period_data(mf_record_input, nper=len(nstp)) + flopy_input['stress_period_data'] = dict(mf_record_input) + if output_fmt == 'mf6': + return mf6_flopy_input + return flopy_input + + +def fill_oc_stress_period_data(stress_period_data, nper): + """For MODFLOW 2005-style models, repeat last entry in stress_period_data + for subsequent stress periods (until another entry is encountered), + as is done by default in MODFLOW 6. + """ + filled_spd = {} + last_period_data = {} + for period in range(nper): + for (kper, kstp), data in stress_period_data.items(): + if kper == period: + last_period_data[(kper, kstp)] = data + last_period_data = {(period, kstp): data for (kper, kstp), data + in last_period_data.items()} + filled_spd.update(last_period_data) + return filled_spd diff --git a/mfsetup/tdis.py b/mfsetup/tdis.py index b461274b..ad2ea2e3 100644 --- a/mfsetup/tdis.py +++ b/mfsetup/tdis.py @@ -169,7 +169,7 @@ def setup_perioddata_group(start_date_time, end_date_time=None, oc_saverecord : dict Dictionary with zero-based stress periods as keys and output control options as values. Similar to MODFLOW-6 input, the information specified for a period will - continue to apply until information for another perior is specified. + continue to apply until information for another period is specified. Returns ------- diff --git a/mfsetup/tests/data/shellmound.yml b/mfsetup/tests/data/shellmound.yml index 8411d1a4..04a042e6 100644 --- a/mfsetup/tests/data/shellmound.yml +++ b/mfsetup/tests/data/shellmound.yml @@ -318,7 +318,7 @@ oc: # MODFLOW 6-style text input can also be used # e.g. # period_options: {0: ['save head last', - # 'save budget last' ] + # 'save budget last' ] obs: source_data: diff --git a/mfsetup/tests/test_mf6_shellmound.py b/mfsetup/tests/test_mf6_shellmound.py index 6806602b..b44f233b 100644 --- a/mfsetup/tests/test_mf6_shellmound.py +++ b/mfsetup/tests/test_mf6_shellmound.py @@ -487,15 +487,19 @@ def test_obs_setup(shellmound_model_with_dis, config): break -@pytest.mark.parametrize('options', [{'saverecord': {0: {'head': 'last', - 'budget': 'last'}}}, - {'period_options': {0: ['save head last', - 'save budget last']}} +@pytest.mark.parametrize('input', [ + # flopy-style input + {'saverecord': {0: {'head': 'last', 'budget': 'last'}}}, + # MODFLOW 6-style input + {'period_options': {0: ['save head last', 'save budget last']}}, + # blank period to skip subsequent periods + {'period_options': {0: ['save head last', 'save budget last'], + 1: []}} ]) -def test_oc_setup(shellmound_model_with_dis, options): +def test_oc_setup(shellmound_model_with_dis, input): cfg = {'head_fileout_fmt': '{}.hds', 'budget_fileout_fmt': '{}.cbc'} - cfg.update(options) + cfg.update(input) m = shellmound_model_with_dis # deepcopy(model) m.cfg['oc'] = cfg oc = m.setup_oc() @@ -506,10 +510,15 @@ def test_oc_setup(shellmound_model_with_dis, options): options = read_mf6_block(ocfile, 'options') options = {k: ' '.join(v).lower() for k, v in options.items()} perioddata = read_mf6_block(ocfile, 'period') + # convert back to zero-based + perioddata = {k-1:v for k, v in perioddata.items()} assert 'fileout' in options['budget'] and '.cbc' in options['budget'] assert 'fileout' in options['head'] and '.hds' in options['head'] - assert 'save head last' in perioddata[1] - assert 'save budget last' in perioddata[1] + if 'saverecord' in input: + assert 'save head last' in perioddata[0] + assert 'save budget last' in perioddata[0] + else: + assert perioddata == input['period_options'] def test_rch_setup(shellmound_model_with_dis): @@ -728,7 +737,7 @@ def test_sfr_inflows_from_csv(model_with_sfr): pd.testing.assert_series_equal(left, right, check_names=False, check_freq=False) -#@pytest.mark.xfail(reason='flopy remove_package() issue') +@pytest.mark.xfail(reason='flopy remove_package() issue') def test_idomain_above_sfr(model_with_sfr): m = model_with_sfr sfr = m.sfr diff --git a/mfsetup/tests/test_mf6_shellmound_rot_grid.py b/mfsetup/tests/test_mf6_shellmound_rot_grid.py index f0b0f582..d0d9c79c 100644 --- a/mfsetup/tests/test_mf6_shellmound_rot_grid.py +++ b/mfsetup/tests/test_mf6_shellmound_rot_grid.py @@ -117,7 +117,7 @@ def test_rotated_grid(shellmound_cfg, shellmound_simulation, mf6_exe): cfg['setup_grid']['yoff'] = yoff cfg['setup_grid']['rotation'] = rotation cfg['dis']['dimensions']['nrow'] = nrow - cfg['dis']['dimensions']['ncol'] = 25 + cfg['dis']['dimensions']['ncol'] = ncol cfg = MF6model._parse_model_kwargs(cfg) kwargs = get_input_arguments(cfg['model'], mf6.ModflowGwf, diff --git a/mfsetup/tests/test_oc.py b/mfsetup/tests/test_oc.py new file mode 100644 index 00000000..4fdda516 --- /dev/null +++ b/mfsetup/tests/test_oc.py @@ -0,0 +1,64 @@ +"""Tests for the oc.py module +""" +import pytest + +from mfsetup.oc import fill_oc_stress_period_data, parse_oc_period_input + + +@pytest.mark.parametrize('input,expected,output_fmt', [ + # dictionary-based flopy-like input to (mf6) flopy-style input + ({'saverecord': {0: {'head': 'last', 'budget': 'last'}}}, + {'saverecord': {0: [('head', 'last'), ('budget', 'last')]}}, + 'mf6'), + # mf6-style input to (mf6) flopy-style input + ({'period_options': {0: ['save head last', 'save budget last']}}, + {'saverecord': {0: [('head', 'last'), ('budget', 'last')]}}, + 'mf6'), + # mf6-style input to flopy-style input + ({'period_options': {0: ['save head first', 'save budget first']}}, + {'stress_period_data': {(0, 0): ['save head', 'save budget']}}, + 'mfnwt'), + ({'period_options': {0: ['save head last', 'save budget last']}}, + {'stress_period_data': {(0, 9): ['save head', 'save budget']}}, + 'mfnwt'), + ({'period_options': {0: ['save head frequency 5', 'save budget frequency 5']}}, + {'stress_period_data': {(0, 0): ['save head', 'save budget'], + (0, 5): ['save head', 'save budget']}}, + 'mfnwt'), + ({'period_options': {0: ['save head steps 2 3', 'save budget steps 2 3']}}, + {'stress_period_data': {(0, 2): ['save head', 'save budget'], + (0, 3): ['save head', 'save budget']}}, + 'mfnwt'), + ({'period_options': {0: ['save head all', 'save budget all']}}, + None, 'mfnwt' + ), + # input already in flopy format + ({'stress_period_data': {(0, 2): ['save head', 'save budget'], + (0, 3): ['save head', 'save budget']}}, + {}, + 'mfnwt') + ], + ) +def test_parse_oc_period_input(input, expected, output_fmt): + results = parse_oc_period_input(input, nstp=[10], output_fmt=output_fmt) + # kludge for testing 'all' + if expected is None: + expected = {'stress_period_data': {(0, i): ['save head', 'save budget'] + for i in range(10)}} + assert results == expected + + +@pytest.mark.parametrize('stress_period_data,nper', + [({(0, 0): ['save head', 'save budget'], + (3, 0): ['save head'] + }, + 5)] + ) +def test_fill_oc_stress_period_data(stress_period_data, nper): + results = fill_oc_stress_period_data(stress_period_data, nper) + expected = {(0, 0): ['save head', 'save budget'], + (1, 0): ['save head', 'save budget'], + (2, 0): ['save head', 'save budget'], + (3, 0): ['save head'], + (4, 0): ['save head']} + assert results == expected diff --git a/mfsetup/tests/test_pfl_mfnwt_inset.py b/mfsetup/tests/test_pfl_mfnwt_inset.py index cbcac328..c2c3fe3b 100644 --- a/mfsetup/tests/test_pfl_mfnwt_inset.py +++ b/mfsetup/tests/test_pfl_mfnwt_inset.py @@ -515,7 +515,6 @@ def test_lak_setup(pfl_nwt_with_dis): def test_nwt_setup(pfl_nwt, project_root_path): - m = pfl_nwt #deepcopy(pfl_nwt) m.cfg['nwt']['use_existing_file'] = project_root_path + '/mfsetup/tests/data/RGN_rjh_3_23_18.NWT' nwt = m.setup_nwt() @@ -525,12 +524,23 @@ def test_nwt_setup(pfl_nwt, project_root_path): nwt.write_file() -def test_oc_setup(pfl_nwt): +@pytest.mark.parametrize('input,expected', [ + # MODFLOW 6-style input + ({'period_options': {0: ['save head last', 'save budget last'], + 1: []}}, + {'stress_period_data': {(0, 0): ['save head', 'save budget'], + (0, 1): []}}), + # MODFLOW 2005-style input + ({'stress_period_data': {(0, 0): ['save head', 'save budget'], + (0, 1): []}}, + {'stress_period_data': {(0, 0): ['save head', 'save budget'], + (0, 1): []}}) +]) +def test_oc_setup(pfl_nwt, input, expected): m = pfl_nwt + m.cfg['oc'].update(input) oc = m.setup_oc() - for (kper, kstp), words in oc.stress_period_data.items(): - assert kstp == m.perioddata.loc[kper, 'nstp'] - 1 - assert words == m.perioddata.loc[kper, 'oc'] + assert oc.stress_period_data == expected['stress_period_data'] # TODO: add datetime comments to OC file diff --git a/mfsetup/tests/test_pleasant_mfnwt_inset.py b/mfsetup/tests/test_pleasant_mfnwt_inset.py index d4f1a1d7..2bbdec52 100644 --- a/mfsetup/tests/test_pleasant_mfnwt_inset.py +++ b/mfsetup/tests/test_pleasant_mfnwt_inset.py @@ -103,6 +103,12 @@ def test_wel_setup(get_pleasant_nwt_with_dis_bas6): assert len(spd) >= nwells0 + n_added_wels +def test_oc_setup(get_pleasant_nwt_with_dis_bas6): + m = get_pleasant_nwt_with_dis_bas6 # deepcopy(model) + oc = m.setup_oc() + # oc stress period data should be filled + assert len(oc.stress_period_data) == m.nper + def test_model_setup(full_pleasant_nwt): m = full_pleasant_nwt assert isinstance(m, MFnwtModel)