Skip to content

Commit 5153cc6

Browse files
committed
feat: allow/warn boundary conditions on inactive cells
1 parent f10fbb1 commit 5153cc6

12 files changed

Lines changed: 229 additions & 20 deletions

File tree

Lines changed: 167 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,167 @@
1+
"""
2+
Test that stress packages issue a warning (not an error) when a boundary
3+
condition entry references a cell that is not in the active grid domain.
4+
The model should run to completion and the listing file should contain a
5+
warning about the inactive cell.
6+
7+
Packages tested: RCH, WEL, GHB, DRN, RIV, CHD
8+
"""
9+
10+
import os
11+
12+
import flopy
13+
import numpy as np
14+
import pytest
15+
from framework import TestFramework
16+
17+
cases = [
18+
"inact_rch",
19+
"inact_wel",
20+
"inact_ghb",
21+
"inact_drn",
22+
"inact_riv",
23+
"inact_chd",
24+
]
25+
26+
nlay, nrow, ncol = 1, 3, 3
27+
inactive_cell = (0, 0, 0)
28+
chd_cell = (0, 2, 2)
29+
top = 10.0
30+
botm = 0.0
31+
strt = 5.0
32+
idomain = np.ones((nlay, nrow, ncol), dtype=int)
33+
idomain[inactive_cell] = 0
34+
35+
36+
def _base_sim(test, name):
37+
ws = test.workspace
38+
sim = flopy.mf6.MFSimulation(
39+
sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws
40+
)
41+
flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(1.0, 1, 1.0)])
42+
flopy.mf6.ModflowIms(
43+
sim,
44+
outer_dvclose=1e-6,
45+
inner_dvclose=1e-6,
46+
)
47+
gwf = flopy.mf6.ModflowGwf(sim, modelname=name)
48+
flopy.mf6.ModflowGwfdis(
49+
gwf,
50+
nlay=nlay,
51+
nrow=nrow,
52+
ncol=ncol,
53+
top=top,
54+
botm=botm,
55+
idomain=idomain,
56+
)
57+
flopy.mf6.ModflowGwfnpf(gwf, k=1.0)
58+
flopy.mf6.ModflowGwfic(gwf, strt=strt)
59+
flopy.mf6.ModflowGwfchd(
60+
gwf,
61+
stress_period_data={0: [[chd_cell, strt]]},
62+
)
63+
64+
return sim, gwf
65+
66+
67+
def build_models(idx, test):
68+
name = cases[idx]
69+
sim, gwf = _base_sim(test, name)
70+
71+
if idx == 0:
72+
# RCH: one entry on the inactive cell, one on an active cell
73+
flopy.mf6.ModflowGwfrch(
74+
gwf,
75+
stress_period_data={
76+
0: [
77+
[inactive_cell, 1e-3],
78+
[(0, 1, 1), 1e-3],
79+
]
80+
},
81+
)
82+
83+
elif idx == 1:
84+
# WEL: one entry on the inactive cell
85+
flopy.mf6.ModflowGwfwel(
86+
gwf,
87+
stress_period_data={
88+
0: [
89+
[inactive_cell, -1.0],
90+
[(0, 1, 1), -1.0],
91+
]
92+
},
93+
)
94+
95+
elif idx == 2:
96+
# GHB: one entry on the inactive cell
97+
flopy.mf6.ModflowGwfghb(
98+
gwf,
99+
stress_period_data={
100+
0: [
101+
[inactive_cell, strt, 1.0],
102+
[(0, 1, 1), strt, 1.0],
103+
]
104+
},
105+
)
106+
107+
elif idx == 3:
108+
# DRN: one entry on the inactive cell
109+
flopy.mf6.ModflowGwfdrn(
110+
gwf,
111+
stress_period_data={
112+
0: [
113+
[inactive_cell, 3.0, 1.0],
114+
[(0, 1, 1), 3.0, 1.0],
115+
]
116+
},
117+
)
118+
119+
elif idx == 4:
120+
# RIV: one entry on the inactive cell
121+
flopy.mf6.ModflowGwfriv(
122+
gwf,
123+
stress_period_data={
124+
0: [
125+
[inactive_cell, strt, 1.0, botm],
126+
[(0, 1, 1), strt, 1.0, botm],
127+
]
128+
},
129+
)
130+
131+
elif idx == 5:
132+
# CHD: one entry on the inactive cell (the base sim already has a
133+
# CHD on an active cell; add a separate package with the bad entry)
134+
flopy.mf6.ModflowGwfchd(
135+
gwf,
136+
stress_period_data={
137+
0: [
138+
[inactive_cell, strt],
139+
[(0, 0, 1), strt],
140+
]
141+
},
142+
pname="chd_bad",
143+
filename=f"{name}_bad.chd",
144+
)
145+
146+
return sim, None
147+
148+
149+
def check_output(idx, test):
150+
mfsim_lst = os.path.join(test.workspace, "mfsim.lst")
151+
assert os.path.isfile(mfsim_lst)
152+
lst_text = open(mfsim_lst).read()
153+
assert "Normal termination" in lst_text
154+
warning_phrase = "outside active grid domain"
155+
assert warning_phrase.lower() in lst_text.lower()
156+
157+
158+
@pytest.mark.parametrize("idx, name", enumerate(cases))
159+
def test_mf6model(idx, name, function_tmpdir, targets):
160+
test = TestFramework(
161+
name=name,
162+
workspace=function_tmpdir,
163+
build=lambda t: build_models(idx, t),
164+
check=lambda t: check_output(idx, t),
165+
targets=targets,
166+
)
167+
test.run()

