From 9f53a9c10269a23995e21c1e2d8abac9d5598665 Mon Sep 17 00:00:00 2001 From: "Leaf, Andrew T" Date: Mon, 23 Sep 2024 10:30:45 -0500 Subject: [PATCH] refactor(discretization.fill_cells_vertically): only return bottom array, 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 | 6 ++-- mfsetup/sourcedata.py | 51 ++++++++++++---------------- mfsetup/tests/test_discretization.py | 10 ++++-- 3 files changed, 31 insertions(+), 36 deletions(-) diff --git a/mfsetup/discretization.py b/mfsetup/discretization.py index 7b95bede..666c5a55 100644 --- a/mfsetup/discretization.py +++ b/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 diff --git a/mfsetup/sourcedata.py b/mfsetup/sourcedata.py index 4d1fbd78..6727bde2 100644 --- a/mfsetup/sourcedata.py +++ b/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: diff --git a/mfsetup/tests/test_discretization.py b/mfsetup/tests/test_discretization.py index 73d55446..fb6fc359 100644 --- a/mfsetup/tests/test_discretization.py +++ b/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,