gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

SCHEMES: Fix adiabatic HLLD solver for weak Bx.

Browse Source 
If the parallel component of magnetic field is too small, the
intermediate states might produce numerical instabilities, since they
are obtained by the division by a small factor.

In order to remove this situation, we apply full HLLD solver for strong
enough Alfvén wave. If it is weak, the HLLC solver is applied.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2017-04-10 08:03:50 -03:00
parent 8dc157620c
commit da2f2a1bf0
1 changed files with 88 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

125
src/schemes.F90
View File

@ -3668,20 +3668,6 @@ module schemes
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)
! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to
! the HLLC solver
!
if (b2 > 0.0d+00) then ! Bₓ² > 0
! the total pressure, constant across the contact discontinuity and Alfvén waves
!
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
! speed differences
!
slmm = sl - sm
@ -3692,11 +3678,25 @@ module schemes
dnl = wl(idn) / slmm
dnr = wr(idn) / srmm
! square of Bₓ, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! the total pressure, constant across the contact discontinuity and Alfvén waves
!
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
! if Alfvén wave is strong enough, apply the full HLLD solver, otherwise
! revert to the HLLC one
!
if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε
! left and right Alfvén speeds
!
cal = - sqrt(b2 / dnl)
cal = sqrt(b2 / dnl)
car = sqrt(b2 / dnr)
sml = sm + cal
sml = sm - cal
smr = sm + car
! calculate division factors (also used to determine degeneracies)
@ -3814,10 +3814,10 @@ module schemes
! prepare constant primitive variables of the intermediate states
!
dv = car * dnr - cal * dnl
dv = car * dnr + cal * dnl
vy = (wcr(imy) - wcl(imy)) / dv
vz = (wcr(imz) - wcl(imz)) / dv
dv = car - cal
dv = car + cal
by = (wcr(iby) - wcl(iby)) / dv
bz = (wcr(ibz) - wcl(ibz)) / dv
vb = sm * bx + vy * by + vz * bz
@ -3837,7 +3837,7 @@ module schemes
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal
ui(ien) = (wcl(ien) + sm * pt - bx * vb) / (sml - sm)
! the inmost left intermediate flux
!
@ -3855,7 +3855,7 @@ module schemes
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / (smr - sm)
! the inmost right intermediate flux
!
@ -3914,7 +3914,7 @@ module schemes
! primitive variables for the inner left intermediate state
!
dv = srmm * dnr - cal * dnl
dv = srmm * dnr + cal * dnl
vy = (wr(imy) - wcl(imy)) / dv
vz = (wr(imz) - wcl(imz)) / dv
dv = sr - sml
@ -3932,7 +3932,7 @@ module schemes
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal
ui(ien) = (wcl(ien) + sm * pt - bx * vb) / (sml - sm)
! choose the correct state depending on the sign of contact discontinuity
! advection speed
@ -4013,7 +4013,7 @@ module schemes
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / (smr - sm)
! choose the correct state depending on the sign of contact discontinuity
! advection speed
@ -4048,24 +4048,77 @@ module schemes
end if ! one degeneracy
else ! Bₓ² = 0
else if (b2 > 0.0d+00) then ! Bₓ² > 0
! the total pressure, constant across the contact discontinuity
! constant intermediate state tangential components of velocity and
! magnetic field
!
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn
vy = (wr(imy) - wl(imy)) / dn
vz = (wr(imz) - wl(imz)) / dn
by = (wr(iby) - wl(iby)) / srml
bz = (wr(ibz) - wl(ibz)) / srml
! scalar product of velocity and magnetic field for the intermediate states
!
vb = sm * bx + vy * by + vz * bz
! separate intermediate states depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! the left speed difference
! conservative variables for the left intermediate state
!
slmm = sl - sm
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * vy
ui(imz) = dnl * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
! 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) = dnr * vy
ui(imz) = dnr * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
else ! Bₓ² = 0
! separate intermediate states 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) = wl(idn) / slmm
ui(imx) = ui(idn) * sm
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = wl(imy) / slmm
ui(imz) = wl(imz) / slmm
ui(ibx) = 0.0d+00
@ -4080,14 +4133,10 @@ module schemes
else if (sm < 0.0d+00) then ! sm < 0
! the right speed difference
!
srmm = sr - sm
! conservative variables for the right intermediate state
!
ui(idn) = wr(idn) / srmm
ui(imx) = ui(idn) * sm
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = wr(imy) / srmm
ui(imz) = wr(imz) / srmm
ui(ibx) = 0.0d+00
@ -4102,8 +4151,8 @@ module schemes
else ! sm = 0
! intermediate flux is constant across the contact discontinuity, so we end up
! having only one intermediate state, which is equivalent to the HLL solver
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml