mod_limiter.t

Go to the documentation of this file.

!> Module with slope/flux limiters
module mod_limiter
  use mod_ppm
  use mod_mp5
  use mod_weno
  use mod_venk
 
  implicit none
  public
 
  !> radius of the asymptotic region [0.001, 10], larger means more accurate in smooth
  !> region but more overshooting at discontinuities
  double precision :: cada3_radius
  double precision :: schmid_rad^d
  integer, parameter :: limiter_minmod = 1
  integer, parameter :: limiter_woodward = 2
  integer, parameter :: limiter_mcbeta = 3
  integer, parameter :: limiter_superbee = 4
  integer, parameter :: limiter_vanleer = 5
  integer, parameter :: limiter_albada = 6
  integer, parameter :: limiter_koren = 7
  integer, parameter :: limiter_cada = 8
  integer, parameter :: limiter_cada3 = 9
  integer, parameter :: limiter_schmid = 10
  integer, parameter :: limiter_venk = 11
  ! Special cases
  integer, parameter :: limiter_ppm = 12
  integer, parameter :: limiter_mp5 = 13
  integer, parameter :: limiter_weno3  = 14
  integer, parameter :: limiter_wenoyc3  = 15
  integer, parameter :: limiter_weno5 = 16
  integer, parameter :: limiter_weno5nm = 17
  integer, parameter :: limiter_wenoz5  = 18
  integer, parameter :: limiter_wenoz5nm = 19
  integer, parameter :: limiter_wenozp5  = 20
  integer, parameter :: limiter_wenozp5nm = 21
  integer, parameter :: limiter_weno5cu6 = 22
  integer, parameter :: limiter_teno5ad = 23
  integer, parameter :: limiter_weno7 = 24
  integer, parameter :: limiter_mpweno7 = 25
 
