fix(ic.py): handle cases where there is no strt source data configura…

…tion or default parent model source data. Add error trap for invalid strt source data and additional testing of strt input options.

mfsetup/ic.py

@@ -19,6 +19,13 @@ def setup_strt(model, package, strt=None, source_data_config=None,
                write_nodata=None,
                **kwargs):
 
+    # default arguments to ArraySourceData
+    default_kwargs = {
+        'resample_method': 'linear'
+    }
+    for k, v in default_kwargs.items():
+        if k not in kwargs:
+            kwargs[k] = v
     var = 'strt'
     datatype = 'array3d'
     # model that strt values could come from
@@ -53,8 +60,6 @@ def setup_strt(model, package, strt=None, source_data_config=None,
     # data read from binary file with parent model head solution
     elif binary_file:
         kwargs = get_input_arguments(kwargs, MFBinaryArraySourceData)
-        if 'resample_method' not in kwargs:
-            kwargs['resample_method'] = 'linear'
         sd = MFBinaryArraySourceData(variable='strt', filename=binary_file,
                                      datatype=datatype,
                                      dest_model=model,
@@ -79,14 +84,19 @@ def setup_strt(model, package, strt=None, source_data_config=None,
                              time_units=model.time_units,
                              **kwargs)
     # data from files
-    elif source_data_config[var] is not None:
-        #ext = get_source_data_file_ext(source_data_config, package, var)
-        kwargs = get_input_arguments(kwargs, ArraySourceData)
-        sd = ArraySourceData.from_config(source_data_config,
-                                         datatype=datatype,
-                                         variable=var,
-                                         dest_model=model,
-                                         **kwargs)
+    elif source_data_config:
+        if source_data_config.get(var) is not None:
+            #ext = get_source_data_file_ext(source_data_config, package, var)
+            source_data_files = source_data_config.get(var)
+            kwargs = get_input_arguments(kwargs, ArraySourceData)
+            sd = ArraySourceData.from_config(source_data_files,
+                                            datatype=datatype,
+                                            variable=var,
+                                            dest_model=model,
+                                            **kwargs)
+        else:
+            raise ValueError(f"Invalid configuration input: {package}: source_data:\n"
+                             f"Need a {var}: sub-block")
 
     # default to setting strt from model top
     else:

mfsetup/mf6_defaults.yml

@@ -77,9 +77,6 @@ tdis:
 ic:
   griddata:
     strt:
-  source_data:
-    strt:
-    resample_method: 'linear'
   strt_filename_fmt: "strt_{:03d}.dat"
   write_fmt: '%.2f'
 

mfsetup/tests/data/pfl_nwt_test.yml

@@ -64,11 +64,11 @@ dis:
         4: 3
 
 bas6:
-  source_data:
-    strt:
-      from_parent:
-        binaryfile: 'plainfieldlakes/pfl.hds'
-        stress_period: 0
+    source_data:
+      strt:
+        from_parent:
+          binaryfile: 'plainfieldlakes/pfl.hds'
+          stress_period: 0
 
 upw:
   ipakcb: 53 # unit for writing cell budget output

mfsetup/tests/test_pfl_mfnwt_inset.py

@@ -5,6 +5,7 @@
 import glob
 import os
 import shutil
+from collections import defaultdict
 from copy import deepcopy
 from pathlib import Path
 
@@ -18,7 +19,7 @@
 mf6 = flopy.mf6
 from mfsetup import MFnwtModel
 from mfsetup.checks import check_external_files_for_nans
-from mfsetup.fileio import exe_exists, load_cfg, remove_file_header
+from mfsetup.fileio import exe_exists, load, load_cfg, remove_file_header
 from mfsetup.grid import MFsetupGrid, get_ij
 from mfsetup.units import convert_length_units, convert_time_units
 from mfsetup.utils import get_input_arguments
@@ -226,32 +227,115 @@ def test_dis_setup(pfl_nwt_with_grid):
     assert Path(m.cfg['intermediate_data']['output_folder']).is_dir()
 
 
-def test_bas_setup(pfl_nwt_with_dis):
+@pytest.mark.parametrize('bas6_config,default_parent_source_data', (
+    ('config file',True),  # test whatever is in the configuration file
+    # with default_parent_source_data
+    # starting heads are automatically resampled from parent model
+    # unless a strt array or binary head file is provided
+    #
+    # test case of default configuration (no bas: block on configuration file)
+    # starting heads are set to the model top
+    ('defaults', False),
+    # test case where just bas: (None configuration) is argued in configration file
+    # (this overrides the defaults)
+    # starting heads are set to the model top
+    (None, False),
+    # starting heads from a raster
+    ({'source_data': {
+        'strt': {
+            'filename': 'plainfieldlakes/source_data/dem10m.tif'
+        }
+        }}, False),
+    # need a strt variable
+    pytest.param({'source_data': {
+        'filename': 'plainfieldlakes/source_data/dem10m.tif'
+        }}, False, marks=pytest.mark.xfail(reason='some bug')),
+    # with a layer specified (gets repeated)
+    ({'source_data': {
+        'strt': {
+            'filenames': {
+                0: 'plainfieldlakes/source_data/dem10m.tif'
+            }
+        }
+    }}, False),
+    # starting heads from MODFLOW binary head output
+    ({'source_data': {
+        'strt': {
+            'from_parent': {
+                'binaryfile': 'plainfieldlakes/pfl.hds',
+                'stress_period': 0
+            }
+        }
+        }}, False)
+)
+                         )
+def test_bas_setup(pfl_nwt_cfg, pfl_nwt_with_dis, bas6_config,
+                   default_parent_source_data, project_root_path):
+    """Test setup of the BAS6 package, especially the starting heads."""
+    cfg = pfl_nwt_cfg.copy()
+    project_root_path = Path(project_root_path)
+    # load defaults here
+    default_cfg = project_root_path / 'mfsetup/mfnwt_defaults.yml'
+    defaults = load(default_cfg)
+    if bas6_config == 'config file':
+        pass
+    elif bas6_config == 'defaults':
+        cfg['bas6'] = defaults['bas6']
+    elif bas6_config is None:
+        cfg['bas6'] = defaultdict(dict)
+    else:
+        cfg['bas6'] = bas6_config
+    cfg['parent']['default_source_data'] = default_parent_source_data
 
-    m = pfl_nwt_with_dis  #deepcopy(pfl_nwt_with_dis)
+    m = MFnwtModel(cfg=cfg, **cfg['model'])
+    m.setup_grid()
+    m.setup_dis()
 
-    # test intermediate array creation
     bas = m.setup_bas6()
-    arrayfiles = m.cfg['intermediate_data']['strt'] + \
-                 m.cfg['intermediate_data']['ibound']
-                 #m.cfg['intermediate_data']['lakarr']
-    for f in arrayfiles:
-        assert os.path.exists(f)
 
-    # test using previously made external files as input
-    if m.version == 'mf6':
-        assert m.cfg['bas6']['strt'] == m.cfg['external_files']['strt']
-        assert m.cfg['bas6']['ibound'] == m.cfg['external_files']['ibound']
+    assert bas.strt.array.shape == m.modelgrid.shape
+    assert not np.isnan(bas.strt.array).any(axis=(0, 1, 2))
+
+    # In the absence of a parent model with default_source_data
+    # or input to strt
+    # check that strt was set from model top
+    if bas6_config in ('defaults', None):
+        assert np.allclose(bas.strt.array,
+                           np.array([m.dis.top.array] * m.modelgrid.nlay))
+    # assumes that all rasters being tested
+    # are the same as the dem used to make the model top
+    elif isinstance(bas6_config, dict):
+        if 'filename' in bas6_config.get('source_data', dict()) or\
+            'filenames' in bas6_config.get('source_data', dict()):
+            assert np.allclose(bas.strt.array,
+                            np.array([m.dis.top.array] * m.modelgrid.nlay))
+    # TODO: placeholder for more rigorous test that starting heads
+    # are consistent with parent model head solution
     else:
-        assert m.cfg['bas6']['strt'] == m.cfg['intermediate_data']['strt']
-        assert m.cfg['bas6']['ibound'] == m.cfg['intermediate_data']['ibound']
-    bas = m.setup_bas6()
-    bas.write_file()
-    arrayfiles = m.cfg['bas6']['strt'] + \
-                 m.cfg['bas6']['ibound']
-    for f in arrayfiles:
-        assert os.path.exists(f)
-    assert os.path.exists(bas.fn_path)
+        pass
+
+    if bas6_config == 'defaults':
+        # test intermediate array creation
+        arrayfiles = m.cfg['intermediate_data']['strt'] + \
+                    m.cfg['intermediate_data']['ibound']
+                    #m.cfg['intermediate_data']['lakarr']
+        for f in arrayfiles:
+            assert os.path.exists(f)
+
+        # test using previously made external files as input
+        if m.version == 'mf6':
+            assert m.cfg['bas6']['strt'] == m.cfg['external_files']['strt']
+            assert m.cfg['bas6']['ibound'] == m.cfg['external_files']['ibound']
+        else:
+            assert m.cfg['bas6']['strt'] == m.cfg['intermediate_data']['strt']
+            assert m.cfg['bas6']['ibound'] == m.cfg['intermediate_data']['ibound']
+        bas = m.setup_bas6()
+        bas.write_file()
+        arrayfiles = m.cfg['bas6']['strt'] + \
+                    m.cfg['bas6']['ibound']
+        for f in arrayfiles:
+            assert os.path.exists(f)
+        assert os.path.exists(bas.fn_path)
 
 
 @pytest.mark.parametrize('simulate_high_k_lakes', (False, True))