Merge branch 'master' into reconnection
This commit is contained in:
commit
bd3143c499
779
src/schemes.F90
779
src/schemes.F90
@ -195,6 +195,26 @@ module schemes
|
|||||||
!
|
!
|
||||||
select case(trim(solver))
|
select case(trim(solver))
|
||||||
|
|
||||||
|
case ("hlld", "HLLD")
|
||||||
|
|
||||||
|
! set the solver name
|
||||||
|
!
|
||||||
|
name_sol = "HLLD"
|
||||||
|
|
||||||
|
! set the solver pointer
|
||||||
|
!
|
||||||
|
riemann => riemann_mhd_iso_hlld
|
||||||
|
|
||||||
|
case ("hlldm", "hlld-m", "HLLDM", "HLLD-M")
|
||||||
|
|
||||||
|
! set the solver name
|
||||||
|
!
|
||||||
|
name_sol = "HLLD-M"
|
||||||
|
|
||||||
|
! set the solver pointer
|
||||||
|
!
|
||||||
|
riemann => riemann_mhd_iso_hlldm
|
||||||
|
|
||||||
! in the case of unknown Riemann solver, revert to HLL
|
! in the case of unknown Riemann solver, revert to HLL
|
||||||
!
|
!
|
||||||
case default
|
case default
|
||||||
@ -1593,6 +1613,765 @@ 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:
|
||||||
|
!
|
||||||
|
! n - the length of input vectors;
|
||||||
|
! h - the spatial step;
|
||||||
|
! q - the input array of primitive variables;
|
||||||
|
! 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(n, h, q, f)
|
||||||
|
|
||||||
|
! include external procedures
|
||||||
|
!
|
||||||
|
use equations , only : nv
|
||||||
|
use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp
|
||||||
|
use equations , only : cmax
|
||||||
|
use equations , only : prim2cons, fluxspeed
|
||||||
|
use interpolations, only : reconstruct, fix_positivity
|
||||||
|
|
||||||
|
! local variables are not implicit by default
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! subroutine arguments
|
||||||
|
!
|
||||||
|
integer , intent(in) :: n
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(nv,n), intent(in) :: q
|
||||||
|
real(kind=8), dimension(nv,n), intent(out) :: f
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
!
|
||||||
|
integer :: p, i
|
||||||
|
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
|
||||||
|
real(kind=8) :: bx, b2, dn, dnl, dnr, dvl, dvr
|
||||||
|
|
||||||
|
! local arrays to store the states
|
||||||
|
!
|
||||||
|
real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr
|
||||||
|
real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui
|
||||||
|
real(kind=8), dimension(n) :: cl, cr
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
! reconstruct the left and right states of primitive variables
|
||||||
|
!
|
||||||
|
do p = 1, nv
|
||||||
|
call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:))
|
||||||
|
end do
|
||||||
|
|
||||||
|
! obtain the state values for Bx and Psi for the GLM-MHD equations
|
||||||
|
!
|
||||||
|
cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax)
|
||||||
|
cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax)
|
||||||
|
ql(ibx,:) = cl(:)
|
||||||
|
qr(ibx,:) = cl(:)
|
||||||
|
ql(ibp,:) = cr(:)
|
||||||
|
qr(ibp,:) = cr(:)
|
||||||
|
|
||||||
|
! check if the reconstruction doesn't give the negative density or pressure,
|
||||||
|
! if so, correct the states
|
||||||
|
!
|
||||||
|
call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:))
|
||||||
|
|
||||||
|
! calculate corresponding conserved variables of the left and right states
|
||||||
|
!
|
||||||
|
call prim2cons(n, ql(:,:), ul(:,:))
|
||||||
|
call prim2cons(n, qr(:,:), ur(:,:))
|
||||||
|
|
||||||
|
! calculate the physical fluxes and speeds at the states
|
||||||
|
!
|
||||||
|
call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:))
|
||||||
|
call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:))
|
||||||
|
|
||||||
|
! iterate over all points
|
||||||
|
!
|
||||||
|
do i = 1, n
|
||||||
|
|
||||||
|
! estimate the minimum and maximum speeds
|
||||||
|
!
|
||||||
|
sl = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i))
|
||||||
|
sr = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(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
|
||||||
|
|
||||||
|
! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to
|
||||||
|
! the HLL one
|
||||||
|
!
|
||||||
|
if (b2 > 0.0d+00) then ! Bₓ² > 0
|
||||||
|
|
||||||
|
! 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) = ql(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) = qr(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) = ql(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) = qr(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) = ql(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) = qr(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 ! Bₓ² = 0
|
||||||
|
|
||||||
|
! in the vase of vanishing Bₓ there is no Alfvén wave, density is constant, and
|
||||||
|
! still we can have a discontinuity in the perpendicular components
|
||||||
|
!
|
||||||
|
dn = dn / srml
|
||||||
|
|
||||||
|
! 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) = dn
|
||||||
|
ui(imx) = dn * 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) = ql(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) = dn
|
||||||
|
ui(imx) = dn * 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) = qr(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
|
||||||
|
!
|
||||||
|
f(idn,i) = (sl * wr(idn) - sr * wl(idn)) / srml
|
||||||
|
f(imx,i) = (sl * wr(imx) - sr * wl(imx)) / srml
|
||||||
|
f(imy,i) = 0.0d+00
|
||||||
|
f(imz,i) = 0.0d+00
|
||||||
|
f(ibx,i) = (sl * wr(ibx) - sr * wl(ibx)) / srml
|
||||||
|
f(iby,i) = 0.0d+00
|
||||||
|
f(ibz,i) = 0.0d+00
|
||||||
|
f(ibp,i) = (sl * wr(ibp) - sr * wl(ibp)) / srml
|
||||||
|
|
||||||
|
end if ! sm = 0
|
||||||
|
|
||||||
|
end if ! Bₓ² = 0
|
||||||
|
|
||||||
|
end if ! sl < 0 < sr
|
||||||
|
|
||||||
|
end do ! i = 1, n
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
end subroutine riemann_mhd_iso_hlld
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
! subroutine RIEMANN_MHD_ISO_HLLDM:
|
||||||
|
! --------------------------------
|
||||||
|
!
|
||||||
|
! Subroutine solves one dimensional Riemann problem using the isothermal HLLD
|
||||||
|
! method described by Mignone.
|
||||||
|
!
|
||||||
|
! Arguments:
|
||||||
|
!
|
||||||
|
! n - the length of input vectors;
|
||||||
|
! h - the spatial step;
|
||||||
|
! q - the input array of primitive variables;
|
||||||
|
! f - the output array of fluxes;
|
||||||
|
!
|
||||||
|
! References:
|
||||||
|
!
|
||||||
|
! [1] Mignone, A.,
|
||||||
|
! "A simple and accurate Riemann solver for isothermal MHD",
|
||||||
|
! Journal of Computational Physics, 2007, 225, pp. 1427-1441
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
subroutine riemann_mhd_iso_hlldm(n, h, q, f)
|
||||||
|
|
||||||
|
! include external procedures
|
||||||
|
!
|
||||||
|
use equations , only : nv
|
||||||
|
use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp
|
||||||
|
use equations , only : cmax
|
||||||
|
use equations , only : prim2cons, fluxspeed
|
||||||
|
use interpolations, only : reconstruct, fix_positivity
|
||||||
|
|
||||||
|
! local variables are not implicit by default
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! subroutine arguments
|
||||||
|
!
|
||||||
|
integer , intent(in) :: n
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(nv,n), intent(in) :: q
|
||||||
|
real(kind=8), dimension(nv,n), intent(out) :: f
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
!
|
||||||
|
integer :: p, i
|
||||||
|
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
|
||||||
|
real(kind=8) :: bx, b2, dn, dnl, dnr, dvl, dvr, ca
|
||||||
|
|
||||||
|
! local arrays to store the states
|
||||||
|
!
|
||||||
|
real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr
|
||||||
|
real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui
|
||||||
|
real(kind=8), dimension(n) :: cl, cr
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
! reconstruct the left and right states of primitive variables
|
||||||
|
!
|
||||||
|
do p = 1, nv
|
||||||
|
call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:))
|
||||||
|
end do
|
||||||
|
|
||||||
|
! obtain the state values for Bx and Psi for the GLM-MHD equations
|
||||||
|
!
|
||||||
|
cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax)
|
||||||
|
cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax)
|
||||||
|
ql(ibx,:) = cl(:)
|
||||||
|
qr(ibx,:) = cl(:)
|
||||||
|
ql(ibp,:) = cr(:)
|
||||||
|
qr(ibp,:) = cr(:)
|
||||||
|
|
||||||
|
! check if the reconstruction doesn't give the negative density or pressure,
|
||||||
|
! if so, correct the states
|
||||||
|
!
|
||||||
|
call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:))
|
||||||
|
|
||||||
|
! calculate corresponding conserved variables of the left and right states
|
||||||
|
!
|
||||||
|
call prim2cons(n, ql(:,:), ul(:,:))
|
||||||
|
call prim2cons(n, qr(:,:), ur(:,:))
|
||||||
|
|
||||||
|
! calculate the physical fluxes and speeds at the states
|
||||||
|
!
|
||||||
|
call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:))
|
||||||
|
call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:))
|
||||||
|
|
||||||
|
! iterate over all points
|
||||||
|
!
|
||||||
|
do i = 1, n
|
||||||
|
|
||||||
|
! estimate the minimum and maximum speeds
|
||||||
|
!
|
||||||
|
sl = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i))
|
||||||
|
sr = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(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
|
||||||
|
|
||||||
|
! calculate density of the intermediate states
|
||||||
|
!
|
||||||
|
dn = dn / srml
|
||||||
|
|
||||||
|
! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to
|
||||||
|
! the HLL one
|
||||||
|
!
|
||||||
|
if (b2 > 0.0d+00) then ! Bₓ² > 0
|
||||||
|
|
||||||
|
! left and right Alfvén speeds
|
||||||
|
!
|
||||||
|
ca = sqrt(b2 / dn)
|
||||||
|
sml = sm - ca
|
||||||
|
smr = sm + ca
|
||||||
|
|
||||||
|
! calculate denominators in order to check degeneracy
|
||||||
|
!
|
||||||
|
dvl = slmm * slmm * dn - b2
|
||||||
|
dvr = srmm * srmm * dn - 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) = dn
|
||||||
|
ui(imx) = dn * sm
|
||||||
|
ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||||
|
ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||||
|
ui(ibx) = bx
|
||||||
|
ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
|
||||||
|
ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
|
||||||
|
ui(ibp) = ql(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) = dn
|
||||||
|
ui(imx) = dn * sm
|
||||||
|
ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||||
|
ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||||
|
ui(ibx) = bx
|
||||||
|
ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
|
||||||
|
ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
|
||||||
|
ui(ibp) = qr(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) = dn
|
||||||
|
ui(imx) = dn * sm
|
||||||
|
ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||||
|
ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||||
|
ui(ibx) = bx
|
||||||
|
ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
|
||||||
|
ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
|
||||||
|
ui(ibp) = ql(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) = dn
|
||||||
|
ui(imx) = dn * sm
|
||||||
|
ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||||
|
ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||||
|
ui(ibx) = bx
|
||||||
|
ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
|
||||||
|
ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
|
||||||
|
ui(ibp) = qr(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) = dn
|
||||||
|
ui(imx) = dn * sm
|
||||||
|
ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||||||
|
ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||||||
|
ui(ibx) = bx
|
||||||
|
ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
|
||||||
|
ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
|
||||||
|
ui(ibp) = ql(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) = dn
|
||||||
|
ui(imx) = dn * sm
|
||||||
|
ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||||||
|
ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||||||
|
ui(ibx) = bx
|
||||||
|
ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
|
||||||
|
ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
|
||||||
|
ui(ibp) = qr(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 ! Bₓ² = 0
|
||||||
|
|
||||||
|
! in the vase of vanishing Bₓ there is no Alfvén wave, density is constant, and
|
||||||
|
! still we can have a discontinuity in the perpendicular components
|
||||||
|
!
|
||||||
|
! 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) = dn
|
||||||
|
ui(imx) = dn * 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) = ql(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) = dn
|
||||||
|
ui(imx) = dn * 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) = qr(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
|
||||||
|
!
|
||||||
|
f(idn,i) = (sl * wr(idn) - sr * wl(idn)) / srml
|
||||||
|
f(imx,i) = (sl * wr(imx) - sr * wl(imx)) / srml
|
||||||
|
f(imy,i) = 0.0d+00
|
||||||
|
f(imz,i) = 0.0d+00
|
||||||
|
f(ibx,i) = (sl * wr(ibx) - sr * wl(ibx)) / srml
|
||||||
|
f(iby,i) = 0.0d+00
|
||||||
|
f(ibz,i) = 0.0d+00
|
||||||
|
f(ibp,i) = (sl * wr(ibp) - sr * wl(ibp)) / srml
|
||||||
|
|
||||||
|
end if ! sm = 0
|
||||||
|
|
||||||
|
end if ! Bₓ² = 0
|
||||||
|
|
||||||
|
end if ! sl < 0 < sr
|
||||||
|
|
||||||
|
end do ! i = 1, n
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
end subroutine riemann_mhd_iso_hlldm
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
! subroutine RIEMANN_MHD_ADI_HLL:
|
! subroutine RIEMANN_MHD_ADI_HLL:
|
||||||
! ------------------------------
|
! ------------------------------
|
||||||
!
|
!
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user