diff --git a/examples/01_basic_model.ipynb b/examples/01_basic_model.ipynb index 8ea46219..f67a9e49 100644 --- a/examples/01_basic_model.ipynb +++ b/examples/01_basic_model.ipynb @@ -276,8 +276,8 @@ "\n", "\n", "# add constant head cells at model boundaries\n", - "chd = nlmod.mfpackages.chd_at_model_edge_from_model_ds(\n", - " model_ds, gwf, head='starting_head')" + "model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain'])) \n", + "chd = nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head')" ] }, { @@ -694,7 +694,7 @@ "Attributes: (12/24)\n", " model_name: IJmuiden\n", " mfversion: mf6\n", - " model_dataset_created_on: 20211202_17:29:52\n", + " model_dataset_created_on: 20211203_12:31:15\n", " exe_name: c:\\users\\oebbe\\02_python\\nlmod\\nlmod\\mdims\\..\\...\n", " model_ws: model1\n", " figdir: model1\\figure\n", @@ -704,7 +704,7 @@ " anisotropy: 10\n", " fill_value_kh: 1.0\n", " fill_value_kv: 0.1\n", - " surface_drn_cond: 1000
  • model_name :
    IJmuiden
    mfversion :
    mf6
    model_dataset_created_on :
    20211203_12:31:15
    exe_name :
    c:\\users\\oebbe\\02_python\\nlmod\\nlmod\\mdims\\..\\..\\bin\\mf6.exe
    model_ws :
    model1
    figdir :
    model1\\figure
    cachedir :
    model1\\cache
    time_units :
    DAYS
    start_time :
    2015-01-01 00:00:00
    nper :
    1
    perlen :
    1.0
    nstp :
    1
    tsmult :
    1.0
    steady_start :
    0
    steady_state :
    1
    gridtype :
    structured
    extent :
    [ 95000. 105000. 494000. 500000.]
    delr :
    100.0
    delc :
    100.0
    nodata :
    -999
    anisotropy :
    10
    fill_value_kh :
    1.0
    fill_value_kv :
    0.1
    surface_drn_cond :
    1000
    • " ], "text/plain": [ "\n", @@ -1095,7 +1095,7 @@ "Attributes: (12/24)\n", " model_name: IJmuiden\n", " mfversion: mf6\n", - " model_dataset_created_on: 20211202_17:29:52\n", + " model_dataset_created_on: 20211203_12:31:15\n", " exe_name: c:\\users\\oebbe\\02_python\\nlmod\\nlmod\\mdims\\..\\...\n", " model_ws: model1\n", " figdir: model1\\figure\n", @@ -1136,7 +1136,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "INFO:nlmod.util:write script 2021120201_basic_model.ipynb to model workspace\n", + "INFO:nlmod.util:write script 2021120301_basic_model.ipynb to model workspace\n", "INFO:nlmod.util:write model dataset to cache\n", "INFO:nlmod.util:write modflow files to model workspace\n" ] @@ -1193,15 +1193,15 @@ "and distribution information.\n", "\n", " \n", - " Run start date and time (yyyy/mm/dd hh:mm:ss): 2021/12/02 17:30:04\n", + " Run start date and time (yyyy/mm/dd hh:mm:ss): 2021/12/03 12:31:27\n", " \n", " Writing simulation list file: mfsim.lst\n", " Using Simulation name file: mfsim.nam\n", " \n", " Solving: Stress period: 1 Time step: 1\n", " \n", - " Run end date and time (yyyy/mm/dd hh:mm:ss): 2021/12/02 17:30:05\n", - " Elapsed run time: 1.478 Seconds\n", + " Run end date and time (yyyy/mm/dd hh:mm:ss): 2021/12/03 12:31:29\n", + " Elapsed run time: 1.464 Seconds\n", " \n", " Normal termination of simulation.\n" ] diff --git a/examples/03_local_grid_refinement.ipynb b/examples/03_local_grid_refinement.ipynb index 600ad84b..975f77ec 100644 --- a/examples/03_local_grid_refinement.ipynb +++ b/examples/03_local_grid_refinement.ipynb @@ -265,8 +265,8 @@ "\n", "\n", "# add constant head cells at model boundaries\n", - "chd = nlmod.mfpackages.chd_at_model_edge_from_model_ds(\n", - " model_ds, gwf, head='starting_head')" + "model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain'])) \n", + "chd = nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head')" ] }, { diff --git a/nlmod/mfpackages/__init__.py b/nlmod/mfpackages/__init__.py index 35b02b74..ce589d0a 100644 --- a/nlmod/mfpackages/__init__.py +++ b/nlmod/mfpackages/__init__.py @@ -1,2 +1,2 @@ -from . import recharge, surface_water +from . import recharge, surface_water, constant_head from .mfpackages import * diff --git a/nlmod/mfpackages/constant_head.py b/nlmod/mfpackages/constant_head.py new file mode 100644 index 00000000..551f5abc --- /dev/null +++ b/nlmod/mfpackages/constant_head.py @@ -0,0 +1,60 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Dec 3 11:58:49 2021 + +@author: oebbe +""" + +import xarray as xr +import numpy as np + +from .. import util, cache + +@cache.cache_netcdf +def get_chd_at_model_edge(model_ds, idomain): + """get data array which is 1 at every active cell (defined by idomain) at + the boundaries of the model (xmin, xmax, ymin, ymax). Other cells are 0. + + Parameters + ---------- + model_ds : xarray.Dataset + dataset with model data. + idomain : xarray.DataArray + idomain used to get active cells and shape of DataArray + + Returns + ------- + model_ds_out : xarray.Dataset + dataset with chd array + """ + + # add constant head cells at model boundaries + + # get mask with grid edges + xmin = model_ds['x'] == model_ds['x'].min() + xmax = model_ds['x'] == model_ds['x'].max() + ymin = model_ds['y'] == model_ds['y'].min() + ymax = model_ds['y'] == model_ds['y'].max() + + model_ds_out = util.get_model_ds_empty(model_ds) + + if model_ds.gridtype == 'structured': + mask2d = (ymin | ymax | xmin | xmax) + + # assign 1 to cells that are on the edge and have an active idomain + model_ds_out['chd'] = xr.zeros_like(idomain) + for lay in model_ds.layer: + model_ds_out['chd'].loc[lay] = np.where( + mask2d & (idomain.loc[lay] == 1), 1, 0) + + elif model_ds.gridtype == 'unstructured': + mask = np.where([xmin | xmax | ymin | ymax])[1] + + # assign 1 to cells that are on the edge, have an active idomain + model_ds_out['chd'] = xr.zeros_like(idomain) + model_ds_out['chd'].loc[:, mask] = 1 + model_ds_out['chd'] = xr.where(idomain == 1, + model_ds_out['chd'], 0) + + + return model_ds_out \ No newline at end of file diff --git a/nlmod/mfpackages/mfpackages.py b/nlmod/mfpackages/mfpackages.py index f05dc1ab..c0a34d23 100644 --- a/nlmod/mfpackages/mfpackages.py +++ b/nlmod/mfpackages/mfpackages.py @@ -336,7 +336,8 @@ def sto_from_model_ds(model_ds, gwf, return sto -def chd_at_model_edge_from_model_ds(model_ds, gwf, head='starting_head'): +def chd_from_model_ds(model_ds, gwf, chd='chd', + head='starting_head'): """get constant head boundary at the model's edges from the model dataset. Parameters @@ -345,6 +346,9 @@ def chd_at_model_edge_from_model_ds(model_ds, gwf, head='starting_head'): dataset with model data. gwf : flopy ModflowGwf groundwaterflow object. + chd : str, optional + name of data variable in model_ds that is 1 for cells with a constant + head and zero for all other cells. The default is 'chd'. head : str, optional name of data variable in model_ds that is used as the head in the chd cells. The default is 'starting_head'. @@ -354,43 +358,18 @@ def chd_at_model_edge_from_model_ds(model_ds, gwf, head='starting_head'): chd : flopy ModflowGwfchd chd package """ - # add constant head cells at model boundaries - - # get mask with grid edges - xmin = model_ds['x'] == model_ds['x'].min() - xmax = model_ds['x'] == model_ds['x'].max() - ymin = model_ds['y'] == model_ds['y'].min() - ymax = model_ds['y'] == model_ds['y'].max() - + # get the stress_period_data if model_ds.gridtype == 'structured': - mask2d = (ymin | ymax | xmin | xmax) - - # assign 1 to cells that are on the edge and have an active idomain - model_ds['chd'] = xr.zeros_like(model_ds['idomain']) - for lay in model_ds.layer: - model_ds['chd'].loc[lay] = np.where( - mask2d & (model_ds['idomain'].loc[lay] == 1), 1, 0) - - # get the stress_period_data chd_rec = mdims.data_array_3d_to_rec_list(model_ds, - model_ds['chd'] != 0, + model_ds[chd] != 0, col1=head) elif model_ds.gridtype == 'unstructured': - mask = np.where([xmin | xmax | ymin | ymax])[1] - - # assign 1 to cells that are on the edge, have an active idomain - model_ds['chd'] = xr.zeros_like(model_ds['idomain']) - model_ds['chd'].loc[:, mask] = 1 - model_ds['chd'] = xr.where(model_ds['idomain'] == 1, - model_ds['chd'], 0) - - # get the stress_period_data - cellids = np.where(model_ds['chd']) + cellids = np.where(model_ds[chd]) chd_rec = list(zip(zip(cellids[0], cellids[1]), [1.0] * len(cellids[0]))) - chd = flopy.mf6.ModflowGwfchd(gwf, pname='chd', + chd = flopy.mf6.ModflowGwfchd(gwf, pname=chd, maxbound=len(chd_rec), stress_period_data=chd_rec, save_flows=True) diff --git a/tests/test_001_model.py b/tests/test_001_model.py index ed46452f..a789e3ef 100644 --- a/tests/test_001_model.py +++ b/tests/test_001_model.py @@ -191,8 +191,8 @@ def test_create_sea_model(tmpdir): nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) # add constant head cells at model boundaries - nlmod.mfpackages.chd_at_model_edge_from_model_ds( - model_ds, gwf, head='starting_head') + model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain'])) + nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head') # add knmi recharge to the model datasets model_ds.update(nlmod.read.knmi.add_knmi_to_model_dataset(model_ds)) @@ -262,8 +262,8 @@ def test_create_sea_model_perlen_list(tmpdir): nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) # add constant head cells at model boundaries - nlmod.mfpackages.chd_at_model_edge_from_model_ds( - model_ds, gwf, head='starting_head') + model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain'])) + nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head') # add knmi recharge to the model datasets model_ds.update(nlmod.read.knmi.add_knmi_to_model_dataset(model_ds)) @@ -325,8 +325,8 @@ def test_create_sea_model_perlen_14(tmpdir): nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) # add constant head cells at model boundaries - nlmod.mfpackages.chd_at_model_edge_from_model_ds( - model_ds, gwf, head='starting_head') + model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain'])) + nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head') # add knmi recharge to the model datasets model_ds.update(nlmod.read.knmi.add_knmi_to_model_dataset(model_ds)) diff --git a/tests/test_003_mfpackages.py b/tests/test_003_mfpackages.py index db333d5b..5da6ee0d 100644 --- a/tests/test_003_mfpackages.py +++ b/tests/test_003_mfpackages.py @@ -113,6 +113,8 @@ def chd_from_model_ds(tmpdir): nlmod.mfpackages.ic_from_model_ds(model_ds, gwf, starting_head=1.0) - chd = nlmod.mfpackages.chd_at_model_edge_from_model_ds(model_ds, gwf) + # add constant head cells at model boundaries + model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain'])) + chd = nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head') return chd