Skip to content

River well #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 16 commits into
base: main
Choose a base branch
from
Open

River well #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
5284cca
This commit is a mistake.
uniaofazaforca Jan 26, 2024
aa6f57a
Revert "This commit is a mistake."
uniaofazaforca Jan 26, 2024
b9ef788
Added example 1 data and updated model dimensions
uniaofazaforca Feb 9, 2024
acecdfd
Added functions to calculate the CHD boundaries
uniaofazaforca Apr 5, 2024
61c76e7
Added CHD 'hard-wired' boundaries to both sides of
uniaofazaforca Apr 5, 2024
6686896
Added boundary conditions: CHD for transport package,
uniaofazaforca Apr 5, 2024
33295bd
Added stress_period to wells of example 1.
uniaofazaforca Apr 5, 2024
e87f540
Substituted model name for exe_model.
uniaofazaforca Apr 5, 2024
7c79c0e
Added Example 1 and 2 to test_base.
uniaofazaforca Apr 5, 2024
235aca8
Control-script for example 1.
uniaofazaforca Apr 5, 2024
faa92ec
Example 1 - River and system of wells interacting.
uniaofazaforca Apr 5, 2024
3ee779d
Example 2 - Hydraulic barrier to contain contamination.
uniaofazaforca Apr 5, 2024
263ae38
Example 3 - Contamintion of a urban watershed with
uniaofazaforca Apr 5, 2024
c23a85c
Python file to test the control of the example 1.
uniaofazaforca Apr 5, 2024
886077f
Python script to test control of example 2.
uniaofazaforca Apr 5, 2024
663b313
Organized files for the publication.
uniaofazaforca Apr 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
236 changes: 215 additions & 21 deletions src/pymf6tools/base_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,33 @@
from copy import deepcopy


# Function to calculate the constant heads automatically
# def add_chd_to_dict(base_model_data):
# chd_value_left = base_model_data.get('chd_value_left')
# chd_value_right = base_model_data.get('chd_value_right')
# # Loop to reproduce the matrix for the CHD package
# chd = []
# for col, chd_value in zip([0, base_model_data['ncol']], [chd_value_left, chd_value_right]):
# for lay in range(base_model_data['nlay']):
# for row in range(0, base_model_data['nrow'] - 1):
# chd.append([(lay, row, col), chd_value])
# return chd

# Function to calculate the heads automatically for transport package
# def add_chd_to_dict_transp(base_transport_model_data, base_model_data):
# chd_value_left = base_model_data.get('chd_value_left')
# chd_value_right = base_model_data.get('chd_value_right')
# cnc_left = base_transport_model_data.get('concentration_source_wall')
# cnc_right = 0
# # Loop to reproduce the matrix for the CHD package
# chd = []
# for col, chd_value, cnc in zip([0, base_model_data['ncol']],
# [chd_value_left, chd_value_right],
# [cnc_left, cnc_right]):
# for lay in range(base_model_data['nlay']):
# for row in range(0, base_model_data['nrow'] - 1):
# chd.append([(lay, row, col), chd_value, cnc])

