mod_rhd_hllc.t

Go to the documentation of this file.

!> Hydrodynamics HLLC module
module mod_rhd_hllc
  use mod_rhd_phys
 
  implicit none
  private
 
  public :: rhd_hllc_init
 
contains
 
  subroutine rhd_hllc_init()
    use mod_physics_hllc
 
    phys_diffuse_hllcd => rhd_diffuse_hllcd
    phys_get_lcd => rhd_get_lcd
    phys_get_wcd => rhd_get_wcd
  end subroutine rhd_hllc_init
 
  subroutine rhd_diffuse_hllcd(ixI^L,ixO^L,idim,wLC,wRC,fLC,fRC,patchf)
 
    ! when method is hllcd or hllcd1 then:
 
    ! this subroutine is to enforce regions where we AVOID HLLC
    ! and use TVDLF instead: this is achieved by setting patchf to 4 in
    ! certain regions. An additional input parameter is nxdiffusehllc
    ! which sets the size of the fallback region.
 
    use mod_global_parameters
 
    integer, intent(in)                                      :: ixi^l,ixo^l,idim
    double precision, dimension(ixI^S,1:nw), intent(in)      :: wrc,wlc
    double precision, dimension(ixI^S,1:nwflux),intent(in) :: flc, frc
    integer         , dimension(ixI^S), intent(inout)        :: patchf
 
    integer                                           :: ixoo^d,txoo^l
 
    ! In a user-controlled region around any point with flux sign change between
    ! left and right, ensure fallback to TVDLF
    {do ixoo^d= ixo^lim^d\}
    {
    txoomin^d= max(ixoo^d - nxdiffusehllc*kr(idim,^d), ixomin^d);
    txoomax^d= min(ixoo^d + nxdiffusehllc*kr(idim,^d), ixomax^d);
    \}
    if(abs(patchf(ixoo^d)) == 1 .or. abs(patchf(ixoo^d)) == 4)Then
       if(any(frc(ixoo^d,1:nwflux)*flc(ixoo^d,1:nwflux)<-smalldouble))Then
          where(abs(patchf(txoo^s))==1)
             patchf(txoo^s) = 4
          endwhere
       endif
    endif
    {enddo^d&\}
 
  end subroutine rhd_diffuse_hllcd
 
  subroutine rhd_get_lcd(wLC,wRC,fLC,fRC,cmin,cmax,idim,ixI^L,ixO^L, &
       whll,Fhll,lambdaCD,patchf)
    ! Calculate lambda at CD and set the patchf to know the orientation
    ! of the riemann fan and decide on the flux choice
    ! We also compute here the HLL flux and w value, for fallback strategy
 
    ! In this nul version, we simply compute nothing and ensure TVDLF fallback
    use mod_global_parameters
 
    integer, intent(in)                                      :: ixi^l,ixo^l,idim
    double precision, dimension(ixI^S,1:nw), intent(in)      :: wlc,wrc
    double precision, dimension(ixI^S,1:nwflux), intent(in)  :: flc,frc
    double precision, dimension(ixI^S), intent(in)           :: cmax,cmin
    integer         , dimension(ixI^S), intent(inout)        :: patchf
    double precision, dimension(ixI^S,1:nwflux), intent(out) :: fhll,whll
    double precision, dimension(ixI^S), intent(out)          :: lambdacd
 
    logical         , dimension(ixI^S)     :: cond_patchf
    double precision                       :: epsilon
    integer                                :: iw,ix^d
 
    !--------------------------------------------
    ! on entry, patch is preset to contain values from -2,1,2,4
    !      -2: take left flux, no computation here
    !      +2: take right flux, no computation here
    !      +4: take TVDLF flux, no computation here
    !       1: compute the characteristic speed for the CD
 
    cond_patchf(ixo^s)=(abs(patchf(ixo^s))==1)
    lambdacd=0.d0
 
    do iw=1,nwflux
      where(cond_patchf(ixo^s))
        !============= compute HLL flux ==============!
        fhll(ixo^s,iw)= (cmax(ixo^s)*flc(ixo^s,iw)-cmin(ixo^s)*frc(ixo^s,iw) &
             + cmin(ixo^s)*cmax(ixo^s)*(wrc(ixo^s,iw)-wlc(ixo^s,iw)))&
             /(cmax(ixo^s)-cmin(ixo^s))
        !======== compute intermediate HLL state =======!
        whll(ixo^s,iw) = (cmax(ixo^s)*wrc(ixo^s,iw)-cmin(ixo^s)*wlc(ixo^s,iw)&
             +flc(ixo^s,iw)-frc(ixo^s,iw))/(cmax(ixo^s)-cmin(ixo^s))
      end where
    end do
 
    ! deduce the characteristic speed at the CD
    where(cond_patchf(ixo^s))
       lambdacd(ixo^s)=whll(ixo^s,mom(idim))/whll(ixo^s,rho_)
    end where
 
    {do ix^db=ixomin^db,ixomax^db\}
       if(cond_patchf(ix^d)) then
         ! double check whether obtained speed is in between min and max speeds given
         ! and identify in which part of the Riemann fan the time-axis is
         if(cmin(ix^d) < zero .and. lambdacd(ix^d)>zero&
              .and.lambdacd(ix^d)<cmax(ix^d)) then
            patchf(ix^d) = -1
         else if(cmax(ix^d) > zero .and. lambdacd(ix^d) < zero&
              .and.lambdacd(ix^d)>cmin(ix^d)) then
            patchf(ix^d) =  1
         else if(lambdacd(ix^d) >= cmax(ix^d) .or. &
              lambdacd(ix^d) <= cmin(ix^d)) then
            lambdacd(ix^d) = zero
            ! we will fall back to HLL flux case in this degeneracy
            patchf(ix^d) =  3
         end if
       end if
    {end do\}
 
    where(patchf(ixo^s)== 3)
      cond_patchf(ixo^s)=.false.
    end where
 
    ! handle the specific case where the time axis is exactly on the CD
    if(any(cond_patchf(ixo^s).and.lambdacd(ixo^s)==zero))then
      ! determine which sector (forward or backward) of the Riemann fan is smallest
      ! and select left or right flux accordingly
      {do ix^db=ixomin^db,ixomax^db\}
        if(lambdacd(ix^d)==zero.and.cond_patchf(ix^d)) then
          if(-cmin(ix^d)>=cmax(ix^d)) then
            patchf(ix^d) =  1
          else
            patchf(ix^d) = -1
          end if
        end if
      {end do\}
    end if
 
    ! eigenvalue lambda for contact is near zero: decrease noise by this trick
    if(flathllc)then
      epsilon=1.0d-6
      where(cond_patchf(ixo^s).and. &
          dabs(lambdacd(ixo^s))/max(cmax(ixo^s),epsilon)< epsilon  .and. &
          dabs(lambdacd(ixo^s))/max(dabs(cmin(ixo^s)),epsilon)< epsilon)
        lambdacd(ixo^s) =  zero
      end where
    end if
 
  end subroutine rhd_get_lcd
 
  subroutine rhd_get_wcd(wLC,wRC,whll,fRC,fLC,Fhll,patchf,lambdaCD,cmin,cmax,&
       ixI^L,ixO^L,idim,f)
    ! compute the intermediate state U*
    ! only needed where patchf=-1/1
 
    ! reference Li S., JCP, 203, 2005, 344-357
    ! reference T. Miyoski, Kusano JCP, 2008, 2005.
 
    use mod_global_parameters
    use mod_dust
 
    integer, intent(in)                                      :: ixi^l,ixo^l,idim
    double precision, dimension(ixI^S,1:nw), intent(in)      :: wrc,wlc
    double precision, dimension(ixI^S,1:nwflux), intent(in):: whll, fhll
    double precision, dimension(ixI^S), intent(in)           :: lambdacd
    double precision, dimension(ixI^S), intent(in)           :: cmax,cmin
    double precision, dimension(ixI^S,1:nwflux), intent(in):: frc,flc
    double precision, dimension(ixI^S,1:nwflux),intent(out):: f
    double precision, dimension(ixI^S,1:nw)      :: wcd,wsub
    double precision, dimension(ixI^S,1:nwflux)  :: fsub
    double precision, dimension(ixI^S)           :: vsub,cspeed,pcd
    integer         , dimension(ixI^S), intent(in)           :: patchf
 
    double precision :: csmls
    integer                                      :: n, iw, ix^d
 
    !-------------- auxiliary Speed and array-------------!
    {do ix^db=ixomin^db,ixomax^db\}
       if(patchf(ix^d)==1) then
         cspeed(ix^d)=cmax(ix^d)
         vsub(ix^d)=wrc(ix^d,mom(idim))/wrc(ix^d,rho_)
         wsub(ix^d,:)=wrc(ix^d,:)
         fsub(ix^d,:)=frc(ix^d,:)
       else if(patchf(ix^d)==-1) then
         cspeed(ix^d)=cmin(ix^d)
         vsub(ix^d)=wlc(ix^d,mom(idim))/wlc(ix^d,rho_)
         wsub(ix^d,:)=wlc(ix^d,:)
         fsub(ix^d,:)=flc(ix^d,:)
       end if
    {end do\}
 
    {do ix^db=ixomin^db,ixomax^db\}
      if(abs(patchf(ix^d))==1) then
        csmls=one/(cspeed(ix^d)-lambdacd(ix^d))
        wcd(ix^d,rho_) = wsub(ix^d,rho_)&
                         *(cspeed(ix^d)-vsub(ix^d))*csmls
        do n=1,rhd_n_tracer
          iw = tracer(n)
          wcd(ix^d,iw) = wsub(ix^d,iw)*(cspeed(ix^d)-vsub(ix^d))*csmls
        end do
 
        !------- Momentum ------!
        do iw=1, ndir
          if(iw /= idim)then
            ! eq. 21 22
            wcd(ix^d,mom(iw))=(cspeed(ix^d)*wsub(ix^d,mom(iw))-fsub(ix^d,mom(iw)))*csmls
          else
            ! eq. 20
            wcd(ix^d,mom(iw)) =  wcd(ix^d,rho_) * lambdacd(ix^d)
          endif
        enddo
        if(rhd_energy) then
          pcd(ix^d) = wsub(ix^d,rho_)*(cspeed(ix^d)-vsub(ix^d))&
                        *(lambdacd(ix^d)-vsub(ix^d))&
                        +fsub(ix^d,mom(idim))-wsub(ix^d,mom(idim))*vsub(ix^d)
          ! Eq 31
          wcd(ix^d,e_) = (cspeed(ix^d)*wsub(ix^d,e_) &
                          -fsub(ix^d,e_)+lambdacd(ix^d)*pcd(ix^d))*csmls
        end if
        !if(rhd_dust) then
        !  do n=1,dust_n_species
        !    wCD(ix^D,dust_rho(n)) = wSub(ix^D,dust_rho(n))*(cspeed(ix^D)-vSub(ix^D))*csmls
        !    do iw=1,ndir
        !      if(iw /= idim)then
        !        ! eq. 21 22
        !        wCD(ix^D,dust_mom(iw,n))=wSub(ix^D,dust_mom(iw,n))*(cspeed(ix^D)-vSub(ix^D))*csmls
        !      else
        !        ! eq. 20
        !        wCD(ix^D,dust_mom(iw,n)) =  wCD(ix^D,dust_rho(n)) * lambdaCD(ix^D)
        !      endif
        !    end do
        !  end do
        !end if
 
        if(rhd_dust) then
          do iw=1,nwflux
            if(iw>=dust_rho(1).and.iw<=dust_mom(ndir,dust_n_species)) then
              ! use HLL flux for dust
              f(ix^d,iw)=fhll(ix^d,iw)
            else
              ! f_i=fsub+lambda (wCD-wSub)
              f(ix^d,iw)=fsub(ix^d,iw)+cspeed(ix^d)*(wcd(ix^d,iw)-wsub(ix^d,iw))
            end if
          end do
        else
          do iw=1,nwflux
            ! f_i=fsub+lambda (wCD-wSub)
            f(ix^d,iw)=fsub(ix^d,iw)+cspeed(ix^d)*(wcd(ix^d,iw)-wsub(ix^d,iw))
          end do
        end if
      end if
    {end do\}
 
  end subroutine rhd_get_wcd
 
end module mod_rhd_hllc