mod_mhd_hllc.t

Go to the documentation of this file.

module mod_mhd_hllc
  use mod_mhd_phys
  use mod_functions_bfield, only: mag
 
  implicit none
  private
 
  public :: mhd_hllc_init
 
contains
 
  subroutine mhd_hllc_init()
    use mod_physics_hllc
 
    phys_diffuse_hllcd => mhd_diffuse_hllcd
    phys_get_lcd => mhd_get_lcd
    phys_get_wcd => mhd_get_wcd
 
  end subroutine mhd_hllc_init
 
  subroutine mhd_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 mhd_diffuse_hllcd
 
  subroutine mhd_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
    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(ixO^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(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
      {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 mhd_get_lcd
 
  subroutine mhd_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
    
    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
 
    !-------------- 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
        wcd(ix^d,rho_) = wsub(ix^d,rho_)&
                         *(cspeed(ix^d)-vsub(ix^d))/(cspeed(ix^d)-lambdacd(ix^d))
        do n=1,mhd_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,mom(iw))=(cspeed(ix^d)*wsub(ix^d,mom(iw))-fsub(ix^d,mom(iw))&
                -wcd(ix^d,mag(idim))*wcd(ix^d,mag(iw)))/&
                (cspeed(ix^d)-lambdacd(ix^d))
          else
            ! eq. 20
            wcd(ix^d,mom(iw)) =  wcd(ix^d,rho_) * lambdacd(ix^d)
          endif
        enddo
        if(mhd_energy) then
          vdotbcd(ix^d) = sum(whll(ix^d,mom(:))*whll(ix^d,mag(:)))/whll(ix^d,rho_)
          ! Eq 17
          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)&
                        + wcd(ix^d,mag(idim))**2
          ! Eq 31
          wcd(ix^d,e_) = (cspeed(ix^d)*wsub(ix^d,e_) &
                          -fsub(ix^d,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(mhd_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 mhd_get_wcd
 
end module mod_mhd_hllc