BASE_MODEL_DATA = {
# flopy.mf6.ModflowTdis
'times': (
Expand All @@ -29,11 +56,11 @@
'length_units': 'meters',
'repeat_times': 3, # nper = repeat_times + 1
# flopy.mf6.ModflowGwfdis
'nrow': 15,
'ncol': 10,
'nrow': 30,
'ncol': 20,
'nlay': 3,
'delr': 100.0,
'delc': 100.0,
'delr': 5.0,
'delc': 5.0,
'top': 15.0,
'botm': [-5.0, -10.0, -15.0],
# flopy.mf6.ModflowGwfnpf
Expand All @@ -43,41 +70,174 @@
'sy': 0.2,
'ss': 0.000001,
'initial_head': 10.0,
# flopy.mf6.ModflowGwfchd(
# flopy.mf6.ModflowGwfchd
# default values if function not running
'chd_value_left': 10,
'chd_value_right': 12,
'chd': [
[(0, 0, 0), 10.],
[(0, 14, 9), 10.]
],
[(0, 0, 0), 10.],
[(0, 1, 0), 10.],
[(0, 2, 0), 10.],
[(0, 3, 0), 10.],
[(0, 4, 0), 10.],
[(0, 5, 0), 10.],
[(0, 6, 0), 10.],
[(0, 7, 0), 10.],
[(0, 8, 0), 10.],
[(0, 9, 0), 10.],
[(0, 10, 0), 10.],
[(0, 11, 0), 10.],
[(0, 12, 0), 10.],
[(0, 13, 0), 10.],
[(0, 14, 0), 10.],
[(0, 15, 0), 10.],
[(0, 16, 0), 10.],
[(0, 17, 0), 10.],
[(0, 18, 0), 10.],
[(0, 19, 0), 10.],
[(0, 20, 0), 10.],
[(0, 21, 0), 10.],
[(0, 22, 0), 10.],
[(0, 23, 0), 10.],
[(0, 24, 0), 10.],
[(0, 25, 0), 10.],
[(0, 26, 0), 10.],
[(0, 27, 0), 10.],
[(0, 28, 0), 10.],
[(0, 29, 0), 10.],
[(0, 0, 19), 10.5],
[(0, 1, 19), 10.5],
[(0, 2, 19), 10.5],
[(0, 3, 19), 10.5],
[(0, 4, 19), 10.5],
[(0, 5, 19), 10.5],
[(0, 6, 19), 10.5],
[(0, 7, 19), 10.5],
[(0, 8, 19), 10.5],
[(0, 9, 19), 10.5],
[(0, 10, 19), 10.5],
[(0, 11, 19), 10.5],
[(0, 12, 19), 10.5],
[(0, 13, 19), 10.5],
[(0, 14, 19), 10.5],
[(0, 15, 19), 10.5],
[(0, 16, 19), 10.5],
[(0, 17, 19), 10.5],
[(0, 18, 19), 10.5],
[(0, 19, 19), 10.5],
[(0, 20, 19), 10.5],
[(0, 21, 19), 10.5],
[(0, 22, 19), 10.5],
[(0, 23, 19), 10.5],
[(0, 24, 19), 10.5],
[(0, 25, 19), 10.5],
[(0, 26, 19), 10.5],
[(0, 27, 19), 10.5],
[(0, 28, 19), 10.5],
[(0, 29, 19), 10.5],
],
#,
'transport': False,
'river': False,
}
}

#BASE_MODEL_DATA['chd'] = add_chd_to_dict(BASE_MODEL_DATA)

