mod_tvd.t

Go to the documentation of this file.

!> Subroutines for TVD-MUSCL schemes
module mod_tvd
 
  implicit none
  private
 
  public :: tvdlimit
  public :: tvdlimit2
  public :: entropyfix
 
contains
 
  subroutine tvdlimit(method,qdt,ixI^L,ixO^L,idim^LIM,s,qt,snew,fC,dxs,x)
    use mod_global_parameters
 
    integer, intent(in) :: method
    double precision, intent(in) :: qdt, qt, dxs(ndim)
    integer, intent(in) :: ixi^l, ixo^l, idim^lim
    double precision, dimension(ixI^S,nw) :: w, wnew
    type(state) :: s, snew
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision :: fc(ixi^s,1:nwflux,1:ndim)
 
    integer :: idims, ixic^l, jxic^l
    double precision, dimension(ixI^S,nw) :: wr, wl
 
    associate(w=>s%w,wnew=>snew%w)
    do idims= idim^lim
       ixicmax^d=ixomax^d+kr(idims,^d); ixicmin^d=ixomin^d-2*kr(idims,^d);
       wl(ixic^s,1:nw)=w(ixic^s,1:nw)
       jxic^l=ixic^l+kr(idims,^d);
       wr(ixic^s,1:nw)=w(jxic^s,1:nw)
       call tvdlimit2(method,qdt,ixi^l,ixic^l,ixo^l,idims,wl,wr,wnew,x,fc,dxs)
    end do
    end associate
  end subroutine tvdlimit
 
  subroutine tvdlimit2(method,qdt,ixI^L,ixIC^L,ixO^L,idims,wL,wR,wnew,x,fC,dxs)
 
    ! Limit the flow variables in wnew according to typetvd. 
    ! wroeC is based on wL and wR.
    ! If method=fs_tvd an extra adtdx**2*jumpC is added to phiC for 2nd order
    ! accuracy in time.
 
    use mod_global_parameters
    use mod_physics_roe
 
    integer, intent(in) :: method
    double precision, intent(in) :: qdt, dxs(ndim)
    integer, intent(in) :: ixi^l, ixic^l, ixo^l, idims
    double precision, dimension(ixG^T,nw) :: wl, wr
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision :: wnew(ixi^s,1:nw)
    double precision :: fc(ixi^s,1:nwflux,1:ndim)
 
    double precision:: workroe(ixg^t,1:nworkroe)
    double precision, dimension(ixG^T,nw) :: wroec
    double precision, dimension(ixG^T) :: phic, rphic, jumpc, adtdxc, smallac
    double precision :: dxinv(1:ndim)
    integer :: hxo^l, ixc^l, jxc^l, jxic^l, iw, il
    !-----------------------------------------------------------------------------
 
    hxo^l=ixo^l-kr(idims,^d);
    ixcmax^d=ixomax^d; ixcmin^d=hxomin^d; 
 
    jxc^l=ixc^l+kr(idims,^d);
    jxic^l=ixic^l+kr(idims,^d);
 
    call phys_average(wl,wr,x,ixic^l,idims,wroec,workroe)
 
    dxinv=qdt/dxs
 
    ! A loop on characteristic variables to calculate the dissipative flux phiC.
    do il=1,nwflux
       !Calculate the jump in the il-th characteristic variable: L(wroe)*dw
       call phys_get_eigenjump(wl,wr,wroec,x,ixic^l,il,idims,smallac,adtdxc,jumpc,workroe)
 
       ! Normalize the eigenvalue "a" (and its limit "smalla" if needed):
       if (slab_uniform) then
          adtdxc(ixic^s)=adtdxc(ixic^s)*dxinv(idims)
          if (typeentropy(il)=='harten' .or. typeentropy(il)=='powell')&
            smallac(ixic^s)=smallac(ixic^s)*dxinv(idims)
       else
          adtdxc(ixic^s)=adtdxc(ixic^s)*qdt*block%surfaceC(ixic^s,idims)*&
             2.0d0/(block%dvolume(ixic^s)+block%dvolume(jxic^s))
          if (typeentropy(il)=='harten' .or. typeentropy(il)=='powell')&
            smallac(ixic^s)=smallac(ixic^s)*qdt*block%surfaceC(ixic^s,idims)*&
             2.0d0/(block%dvolume(ixic^s)+block%dvolume(jxic^s))
       endif
 
       ! Calculate the flux limiter function phi
       call getphi(method,jumpc,adtdxc,smallac,ixi^l,ixic^l,ixc^l,il,idims,phic)
 
       !Add R(iw,il)*phiC(il) to each variable iw in wnew
       do iw=1,nwflux
          call phys_rtimes(phic,wroec,ixc^l,iw,il,idims,rphic,workroe)
 
          if (slab_uniform) then
             rphic(ixc^s)=rphic(ixc^s)*half
             fc(ixc^s,iw,idims)=fc(ixc^s,iw,idims)+rphic(ixc^s)
             wnew(ixo^s,iw)=wnew(ixo^s,iw)+rphic(ixo^s)-rphic(hxo^s)
          else
             rphic(ixc^s)=rphic(ixc^s)*quarter* &
                   (block%dvolume(ixc^s)+block%dvolume(jxc^s))
             fc(ixc^s,iw,idims)=fc(ixc^s,iw,idims)+rphic(ixc^s)
             wnew(ixo^s,iw)=wnew(ixo^s,iw)+(rphic(ixo^s)-rphic(hxo^s)) &
                                            /block%dvolume(ixo^s)
          endif
       end do  !iw
    end do     !il
 
  end subroutine tvdlimit2
 
  subroutine getphi(method,jumpC,adtdxC,smallaC,ixI^L,ixIC^L,ixC^L,il,idims,phiC)
 
    ! Calculate the dissipative flux from jumpC=L*dw and adtdx=eigenvalue*dt/dx.
    ! Add Lax-Wendroff type correction if method=fs_tvd.
    ! Limit according to method and typetvd.
    use mod_limiter
    use mod_global_parameters
 
    integer, intent(in) :: method
    integer, intent(in) :: ixI^L, ixIC^L, ixC^L, il, idims
    double precision, dimension(ixG^T) :: jumpC, adtdxC, smallaC, phiC
 
    double precision, dimension(ixG^T) :: ljumpC, tmp
    integer :: jxC^L, ix^L, hx^L, typelimiter
    !-----------------------------------------------------------------------------
 
    typelimiter=type_limiter(block%level)
    if(method==fs_tvdmu)then
       ! In the MUSCL scheme phi=|a|*jump, apply entropy fix to it
       if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
          phic(ixc^s)=abs(adtdxc(ixc^s))*jumpc(ixc^s)
       else
          where(abs(adtdxc(ixc^s))>=smallac(ixc^s))
             phic(ixc^s)=abs(adtdxc(ixc^s))*jumpc(ixc^s)
          elsewhere
             phic(ixc^s)=half*(smallac(ixc^s)+adtdxc(ixc^s)**2/smallac(ixc^s))&
                  *jumpc(ixc^s)
          endwhere
       endif
       ! That's all for the MUSCL scheme
       return
    endif
 
    if(method==fs_tvd)then
       !Entropy fix to |a|-a**2
       if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
          phic(ixic^s)=abs(adtdxc(ixic^s))-adtdxc(ixic^s)**2
       else
          where(abs(adtdxc(ixic^s))>=smallac(ixic^s))
             phic(ixic^s)=abs(adtdxc(ixic^s))-adtdxc(ixic^s)**2
          elsewhere
             phic(ixic^s)=half*smallac(ixic^s)+&
                  (half/smallac(ixic^s)-one)*adtdxc(ixic^s)**2
          endwhere
       endif
    else
       !Entropy fix to |a|
       if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
          phic(ixic^s)=abs(adtdxc(ixic^s))
       else
          where(abs(adtdxc(ixic^s))>=smallac(ixic^s))
             phic(ixic^s)=abs(adtdxc(ixic^s))
          elsewhere
             phic(ixic^s)=half*smallac(ixic^s)+&
                  half/smallac(ixic^s)*adtdxc(ixic^s)**2
          endwhere
       endif
    endif
 
    jxc^l=ixc^l+kr(idims,^d);
    hxmin^d=ixicmin^d; hxmax^d=ixicmax^d-kr(idims,^d);
    ix^l=hx^l+kr(idims,^d);
 
    if (.not. limiter_symmetric(typelimiter)) then
       call mpistop("TVD only supports symmetric limiters")
    end if
 
    select case(typetvd)
    case('roe')
       call dwlimiter2(jumpc,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
       where(adtdxc(ixc^s)<=0)
          phic(ixc^s)=phic(ixc^s)*(jumpc(ixc^s)-ljumpc(jxc^s))
       elsewhere
          phic(ixc^s)=phic(ixc^s)*(jumpc(ixc^s)-ljumpc(ixc^s))
       end where
       !extra (a*lambda)**2*delta
       if(method==fs_tvd)phic(ixc^s)=phic(ixc^s)+adtdxc(ixc^s)**2*jumpc(ixc^s)
    case('sweby')
       !Sweby eqs.4.11-4.15, but no 0.5 ?!
       phic(ixic^s)=phic(ixic^s)*jumpc(ixic^s)
       call dwlimiter2(phic,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
       where(adtdxc(ixc^s)<=0)
          phic(ixc^s)=phic(ixc^s)-ljumpc(jxc^s)
       elsewhere
          phic(ixc^s)=phic(ixc^s)-ljumpc(ixc^s)
       end where
       !extra (a*lambda)**2*delta
       if(method==fs_tvd)phic(ixc^s)=phic(ixc^s)+adtdxc(ixc^s)**2*jumpc(ixc^s)
    case('yee')
       !eq.3.51 with correction
       call dwlimiter2(jumpc,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
 
       !Use phiC as 0.5*(|nu|-nu**2) eq.3.45e for tvd otherwise 0.5*|nu|
       phic(ixc^s)=half*phic(ixc^s)
       !gamma*lambda eq.3.51d, use tmp to store agdtdxC
       where(abs(jumpc(ixc^s))>smalldouble)
          tmp(ixc^s)=adtdxc(ixc^s)+phic(ixc^s)*&
               (ljumpc(jxc^s)-ljumpc(ixc^s))/jumpc(ixc^s)
       elsewhere
          tmp(ixc^s)=adtdxc(ixc^s)
       end where
 
       !eq.3.51a
       if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
          phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
               abs(tmp(ixc^s))*jumpc(ixc^s)
       else
          where(abs(tmp(ixc^s))>=smallac(ixc^s))
             phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
                  abs(tmp(ixc^s))*jumpc(ixc^s)
          elsewhere
             phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
                  (half*smallac(ixc^s)+half/smallac(ixc^s)*&
                  tmp(ixc^s)**2)*jumpc(ixc^s)
          endwhere
       endif
    case('harten')
       !See Ryu, section 2.3
       !Use phiC as 0.5*(|nu|-nu**2)*jumpC eq.3.45b,e
       phic(ixic^s)=half*phic(ixic^s)*jumpc(ixic^s)
       call dwlimiter2(phic,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
 
       !gamma*lambda eq.3.45d, use tmp as agdtdxC
       where(abs(jumpc(ixc^s))>smalldouble)
          tmp(ixc^s)=adtdxc(ixc^s)+&
               (ljumpc(jxc^s)-ljumpc(ixc^s))/jumpc(ixc^s)
       elsewhere
          tmp(ixc^s)=adtdxc(ixc^s)
       end where
       !eq.3.45a with correction
       if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
          phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
               abs(tmp(ixc^s))
       else
          where(abs(tmp(ixc^s))>=smallac(ixc^s))
             phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
                  abs(tmp(ixc^s))
          elsewhere
             phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
                  (half*smallac(ixc^s)+half/smallac(ixc^s)*tmp(ixc^s)**2)
          endwhere
       endif
       !extra -(a*lambda)**2*delta
    case default
       call mpistop("Error in TVDLimit: Unknown TVD type")
    end select
 
  end subroutine getphi
 
  subroutine entropyfix(ix^L,il,aL,aR,a,smalla)
 
    ! Apply entropyfix based on typeentropy(il),aL,aR, and a
    ! Calculate "smalla" (Harten,Powell) or modify "a" (ratio)
 
    use mod_global_parameters
 
    integer, intent(in) :: ix^l, il
    double precision, dimension(ixG^T) :: al, ar, a, smalla
    !-----------------------------------------------------------------------------
 
    select case(typeentropy(il))
    case('harten')
       smalla(ix^s)=max(zero,a(ix^s)-al(ix^s),ar(ix^s)-a(ix^s))
    case('powell')
       smalla(ix^s)=max(zero,two*(ar(ix^s)-al(ix^s)))
       !!case('ratio')
       !!   where(aL(ix^S)<zero .and. aR(ix^S)>zero)&
       !!      a(ix^S)=a(ix^S)-2*aR(ix^S)*aL(ix^S)/(aR(ix^S)-aL(ix^S))
    case('yee')
       ! This has been done in geteigenjump already
    case('nul')
       ! No entropyfix is applied
    case default
       call mpistop("No such type of entropy fix")
    end select
 
  end subroutine entropyfix
 
end module mod_tvd