TVDScheme.f90

Go to the documentation of this file.

module tvdschememodule
  use kindmodule, only: dp, i4b
  use constantsmodule, only: done, dzero, dprec, dhalf, dtwo
  use interpolationschemeinterfacemodule, only: interpolationschemeinterface, &
                                                coefficientstype
  use basedismodule, only: disbasetype
  use tspfmimodule, only: tspfmitype
 
  implicit none
  private
 
  public :: tvdschemetype
 
  type, extends(interpolationschemeinterface) :: tvdschemetype
    private
    class(disbasetype), pointer :: dis
    type(tspfmitype), pointer :: fmi
    integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound
  contains
    procedure :: compute
  end type tvdschemetype
 
  interface tvdschemetype
    module procedure constructor
  end interface tvdschemetype
 
contains
  function constructor(dis, fmi, ibound) result(interpolation_scheme)
    ! -- return
    type(tvdschemetype) :: interpolation_scheme
    ! --dummy
    class(disbasetype), pointer, intent(in) :: dis
    type(tspfmitype), pointer, intent(in) :: fmi
    integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibound
 
    interpolation_scheme%dis => dis
    interpolation_scheme%fmi => fmi
    interpolation_scheme%ibound => ibound
 
  end function constructor
 
  function compute(this, n, m, iposnm, phi) result(phi_face)
    !-- return
    type(coefficientstype), target :: phi_face
    ! -- dummy
    class(tvdschemetype), target :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: m
    integer(I4B), intent(in) :: iposnm
    real(dp), intent(in), dimension(:) :: phi
    ! -- local
    integer(I4B) :: ipos, isympos, iup, idn, i2up, j
    real(dp) :: qnm, qmax, qupj, elupdn, elup2up
    real(dp) :: smooth, cdiff, alimiter
    real(dp), pointer :: coef_up, coef_dn
    !
    ! -- Find upstream node
    isympos = this%dis%con%jas(iposnm)
    qnm = this%fmi%gwfflowja(iposnm)
    if (qnm > dzero) then
      ! -- positive flow into n means m is upstream
      iup = m
      idn = n
 
      coef_up => phi_face%c_m
      coef_dn => phi_face%c_n
    else
      iup = n
      idn = m
 
      coef_up => phi_face%c_n
      coef_dn => phi_face%c_m
    end if
    elupdn = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos)
    !
    ! -- Add low order terms
    coef_up = done
    !
    ! -- Add high order terms
    !
    ! -- Find second node upstream to iup
    i2up = 0
    qmax = dzero
    do ipos = this%dis%con%ia(iup) + 1, this%dis%con%ia(iup + 1) - 1
      j = this%dis%con%ja(ipos)
      if (this%ibound(j) == 0) cycle
      qupj = this%fmi%gwfflowja(ipos)
      isympos = this%dis%con%jas(ipos)
      if (qupj > qmax) then
        qmax = qupj
        i2up = j
        elup2up = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos)
      end if
    end do
    !
    ! -- Calculate flux limiting term
    if (i2up > 0) then
      smooth = dzero
      cdiff = abs(phi(idn) - phi(iup))
      if (cdiff > dprec) then
        smooth = (phi(iup) - phi(i2up)) / elup2up * &
                 elupdn / (phi(idn) - phi(iup))
      end if
      if (smooth > dzero) then
        alimiter = dtwo * smooth / (done + smooth)
        phi_face%rhs = -dhalf * alimiter * (phi(idn) - phi(iup))
      end if
    end if
 
  end function compute
 
end module tvdschememodule