Skip to content

Commit 849578f

Browse files
wpbonelliwpbonelli
authored
fix(exg-gwf*): check for equal idomain (#2149)
Check that models connected to a flow model by an exchange have an idomain identical to the flow model's idomain. This is an easy error to make, and (up to now) hard to detect. I only added a test for PRT for now. The logic is identical in GWT and GWE. #2144 facilitates adding the same check for flow model interface. That will be done in a separate PR since it deserves a separate release note.
1 parent 216b363 commit 849578f

File tree

5 files changed

+140
-4
lines changed

5 files changed

+140
-4
lines changed

‎autotest/test_prt_exg.py‎

Lines changed: 31 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,15 @@
88
99
Results are compared against a MODPATH 7 model.
1010
11-
This test includes two cases, one which gives
12-
boundnames to particles and one which does not.
11+
This test includes four cases: one which gives
12+
boundnames to particles, one which does not, a
13+
third in which the flow model has a uniformly
14+
active idomain array while the tracking model
15+
does not, and a final case in which flow and
16+
tracking model have different IDOMAIN arrays,
17+
both non-uniform, where the active region is
18+
the same size but consists of different cells.
19+
Both latter cases should be caught as errors.
1320
"""
1421

1522
from pathlib import Path
@@ -30,7 +37,7 @@
3037
)
3138

3239
simname = "prtexg01"
33-
cases = [simname, f"{simname}bnms"]
40+
cases = [simname, f"{simname}bnms", f"{simname}idmu", f"{simname}idmn"]
3441

3542

3643
def get_model_name(idx, mdl):
@@ -41,18 +48,34 @@ def build_mf6_sim(idx, test):
4148
# create simulation
4249
name = cases[idx]
4350
sim = FlopyReadmeCase.get_gwf_sim(name, test.workspace, test.targets["mf6"])
51+
gwf = sim.get_model()
4452

4553
# create prt model
4654
prt_name = get_model_name(idx, "prt")
4755
prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name)
4856

4957
# create prt discretization
58+
idomain = np.ones(
59+
(FlopyReadmeCase.nlay, FlopyReadmeCase.nrow, FlopyReadmeCase.ncol)
60+
)
61+
if "idm" in name:
62+
# add an inactive cell to
63+
# tracking model idomain
64+
idomain[-1, -1, -1] = 0
65+
if "idmn" in name:
66+
# add a (different) inactive
67+
# cell to flow model idomain
68+
gwf_idomain = idomain.copy()
69+
gwf_idomain[-1, -1, -1] = 1
70+
gwf_idomain[0, 0, 0] = 0
71+
gwf.dis.idomain = gwf_idomain
5072
flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
5173
prt,
5274
pname="dis",
5375
nlay=FlopyReadmeCase.nlay,
5476
nrow=FlopyReadmeCase.nrow,
5577
ncol=FlopyReadmeCase.ncol,
78+
idomain=idomain,
5679
)
5780

5881
# create mip package
@@ -155,7 +178,7 @@ def build_models(idx, test):
155178
gwf_name = get_model_name(idx, "gwf")
156179
gwf = mf6sim.get_model(gwf_name)
157180
mp7sim = build_mp7_sim(idx, test.workspace / "mp7", test.targets["mp7"], gwf)
158-
return mf6sim, mp7sim
181+
return mf6sim, None if "idm" in test.name else mp7sim
159182

160183

161184
def check_output(idx, test):
@@ -178,6 +201,9 @@ def check_output(idx, test):
178201
# extract model grid
179202
mg = gwf.modelgrid
180203

204+
if "idm" in name:
205+
return
206+
181207
# check mf6 output files exist
182208
gwf_budget_file = f"{gwf_name}.bud"
183209
gwf_head_file = f"{gwf_name}.hds"
@@ -362,5 +388,6 @@ def test_mf6model(idx, name, function_tmpdir, targets, plot):
362388
plot=lambda t: plot_output(idx, t) if plot else None,
363389
targets=targets,
364390
compare=None,
391+
xfail="idm" in name,
365392
)
366393
test.run()

‎doc/ReleaseNotes/develop.tex‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
\textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\
2121
\underline{BASIC FUNCTIONALITY}
2222
\begin{itemize}
23+
\item Previously, GWF Model exchanges with other model types (GWT, GWE, PRT) verified that the flow model and the coupled model had the same number of active nodes, but did not check that the models' IDOMAIN arrays selected precisely the same set. Exchanges will now check that model IDOMAIN arrays are identical.
2324
\item A regression was recently introduced into the PRT model's generalized tracking method, in which a coordinate transformation was carried out prematurely. This could result in incorrect particle positions in and near quad-refined cells. This bug has been fixed.
2425
\item The PRT model previously allowed particles to be released at any time. Release times falling outside the bounds of the simulation's time discretization could produce undefined behavior. Any release times occurring before the simulation begins (i.e. negative times) will now be skipped with a warning message. If EXTEND\_TRACKING is not enabled, release times occurring after the end of the simulation will now be skipped with a warning message as well. If EXTEND\_TRACKING is enabled, release times after the end of the simulation are allowed.
2526
\end{itemize}

‎src/Exchange/exg-gwfgwe.f90‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -181,6 +181,9 @@ end subroutine exg_df
181181
subroutine exg_ar(this)
182182
! -- modules
183183
use MemoryManagerModule, only: mem_checkin
184+
use DisModule, only: DisType
185+
use DisvModule, only: DisvType
186+
use DisuModule, only: DisuType
184187
! -- dummy
185188
class(GwfGweExchangeType) :: this
186189
! -- local
@@ -194,6 +197,11 @@ subroutine exg_ar(this)
194197
& GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
195198
& GWE Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
196199
& Ensure discretization packages, including IDOMAIN, are identical.')"
200+
character(len=*), parameter :: fmtidomerr = &
201+
"('GWF and GWE Models do not have the same discretization for exchange&
202+
& ',a,'.&
203+
& GWF Model and GWE Model have different IDOMAIN arrays.&
204+
& Ensure discretization packages, including IDOMAIN, are identical.')"
197205
!
198206
! -- set gwfmodel
199207
mb => GetBaseModelFromList(basemodellist, this%m1_idx)
@@ -220,6 +228,34 @@ subroutine exg_ar(this)
220228
call store_error(errmsg, terminate=.TRUE.)
221229
end if
222230
!
231+
! -- Make sure idomains are identical
232+
select type (gwfdis => gwfmodel%dis)
233+
type is (DisType)
234+
select type (gwedis => gwemodel%dis)
235+
type is (DisType)
236+
if (.not. all(gwfdis%idomain == gwedis%idomain)) then
237+
write (errmsg, fmtidomerr) trim(this%name)
238+
call store_error(errmsg, terminate=.TRUE.)
239+
end if
240+
end select
241+
type is (DisvType)
242+
select type (gwedis => gwemodel%dis)
243+
type is (DisvType)
244+
if (.not. all(gwfdis%idomain == gwedis%idomain)) then
245+
write (errmsg, fmtidomerr) trim(this%name)
246+
call store_error(errmsg, terminate=.TRUE.)
247+
end if
248+
end select
249+
type is (DisuType)
250+
select type (gwedis => gwemodel%dis)
251+
type is (DisuType)
252+
if (.not. all(gwfdis%idomain == gwedis%idomain)) then
253+
write (errmsg, fmtidomerr) trim(this%name)
254+
call store_error(errmsg, terminate=.TRUE.)
255+
end if
256+
end select
257+
end select
258+
!
223259
! -- setup pointers to gwf variables allocated in gwf_ar
224260
gwemodel%fmi%gwfhead => gwfmodel%x
225261
call mem_checkin(gwemodel%fmi%gwfhead, &

‎src/Exchange/exg-gwfgwt.f90‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -184,6 +184,9 @@ end subroutine exg_df
184184
subroutine exg_ar(this)
185185
! -- modules
186186
use MemoryManagerModule, only: mem_checkin
187+
use DisModule, only: DisType
188+
use DisvModule, only: DisvType
189+
use DisuModule, only: DisuType
187190
! -- dummy
188191
class(GwfGwtExchangeType) :: this
189192
! -- local
@@ -197,6 +200,11 @@ subroutine exg_ar(this)
197200
& GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
198201
& GWT Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
199202
& Ensure discretization packages, including IDOMAIN, are identical.')"
203+
character(len=*), parameter :: fmtidomerr = &
204+
"('GWF and GWT Models do not have the same discretization for exchange&
205+
& ',a,'.&
206+
& GWF Model and GWT Model have different IDOMAIN arrays.&
207+
& Ensure discretization packages, including IDOMAIN, are identical.')"
200208
!
201209
! -- set gwfmodel
202210
mb => GetBaseModelFromList(basemodellist, this%m1_idx)
@@ -223,6 +231,34 @@ subroutine exg_ar(this)
223231
call store_error(errmsg, terminate=.TRUE.)
224232
end if
225233
!
234+
! -- Make sure idomains are identical
235+
select type (gwfdis => gwfmodel%dis)
236+
type is (DisType)
237+
select type (gwtdis => gwtmodel%dis)
238+
type is (DisType)
239+
if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
240+
write (errmsg, fmtidomerr) trim(this%name)
241+
call store_error(errmsg, terminate=.TRUE.)
242+
end if
243+
end select
244+
type is (DisvType)
245+
select type (gwtdis => gwtmodel%dis)
246+
type is (DisvType)
247+
if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
248+
write (errmsg, fmtidomerr) trim(this%name)
249+
call store_error(errmsg, terminate=.TRUE.)
250+
end if
251+
end select
252+
type is (DisuType)
253+
select type (gwtdis => gwtmodel%dis)
254+
type is (DisuType)
255+
if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
256+
write (errmsg, fmtidomerr) trim(this%name)
257+
call store_error(errmsg, terminate=.TRUE.)
258+
end if
259+
end select
260+
end select
261+
!
226262
! -- setup pointers to gwf variables allocated in gwf_ar
227263
gwtmodel%fmi%gwfhead => gwfmodel%x
228264
call mem_checkin(gwtmodel%fmi%gwfhead, &

‎src/Exchange/exg-gwfprt.f90‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -177,6 +177,9 @@ end subroutine exg_df
177177
subroutine exg_ar(this)
178178
! -- modules
179179
use MemoryManagerModule, only: mem_checkin
180+
use DisModule, only: DisType
181+
use DisvModule, only: DisvType
182+
use DisuModule, only: DisuType
180183
! -- dummy
181184
class(GwfPrtExchangeType) :: this
182185
! -- local
@@ -190,6 +193,11 @@ subroutine exg_ar(this)
190193
& GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
191194
& PRT Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
192195
& Ensure discretization packages, including IDOMAIN, are identical.')"
196+
character(len=*), parameter :: fmtidomerr = &
197+
"('GWF and PRT Models do not have the same discretization for exchange&
198+
& ',a,'.&
199+
& GWF Model and PRT Model have different IDOMAIN arrays.&
200+
& Ensure discretization packages, including IDOMAIN, are identical.')"
193201
!
194202
! -- set gwfmodel
195203
mb => GetBaseModelFromList(basemodellist, this%m1id)
@@ -216,6 +224,34 @@ subroutine exg_ar(this)
216224
call store_error(errmsg, terminate=.TRUE.)
217225
end if
218226
!
227+
! -- Make sure idomains are identical
228+
select type (gwfdis => gwfmodel%dis)
229+
type is (DisType)
230+
select type (prtdis => prtmodel%dis)
231+
type is (DisType)
232+
if (.not. all(gwfdis%idomain == prtdis%idomain)) then
233+
write (errmsg, fmtidomerr) trim(this%name)
234+
call store_error(errmsg, terminate=.TRUE.)
235+
end if
236+
end select
237+
type is (DisvType)
238+
select type (prtdis => prtmodel%dis)
239+
type is (DisvType)
240+
if (.not. all(gwfdis%idomain == prtdis%idomain)) then
241+
write (errmsg, fmtidomerr) trim(this%name)
242+
call store_error(errmsg, terminate=.TRUE.)
243+
end if
244+
end select
245+
type is (DisuType)
246+
select type (prtdis => prtmodel%dis)
247+
type is (DisuType)
248+
if (.not. all(gwfdis%idomain == prtdis%idomain)) then
249+
write (errmsg, fmtidomerr) trim(this%name)
250+
call store_error(errmsg, terminate=.TRUE.)
251+
end if
252+
end select
253+
end select
254+
!
219255
! -- setup pointers to gwf variables allocated in gwf_ar
220256
prtmodel%fmi%gwfhead => gwfmodel%x
221257
call mem_checkin(prtmodel%fmi%gwfhead, &

0 commit comments

Comments
 (0)