From 8634f89ce0e4859fd3cb00139a77aada5fbe5d5c Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 20 Aug 2020 11:27:00 -0300 Subject: [PATCH 1/3] SCHEMES: Properly average total pressure in MHD HLLC solver. Signed-off-by: Grzegorz Kowal --- sources/schemes.F90 | 46 ++++++++++++++++----------------------------- 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/sources/schemes.F90 b/sources/schemes.F90 index ed247cb..58ee975 100644 --- a/sources/schemes.F90 +++ b/sources/schemes.F90 @@ -1325,34 +1325,33 @@ 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 ! calculate vectors of the left and right-going waves ! - wl(:) = sl * ul(:,i) - fl(:,i) - wr(:) = sr * ur(:,i) - fr(:,i) + wl(:) = sl * ul(:,i) - fl(:,i) + wr(:) = sr * ur(:,i) - fr(:,i) ! the speed of contact discontinuity ! 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 ! @@ -1371,7 +1370,7 @@ module schemes ! the left speed difference ! - slmm = sl - sm + slmm = sl - sm ! conservative variables for the left intermediate state ! @@ -1393,7 +1392,7 @@ module schemes ! the right speed difference ! - srmm = sr - sm + srmm = sr - sm ! conservative variables for the right intermediate state ! @@ -1422,18 +1421,13 @@ 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 ! the left speed difference ! - slmm = sl - sm + slmm = sl - sm ! conservative variables for the left intermediate state ! @@ -1455,7 +1449,7 @@ module schemes ! the right speed difference ! - srmm = sr - sm + srmm = sr - sm ! conservative variables for the right intermediate state ! @@ -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 From 40619aa119a7163fcf3a819e47406de4b6784a00 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 20 Aug 2020 14:52:43 -0300 Subject: [PATCH 2/3] SCHEMES: Slightly rewrite HLLD solver for isothermal MHD. This fixes some weird behavior when Bx is small. Signed-off-by: Grzegorz Kowal --- sources/schemes.F90 | 248 +++++++++++++++++++------------------------- 1 file changed, 108 insertions(+), 140 deletions(-) diff --git a/sources/schemes.F90 b/sources/schemes.F90 index 58ee975..0b8a74f 100644 --- a/sources/schemes.F90 +++ b/sources/schemes.F90 @@ -2057,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 ! @@ -2086,25 +2086,25 @@ 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 ! calculate vectors of the left and right-going waves ! - wl(:) = sl * ul(:,i) - fl(:,i) - wr(:) = sr * ur(:,i) - fr(:,i) + 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 @@ -2114,76 +2114,72 @@ 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) - sml = sm - ca - smr = sm + ca + 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 ! - 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) = ul(ibp,i) + 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) = (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(:) + 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 ! - 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) = ur(ibp,i) + 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) = (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(:) + 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 ! ui(idn) = dn - ui(imx) = dn * sm + 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(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 @@ -2193,12 +2189,12 @@ module schemes ! conservative variables for the outer right intermediate state ! ui(idn) = dn - ui(imx) = dn * sm + 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(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 @@ -2209,34 +2205,46 @@ module schemes ! f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml) - end if ! sl* < 0 < sr* + else ! Bx = 0 -> Sₘ = 0 - else ! one degeneracy +! 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 ! 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 ! - 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) = ul(ibp,i) + 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) = (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(:) + 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 ! @@ -2246,30 +2254,42 @@ module schemes ! f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml) - end if ! sl* < 0 + else ! Bx = 0 -> Sₘ = 0 - else if (dvr > 0.0d+00) then ! sl* < sl +! 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 (smr < sr) 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) = ur(ibp,i) + 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) = (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(:) + 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 ! @@ -2279,77 +2299,25 @@ module schemes ! f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl) - end if ! sr* > 0 + else ! Bx = 0 -> Sₘ = 0 - else ! sl* < sl & sr* > sr +! 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 ! 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 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 ! sl* < sl & sr* > sr - end if ! Bₓ² = 0 + end if ! degeneracies end if ! sl < 0 < sr From cac043487c60e400e5c9a4b611d4cc739ec3ff09 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 20 Aug 2020 17:39:03 -0300 Subject: [PATCH 3/3] SCHEMES: Slightly rewrite HLLD solver for adiabatic MHD. This fixes some weird behavior when Bx is small too. Signed-off-by: Grzegorz Kowal --- sources/schemes.F90 | 481 +++++++++++++++++++++----------------------- 1 file changed, 228 insertions(+), 253 deletions(-) diff --git a/sources/schemes.F90 b/sources/schemes.F90 index 0b8a74f..5d7d646 100644 --- a/sources/schemes.F90 +++ b/sources/schemes.F90 @@ -2401,14 +2401,19 @@ 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 ! calculate vectors of the left and right-going waves ! - wl(:) = sl * ul(:,i) - fl(:,i) - wr(:) = sr * ur(:,i) - fr(:,i) + wl(:) = sl * ul(:,i) - fl(:,i) + wr(:) = sr * ur(:,i) - fr(:,i) ! the speed of contact discontinuity ! @@ -2425,99 +2430,93 @@ 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) - sml = sm - cal - smr = sm + car + cal = abs(bx) / sqrt(dnl) + car = abs(bx) / sqrt(dnr) + sml = sm - cal + smr = sm + car ! calculate division factors (also used to determine degeneracies) ! - dvl = slmm * wl(idn) - b2 - dvr = srmm * wr(idn) - b2 + dvl = slmm * wl(idn) - b2 + dvr = srmm * wr(idn) - b2 ! 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 ! - vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl - vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl - by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl - bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl - vb = sm * bx + vy * by + vz * bz + vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl + vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl + by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl + bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl + vb = sm * bx + 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) = bx - ui(iby) = by - ui(ibz) = bz - ui(ibp) = ul(ibp,i) - ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm + 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 outer left intermediate flux ! - f(:,i) = sl * ui(:) - wl(:) + 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 ! - vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr - vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr - by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr - bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr - vb = sm * bx + vy * by + vz * bz + vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr + vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr + by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr + bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr + vb = sm * bx + 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) = bx - ui(iby) = by - ui(ibz) = bz - ui(ibp) = ur(ibp,i) - ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm + 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 outer right intermediate flux ! - f(:,i) = sr * ui(:) - wr(:) + f(:,i) = sr * ui(:) - wr(:) - else ! sl* < 0 < sr* + 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 - vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl - by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl - bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl - vb = sm * bx + vy * by + vz * bz + vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl + vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl + by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl + bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl + vb = sm * bx + vy * by + vz * bz ! conservative variables for the outer left intermediate state ! @@ -2533,15 +2532,15 @@ module schemes ! vector of the left-going Alfvén wave ! - wcl(:) = (sml - sl) * ui(:) + wl(:) + wcl(:) = (sml - sl) * ui(:) + wl(:) ! primitive variables for the outer right intermediate state ! - vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr - vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr - by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr - bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr - vb = sm * bx + vy * by + vz * bz + vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr + vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr + by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr + bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr + vb = sm * bx + vy * by + vz * bz ! conservative variables for the outer right intermediate state ! @@ -2557,17 +2556,17 @@ module schemes ! vector of the right-going Alfvén wave ! - wcr(:) = (smr - sr) * ui(:) + wr(:) + wcr(:) = (smr - sr) * ui(:) + wr(:) ! prepare constant primitive variables of the intermediate states ! - dv = car * dnr + cal * dnl - vy = (wcr(imy) - wcl(imy)) / dv - vz = (wcr(imz) - wcl(imz)) / dv - dv = car + cal - by = (wcr(iby) - wcl(iby)) / dv - bz = (wcr(ibz) - wcl(ibz)) / dv - vb = sm * bx + vy * by + vz * bz + dv = car * dnr + cal * dnl + vy = (wcr(imy) - wcl(imy)) / dv + vz = (wcr(imz) - wcl(imz)) / dv + dv = car + cal + by = (wcr(iby) - wcl(iby)) / dv + bz = (wcr(ibz) - wcl(ibz)) / dv + vb = sm * bx + vy * by + vz * bz ! choose the correct state depending on the sign of contact discontinuity ! advection speed @@ -2588,7 +2587,7 @@ module schemes ! the inmost left intermediate flux ! - f(:,i) = sml * ui(:) - wcl(:) + f(:,i) = sml * ui(:) - wcl(:) else if (sm < 0.0d+00) then ! sm < 0 @@ -2606,54 +2605,122 @@ module schemes ! the inmost right intermediate flux ! - f(:,i) = smr * ui(:) - wcr(:) + f(:,i) = smr * ui(:) - wcr(:) 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 - end if ! sl* < 0 < sr* + else ! Bx = 0 - else ! one degeneracy - -! separate degeneracies -! - if (dvl > 0.0d+00) then ! sr* > sr + if (sm > 0.0d+00) then ! primitive variables for the outer left intermediate state ! - vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl - vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl - by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl - bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl - vb = sm * bx + vy * by + vz * bz + 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) = bx - ui(iby) = by - ui(ibz) = bz - ui(ibp) = ul(ibp,i) - ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm - -! choose the correct state depending on the speed signs -! - if (sml >= 0.0d+00) then ! sl* ≥ 0 + 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(:) + f(:,i) = sl * ui(:) - wl(:) - else ! sl* < 0 + 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 ! some degeneracies + +! separate degeneracies +! + if (sml > sl) then ! sr* > sr + +! primitive variables for the outer left intermediate state +! + vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl + vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl + by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl + bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl + vb = sm * bx + 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) = bx + ui(iby) = by + ui(ibz) = bz + ui(ibp) = ul(ibp,i) + ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm + +! 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 + +! separate cases with non-zero and zero Bx +! + if (b2 > 0.0d+00) then ! vector of the left-going Alfvén wave ! @@ -2661,13 +2728,13 @@ module schemes ! primitive variables for the inner left intermediate state ! - dv = srmm * dnr + cal * dnl - vy = (wr(imy) - wcl(imy)) / dv - vz = (wr(imz) - wcl(imz)) / dv - dv = sr - sml - by = (wr(iby) - wcl(iby)) / dv - bz = (wr(ibz) - wcl(ibz)) / dv - vb = sm * bx + vy * by + vz * bz + dv = srmm * dnr + cal * dnl + vy = (wr(imy) - wcl(imy)) / dv + vz = (wr(imz) - wcl(imz)) / dv + dv = sr - sml + by = (wr(iby) - wcl(iby)) / dv + bz = (wr(ibz) - wcl(ibz)) / dv + vb = sm * bx + vy * by + vz * bz ! conservative variables for the inner left intermediate state ! @@ -2684,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 ! @@ -2702,39 +2769,51 @@ module schemes end if ! sm < 0 - end if ! sl* < 0 + else ! bx = 0 - else if (dvr > 0.0d+00) then ! sl* < sl +! 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 (smr < sr) then ! sl* < sl ! primitive variables for the outer right intermediate state ! - vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr - vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr - by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr - bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr - vb = sm * bx + vy * by + vz * bz + vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr + vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr + by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr + bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr + vb = sm * bx + 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) = bx - ui(iby) = by - ui(ibz) = bz - ui(ibp) = ur(ibp,i) - ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm + 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 ! 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(:) + f(:,i) = sr * ui(:) - wr(:) - else ! sr* > 0 + 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 ! @@ -2742,13 +2821,13 @@ module schemes ! primitive variables for the inner right intermediate state ! - dv = slmm * dnl - car * dnr - vy = (wl(imy) - wcr(imy)) / dv - vz = (wl(imz) - wcr(imz)) / dv - dv = sl - smr - by = (wl(iby) - wcr(iby)) / dv - bz = (wl(ibz) - wcr(ibz)) / dv - vb = sm * bx + vy * by + vz * bz + dv = slmm * dnl - car * dnr + vy = (wl(imy) - wcr(imy)) / dv + vz = (wl(imz) - wcr(imz)) / dv + dv = sl - smr + by = (wl(iby) - wcr(iby)) / dv + bz = (wl(ibz) - wcr(ibz)) / dv + vb = sm * bx + vy * by + vz * bz ! conservative variables for the inner left intermediate state ! @@ -2765,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 ! @@ -2775,137 +2854,33 @@ module schemes ! vector of the right-going Alfvén wave ! - wcl(:) = (sm - smr) * ui(:) + wcr(:) + wcl(:) = (sm - smr) * ui(:) + wcr(:) ! calculate the average flux over the left inner intermediate state ! - f(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm + f(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm end if ! sm > 0 - end if ! sr* > 0 + else ! bx = 0 - else ! sl* < sl & sr* > sr +! 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 ! so far we revert to HLL flux in the case of degeneracies -! - 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 - -! 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 + end if ! sl* < sl & sr* > sr - 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