src/Model/GroundWaterEnergy/gwe-ctp.f90

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,7 @@ subroutine ctp_rp(this)
142142
! -- Reset previous CTPs to active cell
143143
do i = 1, this%nbound
144144
node = this%nodelist(i)
145+
if (node <= 0) cycle
145146
this%ibound(node) = this%ibcnum
146147
end do
147148
!
@@ -152,6 +153,7 @@ subroutine ctp_rp(this)
152153
ierr = 0
153154
do i = 1, this%nbound
154155
node = this%nodelist(i)
156+
if (node <= 0) cycle
155157
ibd = this%ibound(node)
156158
if (ibd < 0) then
157159
call this%dis%noder_to_string(node, nodestr)
@@ -188,6 +190,7 @@ subroutine ctp_ad(this)
188190
! -- Process each entry in the constant temperature cell list
189191
do i = 1, this%nbound
190192
node = this%nodelist(i)
193+
if (node <= 0) cycle
191194
cb = this%temp_mult(i)
192195
!
193196
this%xnew(node) = cb
@@ -217,6 +220,7 @@ subroutine ctp_ck(this)
217220
! -- check stress period data
218221
do i = 1, this%nbound
219222
node = this%nodelist(i)
223+
if (node <= 0) cycle
220224
! -- accumulate errors
221225
if (this%temp_mult(i) < DZERO) then
222226
call this%dis%noder_to_string(node, nodestr)
@@ -271,6 +275,7 @@ subroutine ctp_cq(this, x, flowja, iadv)
271275
! -- Loop through each boundary calculating flow.
272276
do i = 1, this%nbound
273277
node = this%nodelist(i)
278+
if (node <= 0) cycle
274279
idiag = this%dis%con%ia(node)
275280
rate = DZERO
276281
ratein = DZERO

src/Model/GroundWaterEnergy/gwe-esl.f90

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -190,6 +190,10 @@ subroutine esl_cf(this)
190190
do i = 1, this%nbound
191191
node = this%nodelist(i)
192192
this%hcof(i) = DZERO
193+
if (node <= 0) then
194+
this%rhs(i) = DZERO
195+
cycle
196+
end if
193197
if (this%ibound(node) <= 0) then
194198
this%rhs(i) = DZERO
195199
cycle
@@ -224,6 +228,7 @@ subroutine esl_fc(this, rhs, ia, idxglo, matrix_sln)
224228
! -- Copy package rhs and hcof into solution rhs and amat
225229
do i = 1, this%nbound
226230
n = this%nodelist(i)
231+
if (n <= 0) cycle
227232
rhs(n) = rhs(n) + this%rhs(i)
228233
ipos = ia(n)
229234
call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i))

