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

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2020-08-20 17:40:07 -03:00
commit f6b2a61232
1 changed files with 352 additions and 423 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

453
sources/schemes.F90
View File

@ -1325,6 +1325,11 @@ module schemes
else ! sl < 0 < sr
! Bx and square of Bx, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! speed difference
!
srml = sr - sl
@ -1339,20 +1344,14 @@ module schemes
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
! square of Bx, i.e. Bₓ²
! the total pressure, constant across the contact discontinuity
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
! separate the cases when Bₓ = 0 and Bₓ ≠ 0
!
if (b2 > 0.0d+00) then
! the total pressure, constant across the contact discontinuity
!
pt = 0.5d+00 * ((wl(idn) * sm - wl(imx)) &
+ (wr(idn) * sm - wr(imx))) + b2
! constant intermediate state tangential components of velocity and
! magnetic field
!
@ -1422,11 +1421,6 @@ module schemes
else ! Bₓ = 0
! the total pressure, constant across the contact discontinuity
!
pt = 0.5d+00 * ((wl(idn) * sm - wl(imx)) &
+ (wr(idn) * sm - wr(imx)))
! separate intermediate states depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
@ -1475,18 +1469,10 @@ module schemes
else ! sm = 0
! intermediate flux is constant across the contact discontinuity and all except
! the parallel momentum flux are zero
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(idn,i) = 0.0d+00
f(imx,i) = - 0.5d+00 * (wl(imx) + wr(imx))
f(imy,i) = 0.0d+00
f(imz,i) = 0.0d+00
f(ibx,i) = fl(ibx,i)
f(iby,i) = 0.0d+00
f(ibz,i) = 0.0d+00
f(ibp,i) = fl(ibp,i)
f(ien,i) = 0.0d+00
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
@ -2071,7 +2057,7 @@ module schemes
!
integer :: i
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
real(kind=8) :: bx, b2, dn, dvl, dvr, ca, ca2
real(kind=8) :: bx, b2, dn, dvl, dvr, ca
! local arrays to store the states
!
@ -2100,6 +2086,11 @@ module schemes
else ! sl < 0 < sr
! square of Bₓ, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! calculate speed difference
!
srml = sr - sl
@ -2114,11 +2105,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)
! speed differences
!
slmm = sl - sm
@ -2128,32 +2114,24 @@ module schemes
!
dn = dn / srml
! get the square of the Alfvén speed
!
ca2 = b2 / dn
! apply full HLLD solver only if the Alfvén wave is strong enough
!
if (ca2 > 1.0d-08 * max(sl**2, sr**2)) then ! Bₓ² > ε
! left and right Alfvén speeds
!
ca = sqrt(ca2)
ca = abs(bx) / sqrt(dn)
sml = sm - ca
smr = sm + ca
! calculate denominators in order to check degeneracy
!
dvl = slmm * slmm * dn - b2
dvr = srmm * srmm * dn - b2
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
if (sml > sl .and. smr < sr) then ! no degeneracies
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
if (sml > 0.0d+00) then ! sl* > 0
! conservative variables for the outer left intermediate state
!
@ -2162,15 +2140,15 @@ module schemes
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(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
else if (smr < 0.0d+00) then ! sr* < 0
! conservative variables for the outer right intermediate state
!
@ -2179,15 +2157,19 @@ module schemes
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(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*
else ! sl* <= 0 <= sr*
! separate cases for non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! conservative variables for the outer left intermediate state
!
@ -2196,8 +2178,8 @@ module schemes
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(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
@ -2211,8 +2193,8 @@ module schemes
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(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
@ -2223,13 +2205,21 @@ module schemes
!
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if
end if ! sl* < 0 < sr*
else ! one degeneracy
else ! some degeneracies
! separate degeneracies
!
if (dvl > 0.0d+00) then ! sr* > sr
if (sml > sl) then ! sr* > sr
! conservative variables for the outer left intermediate state
!
@ -2238,19 +2228,23 @@ module schemes
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(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
if (sml > 0.0d+00) then ! sl* > 0
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else ! sl* < 0
else ! sl* <= 0
! separate cases for non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! vector of the left-going Alfvén wave
!
@ -2260,9 +2254,17 @@ module schemes
!
f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml)
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if
end if ! sl* < 0
else if (dvr > 0.0d+00) then ! sl* < sl
else if (smr < sr) then ! sl* < sl
! conservative variables for the outer right intermediate state
!
@ -2271,19 +2273,23 @@ module schemes
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(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
if (smr < 0.0d+00) then ! sr* < 0
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sr* > 0
else ! sr* >= 0
! separate cases for non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! vector of the right-going Alfvén wave
!
@ -2293,6 +2299,14 @@ module schemes
!
f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl)
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if
end if ! sr* > 0
else ! sl* < sl & sr* > sr
@ -2303,67 +2317,7 @@ module schemes
end if ! sl* < sl & sr* > sr
end if ! one degeneracy
else if (b2 > 0.0d+00) then ! Bₓ² > 0
! the Alfvén wave is very 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
! 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) = 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) = 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) = ur(ibp,i)
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! both states are equal so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
end if ! Bₓ² = 0
end if ! degeneracies
end if ! sl < 0 < sr
@ -2447,6 +2401,11 @@ module schemes
else ! sl < 0 < sr
! square of Bₓ, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! speed difference
!
srml = sr - sl
@ -2471,24 +2430,14 @@ 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)
car = sqrt(b2 / dnr)
cal = abs(bx) / sqrt(dnl)
car = abs(bx) / sqrt(dnr)
sml = sm - cal
smr = sm + car
@ -2499,11 +2448,11 @@ module schemes
! check degeneracy Sl* -> Sl or Sr* -> Sr
!
if (dvl > 0.0d+00 .and. dvr > 0.0d+00) then ! no degeneracy
if (sml > sl .and. smr < sr) then ! no degeneracy
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
if (sml > 0.0d+00) then ! sl* > 0
! primitive variables for the outer left intermediate state
!
@ -2529,7 +2478,7 @@ module schemes
!
f(:,i) = sl * ui(:) - wl(:)
else if (smr <= 0.0d+00) then ! sr* ≤ 0
else if (smr < 0.0d+00) then ! sr* < 0
! primitive variables for the outer right intermediate state
!
@ -2557,6 +2506,10 @@ module schemes
else ! sl* < 0 < sr*
! separate cases with non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! primitive variables for the outer left intermediate state
!
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
@ -2657,19 +2610,83 @@ module schemes
else ! sm = 0
! in the case when Sₘ = 0 and Bₓ² > 0, all variables are continuous, therefore
! the flux can be averaged from the Alfvén waves using the simple HLL formula
! the flux can be averaged from the Alfvén waves using a simple HLL formula;
!
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
end if ! sm = 0
else ! Bx = 0
if (sm > 0.0d+00) then
! primitive variables for the outer left intermediate state
!
vy = slmm * wl(imy) / dvl
vz = slmm * wl(imz) / dvl
by = wl(idn) * wl(iby) / dvl
bz = wl(idn) * wl(ibz) / dvl
vb = vy * by + vz * bz
! conservative variables for the outer left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * vy
ui(imz) = dnl * vz
ui(ibx) = 0.0d+00
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt) / slmm
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then
! primitive variables for the outer right intermediate state
!
vy = ( srmm * wr(imy)) / dvr
vz = ( srmm * wr(imz)) / dvr
by = (wr(idn) * wr(iby)) / dvr
bz = (wr(idn) * wr(ibz)) / dvr
vb = vy * by + vz * bz
! conservative variables for the outer right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * vy
ui(imz) = dnr * vz
ui(ibx) = 0.0d+00
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt) / srmm
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! Sm = 0
! since Bx = 0 and Sm = 0, then revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if
end if ! Bx = 0
end if ! sl* < 0 < sr*
else ! one degeneracy
else ! some degeneracies
! separate degeneracies
!
if (dvl > 0.0d+00) then ! sr* > sr
if (sml > sl) then ! sr* > sr
! primitive variables for the outer left intermediate state
!
@ -2693,13 +2710,17 @@ module schemes
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
if (sml > 0.0d+00) then ! sl* > 0
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else ! sl* < 0
else ! sl* <= 0
! separate cases with non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! vector of the left-going Alfvén wave
!
@ -2730,7 +2751,7 @@ module schemes
! choose the correct state depending on the sign of contact discontinuity
! advection speed
!
if (sm >= 0.0d+00) then ! sm ≥ 0
if (sm >= 0.0d+00) then ! sm >= 0
! the inner left intermediate flux
!
@ -2748,9 +2769,17 @@ module schemes
end if ! sm < 0
else ! bx = 0
! no Alfvén wave, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if
end if ! sl* < 0
else if (dvr > 0.0d+00) then ! sl* < sl
else if (smr < sr) then ! sl* < sl
! primitive variables for the outer right intermediate state
!
@ -2774,7 +2803,7 @@ module schemes
! choose the correct state depending on the speed signs
!
if (smr <= 0.0d+00) then ! sr* ≤ 0
if (smr < 0.0d+00) then ! sr* < 0
! the outer right intermediate flux
!
@ -2782,6 +2811,10 @@ module schemes
else ! sr* > 0
! separate cases with non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! vector of the right-going Alfvén wave
!
wcr(:) = (smr - sr) * ui(:) + wr(:)
@ -2811,7 +2844,7 @@ module schemes
! choose the correct state depending on the sign of contact discontinuity
! advection speed
!
if (sm <= 0.0d+00) then ! sm ≤ 0
if (sm <= 0.0d+00) then ! sm <= 0
! the inner right intermediate flux
!
@ -2829,6 +2862,14 @@ module schemes
end if ! sm > 0
else ! bx = 0
! no Alfvén wave, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if
end if ! sr* > 0
else ! sl* < sl & sr* > sr
@ -2839,119 +2880,7 @@ module schemes
end if ! sl* < sl & sr* > sr
end if ! one degeneracy
else if (b2 > 0.0d+00) then ! Bₓ² > 0
! constant intermediate state tangential components of velocity and
! magnetic field
!
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
! conservative variables for the left intermediate state
!
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) = 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)
ui(ien) = (wl(ien) + sm * pt) / 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) = 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)
ui(ien) = (wr(ien) + sm * pt) / 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
end if ! Bₓ² = 0
end if ! deneneracies
end if ! sl < 0 < sr