refactor(dis): flatten idomain and move to base type

Store idomain as a 1D array in DisBaseType, rather than the proper grid shape in the concrete types. This deduplicates the variable and simplifes grid equivalence checks.
MODFLOW-USGS · Jan 22, 2025 · 5dc2505 · 5dc2505
1 parent b3d9bc9
commit 5dc2505
Show file tree

Hide file tree

Showing 12 changed files with 115 additions and 179 deletions.
diff --git a/src/Exchange/exg-gwfgwe.f90 b/src/Exchange/exg-gwfgwe.f90
@@ -229,32 +229,10 @@ subroutine exg_ar(this)
     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
+    if (.not. all(gwfmodel%dis%idomain == gwemodel%dis%idomain)) then
+      write (errmsg, fmtidomerr) trim(this%name)
+      call store_error(errmsg, terminate=.TRUE.)
+    end if
     !
     ! -- setup pointers to gwf variables allocated in gwf_ar
     gwemodel%fmi%gwfhead => gwfmodel%x

diff --git a/src/Exchange/exg-gwfgwt.f90 b/src/Exchange/exg-gwfgwt.f90
@@ -232,32 +232,10 @@ subroutine exg_ar(this)
     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
+    if (.not. all(gwfmodel%dis%idomain == gwtmodel%dis%idomain)) then
+      write (errmsg, fmtidomerr) trim(this%name)
+      call store_error(errmsg, terminate=.TRUE.)
+    end if
     !
     ! -- setup pointers to gwf variables allocated in gwf_ar
     gwtmodel%fmi%gwfhead => gwfmodel%x

diff --git a/src/Exchange/exg-gwfprt.f90 b/src/Exchange/exg-gwfprt.f90
@@ -225,32 +225,10 @@ subroutine exg_ar(this)
     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
+    if (.not. all(gwfmodel%dis%idomain == prtmodel%dis%idomain)) then
+      write (errmsg, fmtidomerr) trim(this%name)
+      call store_error(errmsg, terminate=.TRUE.)
+    end if
     !
     ! -- setup pointers to gwf variables allocated in gwf_ar
     prtmodel%fmi%gwfhead => gwfmodel%x

diff --git a/src/Model/Discretization/Dis.f90 b/src/Model/Discretization/Dis.f90
@@ -28,11 +28,9 @@ module DisModule
     real(DP), dimension(:), pointer, contiguous :: delc => null() !< spacing along a column
     real(DP), dimension(:, :), pointer, contiguous :: top2d => null() !< top elevations for each cell at top of model (ncol, nrow)
     real(DP), dimension(:, :, :), pointer, contiguous :: bot3d => null() !< bottom elevations for each cell (ncol, nrow, nlay)
-    integer(I4B), dimension(:, :, :), pointer, contiguous :: idomain => null() !< idomain (ncol, nrow, nlay)
     real(DP), dimension(:, :, :), pointer :: botm => null() !< top and bottom elevations for each cell (ncol, nrow, nlay)
     real(DP), dimension(:), pointer, contiguous :: cellx => null() !< cell center x coordinate for column j
     real(DP), dimension(:), pointer, contiguous :: celly => null() !< cell center y coordinate for row i
-
   contains
 
     procedure :: dis_df => dis3d_df
@@ -244,7 +242,6 @@ subroutine source_dimensions(this)
     ! -- dummy
     class(DisType) :: this
     ! -- locals
-    integer(I4B) :: i, j, k
     type(DisFoundType) :: found
     !
     ! -- update defaults with idm sourced values
@@ -280,23 +277,13 @@ subroutine source_dimensions(this)
     ! -- Allocate delr, delc, and non-reduced vectors for dis
     call mem_allocate(this%delr, this%ncol, 'DELR', this%memoryPath)
     call mem_allocate(this%delc, this%nrow, 'DELC', this%memoryPath)
-    call mem_allocate(this%idomain, this%ncol, this%nrow, this%nlay, 'IDOMAIN', &
-                      this%memoryPath)
+    call mem_allocate(this%idomain, this%nodesuser, 'IDOMAIN', this%memoryPath)
     call mem_allocate(this%top2d, this%ncol, this%nrow, 'TOP2D', this%memoryPath)
     call mem_allocate(this%bot3d, this%ncol, this%nrow, this%nlay, 'BOT3D', &
                       this%memoryPath)
     call mem_allocate(this%cellx, this%ncol, 'CELLX', this%memoryPath)
     call mem_allocate(this%celly, this%nrow, 'CELLY', this%memoryPath)
     !
