utvdschememodule Module Reference

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

Functions/Subroutines
type(utvdschemetype) function	constructor (dis, fmi, gradient)
subroutine	destructor (this)
subroutine	set_field (this, phi)
	Set the scalar field for which interpolation will be computed. More...
type(coefficientstype) function, target	compute (this, n, m, iposnm)
real(dp) function	limiter (this, r)
subroutine	compute_node_distance (this)

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  )

private

Definition at line 106 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
    ! -- 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), pointer :: 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
    !
    ! Determine direction of distance vector from upwind to downwind cell
    ! The cached_node_distance always stores vector from lower-numbered node to higher-numbered node.
    ! Since we need dnm to point from upwind (iup) to downwind (idn), we must adjust the sign:
    ! - If iup > idn: the cached vector points from idn to iup, so we negate it to get iup to idn
    ! - If iup < idn: the cached vector already points from iup to idn, so use it as-is
    if (iup > idn) then
      dnm = -this%cached_node_distance(isympos, :)
    else
      dnm = this%cached_node_distance(isympos, :)
    end if
    !
    ! -- Add low order terms
    coef_up = done
    !
    ! -- Add high order terms
    !
    ! -- Return if straddled cells have same value
    if (abs(this%phi(idn) - this%phi(iup)) < dsame) return
    !
    ! -- Compute cell concentration gradient
    grad_c = this%gradient%get(iup)
    !
    ! Darwish's method to compute virtual node concentration
    c_virtual = this%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.
    min_phi => this%min_max_phi%get_min(iup)
    max_phi => this%min_max_phi%get_max(iup)
 
    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 = (this%phi(iup) - c_virtual) / (this%phi(idn) - this%phi(iup))
    !
    ! -- Compute limiter
    alimiter = this%limiter(smooth)
 
    ! High order term is:
    relative_distance = cl1 / (cl1 + cl2)
    phi_face%rhs = -relative_distance * alimiter * (this%phi(idn) - this%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
 

compute_node_distance()

subroutine utvdschememodule::compute_node_distance ( class(utvdschemetype), target  this )

private

Definition at line 235 of file UTVDScheme.f90.

    ! -- dummy
    class(UTVDSchemeType), target :: this
    ! -- local
    integer(I4B) :: n, m, ipos, isympos
 
    this%cached_node_distance = 0.0_dp
    do n = 1, this%dis%nodes
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        m = this%dis%con%ja(ipos)
        if (m <= n) cycle
 
        isympos = this%dis%con%jas(ipos)
        this%cached_node_distance(isympos, :) = node_distance(this%dis, n, m)
      end do
    end do
 

Here is the call graph for this function:

Here is the caller 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 62 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, intent(in), target :: gradient
 
    interpolation_scheme%dis => dis
    interpolation_scheme%fmi => fmi
 
    interpolation_scheme%gradient => gradient
    interpolation_scheme%min_max_phi = localcellextrematype(dis)
 
    allocate (interpolation_scheme%cached_node_distance(dis%njas, 3))
    call compute_node_distance(interpolation_scheme)
 

Here is the caller graph for this function:

destructor()

subroutine utvdschememodule::destructor ( type(utvdschemetype), intent(inout)  this )

private

Definition at line 82 of file UTVDScheme.f90.

    ! -- dummy
    type(UTVDSchemeType), intent(inout) :: this
 
    deallocate (this%cached_node_distance)

Here is the caller graph for this function:

limiter()

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

private

Definition at line 210 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

set_field()

subroutine utvdschememodule::set_field ( class(utvdschemetype), target  this, real(dp), dimension(:), intent(in), pointer  phi  )

private

This method establishes a pointer to the scalar field data and updates any dependent cached data (gradients and local extrema) to ensure subsequent interpolation computations use the current field values.

Definition at line 96 of file UTVDScheme.f90.

    ! -- dummy
    class(UTVDSchemeType), target :: this
    real(DP), intent(in), dimension(:), pointer :: phi
 
    this%phi => phi
    call this%gradient%set_field(phi)
    call this%min_max_phi%set_field(phi)