src/Model/GroundWaterFlow/gwf-chd.f90

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -136,6 +136,7 @@ subroutine chd_rp(this)
136136
! -- Reset previous CHDs to active cell
137137
do i = 1, this%nbound
138138
node = this%nodelist(i)
139+
if (node <= 0) cycle
139140
this%ibound(node) = this%ibcnum
140141
end do
141142
!
@@ -146,6 +147,7 @@ subroutine chd_rp(this)
146147
ierr = 0
147148
do i = 1, this%nbound
148149
node = this%nodelist(i)
150+
if (node <= 0) cycle
149151
ibd = this%ibound(node)
150152
if (ibd < 0) then
151153
call this%dis%noder_to_string(node, nodestr)
@@ -180,6 +182,7 @@ subroutine chd_ad(this)
180182
! -- Process each entry in the specified-head cell list
181183
do i = 1, this%nbound
182184
node = this%nodelist(i)
185+
if (node <= 0) cycle
183186
hb = this%head_mult(i)
184187
!
185188
this%xnew(node) = hb
@@ -211,6 +214,7 @@ subroutine chd_ck(this)
211214
! -- check stress period data
212215
do i = 1, this%nbound
213216
node = this%nodelist(i)
217+
if (node <= 0) cycle
214218
bt = this%dis%bot(node)
215219
! -- accumulate errors
216220
if (this%head_mult(i) < bt .and. this%icelltype(node) /= 0) then
@@ -277,6 +281,7 @@ subroutine calc_chd_rate(this)
277281
! -- Loop through each boundary calculating flow.
278282
do i = 1, this%nbound
279283
node = this%nodelist(i)
284+
if (node <= 0) cycle
280285
idiag = this%dis%con%ia(node)
281286
rate = DZERO
282287
ratein = DZERO

src/Model/GroundWaterFlow/gwf-drn.f90

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -349,6 +349,11 @@ subroutine drn_cf(this)
349349
! -- Calculate hcof and rhs for each drn entry
350350
do i = 1, this%nbound
351351
node = this%nodelist(i)
352+
if (node <= 0) then
353+
this%hcof(i) = DZERO
354+
this%rhs(i) = DZERO
355+
cycle
356+
end if
352357
if (this%ibound(node) <= 0) then
353358
this%hcof(i) = DZERO
354359
this%rhs(i) = DZERO
@@ -394,6 +399,7 @@ subroutine drn_fc(this, rhs, ia, idxglo, matrix_sln)
394399
! -- Copy package rhs and hcof into solution rhs and amat
395400
do i = 1, this%nbound
396401
n = this%nodelist(i)
402+
if (n <= 0) cycle
397403
rhs(n) = rhs(n) + this%rhs(i)
398404
ipos = ia(n)
399405
call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i))
@@ -436,6 +442,7 @@ subroutine drn_fn(this, rhs, ia, idxglo, matrix_sln)
436442
if (this%iauxddrncol /= 0) then
437443
do i = 1, this%nbound
438444
node = this%nodelist(i)
445+
if (node <= 0) cycle
439446
!
440447
! -- test if node is constant or inactive
441448
if (this%ibound(node) <= 0) then

