tspadvmodule Module Reference

Data Types
type	tspadvtype

Functions/Subroutines
subroutine, public	adv_cr (advobj, name_model, input_mempath, inunit, iout, fmi, eqnsclfac)
	@ brief Create a new ADV object More...
subroutine	adv_df (this, adv_options)
	Define ADV object. More...
subroutine	adv_ar (this, dis, ibound)
	Allocate and read method for package. More...
subroutine	adv_dt (this, dtmax, msg, thetam)
	Calculate maximum time step length. More...
subroutine	adv_fc (this, nodes, matrix_sln, idxglo, cnew, rhs)
	Fill coefficient method for ADV package. More...
subroutine	adv_cq (this, cnew, flowja)
	Calculate advection contribution to flowja. More...
subroutine	adv_da (this)
	Deallocate memory. More...
subroutine	allocate_scalars (this)
	Allocate scalars specific to the streamflow energy transport (SFE) package. More...
subroutine	source_options (this)
	Source input options. More...

Function/Subroutine Documentation

adv_ar()

subroutine tspadvmodule::adv_ar ( class(tspadvtype)  this, class(disbasetype), intent(in), pointer  dis, integer(i4b), dimension(:), intent(in), pointer, contiguous  ibound  )

private

Method to allocate and read static data for the ADV package.

Definition at line 117 of file tsp-adv.f90.

    ! -- modules
    use simmodule, only: store_error
    ! -- dummy
    class(TspAdvType) :: this
    class(DisBaseType), pointer, intent(in) :: dis
    integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibound
    ! -- local
    integer(I4B) :: iadvwt_value
    class(IGradientType), allocatable :: gradient
    ! -- adv pointers to arguments that were passed in
    this%dis => dis
    this%ibound => ibound
    !
    ! -- Create interpolation scheme
    iadvwt_value = this%iadvwt ! Dereference iadvwt to work with case statement
    select case (iadvwt_value)
    case (adv_scheme_upstream)
      this%face_interpolation = &
        upstreamschemetype(this%dis, this%fmi)
    case (adv_scheme_central)
      this%face_interpolation = &
        centraldifferenceschemetype(this%dis, this%fmi)
    case (adv_scheme_tvd)
      this%face_interpolation = &
        tvdschemetype(this%dis, this%fmi, this%ibound)
    case (adv_scheme_utvd)
      gradient = leastsquaresgradienttype(this%dis)
      this%gradient = cachedgradienttype(gradient, this%dis)
      this%face_interpolation = &
        utvdschemetype(this%dis, this%fmi, this%gradient)
    case default
      call store_error("Unknown advection scheme", terminate=.true.)
    end select

Here is the call graph for this function:

adv_cq()

subroutine tspadvmodule::adv_cq ( class(tspadvtype)  this, real(dp), dimension(:), intent(in), target  cnew, real(dp), dimension(:), intent(inout)  flowja  )

private

Definition at line 258 of file tsp-adv.f90.

    ! -- modules
    ! -- dummy
    class(TspAdvType) :: this
    real(DP), intent(in), dimension(:), target :: cnew
    real(DP), intent(inout), dimension(:) :: flowja
    ! -- local
    integer(I4B) :: nodes
    integer(I4B) :: n, m, ipos
    real(DP) :: qnm
    type(CoefficientsType) :: coefficients
    !
    ! -- Calculate advection and add to flowja. qnm is the volumetric flow
    !    rate and has dimensions of L^/T.
    nodes = this%dis%nodes
 
    call this%face_interpolation%set_field(cnew)
    do n = 1, nodes
      if (this%ibound(n) == 0) cycle
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        m = this%dis%con%ja(ipos)
        if (this%ibound(m) == 0) cycle
        qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
 
        coefficients = this%face_interpolation%compute(n, m, ipos)
        flowja(ipos) = flowja(ipos) &
                       + qnm * coefficients%c_n * cnew(n) &
                       + qnm * coefficients%c_m * cnew(m) &
                       - qnm * coefficients%rhs
      end do
    end do
 

adv_cr()

subroutine, public tspadvmodule::adv_cr	(	type(tspadvtype), pointer	advobj,
		character(len=*), intent(in)	name_model,
		character(len=*), intent(in)	input_mempath,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout,
		type(tspfmitype), intent(in), target	fmi,
		real(dp), intent(in), pointer	eqnsclfac
	)

Create a new ADV package

Parameters

[in]

eqnsclfac

governing equation scale factor

Definition at line 58 of file tsp-adv.f90.

    ! -- dummy
    type(TspAdvType), pointer :: advobj
    character(len=*), intent(in) :: name_model
    character(len=*), intent(in) :: input_mempath
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    type(TspFmiType), intent(in), target :: fmi
    real(DP), intent(in), pointer :: eqnsclfac !< governing equation scale factor
    !
    ! -- Create the object
    allocate (advobj)
    !
    ! -- create name and memory path
    call advobj%set_names(1, name_model, 'ADV', 'ADV', input_mempath)
    !
    ! -- Allocate scalars
    call advobj%allocate_scalars()
    !
    ! -- Set variables
    advobj%inunit = inunit
    advobj%iout = iout
    advobj%fmi => fmi
    advobj%eqnsclfac => eqnsclfac