-    ! -- initialize all cells to be active (idomain = 1)
-    do k = 1, this%nlay
-      do i = 1, this%nrow
-        do j = 1, this%ncol
-          this%idomain(j, i, k) = 1
-        end do
-      end do
-    end do
-    !
   end subroutine source_dimensions
 
   !> @brief Write dimensions to list file
@@ -329,14 +316,29 @@ end subroutine log_dimensions
   subroutine source_griddata(this)
     ! -- dummy
     class(DisType) :: this
+    ! -- local
     type(DisFoundType) :: found
+    integer(I4B) :: j, i, k
+    integer(I4B), contiguous, pointer :: idomain(:, :, :)
+    !
+    ! -- initialize all cells to be active (idomain = 1)
+    allocate (idomain(this%ncol, this%nrow, this%nlay))
+    do k = 1, this%nlay
+      do i = 1, this%nrow
+        do j = 1, this%ncol
+          idomain(j, i, k) = 1
+        end do
+      end do
+    end do
     !
     ! -- update defaults with idm sourced values
     call mem_set_value(this%delr, 'DELR', this%input_mempath, found%delr)
     call mem_set_value(this%delc, 'DELC', this%input_mempath, found%delc)
     call mem_set_value(this%top2d, 'TOP', this%input_mempath, found%top)
     call mem_set_value(this%bot3d, 'BOTM', this%input_mempath, found%botm)
