refactor(exg-gwf*): check for equal idomain

Check that models connected to a flow model via an exchange have an identical idomain array to the flow model's.

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,31 @@ 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:
+        idomain[-1, -1, -1] = 0
+    if "idmn" in name:
+        # TODO alter flow model's 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 +175,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 +198,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 +385,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()

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

src/Model/ModelUtilities/FlowModelInterface.f90

@@ -23,7 +23,7 @@ module FlowModelInterfaceModule
     logical, pointer :: flows_from_file => null() !< if .false., then flows come from GWF through GWF-Model exg
     type(ListType), pointer :: gwfbndlist => null() !< list of gwf stress packages
     integer(I4B), pointer :: iflowsupdated => null() !< flows were updated for this time step
-    integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to Model ibound
+    integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound
     real(DP), dimension(:), pointer, contiguous :: gwfflowja => null() !< pointer to the GWF flowja array
     real(DP), dimension(:, :), pointer, contiguous :: gwfspdis => null() !< pointer to npf specific discharge array
     real(DP), dimension(:), pointer, contiguous :: gwfhead => null() !< pointer to the GWF head array