Here is the caller graph for this function:

adv_da()

subroutine tspadvmodule::adv_da ( class(tspadvtype)  this )

private

Definition at line 294 of file tsp-adv.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(TspAdvType) :: this
    !
    ! -- Deallocate arrays if package was active
    if (this%inunit > 0) then
    end if
    !
    ! -- nullify pointers
    this%ibound => null()
    !
    ! -- Scalars
    call mem_deallocate(this%iadvwt)
    call mem_deallocate(this%ats_percel)
    !
    ! -- deallocate parent
    call this%NumericalPackageType%da()

adv_df()

subroutine tspadvmodule::adv_df ( class(tspadvtype)  this, type(tspadvoptionstype), intent(in), optional  adv_options  )

private

Define the ADV package

Parameters

[in]

adv_options

the optional options, for when not constructing from file

Definition at line 89 of file tsp-adv.f90.

    ! -- dummy
    class(TspAdvType) :: this
    type(TspAdvOptionsType), optional, intent(in) :: adv_options !< the optional options, for when not constructing from file
    ! -- local
    character(len=*), parameter :: fmtadv = &
      "(1x,/1x,'ADV -- ADVECTION PACKAGE, VERSION 1, 8/25/2017', &
      &' INPUT READ FROM MEMPATH: ', A, //)"
    !
    ! -- Read or set advection options
    if (.not. present(adv_options)) then
      !
      ! --print a message identifying the advection package.
      write (this%iout, fmtadv) this%input_mempath
      !
      ! --read options from file
      call this%source_options()
    else
      !
      ! --set options from input arg
      this%iadvwt = adv_options%iAdvScheme
    end if

adv_dt()

subroutine tspadvmodule::adv_dt	(	class(tspadvtype)	this,
		real(dp), intent(out)	dtmax,
		character(len=*), intent(inout)	msg,
		real(dp), dimension(:), intent(in)	thetam
	)

Return the largest time step that meets stability constraints

Parameters

	this	this instance
[out]	dtmax	maximum allowable dt subject to stability constraint
[in,out]	msg	package/cell dt constraint message
[in]	thetam	porosity

Definition at line 157 of file tsp-adv.f90.

    ! dummy
    class(TspAdvType) :: this !< this instance
    real(DP), intent(out) :: dtmax !< maximum allowable dt subject to stability constraint
    character(len=*), intent(inout) :: msg !< package/cell dt constraint message
    real(DP), dimension(:), intent(in) :: thetam !< porosity
    ! local
    integer(I4B) :: n
    integer(I4B) :: m
    integer(I4B) :: ipos
    integer(I4B) :: nrmax
    character(len=LINELENGTH) :: cellstr
    real(DP) :: dt
    real(DP) :: flowmax
    real(DP) :: flowsumpos
    real(DP) :: flowsumneg
    real(DP) :: flownm
    real(DP) :: cell_volume
    dtmax = dnodata
    nrmax = 0
    msg = ''
 
    ! If ats_percel not specified by user, then return without making
    ! the courant time step calculation
    if (this%ats_percel == dnodata) then
      return
    end if
 
    ! Calculate time step lengths based on stability constraint for each cell
    ! and store the smallest one
    do n = 1, this%dis%nodes
      if (this%ibound(n) == 0) cycle
      flowsumneg = dzero
      flowsumpos = dzero
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        if (this%dis%con%mask(ipos) == 0) cycle
        m = this%dis%con%ja(ipos)
        if (this%ibound(m) == 0) cycle
        flownm = this%fmi%gwfflowja(ipos)
        if (flownm < dzero) then
          flowsumneg = flowsumneg - flownm
        else
          flowsumpos = flowsumpos + flownm
        end if
      end do
      flowmax = max(flowsumneg, flowsumpos)
      if (flowmax < dprec) cycle
      cell_volume = this%dis%get_cell_volume(n, this%dis%top(n))
      dt = cell_volume * this%fmi%gwfsat(n) * thetam(n) / flowmax
      dt = dt * this%ats_percel
      if (dt < dtmax) then
        dtmax = dt
        nrmax = n
      end if
    end do
    if (nrmax > 0) then
      call this%dis%noder_to_string(nrmax, cellstr)
      write (msg, *) adjustl(trim(this%memoryPath))//'-'//trim(cellstr)
    end if

adv_fc()

subroutine tspadvmodule::adv_fc ( class(tspadvtype)  this, integer(i4b), intent(in)  nodes, class(matrixbasetype), pointer  matrix_sln, integer(i4b), dimension(:), intent(in)  idxglo, real(dp), dimension(:), intent(in), target  cnew, real(dp), dimension(:), intent(inout)  rhs  )

