MODFLOW 6  version 6.7.0.dev3
USGS Modular Hydrologic Model
tsp-adv.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b
7  use basedismodule, only: disbasetype
8  use tspfmimodule, only: tspfmitype
11  ! -- Gradient schemes
12  use igradient, only: igradienttype
15  ! -- Interpolation schemes
23 
24  implicit none
25  private
26  public :: tspadvtype
27  public :: adv_cr
28 
29  type, extends(numericalpackagetype) :: tspadvtype
30  integer(I4B), pointer :: iadvwt => null() !< advection scheme. See ADV_SCHEME_* constants
31  real(dp), pointer :: ats_percel => null() !< user-specified fractional number of cells advection can move a particle during one time step
32  integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound
33  type(tspfmitype), pointer :: fmi => null() !< pointer to fmi object
34  real(dp), pointer :: eqnsclfac => null() !< governing equation scale factor; =1. for solute; =rhow*cpw for energy
35 
36  class(interpolationschemeinterface), allocatable :: face_interpolation !< interpolation scheme for face values
37  class(igradienttype), allocatable :: gradient !< cell centered gradient
38  contains
39 
40  procedure :: adv_df
41  procedure :: adv_ar
42  procedure :: adv_dt
43  procedure :: adv_fc
44  procedure :: adv_cq
45  procedure :: adv_da
46 
47  procedure :: allocate_scalars
48  procedure, private :: source_options
49 
50  end type tspadvtype
51 
52 contains
53 
54  !> @ brief Create a new ADV object
55  !!
56  !! Create a new ADV package
57  !<
58  subroutine adv_cr(advobj, name_model, input_mempath, inunit, iout, fmi, &
59  eqnsclfac)
60  ! -- dummy
61  type(tspadvtype), pointer :: advobj
62  character(len=*), intent(in) :: name_model
63  character(len=*), intent(in) :: input_mempath
64  integer(I4B), intent(in) :: inunit
65  integer(I4B), intent(in) :: iout
66  type(tspfmitype), intent(in), target :: fmi
67  real(dp), intent(in), pointer :: eqnsclfac !< governing equation scale factor
68  !
69  ! -- Create the object
70  allocate (advobj)
71  !
72  ! -- create name and memory path
73  call advobj%set_names(1, name_model, 'ADV', 'ADV', input_mempath)
74  !
75  ! -- Allocate scalars
76  call advobj%allocate_scalars()
77  !
78  ! -- Set variables
79  advobj%inunit = inunit
80  advobj%iout = iout
81  advobj%fmi => fmi
82  advobj%eqnsclfac => eqnsclfac
83  end subroutine adv_cr
84 
85  !> @brief Define ADV object
86  !!
87  !! Define the ADV package
88  !<
89  subroutine adv_df(this, adv_options)
90  ! -- dummy
91  class(tspadvtype) :: this
92  type(tspadvoptionstype), optional, intent(in) :: adv_options !< the optional options, for when not constructing from file
93  ! -- local
94  character(len=*), parameter :: fmtadv = &
95  "(1x,/1x,'ADV -- ADVECTION PACKAGE, VERSION 1, 8/25/2017', &
96  &' INPUT READ FROM MEMPATH: ', A, //)"
97  !
98  ! -- Read or set advection options
99  if (.not. present(adv_options)) then
100  !
101  ! --print a message identifying the advection package.
102  write (this%iout, fmtadv) this%input_mempath
103  !
104  ! --read options from file
105  call this%source_options()
106  else
107  !
108  ! --set options from input arg
109  this%iadvwt = adv_options%iAdvScheme
110  end if
111  end subroutine adv_df
112 
113  !> @brief Allocate and read method for package
114  !!
115  !! Method to allocate and read static data for the ADV package.
116  !<
117  subroutine adv_ar(this, dis, ibound)
118  ! -- modules
119  use simmodule, only: store_error
120  ! -- dummy
121  class(tspadvtype) :: this
122  class(disbasetype), pointer, intent(in) :: dis
123  integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibound
124  ! -- local
125  integer(I4B) :: iadvwt_value
126  class(igradienttype), allocatable :: gradient
127  ! -- adv pointers to arguments that were passed in
128  this%dis => dis
129  this%ibound => ibound
130  !
131  ! -- Create interpolation scheme
132  iadvwt_value = this%iadvwt ! Dereference iadvwt to work with case statement
133  select case (iadvwt_value)
134  case (adv_scheme_upstream)
135  this%face_interpolation = &
136  upstreamschemetype(this%dis, this%fmi)
137  case (adv_scheme_central)
138  this%face_interpolation = &
139  centraldifferenceschemetype(this%dis, this%fmi)
140  case (adv_scheme_tvd)
141  this%face_interpolation = &
142  tvdschemetype(this%dis, this%fmi, this%ibound)
143  case (adv_scheme_utvd)
144  gradient = leastsquaresgradienttype(this%dis)
145  this%gradient = cachedgradienttype(gradient, this%dis)
146  this%face_interpolation = &
147  utvdschemetype(this%dis, this%fmi, this%gradient)
148  case default
149  call store_error("Unknown advection scheme", terminate=.true.)
150  end select
151  end subroutine adv_ar
152 
153  !> @brief Calculate maximum time step length
154  !!
155  !! Return the largest time step that meets stability constraints
156  !<
157  subroutine adv_dt(this, dtmax, msg, thetam)
158  ! dummy
159  class(tspadvtype) :: this !< this instance
160  real(DP), intent(out) :: dtmax !< maximum allowable dt subject to stability constraint
161  character(len=*), intent(inout) :: msg !< package/cell dt constraint message
162  real(DP), dimension(:), intent(in) :: thetam !< porosity
163  ! local
164  integer(I4B) :: n
165  integer(I4B) :: m
166  integer(I4B) :: ipos
167  integer(I4B) :: nrmax
168  character(len=LINELENGTH) :: cellstr
169  real(DP) :: dt
170  real(DP) :: flowmax
171  real(DP) :: flowsumpos
172  real(DP) :: flowsumneg
173  real(DP) :: flownm
174  real(DP) :: cell_volume
175  dtmax = dnodata
176  nrmax = 0
177  msg = ''
178 
179  ! If ats_percel not specified by user, then return without making
180  ! the courant time step calculation
181  if (this%ats_percel == dnodata) then
182  return
183  end if
184 
185  ! Calculate time step lengths based on stability constraint for each cell
186  ! and store the smallest one
187  do n = 1, this%dis%nodes
188  if (this%ibound(n) == 0) cycle
189  flowsumneg = dzero
190  flowsumpos = dzero
191  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
192  if (this%dis%con%mask(ipos) == 0) cycle
193  m = this%dis%con%ja(ipos)
194  if (this%ibound(m) == 0) cycle
195  flownm = this%fmi%gwfflowja(ipos)
196  if (flownm < dzero) then
197  flowsumneg = flowsumneg - flownm
198  else
199  flowsumpos = flowsumpos + flownm
200  end if
201  end do
202  flowmax = max(flowsumneg, flowsumpos)
203  if (flowmax < dprec) cycle
204  cell_volume = this%dis%get_cell_volume(n, this%dis%top(n))
205  dt = cell_volume * this%fmi%gwfsat(n) * thetam(n) / flowmax
206  dt = dt * this%ats_percel
207  if (dt < dtmax) then
208  dtmax = dt
209  nrmax = n
210  end if
211  end do
212  if (nrmax > 0) then
213  call this%dis%noder_to_string(nrmax, cellstr)
214  write (msg, *) adjustl(trim(this%memoryPath))//'-'//trim(cellstr)
215  end if
216  end subroutine adv_dt
217 
218  !> @brief Fill coefficient method for ADV package
219  !!
220  !! Method to calculate coefficients and fill amat and rhs.
221  !<
222  subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs)
223  ! -- modules
224  ! -- dummy
225  class(tspadvtype) :: this !< this instance
226  integer(I4B), intent(in) :: nodes !< number of nodes
227  class(matrixbasetype), pointer :: matrix_sln !< pointer to solution matrix
228  integer(I4B), intent(in), dimension(:) :: idxglo !< global indices for matrix
229  real(DP), intent(in), dimension(:), target :: cnew !< new concentration/temperature values
230  real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector
231  ! -- local
232  integer(I4B) :: n, m, idiag, ipos
233  real(DP) :: qnm !< volumetric flow rate
234  type(coefficientstype) :: coefficients
235 
236  ! Calculate internal domain fluxes and add to matrix_sln and rhs.
237  call this%face_interpolation%set_field(cnew)
238  do n = 1, nodes
239  if (this%ibound(n) == 0) cycle ! skip inactive nodes
240  idiag = this%dis%con%ia(n)
241  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
242  m = this%dis%con%ja(ipos)
243  if (this%ibound(m) == 0) cycle ! skip inactive nodes
244  if (this%dis%con%mask(ipos) == 0) cycle ! skip masked connections
245 
246  qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
247  coefficients = this%face_interpolation%compute(n, m, ipos)
248 
249  call matrix_sln%add_value_pos(idxglo(idiag), qnm * coefficients%c_n)
250  call matrix_sln%add_value_pos(idxglo(ipos), qnm * coefficients%c_m)
251  rhs(n) = rhs(n) + qnm * coefficients%rhs
252  end do
253  end do
254  end subroutine adv_fc
255 
256  !> @brief Calculate advection contribution to flowja
257  !<
258  subroutine adv_cq(this, cnew, flowja)
259  ! -- modules
260  ! -- dummy
261  class(tspadvtype) :: this
262  real(DP), intent(in), dimension(:), target :: cnew
263  real(DP), intent(inout), dimension(:) :: flowja
264  ! -- local
265  integer(I4B) :: nodes
266  integer(I4B) :: n, m, ipos
267  real(DP) :: qnm
268  type(coefficientstype) :: coefficients
269  !
270  ! -- Calculate advection and add to flowja. qnm is the volumetric flow
271  ! rate and has dimensions of L^/T.
272  nodes = this%dis%nodes
273 
274  call this%face_interpolation%set_field(cnew)
275  do n = 1, nodes
276  if (this%ibound(n) == 0) cycle
277  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
278  m = this%dis%con%ja(ipos)
279  if (this%ibound(m) == 0) cycle
280  qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
281 
282  coefficients = this%face_interpolation%compute(n, m, ipos)
283  flowja(ipos) = flowja(ipos) &
284  + qnm * coefficients%c_n * cnew(n) &
285  + qnm * coefficients%c_m * cnew(m) &
286  - qnm * coefficients%rhs
287  end do
288  end do
289 
290  end subroutine adv_cq
291 
292  !> @brief Deallocate memory
293  !<
294  subroutine adv_da(this)
295  ! -- modules
297  ! -- dummy
298  class(tspadvtype) :: this
299  !
300  ! -- Deallocate arrays if package was active
301  if (this%inunit > 0) then
302  end if
303  !
304  ! -- nullify pointers
305  this%ibound => null()
306  !
307  ! -- Scalars
308  call mem_deallocate(this%iadvwt)
309  call mem_deallocate(this%ats_percel)
310  !
311  ! -- deallocate parent
312  call this%NumericalPackageType%da()
313  end subroutine adv_da
314 
315  !> @brief Allocate scalars specific to the streamflow energy transport (SFE)
316  !! package.
317  !<
318  subroutine allocate_scalars(this)
319  ! -- modules
321  ! -- dummy
322  class(tspadvtype) :: this
323  ! -- local
324  !
325  ! -- allocate scalars in NumericalPackageType
326  call this%NumericalPackageType%allocate_scalars()
327  !
328  ! -- Allocate
329  call mem_allocate(this%iadvwt, 'IADVWT', this%memoryPath)
330  call mem_allocate(this%ats_percel, 'ATS_PERCEL', this%memoryPath)
331  !
332  ! -- Initialize
333  this%iadvwt = adv_scheme_upstream
334  this%ats_percel = dnodata
335  !
336  ! -- Advection creates an asymmetric coefficient matrix
337  this%iasym = 1
338  end subroutine allocate_scalars
339 
340  !> @brief Source input options
341  !<
342  subroutine source_options(this)
343  ! -- modules
344  use kindmodule, only: lgp
345  use constantsmodule, only: lenvarname
346  use simvariablesmodule, only: errmsg
349  ! -- dummy
350  class(tspadvtype) :: this
351  ! -- locals
352  character(len=LENVARNAME), dimension(4) :: scheme = &
353  &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD', 'UTVD']
354  logical(lgp) :: found_scheme, found_atspercel
355  ! -- formats
356  character(len=*), parameter :: fmtiadvwt = &
357  &"(4x,'ADVECTION WEIGHTING SCHEME HAS BEEN SET TO: ', a)"
358 
359  ! update defaults with input sourced values
360  call mem_set_value(this%iadvwt, 'SCHEME', this%input_mempath, &
361  scheme, found_scheme)
362  call mem_set_value(this%ats_percel, 'ATS_PERCEL', this%input_mempath, &
363  found_atspercel)
364 
365  if (found_scheme) then
366  ! should currently be set to index of scheme names
367  if (this%iadvwt == 0) then
368  write (errmsg, '(a, a)') &
369  'Unknown scheme, must be "UPSTREAM", "CENTRAL", "TVD" or "UTVD"'
370  call store_error(errmsg)
371  call store_error_filename(this%input_fname)
372  else
373  ! scheme parameters are 0 based
374  this%iadvwt = this%iadvwt - 1
375  end if
376  end if
377 
378  if (this%iadvwt == adv_scheme_utvd) then
379  call dev_feature('UTVD is still under development, install the &
380  &nightly build or compile from source with IDEVELOPMODE = 1.')
381  end if
382 
383  if (found_atspercel) then
384  if (this%ats_percel == dzero) this%ats_percel = dnodata
385  end if
386 
387  ! log options
388  write (this%iout, '(1x,a)') 'PROCESSING ADVECTION OPTIONS'
389  if (found_scheme) then
390  write (this%iout, fmtiadvwt) trim(scheme(this%iadvwt + 1))
391  end if
392  if (found_atspercel) then
393  write (this%iout, '(4x,a,1pg15.6)') &
394  'User-specified fractional cell distance for adaptive time &
395  &steps: ', this%ats_percel
396  end if
397  write (this%iout, '(1x,a)') 'END OF ADVECTION OPTIONS'
398  end subroutine source_options
399 
400 end module tspadvmodule
integer(i4b), parameter adv_scheme_tvd
integer(i4b), parameter adv_scheme_upstream
integer(i4b), parameter adv_scheme_utvd
integer(i4b), parameter adv_scheme_central
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
real(dp), parameter dnodata
real no data constant
Definition: Constants.f90:95
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter dprec
real constant machine precision
Definition: Constants.f90:120
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
Disable development features in release mode.
Definition: DevFeature.f90:2
subroutine, public dev_feature(errmsg, iunit)
Terminate if in release mode (guard development features)
Definition: DevFeature.f90:21
This module defines variable data types.
Definition: kind.f90:8
This module contains the base numerical package type.
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
Definition: Sim.f90:203
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
subroutine, public adv_cr(advobj, name_model, input_mempath, inunit, iout, fmi, eqnsclfac)
@ brief Create a new ADV object
Definition: tsp-adv.f90:60
subroutine source_options(this)
Source input options.
Definition: tsp-adv.f90:343
subroutine adv_df(this, adv_options)
Define ADV object.
Definition: tsp-adv.f90:90
subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs)
Fill coefficient method for ADV package.
Definition: tsp-adv.f90:223
subroutine adv_cq(this, cnew, flowja)
Calculate advection contribution to flowja.
Definition: tsp-adv.f90:259
subroutine adv_da(this)
Deallocate memory.
Definition: tsp-adv.f90:295
subroutine adv_ar(this, dis, ibound)
Allocate and read method for package.
Definition: tsp-adv.f90:118
subroutine adv_dt(this, dtmax, msg, thetam)
Calculate maximum time step length.
Definition: tsp-adv.f90:158
subroutine allocate_scalars(this)
Allocate scalars specific to the streamflow energy transport (SFE) package.
Definition: tsp-adv.f90:319
Decorator that adds caching to any gradient computation implementation.
Abstract interface for cell-based gradient computation.
Definition: IGradient.f90:15
Weighted least-squares gradient method for structured and unstructured grids.
Total Variation Diminishing (TVD) interpolation scheme.
Definition: UTVDScheme.f90:37