src/Model/GroundWaterFlow/gwf-riv.f90

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -276,6 +276,11 @@ subroutine riv_cf(this)
276276
! -- Calculate hcof and rhs for each riv entry
277277
do i = 1, this%nbound
278278
node = this%nodelist(i)
279+
if (node <= 0) then
280+
this%hcof(i) = DZERO
281+
this%rhs(i) = DZERO
282+
cycle
283+
end if
279284
if (this%ibound(node) <= 0) then
280285
this%hcof(i) = DZERO
281286
this%rhs(i) = DZERO
@@ -315,6 +320,7 @@ subroutine riv_fc(this, rhs, ia, idxglo, matrix_sln)
315320
! -- Copy package rhs and hcof into solution rhs and amat
316321
do i = 1, this%nbound
317322
n = this%nodelist(i)
323+
if (n <= 0) cycle
318324
rhs(n) = rhs(n) + this%rhs(i)
319325
ipos = ia(n)
320326
call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i))

src/Model/GroundWaterFlow/gwf-wel.f90

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -444,6 +444,10 @@ subroutine wel_cf(this)
444444
do i = 1, this%nbound
445445
node = this%nodelist(i)
446446
this%hcof(i) = DZERO
447+
if (node <= 0) then
448+
this%rhs(i) = DZERO
449+
cycle
450+
end if
447451
if (this%ibound(node) <= 0) then
448452
this%rhs(i) = DZERO
449453
cycle
@@ -497,6 +501,7 @@ subroutine wel_fc(this, rhs, ia, idxglo, matrix_sln)
497501
! -- Copy package rhs and hcof into solution rhs and amat
498502
do i = 1, this%nbound
499503
n = this%nodelist(i)
504+
if (n <= 0) cycle
500505
rhs(n) = rhs(n) + this%rhs(i)
501506
ipos = ia(n)
502507
call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i))
@@ -536,6 +541,7 @@ subroutine wel_fn(this, rhs, ia, idxglo, matrix_sln)
536541
! -- Copy package rhs and hcof into solution rhs and amat
537542
do i = 1, this%nbound
538543
node = this%nodelist(i)
544+
if (node <= 0) cycle
539545
!
540546
! -- test if node is constant or inactive
541547
if (this%ibound(node) <= 0) then
@@ -605,6 +611,7 @@ subroutine wel_afr_csv_write(this)
605611
! -- format
606612
do i = 1, this%nbound
607613
nodereduced = this%nodelist(i)
614+
if (nodereduced <= 0) cycle
608615
!
609616
! -- test if node is constant or inactive
610617
if (this%ibound(nodereduced) <= 0) then

src/Model/GroundWaterTransport/gwt-cnc.f90

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,7 @@ subroutine cnc_rp(this)
142142
! -- Reset previous CNCs to active cell
143143
do i = 1, this%nbound
144144
node = this%nodelist(i)
145+
if (node <= 0) cycle
145146
this%ibound(node) = this%ibcnum
146147
end do
147148
!
@@ -152,6 +153,7 @@ subroutine cnc_rp(this)
152153
ierr = 0
153154
do i = 1, this%nbound
154155
node = this%nodelist(i)
156+
if (node <= 0) cycle
155157
ibd = this%ibound(node)
156158
if (ibd < 0) then
157159
call this%dis%noder_to_string(node, nodestr)
@@ -190,6 +192,7 @@ subroutine cnc_ad(this)
190192
! -- Process each entry in the constant concentration cell list
191193
do i = 1, this%nbound
192194
node = this%nodelist(i)
195+
if (node <= 0) cycle
193196
cb = this%conc_mult(i)
194197
!
195198
this%xnew(node) = cb
@@ -220,6 +223,7 @@ subroutine cnc_ck(this)
220223
! -- check stress period data
221224
do i = 1, this%nbound
222225
node = this%nodelist(i)
226+
if (node <= 0) cycle
223227
! -- accumulate errors
224228
if (this%conc_mult(i) < DZERO) then
225229
call this%dis%noder_to_string(node, nodestr)
@@ -277,6 +281,7 @@ subroutine cnc_cq(this, x, flowja, iadv)
277281
! -- Loop through each boundary calculating flow.
278282
do i = 1, this%nbound
279283
node = this%nodelist(i)
284+
if (node <= 0) cycle
280285
idiag = this%dis%con%ia(node)
281286
rate = DZERO
282287
ratein = DZERO

0 commit comments

Comments
 (0)