fix(exg-gwf*): check for equal idomain (MODFLOW-USGS#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.

MODFLOW-USGS#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.

autotest/test_prt_exg.py

@@ -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
@@ -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):
@@ -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
@@ -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):
@@ -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"
@@ -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()

doc/ReleaseNotes/develop.tex

@@ -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}

src/Exchange/exg-gwfgwe.f90

@@ -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
@@ -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)
@@ -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)
+    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, &

src/Exchange/exg-gwfgwt.f90

@@ -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
@@ -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)
@@ -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, &

src/Exchange/exg-gwfprt.f90

@@ -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
@@ -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)
@@ -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, &