utvdschememodule Module Reference

Data Types
interface	utvdschemetype
	Total Variation Diminishing (TVD) interpolation scheme. More...

Functions/Subroutines
type(utvdschemetype) function	constructor (dis, fmi, gradient)
type(coefficientstype) function, target	compute (this, n, m, iposnm, phi)
subroutine	find_local_extrema (this, n, phi, min_phi, max_phi)
real(dp) function	limiter (this, r)

Function/Subroutine Documentation

compute()

type(coefficientstype) function, target utvdschememodule::compute ( class(utvdschemetype), target  this, integer(i4b), intent(in)  n, integer(i4b), intent(in)  m, integer(i4b), intent(in)  iposnm, real(dp), dimension(:), intent(in)  phi  )

private

Definition at line 68 of file UTVDScheme.f90.

    !-- return
    type(CoefficientsType), target :: phi_face ! Output: coefficients for the face between cells n and m
    ! -- dummy
    class(UTVDSchemeType), target :: this
    integer(I4B), intent(in) :: n ! Index of the first cell
    integer(I4B), intent(in) :: m ! Index of the second cell
    integer(I4B), intent(in) :: iposnm ! Position in the connectivity array for the n-m connection
    real(DP), intent(in), dimension(:) :: phi ! Array of scalar values (e.g., concentrations) at all cells
    ! -- local
    integer(I4B) :: iup, idn, isympos ! Indices for upwind, downwind, and symmetric position in connectivity
    real(DP) :: qnm ! Flow rate across the face between n and m
    real(DP), pointer :: coef_up, coef_dn ! Pointers to upwind and downwind coefficients in phi_face
    real(DP), dimension(3) :: grad_c ! gradient at upwind cell
    real(DP), dimension(3) :: dnm ! vector from upwind to downwind cell
    real(DP) :: smooth ! Smoothness indicator for limiter i.e. ratio of gradients
    real(DP) :: alimiter ! Value of the TVD limiter
    real(DP) :: cl1, cl2 ! Connection lengths from upwind and downwind cells to the face
    real(DP) :: relative_distance ! Relative distance factor for high-order term
    real(DP) :: c_virtual ! Virtual node concentration (Darwish method)
    real(DP) :: min_phi, max_phi ! Local minimum and maximum among cell and neighbors
 
    isympos = this%dis%con%jas(iposnm)
    qnm = this%fmi%gwfflowja(iposnm)
    !
    ! -- Find upstream node
    if (qnm > dzero) then
      ! -- positive flow into n means m is upstream
      iup = m
      idn = n
 
      cl1 = this%dis%con%cl2(isympos)
      cl2 = this%dis%con%cl1(isympos)
 
      coef_up => phi_face%c_m
      coef_dn => phi_face%c_n
    else
      iup = n
      idn = m
 
      cl1 = this%dis%con%cl1(isympos)
      cl2 = this%dis%con%cl2(isympos)
 
      coef_up => phi_face%c_n
      coef_dn => phi_face%c_m
    end if
    !
    ! -- Add low order terms
    coef_up = done
    !
    ! -- Add high order terms
    !
    ! -- Return if straddled cells have same value
    if (abs(phi(idn) - phi(iup)) < dsame) return
    !
    ! -- Compute cell concentration gradient
    grad_c = this%gradient%get(iup, phi)
    !
    ! Darwish's method to compute virtual node concentration
    dnm = node_distance(this%dis, iup, idn)
    c_virtual = phi(idn) - 2.0_dp * (dot_product(grad_c, dnm))
    !
    ! Enforce local TVD condition.
    ! This is done by limiting the virtual concentration to the range of
    ! the max and min concentration of the neighbouring cells.
    call find_local_extrema(this, iup, phi, min_phi, max_phi)
 
    if (c_virtual > max_phi) then
      c_virtual = max_phi
    end if
 
    if (c_virtual < max(min_phi, dzero)) then
      c_virtual = max(min_phi, dzero)
    end if
    !
    ! -- Compute smoothness factor
    smooth = (phi(iup) - c_virtual) / (phi(idn) - phi(iup))
    !
    ! -- Compute limiter
    alimiter = this%limiter(smooth)
 
    ! High order term is:
    relative_distance = cl1 / (cl1 + cl2)
    phi_face%rhs = -relative_distance * alimiter * (phi(idn) - phi(iup))
 
    ! Alternative way of writing the high order term by adding it to the
    ! coefficients matrix. The equation to be added is:
    ! high_order = cl1 / (cl1 + cl2) * alimiter * qnm * (phi(idn) - phi(iup))
    ! This is split into two parts:
    ! coef_up = coef_up - relative_distance * alimiter
    ! coef_dn = coef_dn + relative_distance * alimiter
 

Here is the call graph for this function:

constructor()

type(utvdschemetype) function utvdschememodule::constructor ( class(disbasetype), intent(in), pointer  dis, type(tspfmitype), intent(in), pointer  fmi, class(igradienttype), intent(in), allocatable, target  gradient  )

private

Definition at line 53 of file UTVDScheme.f90.

    ! -- return
    type(UTVDSchemeType) :: interpolation_scheme
    ! --dummy
    class(DisBaseType), pointer, intent(in) :: dis
    type(TspFmiType), pointer, intent(in) :: fmi
    class(IGradientType), allocatable, target, intent(in) :: gradient
 
    interpolation_scheme%dis => dis
    interpolation_scheme%fmi => fmi
    interpolation_scheme%gradient => gradient
 

Here is the caller graph for this function:

find_local_extrema()

subroutine utvdschememodule::find_local_extrema ( class(utvdschemetype)  this, integer(i4b), intent(in)  n, real(dp), dimension(:), intent(in)  phi, real(dp), intent(out)  min_phi, real(dp), intent(out)  max_phi  )

private

Definition at line 162 of file UTVDScheme.f90.

    ! -- dummy
    class(UTVDSchemeType) :: this
    integer(I4B), intent(in) :: n
    real(DP), intent(in), dimension(:) :: phi
    real(DP), intent(out) :: min_phi, max_phi
    ! -- local
    integer(I4B) :: ipos, m
 
    min_phi = phi(n)
    max_phi = phi(n)
    do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
      m = this%dis%con%ja(ipos)
      min_phi = min(min_phi, phi(m))
      max_phi = max(max_phi, phi(m))
    end do

Here is the caller graph for this function:

limiter()

real(dp) function utvdschememodule::limiter ( class(utvdschemetype)  this, real(dp)  r  )

private

Definition at line 180 of file UTVDScheme.f90.

    ! -- return
    real(DP) :: theta ! limited slope
    ! -- dummy
    class(UTVDSchemeType) :: this
    real(DP) :: r ! ratio of successive gradients
 
    select case (this%limiter_id)
    case (2) ! van Leer
      theta = max(0.0_dp, min((r + dabs(r)) / (1.0_dp + dabs(r)), 2.0_dp))
    case (3) ! Koren
      theta = max(0.0_dp, min(2.0_dp * r, &
                              1.0_dp / 3.0_dp + 2.0_dp / 3.0_dp * r, 2.0_dp))
    case (4) ! Superbee
      theta = max(0.0_dp, min(2.0_dp * r, 1.0_dp), min(r, 2.0_dp))
    case (5) ! van Albada
      theta = max(0.0_dp, (r * r + r) / (r * r + 1.0_dp))
    case (6) ! Koren modified
      theta = max(0.0_dp, min(4.0_dp * r * r + r, &
                              1.0_dp / 3.0_dp + 2.0_dp / 3.0_dp * r, 2.0_dp))
    case default
      theta = dzero
    end select