refactor(discretization.fill_cells_vertically): only return bottom ar…

…ray, since top shouldn't be modified; refactor(sourcedata.setup_array): move resetting of model top to branch where the top modified by supplied lake bathymetry (since this should be the only situation where this is needed)

mfsetup/discretization.py

@@ -271,7 +271,7 @@ def fill_cells_vertically(top, botm):
     filled += np.nanmin(all_surfaces, axis=0)  # botm[-1]
     # append the model bottom elevations
     filled = np.append(filled, [np.nanmin(all_surfaces, axis=0)], axis=0)
-    return filled[0].copy(), filled[1:].copy()
+    return filled[1:].copy()
 
 
 def fix_model_layer_conflicts(top_array, botm_array,
@@ -748,6 +748,6 @@ def voxels_to_layers(voxel_array, z_edges, model_top=None, model_botm=None, no_d
 
     # finally, fill any remaining nans with next layer elevation (going upward)
     # might still have nans in areas where there are no voxel values, but model top and botm values
-    top, botm = fill_cells_vertically(layers[0], layers[1:])
-    layers = np.vstack([np.reshape(top, (1, *top.shape)), botm])
+    botm = fill_cells_vertically(layers[0], layers[1:])
+    layers = np.vstack([np.reshape(layers[0], (1, *layers[0].shape)), botm])
     return layers

mfsetup/sourcedata.py

@@ -1381,7 +1381,6 @@ def setup_array(model, package, var, data=None,
     # special handling of some variables
     # (for lakes)
     simulate_high_k_lakes = model.cfg['high_k_lakes']['simulate_high_k_lakes']
-    apply_lake_bathymetry_to_model_top = False
     if var == 'botm':
         if 'lak' in model.cfg['model']['packages'] and not model._load and\
             np.sum(model.lake_bathymetry) != 0:
@@ -1390,7 +1389,6 @@ def setup_array(model, package, var, data=None,
             # (bathymetry is loaded from the configuration file input
             #  when the model.lake_bathymetry property is called for the first time;
             #  see MFsetupMixin._set_lake_bathymetry)
-            apply_lake_bathymetry_to_model_top = True
             bathy = model.lake_bathymetry
 
             # save a copy of original top elevations
@@ -1427,6 +1425,25 @@ def setup_array(model, package, var, data=None,
             if model.version == 'mf6':
                 # reset the model top to the lake bottom
                 top[bathy != 0] -= bathy[bathy != 0]
+
+                # update the top in the model
+                # todo: refactor this code to be consolidated with _set_idomain and setup_dis
+                model._setup_array('dis', 'top',
+                        data={0: top},
+                        datatype='array2d', resample_method='linear',
+                        write_fmt='%.2f', dtype=float)
+                if hasattr(model, 'dis') and model.dis is not None:
+                    if model.version == 'mf6':
+                        model.dis.top = model.cfg['dis']['griddata']['top']
+                    else:
+                        model.dis.top = model.cfg['dis']['top'][0]
+
+                # write the top array again
+                top_filepath = model.setup_external_filepaths(package, 'top',
+                                                            model.cfg[package]['top_filename_fmt'])[0]
+                save_array(top_filepath, top,
+                        nodata=write_nodata,
+                        fmt=write_fmt)
         # if loading the model; use the model top that was just loaded in
         else:
             top_filename = model.cfg['dis']['griddata'].get('top')
@@ -1502,19 +1519,8 @@ def setup_array(model, package, var, data=None,
                 raise Exception('Model layers less than {} {} thickness'.format(min_thickness,
                                                                             model.length_units))
         # fill any nan values that are above or below active cells to avoid cell thickness errors
-        top, botm = fill_cells_vertically(top, botm)
-        # the top may have been modified by fill_cells_vertically
-        # update the top in the model
-        # todo: refactor this code to be consolidated with _set_idomain and setup_dis
-        model._setup_array('dis', 'top',
-                data={0: top},
-                datatype='array2d', resample_method='linear',
-                write_fmt='%.2f', dtype=float)
-        if hasattr(model, 'dis') and model.dis is not None:
-            if model.version == 'mf6':
-                model.dis.top = model.cfg['dis']['griddata']['top']
-            else:
-                model.dis.top = model.cfg['dis']['top'][0]
+        botm = fill_cells_vertically(top, botm)
+
         data = {i: arr for i, arr in enumerate(botm)}
 
         # special case of LGR models
@@ -1613,21 +1619,6 @@ def setup_array(model, package, var, data=None,
             dst = model.cfg['intermediate_data'][var][i]
             shutil.copy(src, dst)
 
-    # write the top array again, because top was filled
-    # with botm array above
-    # NOTE: that tends to corrupt the top array,
-    # is that only applicable to simulate_high_k_lakes ?)
-    if var == 'botm' and apply_lake_bathymetry_to_model_top:
-        top_filepath = model.setup_external_filepaths(package, 'top',
-                                                      model.cfg[package]['top_filename_fmt'])[0]
-        save_array(top_filepath, top,
-                   nodata=write_nodata,
-                   fmt=write_fmt)
-        if model.version == 'mf6':
-            src = filepaths[i]['filename']
-            dst = model.cfg['intermediate_data'][var][i]
-            shutil.copy(src, dst)
-
 
 def get_source_data_file_ext(cfg_data, package, var):
     if 'filenames' in cfg_data:

mfsetup/tests/test_discretization.py

@@ -118,12 +118,16 @@ def test_fill_na(all_layers):
     botm[-2] = 2
     botm[-2, 2, 2] = np.nan
 
-    top, botm = fill_cells_vertically(top, botm)
+    botm = fill_cells_vertically(top, botm)
     filled = all_layers.copy()
     filled[0] = top
     filled[1:] = botm
     assert filled[:, 2, 2].tolist() == [10., 10.,  8.,  8.,  8.,  5.,  5.,  5.,  5, 1.]
-    assert filled[:, 0, 0].tolist() == [8] * 8 + [2, 1]
+    # only compare layers 1:
+    # because fill_cells_vertically doesn't modify the model top
+    # (assuming that in most cases, there will be value
+    # for the model top elevation at each cell location)
+    assert filled[1:, 0, 0].tolist() == [8.] * 7 + [2., 1.]
 
 
 def test_make_ibound(all_layers):
@@ -134,7 +138,7 @@ def test_make_ibound(all_layers):
     botm[-2] = 2
     botm[-2, 2, 2] = np.nan
     botm[:, 3, 3] = np.arange(1, 10)[::-1] # column of cells all eq. to min thickness
-    filled_top, filled_botm = fill_cells_vertically(top, botm)
+    filled_botm = fill_cells_vertically(top, botm)
     ibound = make_ibound(top, botm, nodata=nodata,
                            minimum_layer_thickness=1,
                            drop_thin_cells=True,