bug(csub): Write better cellid format to strain csv files (#355)

Write integer cellids to csv files instead of string. Add method to discretization classes to return array of user indices. Refactored table output to use table object instead of UWWORD. closes #299
MODFLOW-USGS · Mar 20, 2020 · 09e411d · 09e411d
1 parent fc27013
commit 09e411d
Show file tree

Hide file tree

Showing 9 changed files with 740 additions and 445 deletions.
diff --git a/autotest/test_gwf_csub_subwt01.py b/autotest/test_gwf_csub_subwt01.py
@@ -143,15 +143,18 @@ def get_model(idx, dir):
 
     # build MODFLOW 6 files
     ws = dir
-    sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
+    sim = flopy.mf6.MFSimulation(sim_name=name,
+                                 memory_print_option='all',
+                                 version='mf6',
                                  exe_name='mf6',
                                  sim_ws=ws)
     # create tdis package
     tdis = flopy.mf6.ModflowTdis(sim, time_units='DAYS',
                                  nper=nper, perioddata=tdis_rc)
 
     # create gwf model
-    gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True,
+                               print_input=True)
 
     # create iterative model solution and register the gwf model with it
     ims = flopy.mf6.ModflowIms(sim, print_option='SUMMARY',
@@ -200,8 +203,11 @@ def get_model(idx, dir):
             gg.append([(0, i, j), gs0[idx]])
     sig0 = {0: gg}
     opth = '{}.csub.obs'.format(name)
+    fcgstrain = '{}.csub.strain.cg.csv'.format(name)
+    fibstrain = '{}.csub.strain.ib.csv'.format(name)
     csub = flopy.mf6.ModflowGwfcsub(gwf,
-                                    #interbed_stress_offset=True,
+                                    straincg_filerecord=fcgstrain,
+                                    strainib_filerecord=fibstrain,
                                     boundnames=True,
                                     compression_indices=True,
                                     update_material_properties=ump[idx],

diff --git a/src/Model/GroundWaterFlow/gwf3csub8.f90 b/src/Model/GroundWaterFlow/gwf3csub8.f90
diff --git a/src/Model/GroundWaterFlow/gwf3dis8.f90 b/src/Model/GroundWaterFlow/gwf3dis8.f90
@@ -39,6 +39,7 @@ module GwfDisModule
     procedure :: get_nodenumber_idx1
     procedure :: get_nodenumber_idx3
     procedure :: nodeu_to_string
+    procedure :: nodeu_to_array
     procedure :: nodeu_from_string
     procedure :: nodeu_from_cellid
     procedure :: supports_layers
@@ -862,6 +863,47 @@ subroutine nodeu_to_string(this, nodeu, str)
     return
   end subroutine nodeu_to_string
 
+  subroutine nodeu_to_array(this, nodeu, arr)
+! ******************************************************************************
+! nodeu_to_array -- Convert user node number to cellid and fill array with
+!                   (nodenumber) or (k,j) or (k,i,j)
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    use InputOutputModule, only: get_ijk
+    implicit none
+    class(GwfDisType) :: this
+    integer(I4B), intent(in) :: nodeu
+    integer(I4B), dimension(:), intent(inout) :: arr
+    ! -- local
+    character(len=LINELENGTH) :: errmsg
+    integer(I4B) :: isize
+    integer(I4B) :: i, j, k
+! ------------------------------------------------------------------------------
+    !
+    ! -- check the size of arr
+    isize = size(arr)
+    if (isize /= this%ndim) then
+      write(errmsg,'(a,i0,a,i0,a)')                                              &
+        'Program error: nodeu_to_array size of array (', isize,                  &
+        ') is not equal to the discretization dimension (', this%ndim, ')'
+      call store_error(errmsg)
+      call ustop()
+    end if
+    !
+    ! -- get k, i, j
+    call get_ijk(nodeu, this%nrow, this%ncol, this%nlay, i, j, k)
+    !
+    ! -- fill array
+    arr(1) = k
+    arr(2) = i
+    arr(3) = j
+    !
+    ! -- return
+    return
+  end subroutine nodeu_to_array
+
   function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber)
 ! ******************************************************************************
 ! get_nodenumber -- Return a nodenumber from the user specified node number

diff --git a/src/Model/GroundWaterFlow/gwf3disu8.f90 b/src/Model/GroundWaterFlow/gwf3disu8.f90
@@ -43,6 +43,7 @@ module GwfDisuModule
     procedure :: grid_finalize
     procedure :: get_nodenumber_idx1
     procedure :: nodeu_to_string
+    procedure :: nodeu_to_array
     procedure :: nodeu_from_string
     procedure :: nodeu_from_cellid
     procedure :: connection_normal
@@ -513,6 +514,42 @@ subroutine nodeu_to_string(this, nodeu, str)
     return
   end subroutine nodeu_to_string
 
+  subroutine nodeu_to_array(this, nodeu, arr)
+! ******************************************************************************
+! nodeu_to_array -- Convert user node number to cellid and fill array with
+!                   (nodenumber) or (k,j) or (k,i,j)
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    use InputOutputModule, only: get_ijk
+    implicit none
+    class(GwfDisuType) :: this
+    integer(I4B), intent(in) :: nodeu
+    integer(I4B), dimension(:), intent(inout) :: arr
+    ! -- local
+    character(len=LINELENGTH) :: errmsg
+    integer(I4B) :: isize
+    integer(I4B) :: i, j, k
+! ------------------------------------------------------------------------------
+    !
+    ! -- check the size of arr
+    isize = size(arr)
+    if (isize /= this%ndim) then
+      write(errmsg,'(a,i0,a,i0,a)')                                              &
+        'Program error: nodeu_to_array size of array (', isize,                  &
+        ') is not equal to the discretization dimension (', this%ndim, ')'
+      call store_error(errmsg)
+      call ustop()
+    end if
+    !
+    ! -- fill array
+    arr(1) = nodeu
+    !
+    ! -- return
+    return
+  end subroutine nodeu_to_array
+
   subroutine read_options(this)
 ! ******************************************************************************
 ! read_options -- Read discretization options

diff --git a/src/Model/GroundWaterFlow/gwf3disv8.f90 b/src/Model/GroundWaterFlow/gwf3disv8.f90
@@ -41,6 +41,7 @@ module GwfDisvModule
     procedure :: get_nodenumber_idx1
     procedure :: get_nodenumber_idx2
     procedure :: nodeu_to_string
+    procedure :: nodeu_to_array
     procedure :: nodeu_from_string
     procedure :: nodeu_from_cellid
     procedure :: connection_normal
@@ -1134,6 +1135,46 @@ subroutine nodeu_to_string(this, nodeu, str)
     return
   end subroutine nodeu_to_string
 
+  subroutine nodeu_to_array(this, nodeu, arr)
+! ******************************************************************************
+! nodeu_to_array -- Convert user node number to cellid and fill array with
+!                   (nodenumber) or (k, j) or (k,i,j)
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    use InputOutputModule, only: get_ijk
+    implicit none
+    class(GwfDisvType) :: this
+    integer(I4B), intent(in) :: nodeu
+    integer(I4B), dimension(:), intent(inout) :: arr
+    ! -- local
+    character(len=LINELENGTH) :: errmsg
+    integer(I4B) :: isize
+    integer(I4B) :: i, j, k
+! ------------------------------------------------------------------------------
+    !
+    ! -- check the size of arr
+    isize = size(arr)
+    if (isize /= this%ndim) then
+      write(errmsg,'(a,i0,a,i0,a)')                                              &
+        'Program error: nodeu_to_array size of array (', isize,                  &
+        ') is not equal to the discretization dimension (', this%ndim, ')'
+      call store_error(errmsg)
+      call ustop()
+    end if
+    !
+    ! -- get k, i, j
+    call get_ijk(nodeu, 1, this%ncpl, this%nlay, i, j, k)
+    !
+    ! -- fill array
+    arr(1) = k
+    arr(2) = j
+    !
+    ! -- return
+    return
+  end subroutine nodeu_to_array
+
   function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber)
 ! ******************************************************************************
 ! get_nodenumber -- Return a nodenumber from the user specified node number