-    call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    call mem_set_value(idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    this%idomain = reshape(idomain, [this%nodesuser])
+    deallocate (idomain)
     !
     ! -- log simulation values
     if (this%iout > 0) then
@@ -386,7 +388,7 @@ subroutine grid_finalize(this)
     ! -- dummy
     class(DisType) :: this
     ! -- locals
-    integer(I4B) :: n, i, j, k
+    integer(I4B) :: n, nn, k, i, j
     integer(I4B) :: node
     integer(I4B) :: noder
     integer(I4B) :: nrsize
@@ -403,12 +405,8 @@ subroutine grid_finalize(this)
     !
     ! -- count active cells
     this%nodes = 0
-    do k = 1, this%nlay
-      do i = 1, this%nrow
-        do j = 1, this%ncol
-          if (this%idomain(j, i, k) > 0) this%nodes = this%nodes + 1
-        end do
-      end do
+    do i = 1, this%nodesuser
+      if (this%idomain(i) > 0) this%nodes = this%nodes + 1
     end do
     !
     ! -- Check to make sure nodes is a valid number
@@ -424,7 +422,8 @@ subroutine grid_finalize(this)
     do k = 1, this%nlay
       do i = 1, this%nrow
         do j = 1, this%ncol
-          if (this%idomain(j, i, k) < 1) cycle
+          nn = get_node(k, i, j, this%nlay, this%nrow, this%ncol)
+          if (this%idomain(nn) < 1) cycle
           if (k > 1) then
             top = this%bot3d(j, i, k - 1)
           else
@@ -461,10 +460,11 @@ subroutine grid_finalize(this)
       do k = 1, this%nlay
         do i = 1, this%nrow
           do j = 1, this%ncol
-            if (this%idomain(j, i, k) > 0) then
+            nn = get_node(k, i, j, this%nlay, this%nrow, this%ncol)
+            if (this%idomain(nn) > 0) then
               this%nodereduced(node) = noder
               noder = noder + 1
-            elseif (this%idomain(j, i, k) < 0) then
+            elseif (this%idomain(nn) < 0) then
               this%nodereduced(node) = -1
             else
               this%nodereduced(node) = 0
@@ -482,7 +482,8 @@ subroutine grid_finalize(this)
       do k = 1, this%nlay
         do i = 1, this%nrow
           do j = 1, this%ncol
-            if (this%idomain(j, i, k) > 0) then
+            nn = get_node(k, i, j, this%nlay, this%nrow, this%ncol)
+            if (this%idomain(nn) > 0) then
               this%nodeuser(noder) = node
               noder = noder + 1
             end if

diff --git a/src/Model/Discretization/Dis2d.f90 b/src/Model/Discretization/Dis2d.f90
@@ -26,7 +26,6 @@ module Dis2dModule
     real(DP), dimension(:), pointer, contiguous :: delr => null() !< spacing along a row
     real(DP), dimension(:), pointer, contiguous :: delc => null() !< spacing along a column
     real(DP), dimension(:, :), pointer, contiguous :: bottom => null() !< bottom elevations for each cell (ncol, nrow)
-    integer(I4B), dimension(:, :), pointer, contiguous :: idomain => null() !< idomain (ncol, nrow)
     real(DP), dimension(:), pointer, contiguous :: cellx => null() !< cell center x coordinate for column j
     real(DP), dimension(:), pointer, contiguous :: celly => null() !< cell center y coordinate for row i
 
@@ -237,7 +236,6 @@ subroutine source_dimensions(this)
     ! -- dummy
     class(Dis2dType) :: this
     ! -- locals
-    integer(I4B) :: i, j
     type(DisFoundType) :: found
     !
     ! -- update defaults with idm sourced values
@@ -267,20 +265,12 @@ subroutine source_dimensions(this)
     ! -- Allocate delr, delc, and non-reduced vectors for dis
     call mem_allocate(this%delr, this%ncol, 'DELR', this%memoryPath)
     call mem_allocate(this%delc, this%nrow, 'DELC', this%memoryPath)
-    call mem_allocate(this%idomain, this%ncol, this%nrow, 'IDOMAIN', &
-                      this%memoryPath)
+    call mem_allocate(this%idomain, this%nodesuser, 'IDOMAIN', this%memoryPath)
     call mem_allocate(this%bottom, this%ncol, this%nrow, 'BOTTOM', &
                       this%memoryPath)
     call mem_allocate(this%cellx, this%ncol, 'CELLX', this%memoryPath)
     call mem_allocate(this%celly, this%nrow, 'CELLY', this%memoryPath)
     !
-    ! -- initialize all cells to be active (idomain = 1)
-    do i = 1, this%nrow
-      do j = 1, this%ncol
-        this%idomain(j, i) = 1
-      end do
-    end do
-    !
   end subroutine source_dimensions
 
   !> @brief Write dimensions to list file
@@ -309,13 +299,26 @@ end subroutine log_dimensions
   subroutine source_griddata(this)
     ! -- dummy
     class(Dis2dType) :: this
+    ! -- local
     type(DisFoundType) :: found
+    integer(I4B) :: i, j
+    integer(I4B), contiguous, pointer :: idomain(:, :)
+    !
+    ! -- initialize all cells to be active (idomain = 1)
+    allocate (idomain(this%ncol, this%nrow))
+    do i = 1, this%nrow
+      do j = 1, this%ncol
+        idomain(j, i) = 1
+      end do
+    end do
     !
     ! -- update defaults with idm sourced values
     call mem_set_value(this%delr, 'DELR', this%input_mempath, found%delr)
     call mem_set_value(this%delc, 'DELC', this%input_mempath, found%delc)
     call mem_set_value(this%bottom, 'BOTTOM', this%input_mempath, found%bottom)
-    call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    call mem_set_value(idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    this%idomain = reshape(idomain, [this%nodesuser])
+    deallocate (idomain)
     !
     ! -- log simulation values
     if (this%iout > 0) then
@@ -376,10 +379,8 @@ subroutine grid_finalize(this)
     !
     ! -- count active cells
     this%nodes = 0
-    do i = 1, this%nrow
-      do j = 1, this%ncol
-        if (this%idomain(j, i) > 0) this%nodes = this%nodes + 1
-      end do
+    do i = 1, this%nodesuser
+      if (this%idomain(i) > 0) this%nodes = this%nodes + 1
     end do
     !
     ! -- Check to make sure nodes is a valid number
@@ -407,10 +408,10 @@ subroutine grid_finalize(this)
       noder = 1
       do i = 1, this%nrow
         do j = 1, this%ncol
-          if (this%idomain(j, i) > 0) then
+          if (this%idomain(node) > 0) then
             this%nodereduced(node) = noder
             noder = noder + 1
-          elseif (this%idomain(j, i) < 0) then
+          elseif (this%idomain(node) < 0) then
             this%nodereduced(node) = -1
           else
             this%nodereduced(node) = 0
@@ -426,7 +427,7 @@ subroutine grid_finalize(this)
       noder = 1
       do i = 1, this%nrow
         do j = 1, this%ncol
-          if (this%idomain(j, i) > 0) then
+          if (this%idomain(node) > 0) then
             this%nodeuser(noder) = node
             noder = noder + 1
           end if

diff --git a/src/Model/Discretization/Disu.f90 b/src/Model/Discretization/Disu.f90
@@ -43,9 +43,7 @@ module DisuModule
     integer(I4B), pointer :: iangledegx => null() ! =1 when angle information was present in input, 0 otherwise
     integer(I4B), dimension(:), pointer, contiguous :: iavert => null() ! cell vertex pointer ia array
     integer(I4B), dimension(:), pointer, contiguous :: javert => null() ! cell vertex pointer ja array
-    integer(I4B), dimension(:), pointer, contiguous :: idomain => null() ! idomain (nodes)
     logical(LGP) :: readFromFile ! True, when DIS is read from file (almost always)
-
   contains
 
     procedure :: dis_df => disu_df