BASE_TRANSPORT_MODEL_DATA = {
'wells':{},
'initial_concentration': 1,
'concentration_source_wall':0,
'cnc': [
[(0, 5, 1), 10.],
[(0, 6, 1), 10.] # cell_id, conc (const)
[(0, 15, 4), 100.],
[(0, 16, 4), 100.],
[(1, 15, 4), 100.],
[(1, 16, 4), 100.] # cell_id, conc (const)
],
'scheme': 'UPSTREAM', #'TVD', # or 'UPSTREAM'
'longitudinal_dispersivity': 1.0,
# Ratio of transverse to longitudinal dispersitivity
'dispersivity_ratio': 1.0,
'porosity': 0.35,
'obs': None,
'nrow': 30,
'ncol': 20,
'nlay': 3,
# default values if function not running
'chd': [
[(0, 0, 0), 10.0, 10.0],
[(0, 14, 9), 10.0, 10.0]
],
[(0, 0, 0), 10., 1.],
[(0, 1, 0), 10., 1.],
[(0, 2, 0), 10., 1.],
[(0, 3, 0), 10., 1.],
[(0, 4, 0), 10., 1.],
[(0, 5, 0), 10., 1.],
[(0, 6, 0), 10., 1.],
[(0, 7, 0), 10., 1.],
[(0, 8, 0), 10., 1.],
[(0, 9, 0), 10., 1.],
[(0, 10, 0), 10., 1.],
[(0, 11, 0), 10., 1.],
[(0, 12, 0), 10., 1.],
[(0, 13, 0), 10., 1.],
[(0, 14, 0), 10., 1.],
[(0, 15, 0), 10., 1.],
[(0, 16, 0), 10., 1.],
[(0, 17, 0), 10., 1.],
[(0, 18, 0), 10., 1.],
[(0, 19, 0), 10., 1.],
[(0, 20, 0), 10., 1.],
[(0, 21, 0), 10., 1.],
[(0, 22, 0), 10., 1.],
[(0, 23, 0), 10., 1.],
[(0, 24, 0), 10., 1.],
[(0, 25, 0), 10., 1.],
[(0, 26, 0), 10., 1.],
[(0, 27, 0), 10., 1.],
[(0, 28, 0), 10., 1.],
[(0, 29, 0), 10., 1.],
[(0, 0, 19), 10.5, 1.],
[(0, 1, 19), 10.5, 1.],
[(0, 2, 19), 10.5, 1.],
[(0, 3, 19), 10.5, 1.],
[(0, 4, 19), 10.5, 1.],
[(0, 5, 19), 10.5, 1.],
[(0, 6, 19), 10.5, 1.],
[(0, 7, 19), 10.5, 1.],
[(0, 8, 19), 10.5, 1.],
[(0, 9, 19), 10.5, 1.],
[(0, 10, 19), 10.5, 1.],
[(0, 11, 19), 10.5, 1.],
[(0, 12, 19), 10.5, 1.],
[(0, 13, 19), 10.5, 1.],
[(0, 14, 19), 10.5, 1.],
[(0, 15, 19), 10.5, 1.],
[(0, 16, 19), 10.5, 1.],
[(0, 17, 19), 10.5, 1.],
[(0, 18, 19), 10.5, 1.],
[(0, 19, 19), 10.5, 1.],
[(0, 20, 19), 10.5, 1.],
[(0, 21, 19), 10.5, 1.],
[(0, 22, 19), 10.5, 1.],
[(0, 23, 19), 10.5, 1.],
[(0, 24, 19), 10.5, 1.],
[(0, 25, 19), 10.5, 1.],
[(0, 26, 19), 10.5, 1.],
[(0, 27, 19), 10.5, 1.],
[(0, 28, 19), 10.5, 1.],
[(0, 29, 19), 10.5, 1.],
],

}

NRIV = 12
#BASE_TRANSPORT_MODEL_DATA['chd'] = add_chd_to_dict_transp(BASE_TRANSPORT_MODEL_DATA, BASE_MODEL_DATA)

# Set the number of river cells in the model
NRIV = 27

BASE_RIVER_MODEL_DATA = {
'river_spd': {
'rivlay': [0] * NRIV,
'rivrow': [1, 1, 2, 3, 4, 5, 4, 5, 6, 7, 8, 7],
'rivcol': [0, 2, 1, 2, 3, 3, 4, 5, 6, 7, 7, 8],
'rivrow': [1, 1, 2, 3, 4, 5, 4, 5, 6, 7, 8, 7, 8, 8, 9, 10, 11, 11, 12, 13, 11, 14, 14, 14, 15, 16, 17],
'rivcol': [0, 2, 1, 2, 3, 3, 4, 5, 6, 7, 7, 8, 9, 9, 10, 11, 12, 11, 13, 14, 12, 15, 16, 17, 18, 19, 19],
'rivstg': np.linspace(13, 14, NRIV),
'rivbot': np.linspace(7, 10, NRIV),
'rivcnd': [0.05] * NRIV
Expand All @@ -91,9 +251,30 @@

BASE_WELL_MODEL_DATA = {
'wells': {
'wel_out': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 4, 4)},
'wel_out1': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 6, 4)},
'wel_out2': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 8, 4)},
'wel_out': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 10, 4)},
'wel_out1': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 12, 5)},
'wel_out2': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 14, 6)},
},

}

