mod_twofl_hllc.t

Go to the documentation of this file.

module mod_twofl_hllc
#include "amrvac.h"
  use mod_twofl_phys
  use mod_physics_hllc
  use mod_functions_bfield, only: mag
 
  implicit none
  private
 
 
  public :: twofl_hllc_init
 
  integer :: i_rho, i_e
  integer, allocatable :: i_mom(:)
 
 
contains
 
  subroutine twofl_hllc_init()
    use mod_physics_hllc
    use mod_global_parameters
 
    allocate(i_mom(1:ndir))
 
    phys_hllc_init_species => twofl_hllc_init_species
 
  end subroutine twofl_hllc_init
 
  subroutine twofl_hllc_init_species(ii, rho_, mom, e_)
    use mod_global_parameters
    integer, intent(in) :: ii
    integer, intent(out) :: rho_, mom(1:ndir),e_
 
    if(ii==1) then
      phys_diffuse_hllcd => twofl_diffuse_hllcd_c
      phys_get_lcd => twofl_get_lcd_c
      phys_get_wcd => twofl_get_wcd_c
      i_rho = rho_c_
      i_mom(1:ndir) = mom_c(1:ndir)
      i_e = e_c_
    else
      phys_diffuse_hllcd => twofl_diffuse_hllcd_n
      phys_get_lcd => twofl_get_lcd_n
      phys_get_wcd => twofl_get_wcd_n
      i_rho = rho_n_
      i_mom(1:ndir) = mom_n(1:ndir)
      i_e = e_n_
    end if
 
    rho_ = i_rho
    mom(:) = i_mom(:)
    e_ = i_e
 
  end subroutine twofl_hllc_init_species
 
  subroutine twofl_diffuse_hllcd_n(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 twofl_diffuse_hllcd_n
 
  subroutine twofl_get_lcd_n(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,i_mom(idim))/whll(ixo^s,i_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 twofl_get_lcd_n
 
  subroutine twofl_get_wcd_n(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_physics
 
    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,i_mom(idim))/wrc(ix^d,i_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,i_mom(idim))/wlc(ix^d,i_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,i_rho) = wsub(ix^d,i_rho)&
                         *(cspeed(ix^d)-vsub(ix^d))*csmls
 
        !------- Momentum ------!
        do iw=1, ndir
          if(iw /= idim)then
            ! eq. 21 22
            wcd(ix^d,i_mom(iw))=(cspeed(ix^d)*wsub(ix^d,i_mom(iw))-fsub(ix^d,i_mom(iw)))*csmls
          else
            ! eq. 20
            wcd(ix^d,i_mom(iw)) =  wcd(ix^d,i_rho) * lambdacd(ix^d)
          endif
        enddo
        if(phys_energy) then
          pcd(ix^d) = wsub(ix^d,i_rho)*(cspeed(ix^d)-vsub(ix^d))&
                        *(lambdacd(ix^d)-vsub(ix^d))&
                        +fsub(ix^d,i_mom(idim))-wsub(ix^d,i_mom(idim))*vsub(ix^d)
          ! Eq 31
          wcd(ix^d,i_e) = (cspeed(ix^d)*wsub(ix^d,i_e) &
                          -fsub(ix^d,i_e)+lambdacd(ix^d)*pcd(ix^d))*csmls
        end if
        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 do\}
 
  end subroutine twofl_get_wcd_n
 
  subroutine twofl_diffuse_hllcd_c(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 twofl_diffuse_hllcd_c
 
  subroutine twofl_get_lcd_c(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
  
    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
 
    ! 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)
    
    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))
      endwhere
    enddo
    
    ! deduce the characteristic speed at the CD
    where(cond_patchf(ixo^s))
      lambdacd(ixo^s)=whll(ixo^s,i_mom(idim))/whll(ixo^s,i_rho)
    end where
    
    
    where(cond_patchf(ixo^s))
      ! 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
      where(cmin(ixo^s)<zero.and.lambdacd(ixo^s)>zero&
            .and.lambdacd(ixo^s)<cmax(ixo^s))
        patchf(ixo^s) = -1
      elsewhere(cmax(ixo^s)>zero.and.lambdacd(ixo^s)<zero&
                .and.lambdacd(ixo^s)>cmin(ixo^s))
        patchf(ixo^s) =  1
      elsewhere(lambdacd(ixo^s)>=cmax(ixo^s).or.lambdacd(ixo^s) <= cmin(ixo^s))
        lambdacd(ixo^s) = zero
        ! we will fall back to HLL flux case in this degeneracy
        patchf(ixo^s) =  3
      endwhere
    endwhere
    
    
    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(lambdacd(ixo^s)==zero.and.cond_patchf(ixo^s)))then
      ! determine which sector (forward or backward) of the Riemann fan is smallest
      ! and select left or right flux accordingly
      where(lambdacd(ixo^s)==zero.and.cond_patchf(ixo^s))
        where(-cmin(ixo^s)>=cmax(ixo^s))
          patchf(ixo^s) =  1
        elsewhere
          patchf(ixo^s) = -1
        endwhere
      endwhere
    endif
    
    ! 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 twofl_get_lcd_c
 
  subroutine twofl_get_wcd_c(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_physics
    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,vdotbcd
    integer         , dimension(ixI^S), intent(in)           :: patchf
 
    integer                                        :: n, iw, idir,ix^d
    double precision, dimension(ixI^S)           :: rho
    
 
    !-------------- 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,i_mom(idim))/wrc(ix^d,i_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,i_mom(idim))/wlc(ix^d,i_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
        wcd(ix^d,i_rho) = wsub(ix^d,i_rho)&
                         *(cspeed(ix^d)-vsub(ix^d))/(cspeed(ix^d)-lambdacd(ix^d))
         !TODO tracer 
!        do n=1,twofl_n_tracer
!          iw = tracer(n)
!          wCD(ix^D,iw) = wSub(ix^D,iw)*(cspeed(ix^D)-vSub(ix^D))&
!             /(cspeed(ix^D)-lambdaCD(ix^D))
!        end do
        !==== Magnetic field ====!
        do idir=1,ndir
          ! case from eq 31
          wcd(ix^d,mag(idir)) = whll(ix^d,mag(idir))
        end do
        !------- Momentum ------!
        do iw=1, ndir
          if(iw /= idim)then
            ! eq. 21 22
            wcd(ix^d,i_mom(iw))=(cspeed(ix^d)*wsub(ix^d,i_mom(iw))-fsub(ix^d,i_mom(iw))&
                -wcd(ix^d,mag(idim))*wcd(ix^d,mag(iw)))/&
                (cspeed(ix^d)-lambdacd(ix^d))
          else
            ! eq. 20
            wcd(ix^d,i_mom(iw)) =  wcd(ix^d,i_rho) * lambdacd(ix^d)
          endif
        enddo
        if(phys_energy) then
          vdotbcd(ix^d) = sum(whll(ix^d,i_mom(:))*whll(ix^d,mag(:)))/whll(ix^d,i_rho)
          ! Eq 17
          pcd(ix^d)  = wsub(ix^d,i_rho)*(cspeed(ix^d)-vsub(ix^d))&
                        *(lambdacd(ix^d)-vsub(ix^d))&
                        +fsub(ix^d,i_mom(idim))-wsub(ix^d,i_mom(idim))*vsub(ix^d)&
                        + wcd(ix^d,mag(idim))**2
          ! Eq 31
          wcd(ix^d,i_e) = (cspeed(ix^d)*wsub(ix^d,i_e) &
                          -fsub(ix^d,i_e)+lambdacd(ix^d)*pcd(ix^d)&
                          -vdotbcd(ix^d)*wcd(ix^d,mag(idim)))&
                          /(cspeed(ix^d)-lambdacd(ix^d))
        end if
      end if
    {end do\}
 
    do iw=1,nwflux
     if(iw == mag(idim)) then
       f(ixo^s,iw)=zero
     else if(twofl_glm .and. iw == psi_) then
       f(ixo^s,iw)=zero
     else
       where(abs(patchf(ixo^s))==1)
         ! f_i=fsub+lambda (wCD-wSub)
         f(ixo^s,iw)=fsub(ixo^s,iw)+cspeed(ixo^s)*(wcd(ixo^s,iw)-wsub(ixo^s,iw))
       endwhere
     end if
    end do
 
  end subroutine twofl_get_wcd_c
 
end module mod_twofl_hllc