diff --git a/src/Model/GroundWaterFlow/gwf3sto8.f90 b/src/Model/GroundWaterFlow/gwf3sto8.f90
@@ -606,7 +606,7 @@ subroutine sto_da(this)
     class(GwfStoType) :: this
 ! ------------------------------------------------------------------------------
     !
-    ! -- Deallocate arrays if package was active
+    ! -- Deallocate arrays if package is active
     if(this%inunit > 0) then
       call mem_deallocate(this%iconvert)
       call mem_deallocate(this%sc1)
@@ -615,7 +615,7 @@ subroutine sto_da(this)
       call mem_deallocate(this%strgsy)
     endif
     !
-    ! -- Scalars
+    ! -- Deallocate scalars
     call mem_deallocate(this%isfac)
     call mem_deallocate(this%isseg)
     call mem_deallocate(this%satomega)

diff --git a/src/Model/ModelUtilities/DiscretizationBase.f90 b/src/Model/ModelUtilities/DiscretizationBase.f90
@@ -64,6 +64,7 @@ module BaseDisModule
     procedure :: get_nodenumber_idx3
     procedure :: get_nodeuser
     procedure :: nodeu_to_string
+    procedure :: nodeu_to_array
     procedure :: nodeu_from_string
     procedure :: nodeu_from_cellid
     procedure :: noder_from_string
