Skip to content
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

fix(exg-gwf*): check for equal idomain #2149

Merged
merged 4 commits into from
Jan 18, 2025
Merged

fix(exg-gwf*): check for equal idomain #2149

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2e3a30d
refactor(exg-gwf*): check for equal idomain
wpbonelli Jan 17, 2025
516c0a3
release note
wpbonelli Jan 17, 2025
c3a5be7
reset fmi
wpbonelli Jan 17, 2025
4eace5d
fix comment
wpbonelli Jan 18, 2025
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
35 changes: 31 additions & 4 deletions autotest/test_prt_exg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,15 @@

Results are compared against a MODPATH 7 model.

This test includes two cases, one which gives
boundnames to particles and one which does not.
This test includes four cases: one which gives
boundnames to particles, one which does not, a
third in which the flow model has a uniformly
active idomain array while the tracking model
does not, and a final case in which flow and
tracking model have different IDOMAIN arrays,
both non-uniform, where the active region is
the same size but consists of different cells.
Both latter cases should be caught as errors.
"""

from pathlib import Path
Expand All @@ -30,7 +37,7 @@
)

simname = "prtexg01"
cases = [simname, f"{simname}bnms"]
cases = [simname, f"{simname}bnms", f"{simname}idmu", f"{simname}idmn"]


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

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

# create prt discretization
idomain = np.ones(
(FlopyReadmeCase.nlay, FlopyReadmeCase.nrow, FlopyReadmeCase.ncol)
)
if "idm" in name:
# add an inactive cell to
# tracking model idomain
idomain[-1, -1, -1] = 0
if "idmn" in name:
# add a (different) inactive
# cell to flow model idomain
gwf_idomain = idomain.copy()
gwf_idomain[-1, -1, -1] = 1
gwf_idomain[0, 0, 0] = 0
gwf.dis.idomain = gwf_idomain
flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
prt,
pname="dis",
nlay=FlopyReadmeCase.nlay,
nrow=FlopyReadmeCase.nrow,
ncol=FlopyReadmeCase.ncol,
idomain=idomain,
)

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


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

if "idm" in name:
return

# check mf6 output files exist
gwf_budget_file = f"{gwf_name}.bud"
gwf_head_file = f"{gwf_name}.hds"
Expand Down Expand Up @@ -362,5 +388,6 @@ def test_mf6model(idx, name, function_tmpdir, targets, plot):
plot=lambda t: plot_output(idx, t) if plot else None,
targets=targets,
compare=None,
xfail="idm" in name,
)
test.run()
1 change: 1 addition & 0 deletions doc/ReleaseNotes/develop.tex
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
\textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\
\underline{BASIC FUNCTIONALITY}
\begin{itemize}
\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.
\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.
\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.
\end{itemize}
Expand Down
36 changes: 36 additions & 0 deletions src/Exchange/exg-gwfgwe.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -181,6 +181,9 @@ end subroutine exg_df
subroutine exg_ar(this)
! -- modules
use MemoryManagerModule, only: mem_checkin
use DisModule, only: DisType
use DisvModule, only: DisvType
use DisuModule, only: DisuType
! -- dummy
class(GwfGweExchangeType) :: this
! -- local
Expand All @@ -194,6 +197,11 @@ subroutine exg_ar(this)
& GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
& GWE Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
& Ensure discretization packages, including IDOMAIN, are identical.')"
character(len=*), parameter :: fmtidomerr = &
"('GWF and GWE Models do not have the same discretization for exchange&
& ',a,'.&
& GWF Model and GWE Model have different IDOMAIN arrays.&
& Ensure discretization packages, including IDOMAIN, are identical.')"
!
! -- set gwfmodel
mb => GetBaseModelFromList(basemodellist, this%m1_idx)
Expand All @@ -220,6 +228,34 @@ subroutine exg_ar(this)
call store_error(errmsg, terminate=.TRUE.)
end if
!
! -- Make sure idomains are identical
select type (gwfdis => gwfmodel%dis)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All the select type statements are unfortunate, but idomain is defined separately for each dis type since it's the real grid shape instead of flat

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ahh, bummer. I guess now that all discretizations have IDOMAIN, we could move it to disbasetype and make it a 1D int vector, but that is for another time, I think.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I am understanding correctly, we are trying to assess the equality of grid. Maybe we can have a function is_equal(this, other_dis) result(are_equal) to delegate the check to the particular Discretization module?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PS: happy to help and chat more about that option

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mjr-deltares thanks we should definitely chat on this with @langevin-usgs

It may be relevant to #2153 as well. In that PR, I first planned to load the full flow model dis into fmi, then an is_equal function would be quite natural to use. But the dis types have some strong expectations about how they are initialized and where they get their data from. So in the near term, I was thinking of loading only the minimum necessary info to do the same checks. Hopefully we can come up with a nicer long term solution.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a good point. In this new check, we are saying the grids are the same if the number of nodes and idomain is the same; it is easy to accidentally leave out idomain when using flopy to construct dis. But certainly, if tops and bottoms were also different, then that would not be detected here, and problems would be very difficult to find. A discretization-based is_equal() (and maybe is_same()?) seems like a nice way to ensure grid equality.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, let's talk tomorrow and plan a meeting. Btw I had no intention stalling this PR. This should probably go in and then we refactor when we have a plan.

type is (DisType)
select type (gwedis => gwemodel%dis)
type is (DisType)
if (.not. all(gwfdis%idomain == gwedis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
type is (DisvType)
select type (gwedis => gwemodel%dis)
type is (DisvType)
if (.not. all(gwfdis%idomain == gwedis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
type is (DisuType)
select type (gwedis => gwemodel%dis)
type is (DisuType)
if (.not. all(gwfdis%idomain == gwedis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
end select
!
! -- setup pointers to gwf variables allocated in gwf_ar
gwemodel%fmi%gwfhead => gwfmodel%x
call mem_checkin(gwemodel%fmi%gwfhead, &
Expand Down
36 changes: 36 additions & 0 deletions src/Exchange/exg-gwfgwt.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,9 @@ end subroutine exg_df
subroutine exg_ar(this)
! -- modules
use MemoryManagerModule, only: mem_checkin
use DisModule, only: DisType
use DisvModule, only: DisvType
use DisuModule, only: DisuType
! -- dummy
class(GwfGwtExchangeType) :: this
! -- local
Expand All @@ -197,6 +200,11 @@ subroutine exg_ar(this)
& GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
& GWT Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
& Ensure discretization packages, including IDOMAIN, are identical.')"
character(len=*), parameter :: fmtidomerr = &
"('GWF and GWT Models do not have the same discretization for exchange&
& ',a,'.&
& GWF Model and GWT Model have different IDOMAIN arrays.&
& Ensure discretization packages, including IDOMAIN, are identical.')"
!
! -- set gwfmodel
mb => GetBaseModelFromList(basemodellist, this%m1_idx)
Expand All @@ -223,6 +231,34 @@ subroutine exg_ar(this)
call store_error(errmsg, terminate=.TRUE.)
end if
!
! -- Make sure idomains are identical
select type (gwfdis => gwfmodel%dis)
type is (DisType)
select type (gwtdis => gwtmodel%dis)
type is (DisType)
if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
type is (DisvType)
select type (gwtdis => gwtmodel%dis)
type is (DisvType)
if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
type is (DisuType)
select type (gwtdis => gwtmodel%dis)
type is (DisuType)
if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
end select
!
! -- setup pointers to gwf variables allocated in gwf_ar
gwtmodel%fmi%gwfhead => gwfmodel%x
call mem_checkin(gwtmodel%fmi%gwfhead, &
Expand Down
36 changes: 36 additions & 0 deletions src/Exchange/exg-gwfprt.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -177,6 +177,9 @@ end subroutine exg_df
subroutine exg_ar(this)
! -- modules
use MemoryManagerModule, only: mem_checkin
use DisModule, only: DisType
use DisvModule, only: DisvType
use DisuModule, only: DisuType
! -- dummy
class(GwfPrtExchangeType) :: this
! -- local
Expand All @@ -190,6 +193,11 @@ subroutine exg_ar(this)
& GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
& PRT Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
& Ensure discretization packages, including IDOMAIN, are identical.')"
character(len=*), parameter :: fmtidomerr = &
"('GWF and PRT Models do not have the same discretization for exchange&
& ',a,'.&
& GWF Model and PRT Model have different IDOMAIN arrays.&
& Ensure discretization packages, including IDOMAIN, are identical.')"
!
! -- set gwfmodel
mb => GetBaseModelFromList(basemodellist, this%m1id)
Expand All @@ -216,6 +224,34 @@ subroutine exg_ar(this)
call store_error(errmsg, terminate=.TRUE.)
end if
!
! -- Make sure idomains are identical
select type (gwfdis => gwfmodel%dis)
type is (DisType)
select type (prtdis => prtmodel%dis)
type is (DisType)
if (.not. all(gwfdis%idomain == prtdis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
type is (DisvType)
select type (prtdis => prtmodel%dis)
type is (DisvType)
if (.not. all(gwfdis%idomain == prtdis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
type is (DisuType)
select type (prtdis => prtmodel%dis)
type is (DisuType)
if (.not. all(gwfdis%idomain == prtdis%idomain)) then
write (errmsg, fmtidomerr) trim(this%name)
call store_error(errmsg, terminate=.TRUE.)
end if
end select
end select
!
! -- setup pointers to gwf variables allocated in gwf_ar
prtmodel%fmi%gwfhead => gwfmodel%x
call mem_checkin(prtmodel%fmi%gwfhead, &
Expand Down
Loading