CentralDifferenceScheme.f90

Go to the documentation of this file.

module centraldifferenceschememodule
  use kindmodule, only: dp, i4b
  use constantsmodule, only: dhalf, done
  use interpolationschemeinterfacemodule, only: interpolationschemeinterface, &
                                                coefficientstype
  use basedismodule, only: disbasetype
  use tspfmimodule, only: tspfmitype
 
  implicit none
  private
 
  public :: centraldifferenceschemetype
 
  type, extends(interpolationschemeinterface) :: centraldifferenceschemetype
    private
    class(disbasetype), pointer :: dis
    type(tspfmitype), pointer :: fmi
  contains
    procedure :: compute
  end type centraldifferenceschemetype
 
  interface centraldifferenceschemetype
    module procedure constructor
  end interface centraldifferenceschemetype
 
contains
  function constructor(dis, fmi) result(interpolation_scheme)
    ! -- return
    type(centraldifferenceschemetype) :: interpolation_scheme
    ! --dummy
    class(disbasetype), pointer, intent(in) :: dis
    type(tspfmitype), pointer, intent(in) :: fmi
 
    interpolation_scheme%dis => dis
    interpolation_scheme%fmi => fmi
 
  end function constructor
 
  function compute(this, n, m, iposnm, phi) result(phi_face)
    !-- return
    type(coefficientstype) :: phi_face
    ! -- dummy
    class(centraldifferenceschemetype), target :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: m
    integer(I4B), intent(in) :: iposnm
    real(dp), intent(in), dimension(:) :: phi
    ! -- local
    real(dp) :: lnm, lmn, omega
 
    ! -- calculate weight based on distances between nodes and the shared
    !    face of the connection
    if (this%dis%con%ihc(this%dis%con%jas(iposnm)) == 0) then
      ! -- vertical connection; assume cell is fully saturated
      lnm = dhalf * (this%dis%top(n) - this%dis%bot(n))
      lmn = dhalf * (this%dis%top(m) - this%dis%bot(m))
    else
      ! -- horizontal connection
      lnm = this%dis%con%cl1(this%dis%con%jas(iposnm))
      lmn = this%dis%con%cl2(this%dis%con%jas(iposnm))
    end if
 
    omega = lmn / (lnm + lmn)
 
    phi_face%c_n = omega
    phi_face%c_m = done - omega
 
  end function compute
end module centraldifferenceschememodule