private

Method to calculate coefficients and fill amat and rhs.

Parameters

	this	this instance
[in]	nodes	number of nodes
	matrix_sln	pointer to solution matrix
[in]	idxglo	global indices for matrix
[in]	cnew	new concentration/temperature values
[in,out]	rhs	right-hand side vector

Definition at line 222 of file tsp-adv.f90.

    ! -- modules
    ! -- dummy
    class(TspAdvType) :: this !< this instance
    integer(I4B), intent(in) :: nodes !< number of nodes
    class(MatrixBaseType), pointer :: matrix_sln !< pointer to solution matrix
    integer(I4B), intent(in), dimension(:) :: idxglo !< global indices for matrix
    real(DP), intent(in), dimension(:), target :: cnew !< new concentration/temperature values
    real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector
    ! -- local
    integer(I4B) :: n, m, idiag, ipos
    real(DP) :: qnm !< volumetric flow rate
    type(CoefficientsType) :: coefficients
 
    ! Calculate internal domain fluxes and add to matrix_sln and rhs.
    call this%face_interpolation%set_field(cnew)
    do n = 1, nodes
      if (this%ibound(n) == 0) cycle ! skip inactive nodes
      idiag = this%dis%con%ia(n)
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        m = this%dis%con%ja(ipos)
        if (this%ibound(m) == 0) cycle ! skip inactive nodes
        if (this%dis%con%mask(ipos) == 0) cycle ! skip masked connections
 
        qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
        coefficients = this%face_interpolation%compute(n, m, ipos)
 
        call matrix_sln%add_value_pos(idxglo(idiag), qnm * coefficients%c_n)
        call matrix_sln%add_value_pos(idxglo(ipos), qnm * coefficients%c_m)
        rhs(n) = rhs(n) + qnm * coefficients%rhs
      end do
    end do

allocate_scalars()

subroutine tspadvmodule::allocate_scalars

(

class(tspadvtype)

this

)

Definition at line 318 of file tsp-adv.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    ! -- dummy
    class(TspAdvType) :: this
    ! -- local
    !
    ! -- allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- Allocate
    call mem_allocate(this%iadvwt, 'IADVWT', this%memoryPath)
    call mem_allocate(this%ats_percel, 'ATS_PERCEL', this%memoryPath)
    !
    ! -- Initialize
    this%iadvwt = adv_scheme_upstream
    this%ats_percel = dnodata
    !
    ! -- Advection creates an asymmetric coefficient matrix
    this%iasym = 1

source_options()

subroutine tspadvmodule::source_options

(

class(tspadvtype)

this

)

Definition at line 342 of file tsp-adv.f90.

    ! -- modules
    use kindmodule, only: lgp
    use constantsmodule, only: lenvarname
    use simvariablesmodule, only: errmsg
    use simmodule, only: store_error, store_error_filename
    use memorymanagerextmodule, only: mem_set_value
    ! -- dummy
    class(TspAdvType) :: this
    ! -- locals
    character(len=LENVARNAME), dimension(4) :: scheme = &
      &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD', 'UTVD']
    logical(lgp) :: found_scheme, found_atspercel
    ! -- formats
    character(len=*), parameter :: fmtiadvwt = &
      &"(4x,'ADVECTION WEIGHTING SCHEME HAS BEEN SET TO: ', a)"
 
    ! update defaults with input sourced values
    call mem_set_value(this%iadvwt, 'SCHEME', this%input_mempath, &
                       scheme, found_scheme)
    call mem_set_value(this%ats_percel, 'ATS_PERCEL', this%input_mempath, &
                       found_atspercel)
 
    if (found_scheme) then
      ! should currently be set to index of scheme names
      if (this%iadvwt == 0) then
        write (errmsg, '(a, a)') &
          'Unknown scheme, must be "UPSTREAM", "CENTRAL", "TVD" or "UTVD"'
        call store_error(errmsg)
        call store_error_filename(this%input_fname)
      else
        ! scheme parameters are 0 based
        this%iadvwt = this%iadvwt - 1
      end if
    end if
 
    if (this%iadvwt == adv_scheme_utvd) then
      call dev_feature('UTVD is still under development, install the &
          &nightly build or compile from source with IDEVELOPMODE = 1.')
    end if
 
    if (found_atspercel) then
      if (this%ats_percel == dzero) this%ats_percel = dnodata
    end if
 
    ! log options
    write (this%iout, '(1x,a)') 'PROCESSING ADVECTION OPTIONS'
    if (found_scheme) then
      write (this%iout, fmtiadvwt) trim(scheme(this%iadvwt + 1))
    end if
    if (found_atspercel) then
      write (this%iout, '(4x,a,1pg15.6)') &
        'User-specified fractional cell distance for adaptive time &
        &steps: ', this%ats_percel
    end if
    write (this%iout, '(1x,a)') 'END OF ADVECTION OPTIONS'

Here is the call graph for this function: