feat(cmask): add mask for calculation of amat terms (#215)

* feat(cmask): skip calculation of amat terms in *_fc and *_fn for masked connections

The purpose of this PR is to allow another package or module to shutoff the flow connection between two cells.  A connection can be masked by calling set_mask with the position within the ja and amat arrays of the masked connection.  A value of zero means the connection is masked.  If not used, the masked array points to ja, which has non-zero values for all connections.  The first call to set_mask allocates new memory for the masked array of size nja and initializes all values to one.

src/Model/GroundWaterFlow/gwf3npf8.f90

@@ -83,6 +83,7 @@ module GwfNpfModule
     procedure                               :: npf_ac
     procedure                               :: npf_mc
     procedure                               :: npf_ar
+    procedure                               :: npf_init_mem
     procedure                               :: npf_ad
     procedure                               :: npf_cf
     procedure                               :: npf_fc
@@ -469,6 +470,8 @@ subroutine npf_fc(this, kiter, njasln, amat, idxglo, rhs, hnew)
     !
     do n = 1, this%dis%nodes
       do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
+        if (this%dis%con%mask(ii) == 0) cycle
+
         m = this%dis%con%ja(ii)
         !
         ! -- Calculate conductance only for upper triangle but insert into
@@ -596,6 +599,8 @@ subroutine npf_fn(this, kiter, njasln, amat, idxglo, rhs, hnew)
     do n=1, nodes
       idiag=this%dis%con%ia(n)
       do ii=this%dis%con%ia(n)+1,this%dis%con%ia(n+1)-1
+        if (this%dis%con%mask(ii) == 0) cycle
+
         m=this%dis%con%ja(ii)
         isymcon = this%dis%con%isym(ii)
         ! work on upper triangle

src/Model/ModelUtilities/Connections.f90

@@ -18,6 +18,7 @@ module ConnectionsModule
     integer(I4B), pointer                           :: ianglex    => null()      !indicates whether or not anglex was read
     integer(I4B), dimension(:), pointer, contiguous :: ia         => null()      !(size:nodes+1) csr index array
     integer(I4B), dimension(:), pointer, contiguous :: ja         => null()      !(size:nja) csr pointer array
+    integer(I4B), dimension(:), pointer, contiguous :: mask       => null()      !(size:nja) to mask certain connections: ==0 means masked. Do not set the mask directly, use set_mask instead!    
     real(DP), dimension(:), pointer, contiguous     :: cl1        => null()      !(size:njas) connection length between node n and shared face with node m
     real(DP), dimension(:), pointer, contiguous     :: cl2        => null()      !(size:njas) connection length between node m and shared face with node n
     real(DP), dimension(:), pointer, contiguous     :: hwva       => null()      !(size:njas) horizontal perpendicular width (ihc>0) or vertical flow area (ihc=0)
@@ -40,6 +41,7 @@ module ConnectionsModule
     procedure :: disvconnections
     procedure :: iajausr
     procedure :: getjaindex
+    procedure :: set_mask
   end type ConnectionsType
 
   contains
@@ -78,17 +80,23 @@ subroutine con_da(this)
     else
       call mem_deallocate(this%jausr)
     endif
+    ! -- mask
+    if (associated(this%mask, this%ja)) then
+      nullify(this%mask)
+    else
+      call mem_deallocate(this%mask)
+    end if 
     !
     ! -- Arrays
     call mem_deallocate(this%ia)
-    call mem_deallocate(this%ja)
+    call mem_deallocate(this%ja)  
     call mem_deallocate(this%isym)
     call mem_deallocate(this%jas)
     call mem_deallocate(this%hwva)
     call mem_deallocate(this%anglex)
     call mem_deallocate(this%ihc)
     call mem_deallocate(this%cl1)
-    call mem_deallocate(this%cl2)
+    call mem_deallocate(this%cl2)     
     !
     ! -- return
     return
@@ -150,11 +158,15 @@ subroutine allocate_arrays(this)
     call mem_allocate(this%cl2, this%njas, 'CL2', this%cid)
     call mem_allocate(this%iausr, 1, 'IAUSR', this%cid)
     call mem_allocate(this%jausr, 1, 'JAUSR', this%cid)
+    ! 
+    ! let mask point to ja, which is always nonzero, 
+    ! until someone decides to do a 'set_mask'
+    this%mask => this%ja
     !
     ! -- Return
     return
   end subroutine allocate_arrays
-
+  
   subroutine read_from_block(this, name_model, nodes, nja, inunit, iout)
 ! ******************************************************************************
 ! read_from_block -- Read connection information from input block
@@ -1289,5 +1301,32 @@ subroutine vertexconnect(nodes, nrsize, maxnnz, nlay, ncpl, sparse,          &
     return
   end subroutine vertexconnect
 
+  subroutine set_mask(this, ipos, maskval)
+! ******************************************************************************
+! set_mask -- routine to set a value in the mask array 
+! (which has the same shape as this%ja)
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------ 
+  use MemoryManagerModule, only: mem_allocate
+  class(ConnectionsType) :: this
+  integer(I4B), intent(in) :: ipos
+  integer(I4B), intent(in) :: maskval
+  ! local
+  integer(I4B) :: i
+
+  ! if we still point to this%ja, we first need to allocate space
+  if (associated(this%mask, this%ja)) then
+    call mem_allocate(this%mask, this%nja, 'MASK', this%cid)
+    ! and initialize with unmasked
+    do i = 1, this%nja
+      this%mask(i) = 1 
+    end do
+  end if
+
+  this%mask(ipos) = maskVal
 
+  end subroutine set_mask
+
 end module ConnectionsModule

src/Model/ModelUtilities/Xt3dInterface.f90

@@ -381,7 +381,7 @@ subroutine xt3d_fc(this, kiter, njasln, amat, idxglo, rhs, hnew)
     real(DP),intent(inout),dimension(:) :: hnew
     ! -- local
     integer(I4B) :: nodes, nja
-    integer(I4B) :: n, m
+    integer(I4B) :: n, m, ipos
     !
     logical :: allhc0, allhc1
     integer(I4B) :: nnbr0, nnbr1
@@ -425,6 +425,9 @@ subroutine xt3d_fc(this, kiter, njasln, amat, idxglo, rhs, hnew)
       ! -- Loop over active neighbors of cell 0 that have a higher
       ! -- cell number (taking advantage of reciprocity).
       do il0 = 1,nnbr0
+        ipos = this%dis%con%ia(n) + il0
+        if (this%dis%con%mask(ipos) == 0) cycle
+
         m = inbr0(il0)
         ! -- Skip if neighbor is inactive or has lower cell number.
         if ((m.eq.0).or.(m.lt.n)) cycle
@@ -512,7 +515,7 @@ subroutine xt3d_fcpc(this, nodes)
     class(Xt3dType) :: this
     integer(I4B) :: nodes
     ! -- local
-    integer(I4B) :: n, m
+    integer(I4B) :: n, m, ipos
     !
     logical :: allhc0, allhc1
     integer(I4B) :: nnbr0, nnbr1
@@ -538,6 +541,9 @@ subroutine xt3d_fcpc(this, nodes)
       ! -- Loop over active neighbors of cell 0 that have a higher
       ! -- cell number (taking advantage of reciprocity).
       do il0 = 1,nnbr0
+        ipos = this%dis%con%ia(n) + il0
+        if (this%dis%con%mask(ipos) == 0) cycle
+
         m = inbr0(il0)
         ! -- Skip if neighbor has lower cell number.
         if (m.lt.n) cycle
@@ -719,7 +725,7 @@ subroutine xt3d_fn(this, kiter, nodes, nja, njasln, amat, idxglo, rhs, hnew)
     real(DP),intent(inout),dimension(nodes) :: rhs
     real(DP),intent(inout),dimension(nodes) :: hnew
     ! -- local
-    integer(I4B) :: n, m
+    integer(I4B) :: n, m, ipos
     !
     integer(I4B) :: nnbr0
     integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
@@ -744,6 +750,9 @@ subroutine xt3d_fn(this, kiter, nodes, nja, njasln, amat, idxglo, rhs, hnew)
       ! -- Loop over active neighbors of cell 0 that have a higher
       ! -- cell number (taking advantage of reciprocity).
       do il0 = 1,nnbr0
+        ipos = this%dis%con%ia(n) + il0
+        if (this%dis%con%mask(ipos) == 0) cycle
+
         m = inbr0(il0)
         ! -- Skip if neighbor is inactive or has lower cell number.
         if ((inbr0(il0).eq.0).or.(m.lt.n)) cycle