SCHEMES: Separate HLLD Riemann solver for isothermal MHD.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
5f366ab2eb
commit
1bddacc2d8
@ -292,7 +292,7 @@ module schemes
|
||||
|
||||
! set the solver pointer
|
||||
!
|
||||
numerical_flux => riemann_mhd_iso_hlld
|
||||
numerical_flux => numerical_flux_mhd_iso_hlld
|
||||
|
||||
case ("hlldm", "hlld-m", "HLLDM", "HLLD-M")
|
||||
|
||||
@ -1223,6 +1223,113 @@ module schemes
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine NUMERICAL_FLUX_MHD_ISO_HLLD:
|
||||
! --------------------------------------
|
||||
!
|
||||
! Subroutine solves one dimensional Riemann problem using the isothermal HLLD
|
||||
! method with correct state averaging and degeneracies treatement.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - the array of primitive variables at the Riemann states;
|
||||
! f - the output array of fluxes;
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Kowal, G.,
|
||||
! "An isothermal MHD Riemann solver with correct state averaging and
|
||||
! degeneracies treatement",
|
||||
! Journal of Computational Physics, in preparation
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine numerical_flux_mhd_iso_hlld(q, f)
|
||||
|
||||
! include external procedures
|
||||
!
|
||||
use equations, only : ibx, ibp
|
||||
use equations, only : cmax
|
||||
use equations, only : prim2cons_mhd_iso, fluxspeed_mhd_iso
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:,:), intent(inout) :: q
|
||||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i
|
||||
real(kind=8) :: bx, bp
|
||||
|
||||
! local arrays to store the states
|
||||
!
|
||||
real(kind=8), dimension(size(q,1),size(q,2)) :: ql, qr, ul, ur, fl, fr
|
||||
real(kind=8), dimension( 2,size(q,2)) :: cl, cr
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for Riemann solver
|
||||
!
|
||||
call start_timer(imr)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! copy state vectors
|
||||
!
|
||||
ql(:,:) = q(:,:,1)
|
||||
qr(:,:) = q(:,:,2)
|
||||
|
||||
! obtain the state values for Bx and Psi for the GLM-MHD equations
|
||||
!
|
||||
do i = 1, size(q, 2)
|
||||
|
||||
bx = 0.5d+00 * ((qr(ibx,i) + ql(ibx,i)) &
|
||||
- (qr(ibp,i) - ql(ibp,i)) / cmax)
|
||||
bp = 0.5d+00 * ((qr(ibp,i) + ql(ibp,i)) &
|
||||
- (qr(ibx,i) - ql(ibx,i)) * cmax)
|
||||
|
||||
ql(ibx,i) = bx
|
||||
qr(ibx,i) = bx
|
||||
ql(ibp,i) = bp
|
||||
qr(ibp,i) = bp
|
||||
|
||||
end do ! i = 1, n
|
||||
|
||||
! calculate corresponding conserved variables of the left and right states
|
||||
!
|
||||
call prim2cons_mhd_iso(ql(:,:), ul(:,:))
|
||||
call prim2cons_mhd_iso(qr(:,:), ur(:,:))
|
||||
|
||||
! calculate the physical fluxes and speeds at the states
|
||||
!
|
||||
call fluxspeed_mhd_iso(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
|
||||
call fluxspeed_mhd_iso(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
|
||||
|
||||
! get the HLLC fluxes
|
||||
!
|
||||
call riemann_mhd_iso_hlld(ql(:,:), qr(:,:), ul(:,:), ur(:,:), &
|
||||
fl(:,:), fr(:,:), cl(:,:), cr(:,:), f(:,:))
|
||||
|
||||
! higher order flux corrections
|
||||
!
|
||||
call higher_order_flux_correction(f)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for Riemann solver
|
||||
!
|
||||
call stop_timer(imr)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine numerical_flux_mhd_iso_hlld
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RIEMANN_HLL:
|
||||
! ----------------------
|
||||
!
|
||||
@ -1714,6 +1821,345 @@ module schemes
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RIEMANN_ISO_MHD_HLLD:
|
||||
! -------------------------------
|
||||
!
|
||||
! Subroutine solves one dimensional Riemann problem using the isothermal HLLD
|
||||
! method with correct state averaging and degeneracies treatement.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! ql, qr - the primitive variables of the Riemann states;
|
||||
! ul, ur - the conservative variables of the Riemann states;
|
||||
! fl, fr - the physical fluxes of the Riemann states;
|
||||
! f - the Riemann average flux;
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Kowal, G.,
|
||||
! "An isothermal MHD Riemann solver with correct state averaging and
|
||||
! degeneracies treatement",
|
||||
! Journal of Computational Physics, in preparation
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine riemann_mhd_iso_hlld(ql, qr, ul, ur, fl, fr, cl, cr, f)
|
||||
|
||||
! include external procedures
|
||||
!
|
||||
use equations, only : idn, imx, imy, imz, ibx, iby, ibz, ibp
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
|
||||
real(kind=8), dimension(:,:), intent(out) :: f
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i
|
||||
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
|
||||
real(kind=8) :: bx, bp, b2, dn, dnl, dnr, dvl, dvr
|
||||
|
||||
! local arrays to store the states
|
||||
!
|
||||
real(kind=8), dimension(size(ql,1)) :: wl, wr, wcl, wcr, ui
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! iterate over all points
|
||||
!
|
||||
do i = 1, size(ql, 2)
|
||||
|
||||
! estimate the minimum and maximum speeds
|
||||
!
|
||||
sl = min(cl(1,i), cr(1,i))
|
||||
sr = max(cl(2,i), cr(2,i))
|
||||
|
||||
! calculate the HLL flux
|
||||
!
|
||||
if (sl >= 0.0d+00) then
|
||||
|
||||
f(:,i) = fl(:,i)
|
||||
|
||||
else if (sr <= 0.0d+00) then
|
||||
|
||||
f(:,i) = fr(:,i)
|
||||
|
||||
else ! sl < 0 < sr
|
||||
|
||||
! calculate speed difference
|
||||
!
|
||||
srml = sr - sl
|
||||
|
||||
! calculate vectors of the left and right-going waves
|
||||
!
|
||||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||||
|
||||
! the advection speed in the intermediate states
|
||||
!
|
||||
dn = wr(idn) - wl(idn)
|
||||
sm = (wr(imx) - wl(imx)) / dn
|
||||
|
||||
! square of Bₓ, i.e. Bₓ²
|
||||
!
|
||||
bx = ql(ibx,i)
|
||||
b2 = ql(ibx,i) * qr(ibx,i)
|
||||
|
||||
! speed differences
|
||||
!
|
||||
slmm = sl - sm
|
||||
srmm = sr - sm
|
||||
|
||||
! left and right state densities
|
||||
!
|
||||
dnl = wl(idn) / slmm
|
||||
dnr = wr(idn) / srmm
|
||||
|
||||
! apply full HLLD solver only if the Alfvén wave is strong enough
|
||||
!
|
||||
if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε
|
||||
|
||||
! left and right Alfvén speeds
|
||||
!
|
||||
sml = sm - sqrt(b2 / dnl)
|
||||
smr = sm + sqrt(b2 / dnr)
|
||||
|
||||
! calculate denominators in order to check degeneracy
|
||||
!
|
||||
dvl = slmm * wl(idn) - b2
|
||||
dvr = srmm * wr(idn) - b2
|
||||
|
||||
! check degeneracy Sl* -> Sl or Sr* -> Sr
|
||||
!
|
||||
if (min(dvl, dvr) > 0.0d+00) then ! no degeneracy
|
||||
|
||||
! choose the correct state depending on the speed signs
|
||||
!
|
||||
if (sml >= 0.0d+00) then ! sl* ≥ 0
|
||||
|
||||
! conservative variables for the outer left intermediate state
|
||||
!
|
||||
ui(idn) = dnl
|
||||
ui(imx) = dnl * sm
|
||||
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||||
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||||
ui(ibp) = ul(ibp,i)
|
||||
|
||||
! the outer left intermediate flux
|
||||
!
|
||||
f(:,i) = sl * ui(:) - wl(:)
|
||||
|
||||
else if (smr <= 0.0d+00) then ! sr* ≤ 0
|
||||
|
||||
! conservative variables for the outer right intermediate state
|
||||
!
|
||||
ui(idn) = dnr
|
||||
ui(imx) = dnr * sm
|
||||
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||||
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||||
ui(ibp) = ur(ibp,i)
|
||||
|
||||
! the outer right intermediate flux
|
||||
!
|
||||
f(:,i) = sr * ui(:) - wr(:)
|
||||
|
||||
else ! sl* < 0 < sr*
|
||||
|
||||
! conservative variables for the outer left intermediate state
|
||||
!
|
||||
ui(idn) = dnl
|
||||
ui(imx) = dnl * sm
|
||||
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||||
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||||
ui(ibp) = ul(ibp,i)
|
||||
|
||||
! vector of the left-going Alfvén wave
|
||||
!
|
||||
wcl(:) = (sml - sl) * ui(:) + wl(:)
|
||||
|
||||
! conservative variables for the outer right intermediate state
|
||||
!
|
||||
ui(idn) = dnr
|
||||
ui(imx) = dnr * sm
|
||||
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||||
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||||
ui(ibp) = ur(ibp,i)
|
||||
|
||||
! vector of the right-going Alfvén wave
|
||||
!
|
||||
wcr(:) = (smr - sr) * ui(:) + wr(:)
|
||||
|
||||
! the flux corresponding to the middle intermediate state
|
||||
!
|
||||
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
|
||||
|
||||
end if ! sl* < 0 < sr*
|
||||
|
||||
else ! one degeneracy
|
||||
|
||||
! separate degeneracies
|
||||
!
|
||||
if (dvl > 0.0d+00) then ! sr* > sr
|
||||
|
||||
! conservative variables for the outer left intermediate state
|
||||
!
|
||||
ui(idn) = dnl
|
||||
ui(imx) = dnl * sm
|
||||
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||||
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||||
ui(ibp) = ul(ibp,i)
|
||||
|
||||
! choose the correct state depending on the speed signs
|
||||
!
|
||||
if (sml >= 0.0d+00) then ! sl* ≥ 0
|
||||
|
||||
! the outer left intermediate flux
|
||||
!
|
||||
f(:,i) = sl * ui(:) - wl(:)
|
||||
|
||||
else ! sl* < 0
|
||||
|
||||
! vector of the left-going Alfvén wave
|
||||
!
|
||||
wcl(:) = (sml - sl) * ui(:) + wl(:)
|
||||
|
||||
! the flux corresponding to the middle intermediate state
|
||||
!
|
||||
f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml)
|
||||
|
||||
end if ! sl* < 0
|
||||
|
||||
else if (dvr > 0.0d+00) then ! sl* < sl
|
||||
|
||||
! conservative variables for the outer right intermediate state
|
||||
!
|
||||
ui(idn) = dnr
|
||||
ui(imx) = dnr * sm
|
||||
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||||
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||||
ui(ibp) = ur(ibp,i)
|
||||
|
||||
! choose the correct state depending on the speed signs
|
||||
!
|
||||
if (smr <= 0.0d+00) then ! sr* ≤ 0
|
||||
|
||||
! the outer right intermediate flux
|
||||
!
|
||||
f(:,i) = sr * ui(:) - wr(:)
|
||||
|
||||
else ! sr* > 0
|
||||
|
||||
! vector of the right-going Alfvén wave
|
||||
!
|
||||
wcr(:) = (smr - sr) * ui(:) + wr(:)
|
||||
|
||||
! the flux corresponding to the middle intermediate state
|
||||
!
|
||||
f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl)
|
||||
|
||||
end if ! sr* > 0
|
||||
|
||||
else ! sl* < sl & sr* > sr
|
||||
|
||||
! both outer states are degenerate so revert to the HLL flux
|
||||
!
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
|
||||
end if ! sl* < sl & sr* > sr
|
||||
|
||||
end if ! one degeneracy
|
||||
|
||||
else if (b2 > 0.0d+00) then ! Bₓ² > 0
|
||||
|
||||
! the Alfvén wave is too weak, so ignore it in order to not introduce any
|
||||
! numerical instabilities; since Bₓ is not zero, the perpendicular components
|
||||
! of velocity and magnetic field are continuous; we have only one intermediate
|
||||
! state, therefore simply apply the HLL solver
|
||||
!
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
|
||||
else ! Bₓ² = 0
|
||||
|
||||
! take the right flux depending on the sign of the advection speed
|
||||
!
|
||||
if (sm > 0.0d+00) then ! sm > 0
|
||||
|
||||
! conservative variables for the left intermediate state
|
||||
!
|
||||
ui(idn) = dnl
|
||||
ui(imx) = dnl * sm
|
||||
ui(imy) = wl(imy) / slmm
|
||||
ui(imz) = wl(imz) / slmm
|
||||
ui(ibx) = 0.0d+00
|
||||
ui(iby) = wl(iby) / slmm
|
||||
ui(ibz) = wl(ibz) / slmm
|
||||
ui(ibp) = ul(ibp,i)
|
||||
|
||||
! the left intermediate flux
|
||||
!
|
||||
f(:,i) = sl * ui(:) - wl(:)
|
||||
|
||||
else if (sm < 0.0d+00) then ! sm < 0
|
||||
|
||||
! conservative variables for the right intermediate state
|
||||
!
|
||||
ui(idn) = dnr
|
||||
ui(imx) = dnr * sm
|
||||
ui(imy) = wr(imy) / srmm
|
||||
ui(imz) = wr(imz) / srmm
|
||||
ui(ibx) = 0.0d+00
|
||||
ui(iby) = wr(iby) / srmm
|
||||
ui(ibz) = wr(ibz) / srmm
|
||||
ui(ibp) = ur(ibp,i)
|
||||
|
||||
! the right intermediate flux
|
||||
!
|
||||
f(:,i) = sr * ui(:) - wr(:)
|
||||
|
||||
else ! sm = 0
|
||||
|
||||
! the intermediate flux; since the advection speed is zero, perpendicular
|
||||
! components do not change, so revert it to HLL
|
||||
!
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
|
||||
end if ! sm = 0
|
||||
|
||||
end if ! Bₓ² = 0
|
||||
|
||||
end if ! sl < 0 < sr
|
||||
|
||||
end do
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine riemann_mhd_iso_hlld
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine HIGHER_ORDER_FLUX_CORRECTION:
|
||||
! ---------------------------------------
|
||||
!
|
||||
@ -2102,399 +2548,6 @@ module schemes
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RIEMANN_MHD_ISO_HLLD:
|
||||
! -------------------------------
|
||||
!
|
||||
! Subroutine solves one dimensional Riemann problem using the isothermal HLLD
|
||||
! method with correct state averaging and degeneracies treatement.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - the array of primitive variables at the Riemann states;
|
||||
! f - the output array of fluxes;
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Kowal, G.,
|
||||
! "An isothermal MHD Riemann solver with correct state averaging and
|
||||
! degeneracies treatement",
|
||||
! Journal of Computational Physics, in preparation
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine riemann_mhd_iso_hlld(q, f)
|
||||
|
||||
! include external procedures
|
||||
!
|
||||
use equations, only : nf
|
||||
use equations, only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp
|
||||
use equations, only : cmax
|
||||
use equations, only : prim2cons_mhd_iso, fluxspeed_mhd_iso
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:,:), intent(inout) :: q
|
||||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: n, i
|
||||
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
|
||||
real(kind=8) :: bx, bp, b2, dn, dnl, dnr, dvl, dvr
|
||||
|
||||
! local arrays to store the states
|
||||
!
|
||||
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr
|
||||
real(kind=8), dimension(nf) :: wl, wr, wcl, wcr, ui
|
||||
real(kind=8), dimension(2,size(q,2)) :: cl, cr
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for Riemann solver
|
||||
!
|
||||
call start_timer(imr)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! get the length of vector
|
||||
!
|
||||
n = size(q, 2)
|
||||
|
||||
! copy state vectors
|
||||
!
|
||||
ql(:,:) = q(:,:,1)
|
||||
qr(:,:) = q(:,:,2)
|
||||
|
||||
! obtain the state values for Bx and Psi for the GLM-MHD equations
|
||||
!
|
||||
do i = 1, n
|
||||
|
||||
bx = 0.5d+00 * ((qr(ibx,i) + ql(ibx,i)) &
|
||||
- (qr(ibp,i) - ql(ibp,i)) / cmax)
|
||||
bp = 0.5d+00 * ((qr(ibp,i) + ql(ibp,i)) &
|
||||
- (qr(ibx,i) - ql(ibx,i)) * cmax)
|
||||
|
||||
ql(ibx,i) = bx
|
||||
qr(ibx,i) = bx
|
||||
ql(ibp,i) = bp
|
||||
qr(ibp,i) = bp
|
||||
|
||||
end do ! i = 1, n
|
||||
|
||||
! calculate corresponding conserved variables of the left and right states
|
||||
!
|
||||
call prim2cons_mhd_iso(ql(:,:), ul(:,:))
|
||||
call prim2cons_mhd_iso(qr(:,:), ur(:,:))
|
||||
|
||||
! calculate the physical fluxes and speeds at the states
|
||||
!
|
||||
call fluxspeed_mhd_iso(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
|
||||
call fluxspeed_mhd_iso(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
|
||||
|
||||
! iterate over all points
|
||||
!
|
||||
do i = 1, n
|
||||
|
||||
! estimate the minimum and maximum speeds
|
||||
!
|
||||
sl = min(cl(1,i), cr(1,i))
|
||||
sr = max(cl(2,i), cr(2,i))
|
||||
|
||||
! calculate the HLL flux
|
||||
!
|
||||
if (sl >= 0.0d+00) then
|
||||
|
||||
f(:,i) = fl(:,i)
|
||||
|
||||
else if (sr <= 0.0d+00) then
|
||||
|
||||
f(:,i) = fr(:,i)
|
||||
|
||||
else ! sl < 0 < sr
|
||||
|
||||
! calculate speed difference
|
||||
!
|
||||
srml = sr - sl
|
||||
|
||||
! calculate vectors of the left and right-going waves
|
||||
!
|
||||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||||
|
||||
! the advection speed in the intermediate states
|
||||
!
|
||||
dn = wr(idn) - wl(idn)
|
||||
sm = (wr(imx) - wl(imx)) / dn
|
||||
|
||||
! square of Bₓ, i.e. Bₓ²
|
||||
!
|
||||
bx = ql(ibx,i)
|
||||
b2 = ql(ibx,i) * qr(ibx,i)
|
||||
|
||||
! speed differences
|
||||
!
|
||||
slmm = sl - sm
|
||||
srmm = sr - sm
|
||||
|
||||
! left and right state densities
|
||||
!
|
||||
dnl = wl(idn) / slmm
|
||||
dnr = wr(idn) / srmm
|
||||
|
||||
! apply full HLLD solver only if the Alfvén wave is strong enough
|
||||
!
|
||||
if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε
|
||||
|
||||
! left and right Alfvén speeds
|
||||
!
|
||||
sml = sm - sqrt(b2 / dnl)
|
||||
smr = sm + sqrt(b2 / dnr)
|
||||
|
||||
! calculate denominators in order to check degeneracy
|
||||
!
|
||||
dvl = slmm * wl(idn) - b2
|
||||
dvr = srmm * wr(idn) - b2
|
||||
|
||||
! check degeneracy Sl* -> Sl or Sr* -> Sr
|
||||
!
|
||||
if (min(dvl, dvr) > 0.0d+00) then ! no degeneracy
|
||||
|
||||
! choose the correct state depending on the speed signs
|
||||
!
|
||||
if (sml >= 0.0d+00) then ! sl* ≥ 0
|
||||
|
||||
! conservative variables for the outer left intermediate state
|
||||
!
|
||||
ui(idn) = dnl
|
||||
ui(imx) = dnl * sm
|
||||
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||||
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||||
ui(ibp) = ul(ibp,i)
|
||||
|
||||
! the outer left intermediate flux
|
||||
!
|
||||
f(:,i) = sl * ui(:) - wl(:)
|
||||
|
||||
else if (smr <= 0.0d+00) then ! sr* ≤ 0
|
||||
|
||||
! conservative variables for the outer right intermediate state
|
||||
!
|
||||
ui(idn) = dnr
|
||||
ui(imx) = dnr * sm
|
||||
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||||
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||||
ui(ibp) = ur(ibp,i)
|
||||
|
||||
! the outer right intermediate flux
|
||||
!
|
||||
f(:,i) = sr * ui(:) - wr(:)
|
||||
|
||||
else ! sl* < 0 < sr*
|
||||
|
||||
! conservative variables for the outer left intermediate state
|
||||
!
|
||||
ui(idn) = dnl
|
||||
ui(imx) = dnl * sm
|
||||
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||||
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||||
ui(ibp) = ul(ibp,i)
|
||||
|
||||
! vector of the left-going Alfvén wave
|
||||
!
|
||||
wcl(:) = (sml - sl) * ui(:) + wl(:)
|
||||
|
||||
! conservative variables for the outer right intermediate state
|
||||
!
|
||||
ui(idn) = dnr
|
||||
ui(imx) = dnr * sm
|
||||
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||||
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||||
ui(ibp) = ur(ibp,i)
|
||||
|
||||
! vector of the right-going Alfvén wave
|
||||
!
|
||||
wcr(:) = (smr - sr) * ui(:) + wr(:)
|
||||
|
||||
! the flux corresponding to the middle intermediate state
|
||||
!
|
||||
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
|
||||
|
||||
end if ! sl* < 0 < sr*
|
||||
|
||||
else ! one degeneracy
|
||||
|
||||
! separate degeneracies
|
||||
!
|
||||
if (dvl > 0.0d+00) then ! sr* > sr
|
||||
|
||||
! conservative variables for the outer left intermediate state
|
||||
!
|
||||
ui(idn) = dnl
|
||||
ui(imx) = dnl * sm
|
||||
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||||
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||||
ui(ibp) = ul(ibp,i)
|
||||
|
||||
! choose the correct state depending on the speed signs
|
||||
!
|
||||
if (sml >= 0.0d+00) then ! sl* ≥ 0
|
||||
|
||||
! the outer left intermediate flux
|
||||
!
|
||||
f(:,i) = sl * ui(:) - wl(:)
|
||||
|
||||
else ! sl* < 0
|
||||
|
||||
! vector of the left-going Alfvén wave
|
||||
!
|
||||
wcl(:) = (sml - sl) * ui(:) + wl(:)
|
||||
|
||||
! the flux corresponding to the middle intermediate state
|
||||
!
|
||||
f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml)
|
||||
|
||||
end if ! sl* < 0
|
||||
|
||||
else if (dvr > 0.0d+00) then ! sl* < sl
|
||||
|
||||
! conservative variables for the outer right intermediate state
|
||||
!
|
||||
ui(idn) = dnr
|
||||
ui(imx) = dnr * sm
|
||||
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||
ui(ibx) = bx
|
||||
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||||
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||||
ui(ibp) = ur(ibp,i)
|
||||
|
||||
! choose the correct state depending on the speed signs
|
||||
!
|
||||
if (smr <= 0.0d+00) then ! sr* ≤ 0
|
||||
|
||||
! the outer right intermediate flux
|
||||
!
|
||||
f(:,i) = sr * ui(:) - wr(:)
|
||||
|
||||
else ! sr* > 0
|
||||
|
||||
! vector of the right-going Alfvén wave
|
||||
!
|
||||
wcr(:) = (smr - sr) * ui(:) + wr(:)
|
||||
|
||||
! the flux corresponding to the middle intermediate state
|
||||
!
|
||||
f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl)
|
||||
|
||||
end if ! sr* > 0
|
||||
|
||||
else ! sl* < sl & sr* > sr
|
||||
|
||||
! both outer states are degenerate so revert to the HLL flux
|
||||
!
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
|
||||
end if ! sl* < sl & sr* > sr
|
||||
|
||||
end if ! one degeneracy
|
||||
|
||||
else if (b2 > 0.0d+00) then ! Bₓ² > 0
|
||||
|
||||
! the Alfvén wave is too weak, so ignore it in order to not introduce any
|
||||
! numerical instabilities; since Bₓ is not zero, the perpendicular components
|
||||
! of velocity and magnetic field are continuous; we have only one intermediate
|
||||
! state, therefore simply apply the HLL solver
|
||||
!
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
|
||||
else ! Bₓ² = 0
|
||||
|
||||
! take the right flux depending on the sign of the advection speed
|
||||
!
|
||||
if (sm > 0.0d+00) then ! sm > 0
|
||||
|
||||
! conservative variables for the left intermediate state
|
||||
!
|
||||
ui(idn) = dnl
|
||||
ui(imx) = dnl * sm
|
||||
ui(imy) = wl(imy) / slmm
|
||||
ui(imz) = wl(imz) / slmm
|
||||
ui(ibx) = 0.0d+00
|
||||
ui(iby) = wl(iby) / slmm
|
||||
ui(ibz) = wl(ibz) / slmm
|
||||
ui(ibp) = ul(ibp,i)
|
||||
|
||||
! the left intermediate flux
|
||||
!
|
||||
f(:,i) = sl * ui(:) - wl(:)
|
||||
|
||||
else if (sm < 0.0d+00) then ! sm < 0
|
||||
|
||||
! conservative variables for the right intermediate state
|
||||
!
|
||||
ui(idn) = dnr
|
||||
ui(imx) = dnr * sm
|
||||
ui(imy) = wr(imy) / srmm
|
||||
ui(imz) = wr(imz) / srmm
|
||||
ui(ibx) = 0.0d+00
|
||||
ui(iby) = wr(iby) / srmm
|
||||
ui(ibz) = wr(ibz) / srmm
|
||||
ui(ibp) = ur(ibp,i)
|
||||
|
||||
! the right intermediate flux
|
||||
!
|
||||
f(:,i) = sr * ui(:) - wr(:)
|
||||
|
||||
else ! sm = 0
|
||||
|
||||
! the intermediate flux; since the advection speed is zero, perpendicular
|
||||
! components do not change, so revert it to HLL
|
||||
!
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
|
||||
end if ! sm = 0
|
||||
|
||||
end if ! Bₓ² = 0
|
||||
|
||||
end if ! sl < 0 < sr
|
||||
|
||||
end do ! i = 1, n
|
||||
|
||||
! higher order flux corrections
|
||||
!
|
||||
call higher_order_flux_correction(f)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for Riemann solver
|
||||
!
|
||||
call stop_timer(imr)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine riemann_mhd_iso_hlld
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RIEMANN_MHD_ISO_HLLDM:
|
||||
! --------------------------------
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user