EXAMPLE_1_DATA = {
'wells': {
'wel_out': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 8, 4)},
'wel_out1': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 13, 8)},
'wel_out2': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 3, 9)},
'wel_out3': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 9, 9)},
'wel_out4': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 12, 7)},
'wel_out5': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 3, 4)},
},

}

EXAMPLE_3_DATA = {
'wells': {
'wel_out': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 10, 4)},
'wel_out1': {'q': (-0.05, -0.5, -0.05), 'coords': (2, 12, 5)},
'wel_out2': {'q': (-0.05, -0.5, -0.05), 'coords': (0, 14, 6)},
},

}
Expand All @@ -103,7 +284,11 @@ def make_model_data(
base_model_data=BASE_MODEL_DATA,
base_transport_model_data=BASE_TRANSPORT_MODEL_DATA,
base_river_model_data=BASE_RIVER_MODEL_DATA,
base_well_model_data=BASE_WELL_MODEL_DATA ):
base_well_model_data=BASE_WELL_MODEL_DATA,
example_1_data=EXAMPLE_1_DATA,
example_3_data=EXAMPLE_3_DATA
):

"""Make model data.

specific_model_data - dictionary with data specific for the current model
Expand All @@ -116,13 +301,22 @@ def make_model_data(
base_river_model_data=deepcopy(base_river_model_data)
base_transport_model_data=deepcopy(base_transport_model_data)
base_well_model_data=deepcopy(base_well_model_data)
example_1_data=deepcopy(example_1_data)
example_3_data=deepcopy(example_3_data)


if specific_model_data['transport']:
base_model_data.update(base_transport_model_data)
if specific_model_data['river_active']:
base_model_data.update(base_river_model_data)
if specific_model_data['wells_active']:
base_model_data.update(base_well_model_data)
if specific_model_data['example_1_data']:
base_model_data.update(example_1_data)
if specific_model_data['example_3_data']:
base_model_data.update(example_3_data)


# old way up to Python 3.8
if sys.version_info[:2] < (3, 9):
return {**base_model_data, **specific_model_data}
Expand Down
29 changes: 28 additions & 1 deletion src/pymf6tools/make_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,30 @@ def make_input(
)
file_extensions.append('wel')

if model_data['example_1_data']:
# Stress period data for the well
stress_period_data = {}
for index in range(len(times)):
entry = []
for well in model_data['wells'].values():
value = [well['coords'], well['q'][index]]
if model_data['transport']:
value.append(0)
entry.append(tuple(value))
stress_period_data[index + 1] = entry
wel_kwargs = {}
if model_data['transport']:
wel_kwargs.update({
'auxiliary': 'CONCENTRATION',
'pname': 'WEL-1'})
# Instantiating well package
flopy.mf6.ModflowGwfwel(
gwf,
stress_period_data=stress_period_data,
**wel_kwargs,
)
file_extensions.append('wel')

chd_kwargs = {}
if model_data['transport']:
chd_kwargs.update({
Expand Down Expand Up @@ -261,9 +285,12 @@ def make_transport_model(sim, model_data):
sourcerecarray = [
('CHD-1', 'AUX', 'CONCENTRATION'),
]

# for the wells active data
if model_data['wells_active']:
sourcerecarray.append(('WEL-1', 'AUX', 'CONCENTRATION'))
# for the example 1
if model_data['example_1_data']:
sourcerecarray.append(('WEL-1', 'AUX', 'CONCENTRATION'))

flopy.mf6.ModflowGwtssm(
gwt, sources=sourcerecarray, filename=f'{gwtname}.ssm'
Expand Down
Loading