mod_venk.t

Go to the documentation of this file.

module mod_venk
  ! Venkatakrishnan limiter
  !
  ! 2019.10.11 coded up by nanami;
  !
  ! see Venkatakrishnan 1993 for this limiter;
 
 
  implicit none
  private
 
  public :: venklimiter
 
contains
 
  subroutine venklimiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC)
    use mod_global_parameters
 
    integer, intent(in)             :: ixi^l, il^l, idims
    double precision, intent(in)    :: dxdim
    double precision, intent(in)    :: w(ixi^s,1:nw)
    double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw) 
    !> local
    double precision                :: wmax(ixi^s,1:nw),wmin(ixi^s,1:nw)
    double precision                :: westp(ixi^s,1:nw),westm(ixi^s,1:nw)
    double precision                :: phi1(ixi^s,1:nw),phi2(ixi^s,1:nw)
    double precision                :: phi3(ixi^s,1:nw),phi4(ixi^s,1:nw)
    double precision                :: phim(ixi^s,1:nw),phip(ixi^s,1:nw),phi(ixi^s,1:nw)
    double precision                :: deltap(ixi^s,1:nw),deltam(ixi^s,1:nw)
    double precision                :: eps2
    double precision, parameter     :: venk_omega = 1.d-12
    double precision, parameter     :: venk_k = 0.3d0
    integer                         :: im^l, imm^l, imp^l
    integer                         :: ilm^l, ilp^l, ilpp^l
 
    im^l=il^l;
    immax^d=ilmax^d+kr(idims,^d);
    ilm^l=il^l-kr(idims,^d);
    imm^l=im^l-kr(idims,^d);
    ilp^l=il^l+kr(idims,^d);
    imp^l=im^l+kr(idims,^d);
    ilpp^l=ilp^l+kr(idims,^d);
 
    eps2 = (venk_k * dxdim) ** 3
 
    wmax(im^s,1:nw_recon) = max(w(imm^s,1:nw_recon), w(im^s,1:nw_recon), w(imp^s,1:nw_recon))
    wmin(im^s,1:nw_recon) = min(w(imm^s,1:nw_recon), w(im^s,1:nw_recon), w(imp^s,1:nw_recon))
    !> use central difference approximation as (eq.5) and take phi = 1
    westp(im^s,1:nw_recon) = w(im^s,1:nw_recon) + (w(imp^s,1:nw_recon) - w(imm^s,1:nw_recon)) * 0.25d0
    westm(im^s,1:nw_recon) = w(im^s,1:nw_recon) - (w(imp^s,1:nw_recon) - w(imm^s,1:nw_recon)) * 0.25d0
 
    !> (eq.30) & (eq.31)
    deltap = 0
    deltam = 0
    phi1 = 0
    phi2 = 0
    phip = 1
    where(westp(im^s,1:nw_recon) .gt. w(im^s,1:nw_recon))
      deltap(im^s,1:nw_recon) = wmax(im^s,1:nw_recon) - w(im^s,1:nw_recon)
      deltam(im^s,1:nw_recon) = westp(im^s,1:nw_recon) - w(im^s,1:nw_recon)
      deltam(im^s,1:nw_recon) = sign(dabs(deltam(im^s,1:nw_recon)) + venk_omega, deltam(im^s,1:nw_recon))
      phi1(im^s,1:nw_recon) = (deltap(im^s,1:nw_recon)**2 + eps2) + 2.d0 * deltap(im^s,1:nw_recon) * deltam(im^s,1:nw_recon)
      phi2(im^s,1:nw_recon) = deltap(im^s,1:nw_recon)**2 + 2.d0 * deltam(im^s,1:nw_recon)**2 + deltap(im^s,1:nw_recon) * deltam(im^s,1:nw_recon) + eps2
      phip(im^s,1:nw_recon) = phi1(im^s,1:nw_recon) / phi2(im^s, 1:nw_recon)
    elsewhere(westp(im^s,1:nw_recon) .lt. w(im^s,1:nw_recon))
      deltap(im^s,1:nw_recon) = wmin(im^s,1:nw_recon) - w(im^s,1:nw_recon)
      deltam(im^s,1:nw_recon) = westp(im^s,1:nw_recon) - w(im^s,1:nw_recon)
      deltam(im^s,1:nw_recon) = sign(dabs(deltam(im^s,1:nw_recon)) + venk_omega, deltam(im^s,1:nw_recon))
      phi1(im^s,1:nw_recon) = (deltap(im^s,1:nw_recon)**2 + eps2) + 2.d0 * deltap(im^s,1:nw_recon) * deltam(im^s,1:nw_recon)
      phi2(im^s,1:nw_recon) = deltap(im^s,1:nw_recon)**2 + 2.d0 * deltam(im^s,1:nw_recon)**2 + deltap(im^s,1:nw_recon) * deltam(im^s,1:nw_recon) + eps2
      phip(im^s,1:nw_recon) = phi1(im^s,1:nw_recon) / phi2(im^s, 1:nw_recon)
    elsewhere
      phip(im^s,1:nw_recon) = one
   endwhere
 
    deltap = 0
    deltam = 0
    phi3 = 0
    phi4 = 0
    phim = 0
    where(westm(im^s,1:nw_recon) .lt. w(im^s,1:nw_recon))
      deltap(im^s,1:nw_recon) = - (wmax(im^s,1:nw_recon) - w(im^s,1:nw_recon))
      deltam(im^s,1:nw_recon) = westm(im^s,1:nw_recon) - w(im^s,1:nw_recon)
      deltam(im^s,1:nw_recon) = sign(dabs(deltam(im^s,1:nw_recon)) + venk_omega, deltam(im^s,1:nw_recon))
      phi3(im^s,1:nw_recon) = (deltap(im^s,1:nw_recon)**2 + eps2) + 2.d0 * deltap(im^s,1:nw_recon) * deltam(im^s,1:nw_recon)
      phi4(im^s,1:nw_recon) = deltap(im^s,1:nw_recon)**2 + 2.d0 * deltam(im^s,1:nw_recon)**2 + deltap(im^s,1:nw_recon) * deltam(im^s,1:nw_recon) + eps2
      phim(im^s,1:nw_recon) = phi3(im^s,1:nw_recon) / phi4(im^s, 1:nw_recon)
    elsewhere(westm(im^s,1:nw_recon) .gt. w(im^s,1:nw_recon))
      deltap(im^s,1:nw_recon) = - (wmin(im^s,1:nw_recon) - w(im^s,1:nw_recon))
      deltam(im^s,1:nw_recon) = westm(im^s,1:nw_recon) - w(im^s,1:nw_recon)
      deltam(im^s,1:nw_recon) = sign(dabs(deltam(im^s,1:nw_recon)) + venk_omega, deltam(im^s,1:nw_recon))
      phi3(im^s,1:nw_recon) = (deltap(im^s,1:nw_recon)**2 + eps2) + 2.d0 * deltap(im^s,1:nw_recon) * deltam(im^s,1:nw_recon)
      phi4(im^s,1:nw_recon) = deltap(im^s,1:nw_recon)**2 + 2.d0 * deltam(im^s,1:nw_recon)**2 + deltap(im^s,1:nw_recon) * deltam(im^s,1:nw_recon) + eps2
      phim(im^s,1:nw_recon) = phi3(im^s,1:nw_recon) / phi4(im^s, 1:nw_recon)
    elsewhere
      phim(im^s,1:nw_recon) = one
    endwhere
    !> (eq.3)
    phi(im^s,1:nw_recon) = min(phim(im^s,1:nw_recon),phip(im^s,1:nw_recon))
    !> (eq.5)
    wlc(il^s,1:nw_recon) = w(il^s,1:nw_recon) + 0.25d0 * (w(ilp^s,1:nw_recon)-w(ilm^s,1:nw_recon)) * phi(il^s,1:nw_recon)
    wrc(il^s,1:nw_recon) = w(ilp^s,1:nw_recon) - 0.25d0 * (w(ilpp^s,1:nw_recon)-w(il^s,1:nw_recon)) * phi(ilp^s,1:nw_recon)
 
  end subroutine venklimiter
 
end module mod_venk