contains
 
  integer function limiter_type(namelim)
    use mod_comm_lib, only: mpistop
    character(len=*), intent(in) :: namelim
 
    select case (namelim)
    case ('minmod')
       limiter_type = limiter_minmod
    case ('woodward')
       limiter_type = limiter_woodward
    case ('mcbeta')
       limiter_type = limiter_mcbeta
    case ('superbee')
       limiter_type = limiter_superbee
    case ('vanleer')
       limiter_type = limiter_vanleer
    case ('albada')
       limiter_type = limiter_albada
    case ('koren')
       limiter_type = limiter_koren
    case ('cada')
       limiter_type = limiter_cada
    case ('cada3')
       limiter_type = limiter_cada3
    case ('schmid1')
       limiter_type = limiter_schmid
    case ('schmid2')
       limiter_type = limiter_schmid
    case('venk')
       limiter_type = limiter_venk
    case ('ppm')
       limiter_type = limiter_ppm
    case ('mp5')
       limiter_type = limiter_mp5
    case ('weno3')
       limiter_type = limiter_weno3
    case ('wenoyc3')
       limiter_type = limiter_wenoyc3
    case ('weno5')
       limiter_type = limiter_weno5
    case ('weno5nm')
       limiter_type = limiter_weno5nm
    case ('wenoz5')
       limiter_type = limiter_wenoz5
    case ('wenoz5nm')
       limiter_type = limiter_wenoz5nm
    case ('wenozp5')
       limiter_type = limiter_wenozp5
    case ('wenozp5nm')
       limiter_type = limiter_wenozp5nm
    case ('weno5cu6')
       limiter_type = limiter_weno5cu6
    case ('teno5ad')
       limiter_type = limiter_teno5ad
    case ('weno7')
       limiter_type = limiter_weno7
    case ('mpweno7')
       limiter_type = limiter_mpweno7
 
    case default
       limiter_type = -1
       write(*,*) 'Unknown limiter: ', namelim
       call mpistop("No such limiter")
    end select
  end function limiter_type
 
  pure logical function limiter_symmetric(typelim)
    integer, intent(in) :: typelim
 
    select case (typelim)
    case (limiter_koren, limiter_cada, limiter_cada3)
       limiter_symmetric = .false.
    case default
       limiter_symmetric = .true.
    end select
  end function limiter_symmetric
 
  !> Limit the centered dwC differences within ixC for iw in direction idim.
  !> The limiter is chosen according to typelim.
  !>
  !> Note that this subroutine is called from upwindLR (hence from methods
  !> like tvdlf, hancock, hll(c) etc) or directly from tvd.t,
  !> but also from the gradientL and divvectorS subroutines in geometry.t
  !> Accordingly, the typelim here corresponds to one of limiter
  !> or one of gradient_limiter.
  subroutine dwlimiter2(dwC,ixI^L,ixC^L,idims,typelim,ldw,rdw,a2max)
 
    use mod_global_parameters
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: ixI^L, ixC^L, idims
    double precision, intent(in) :: dwC(ixI^S)
    integer, intent(in) :: typelim
    !> Result using left-limiter (same as right for symmetric)
    double precision, intent(out), optional :: ldw(ixI^S)
    !> Result using right-limiter (same as left for symmetric)
    double precision, intent(out), optional :: rdw(ixI^S)
    double precision, intent(in), optional :: a2max
 
    double precision :: tmp(ixI^S), tmp2(ixI^S)
    double precision, parameter :: qsmall=1.d-12, qsmall2=2.d-12
    double precision, parameter :: eps = sqrt(epsilon(1.0d0))
 
    ! mcbeta limiter parameter value
    double precision, parameter :: c_mcbeta=1.4d0
    ! cada limiter parameter values
    double precision, parameter :: cadalfa=0.5d0, cadbeta=2.0d0, cadgamma=1.6d0
    ! full third order cada limiter
    double precision :: rdelinv
    double precision :: ldwA(ixI^S),ldwB(ixI^S),tmpeta(ixI^S)
    double precision, parameter :: cadepsilon=1.d-14, invcadepsilon=1.d14
    integer :: ixO^L, hxO^L, ix^D, hx^D
    !-----------------------------------------------------------------------------
 
    ! Contract indices in idim for output.
    ixomin^d=ixcmin^d+kr(idims,^d); ixomax^d=ixcmax^d;
    hxo^l=ixo^l-kr(idims,^d);
    hx^d=kr(idims,^d);
 
    ! About the notation: the conventional argument theta (the ratio of slopes)
    ! would be given by dwC(ixO^S)/dwC(hxO^S). However, in the end one
    ! multiplies phi(theta) by dwC(hxO^S), which is incorporated in the
    ! equations below. The minmod limiter can for example be written as:
    ! A:
    ! max(0.0d0, min(1.0d0, dwC(ixO^S)/dwC(hxO^S))) * dwC(hxO^S)
    ! B:
    ! tmp(ixO^S)*max(0.0d0,min(abs(dwC(ixO^S)),tmp(ixO^S)*dwC(hxO^S)))
    ! where tmp(ixO^S)=sign(1.0d0,dwC(ixO^S))
 
    select case (typelim)
    case (limiter_minmod)
       ! Minmod limiter eq(3.51e) and (eq.3.38e) with omega=1
       tmp(ixo^s)=sign(one,dwc(ixo^s))
       tmp(ixo^s)=tmp(ixo^s)* &
            max(zero,min(abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s)))
       if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
       if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
    case (limiter_woodward)
       ! Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
       tmp(ixo^s)=sign(one,dwc(ixo^s))
       tmp(ixo^s)=2*tmp(ixo^s)* &
            max(zero,min(abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s),&
            tmp(ixo^s)*quarter*(dwc(hxo^s)+dwc(ixo^s))))
       if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
       if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
    case (limiter_mcbeta)
       ! Woodward and Collela limiter, with factor beta
       tmp(ixo^s)=sign(one,dwc(ixo^s))
       tmp(ixo^s)=tmp(ixo^s)* &
            max(zero,min(c_mcbeta*abs(dwc(ixo^s)),c_mcbeta*tmp(ixo^s)*dwc(hxo^s),&
            tmp(ixo^s)*half*(dwc(hxo^s)+dwc(ixo^s))))
       if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
       if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
    case (limiter_superbee)
       ! Roes superbee limiter (eq.3.51i)
       tmp(ixo^s)=sign(one,dwc(ixo^s))
       tmp(ixo^s)=tmp(ixo^s)* &
            max(zero,min(2*abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s)),&
            min(abs(dwc(ixo^s)),2*tmp(ixo^s)*dwc(hxo^s)))
       if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
       if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
    case (limiter_vanleer)
       ! van Leer limiter (eq 3.51f), but a missing delta2=1.D-12 is added
       tmp(ixo^s)=2*max(dwc(hxo^s)*dwc(ixo^s),zero) &
            /(dwc(ixo^s)+dwc(hxo^s)+qsmall)
       if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
       if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
    case (limiter_albada)
       ! Albada limiter (eq.3.51g) with delta2=1D.-12
       tmp(ixo^s)=(dwc(hxo^s)*(dwc(ixo^s)**2+qsmall)&
            +dwc(ixo^s)*(dwc(hxo^s)**2+qsmall))&
            /(dwc(ixo^s)**2+dwc(hxo^s)**2+qsmall2)
       if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
       if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
    case (limiter_koren)
       tmp(ixo^s)=sign(one,dwc(ixo^s))
       tmp2(ixo^s)=min(2*abs(dwc(ixo^s)),2*tmp(ixo^s)*dwc(hxo^s))
       if (present(ldw)) then
          ldw(ixo^s)=tmp(ixo^s)* &
               max(zero,min(tmp2(ixo^s),&
               (dwc(hxo^s)*tmp(ixo^s)+2*abs(dwc(ixo^s)))*third))
       end if
       if (present(rdw)) then
          rdw(ixo^s)=tmp(ixo^s)* &
                max(zero,min(tmp2(ixo^s),&
               (2*dwc(hxo^s)*tmp(ixo^s)+abs(dwc(ixo^s)))*third))
       end if
    case (limiter_cada)
       ! This limiter has been rewritten in the usual form, and uses a division
       ! of the gradients.
       if (present(ldw)) then
          ! Cada Left variant
          ! Compute theta, but avoid division by zero
          tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s) + sign(eps, dwc(ixo^s)))
          tmp2(ixo^s)=(2+tmp(ixo^s))*third
          ldw(ixo^s)= max(zero,min(tmp2(ixo^s), &
               max(-cadalfa*tmp(ixo^s), &
               min(cadbeta*tmp(ixo^s), tmp2(ixo^s), &
               cadgamma)))) * dwc(ixo^s)
       end if
 
       if (present(rdw)) then
          ! Cada Right variant
          tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s) + sign(eps, dwc(hxo^s)))
          tmp2(ixo^s)=(2+tmp(ixo^s))*third
          rdw(ixo^s)= max(zero,min(tmp2(ixo^s), &
               max(-cadalfa*tmp(ixo^s), &
               min(cadbeta*tmp(ixo^s), tmp2(ixo^s), &
               cadgamma)))) * dwc(hxo^s)
       end if
    case (limiter_cada3)
       rdelinv=one/(cada3_radius*dxlevel(idims))**2
      {!DEC$ VECTOR ALWAYS
       do ix^db=ixomin^db,ixomax^db\}
         tmpeta(ix^d)=(dwc(ix^d)**2+dwc(ix^d-hx^d)**2)*rdelinv
         if (present(ldw)) then
            tmp(ix^d)=dwc(ix^d-hx^d)/(dwc(ix^d)+sign(eps,dwc(ix^d)))
            ldwa(ix^d)=(two+tmp(ix^d))*third
            if(tmpeta(ix^d)<=one-cadepsilon) then
              ldw(ix^d)=ldwa(ix^d)
            else if(tmpeta(ix^d)>=one+cadepsilon) then
              ldw(ix^d)=max(zero,min(ldwa(ix^d),max(-cadalfa*tmp(ix^d),&
                min(cadbeta*tmp(ix^d),ldwa(ix^d),cadgamma))))
            else
              tmp2(ix^d)=(tmpeta(ix^d)-one)*invcadepsilon
              ldw(ix^d)=half*((one-tmp2(ix^d))*ldwa(ix^d)+(one+tmp2(ix^d))*&
                max(zero,min(ldwa(ix^d),max(-cadalfa*tmp(ix^d),&
                  min(cadbeta*tmp(ix^d),ldwa(ix^d),cadgamma)))))
            end if
            ldw(ix^d)=ldw(ix^d)*dwc(ix^d)
         end if
         if (present(rdw)) then
            tmp(ix^d)=dwc(ix^d)/(dwc(ix^d-hx^d)+sign(eps,dwc(ix^d-hx^d)))
            ldwa(ix^d)=(two+tmp(ix^d))*third
            if(tmpeta(ix^d)<=one-cadepsilon) then
              rdw(ix^d)=ldwa(ix^d)
            else if(tmpeta(ix^d)>=one+cadepsilon) then
              rdw(ix^d)=max(zero,min(ldwa(ix^d),max(-cadalfa*tmp(ix^d),&
                min(cadbeta*tmp(ix^d),ldwa(ix^d),cadgamma))))
            else
              tmp2(ix^d)=(tmpeta(ix^d)-one)*invcadepsilon
              rdw(ix^d)=half*((one-tmp2(ix^d))*ldwa(ix^d)+(one+tmp2(ix^d))*&
                max(zero,min(ldwa(ix^d),max(-cadalfa*tmp(ix^d),&
                  min(cadbeta*tmp(ix^d),ldwa(ix^d),cadgamma)))))
            end if
            rdw(ix^d)=rdw(ix^d)*dwc(ix^d-hx^d)
         end if
      {end do\}
    case(limiter_schmid)
      tmpeta(ixo^s)=(sqrt(0.4d0*(dwc(ixo^s)**2+dwc(hxo^s)**2)))&
        /((a2max+cadepsilon)*dxlevel(idims)**2)
      if(present(ldw)) then
        tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s)+sign(eps,dwc(ixo^s)))
        ldwa(ixo^s)=(two+tmp(ixo^s))*third
        where(tmpeta(ixo^s)<=one-cadepsilon)
          ldw(ixo^s)=ldwa(ixo^s)
        else where(tmpeta(ixo^s)>=one+cadepsilon)
          ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
             min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
          ldw(ixo^s)=ldwb(ixo^s)
        else where
          ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
             min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
          tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
          ldw(ixo^s)=half*((one-tmp2(ixo^s))*ldwa(ixo^s)&
            +(one+tmp2(ixo^s))*ldwb(ixo^s))
        end where
        ldw(ixo^s)=ldw(ixo^s)*dwc(ixo^s)
      end if
      if(present(rdw)) then
        tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s)+sign(eps,dwc(hxo^s)))
        ldwa(ixo^s)=(two+tmp(ixo^s))*third
        where(tmpeta(ixo^s)<=one-cadepsilon)
          rdw(ixo^s)=ldwa(ixo^s)
        else where(tmpeta(ixo^s)>=one+cadepsilon)
          ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
             min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
          rdw(ixo^s)=ldwb(ixo^s)
        else where
          ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
            min(cadbeta*tmp(ixo^s), ldwa(ixo^s), 1.5d0))))
          tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
          rdw(ixo^s)=half*((one-tmp2(ixo^s))*ldwa(ixo^s)&
            +(one+tmp2(ixo^s))*ldwb(ixo^s))
        end where
        rdw(ixo^s)=rdw(ixo^s)*dwc(hxo^s)
       end if
    case default
       call mpistop("Error in dwLimiter: unknown limiter")
    end select
 
  end subroutine dwlimiter2
 
end module mod_limiter