@@ -90,6 +91,7 @@ module BaseDisModule
     procedure, public  :: record_array
     procedure, public  :: record_connection_array
     procedure, public  :: noder_to_string
+    procedure, public  :: noder_to_array
     procedure, public  :: record_srcdst_list_header
     procedure, private :: record_srcdst_list_entry
     generic, public    :: record_mf6_list_entry => record_srcdst_list_entry
@@ -295,7 +297,7 @@ end subroutine dis_da
 
   subroutine nodeu_to_string(this, nodeu, str)
 ! ******************************************************************************
-! noder_to_string -- Convert user node number to a string in the form of
+! nodeu_to_string -- Convert user node number to a string in the form of
 ! (nodenumber) or (k,i,j)
 ! ******************************************************************************
 !
@@ -316,6 +318,29 @@ subroutine nodeu_to_string(this, nodeu, str)
     return
   end subroutine nodeu_to_string
 
+  subroutine nodeu_to_array(this, nodeu, arr)
+! ******************************************************************************
+! nodeu_to_array -- Convert user node number to cellid and fill array with
+!                   (nodenumber) or (k,j) or (k,i,j) 
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- dummy
+    class(DisBaseType) :: this
+    integer(I4B), intent(in) :: nodeu
+    integer(I4B), dimension(:), intent(inout) :: arr
+    ! -- local
+! ------------------------------------------------------------------------------
+    !
+    call store_error('Program error: DisBaseType method nodeu_to_array not &
+                     &implemented.')
+    call ustop()
+    !
+    ! -- return
+    return
+  end subroutine nodeu_to_array
+
   function get_nodeuser(this, noder) result(nodenumber)
 ! ******************************************************************************
 ! get_nodeuser -- Return the user nodenumber from the reduced node number
@@ -1250,6 +1275,30 @@ subroutine noder_to_string(this, noder, str)
     return
   end subroutine noder_to_string
 
+  subroutine noder_to_array(this, noder, arr)
+! ******************************************************************************
+! noder_to_array -- Convert reduced node number to cellid and fill array with
+!                   (nodenumber) or (k,j) or (k,i,j) 
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    ! -- dummy
+    class(DisBaseType) :: this
+    integer(I4B), intent(in) :: noder
+    integer(I4B), dimension(:), intent(inout) :: arr
+    ! -- local
+    integer(I4B) :: nodeu
+! ------------------------------------------------------------------------------
+    !
+    nodeu = this%get_nodeuser(noder)
+    call this%nodeu_to_array(nodeu, arr)
+    !
+    ! -- return
+    return
+  end subroutine noder_to_array
+
   subroutine record_srcdst_list_header(this, text, textmodel, textpackage,      &
                                        dstmodel, dstpackage, naux, auxtxt,      &
                                        ibdchn, nlist, iout)

diff --git a/src/Timing/tdis.f90 b/src/Timing/tdis.f90
@@ -247,14 +247,14 @@ subroutine tdis_ot(iout)
 !C
 !C5------PRINT TIME STEP LENGTH AND ELAPSED TIMES IN ALL TIME UNITS.
       WRITE(IOUT,200)
-  200 FORMAT(19X,' SECONDS     MINUTES      HOURS',7X, &
+  200 FORMAT(19X,' SECONDS     MINUTES      HOURS',7X,                           &
      &    'DAYS        YEARS'/20X,59('-'))
       write(IOUT,201) DELSEC,DELMN,DELHR,DELDY,DELYR
   201 FORMAT(1X,'  TIME STEP LENGTH',1P,5G12.5)
       WRITE(IOUT,202) PERSEC,PERMN,PERHR,PERDY,PERYR
   202 FORMAT(1X,'STRESS PERIOD TIME',1P,5G12.5)
       WRITE(IOUT,203) TOTSEC,TOTMN,TOTHR,TOTDY,TOTYR
-  203 FORMAT(1X,'        TOTAL TIME',1P,5G12.5)
+  203 FORMAT(1X,'        TOTAL TIME',1P,5G12.5,/)
 !C
 !C6------RETURN
       RETURN