EQUATIONS: Rewrite a bit esystem_roe_mhd_adi().
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
e0f11b9b21
commit
8843b6f45f
@ -3281,16 +3281,16 @@ module equations
|
||||
real(kind=8), dimension(9,9), save :: lvec, rvec
|
||||
!$omp threadprivate(first, gammam2, lvec, rvec)
|
||||
|
||||
real(kind=8) :: di, vsq, btsq, bt_starsq, casq, hp, twid_asq
|
||||
real(kind=8) :: ct2, tsum, tdif, cf2_cs2, cfsq, cf, cssq, cs, ca, bt
|
||||
real(kind=8) :: bt_star, bet2, bet3, bet2_star, bet3_star, bet_starsq, vbet
|
||||
real(kind=8) :: alpha_f, alpha_s, af_prime, as_prime
|
||||
real(kind=8) :: sqrtd, isqrtd, s, twid_a, qf, qs, afpbb, aspbb
|
||||
real(kind=8) :: qa, qb, qc, qd
|
||||
real(kind=8) :: norm, cff, css, af, as, afpb, aspb
|
||||
real(kind=8) :: q2_star, q3_star, vqstr
|
||||
real(kind=8) :: ca, ca2, cf, cf2, cs, cs2, cf2_cs2, ct2
|
||||
real(kind=8) :: v2, br2, br, br2s, brs, hp
|
||||
real(kind=8) :: tsum, tdif, twid_a2, twid_a
|
||||
real(kind=8) :: bty, btz, btys, btzs, bt2s, vbet
|
||||
real(kind=8) :: alf, als, af_prime, as_prime
|
||||
real(kind=8) :: sqrtd, sgn, qf, qs, afpbb, aspbb
|
||||
real(kind=8) :: qa, qb, qc, qd, q2s, q3s
|
||||
real(kind=8) :: norm, cff, css, af, as, afpb, aspb, vqstr
|
||||
|
||||
real(kind=8), parameter :: eps = epsilon(di)
|
||||
real(kind=8), parameter :: eps = epsilon(1.0d+00)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -3306,22 +3306,21 @@ module equations
|
||||
|
||||
! coefficients
|
||||
!
|
||||
di = 1.0d+00 / q(idn)
|
||||
casq = q(ibx) * q(ibx) * di
|
||||
ca = sqrt(casq)
|
||||
vsq = q(ivx) * q(ivx) + q(ivy) * q(ivy) + q(ivz) * q(ivz)
|
||||
btsq = q(iby) * q(iby) + q(ibz) * q(ibz)
|
||||
bt_starsq = (gammam1 - gammam2 * y) * btsq
|
||||
hp = q(ien) - (casq + btsq * di)
|
||||
twid_asq = max(eps, (gammam1 * (hp - 0.5d+00 * vsq) - gammam2 * x))
|
||||
ct2 = bt_starsq * di
|
||||
tsum = casq + ct2 + twid_asq
|
||||
tdif = casq + ct2 - twid_asq
|
||||
cf2_cs2 = sqrt((tdif * tdif + 4.0d+00 * twid_asq * ct2))
|
||||
cfsq = 0.5d+00 * (tsum + cf2_cs2)
|
||||
cf = sqrt(cfsq)
|
||||
cssq = twid_asq * casq / cfsq
|
||||
cs = sqrt(cssq)
|
||||
ca2 = q(ibx) * q(ibx) / q(idn)
|
||||
ca = sqrt(ca2)
|
||||
v2 = sum(q(ivx:ivz)**2)
|
||||
br2 = sum(q(iby:ibz)**2)
|
||||
br2s = (gammam1 - gammam2 * y) * br2
|
||||
hp = q(ien) - (ca2 + br2 / q(idn))
|
||||
twid_a2 = max(eps, (gammam1 * (hp - 0.5d+00 * v2) - gammam2 * x))
|
||||
ct2 = br2s / q(idn)
|
||||
tsum = ca2 + ct2 + twid_a2
|
||||
tdif = ca2 + ct2 - twid_a2
|
||||
cf2_cs2 = sqrt((tdif * tdif + 4.0d+00 * twid_a2 * ct2))
|
||||
cf2 = 0.5d+00 * (tsum + cf2_cs2)
|
||||
cf = sqrt(cf2)
|
||||
cs2 = twid_a2 * ca2 / cf2
|
||||
cs = sqrt(cs2)
|
||||
|
||||
! eigenvalues
|
||||
!
|
||||
@ -3341,171 +3340,169 @@ module equations
|
||||
|
||||
! remaining coefficients
|
||||
!
|
||||
bt = sqrt(btsq)
|
||||
bt_star = sqrt(bt_starsq)
|
||||
if (abs(bt) > 0.0d+00) then
|
||||
bet2 = q(iby) / bt
|
||||
bet3 = q(ibz) / bt
|
||||
br = sqrt(br2)
|
||||
brs = sqrt(br2s)
|
||||
if (abs(br) > 0.0d+00) then
|
||||
bty = q(iby) / br
|
||||
btz = q(ibz) / br
|
||||
else
|
||||
bet2 = 1.0d+00
|
||||
bet3 = 0.0d+00
|
||||
bty = 1.0d+00
|
||||
btz = 0.0d+00
|
||||
end if
|
||||
bet2_star = bet2 / sqrt(gammam1 - gammam2 * y)
|
||||
bet3_star = bet3 / sqrt(gammam1 - gammam2 * y)
|
||||
bet_starsq = bet2_star * bet2_star + bet3_star * bet3_star
|
||||
vbet = q(ivy) * bet2_star + q(ivz) * bet3_star
|
||||
btys = bty / sqrt(gammam1 - gammam2 * y)
|
||||
btzs = btz / sqrt(gammam1 - gammam2 * y)
|
||||
bt2s = btys * btys + btzs * btzs
|
||||
vbet = q(ivy) * btys + q(ivz) * btzs
|
||||
|
||||
if ( .not. abs(cfsq - cssq) > 0.0d+00 ) then
|
||||
alpha_f = 1.0d+00
|
||||
alpha_s = 0.0d+00
|
||||
else if ( (twid_asq - cssq) <= 0.0d+00 ) then
|
||||
alpha_f = 0.0d+00
|
||||
alpha_s = 1.0d+00
|
||||
else if ( (cfsq - twid_asq) <= 0.0d+00 ) then
|
||||
alpha_f = 1.0d+00
|
||||
alpha_s = 0.0d+00
|
||||
if ( .not. abs(cf2 - cs2) > 0.0d+00 ) then
|
||||
alf = 1.0d+00
|
||||
als = 0.0d+00
|
||||
else if ( (twid_a2 - cs2) <= 0.0d+00 ) then
|
||||
alf = 0.0d+00
|
||||
als = 1.0d+00
|
||||
else if ( (cf2 - twid_a2) <= 0.0d+00 ) then
|
||||
alf = 1.0d+00
|
||||
als = 0.0d+00
|
||||
else
|
||||
alpha_f = sqrt((twid_asq - cssq) / (cfsq - cssq))
|
||||
alpha_s = sqrt((cfsq - twid_asq) / (cfsq - cssq))
|
||||
alf = sqrt((twid_a2 - cs2) / (cf2 - cs2))
|
||||
als = sqrt((cf2 - twid_a2) / (cf2 - cs2))
|
||||
end if
|
||||
|
||||
sqrtd = sqrt(q(idn))
|
||||
isqrtd = 1.0d+00 / sqrtd
|
||||
s = sign(1.0d+00, q(ibx))
|
||||
twid_a = sqrt(twid_asq)
|
||||
qf = cf * alpha_f * s
|
||||
qs = cs * alpha_s * s
|
||||
af_prime = twid_a * alpha_f * isqrtd
|
||||
as_prime = twid_a * alpha_s * isqrtd
|
||||
afpbb = af_prime * bt_star * bet_starsq
|
||||
aspbb = as_prime * bt_star * bet_starsq
|
||||
sgn = sign(1.0d+00, q(ibx))
|
||||
twid_a = sqrt(twid_a2)
|
||||
qf = cf * alf * sgn
|
||||
qs = cs * als * sgn
|
||||
af_prime = twid_a * alf / sqrtd
|
||||
as_prime = twid_a * als / sqrtd
|
||||
afpbb = af_prime * brs * bt2s
|
||||
aspbb = as_prime * brs * bt2s
|
||||
|
||||
! update the varying elements of the matrix of right eigenvectors
|
||||
! === update the varying elements of the right eigenvectors matrix
|
||||
!
|
||||
rvec(1,idn) = alpha_f
|
||||
rvec(3,idn) = alpha_s
|
||||
rvec(6,idn) = alpha_s
|
||||
rvec(8,idn) = alpha_f
|
||||
rvec(1,idn) = alf
|
||||
rvec(3,idn) = als
|
||||
rvec(6,idn) = als
|
||||
rvec(8,idn) = alf
|
||||
|
||||
rvec(1,ivx) = alpha_f * c(1)
|
||||
rvec(3,ivx) = alpha_s * c(3)
|
||||
rvec(1,ivx) = alf * c(1)
|
||||
rvec(3,ivx) = als * c(3)
|
||||
rvec(4,ivx) = q(ivx)
|
||||
rvec(6,ivx) = alpha_s * c(6)
|
||||
rvec(8,ivx) = alpha_f * c(8)
|
||||
rvec(6,ivx) = als * c(6)
|
||||
rvec(8,ivx) = alf * c(8)
|
||||
|
||||
qa = alpha_f * q(ivy)
|
||||
qb = alpha_s * q(ivy)
|
||||
qc = qs * bet2_star
|
||||
qd = qf * bet2_star
|
||||
qa = alf * q(ivy)
|
||||
qb = als * q(ivy)
|
||||
qc = qs * btys
|
||||
qd = qf * btys
|
||||
|
||||
rvec(1,ivy) = qa + qc
|
||||
rvec(2,ivy) = - bet3
|
||||
rvec(2,ivy) = - btz
|
||||
rvec(3,ivy) = qb - qd
|
||||
rvec(4,ivy) = q(ivy)
|
||||
rvec(6,ivy) = qb + qd
|
||||
rvec(7,ivy) = bet3
|
||||
rvec(7,ivy) = btz
|
||||
rvec(8,ivy) = qa - qc
|
||||
|
||||
qa = alpha_f * q(ivz)
|
||||
qb = alpha_s * q(ivz)
|
||||
qc = qs * bet3_star
|
||||
qd = qf * bet3_star
|
||||
qa = alf * q(ivz)
|
||||
qb = als * q(ivz)
|
||||
qc = qs * btzs
|
||||
qd = qf * btzs
|
||||
|
||||
rvec(1,ivz) = qa + qc
|
||||
rvec(2,ivz) = bet2
|
||||
rvec(2,ivz) = bty
|
||||
rvec(3,ivz) = qb - qd
|
||||
rvec(4,ivz) = q(ivz)
|
||||
rvec(6,ivz) = qb + qd
|
||||
rvec(7,ivz) = - bet2
|
||||
rvec(7,ivz) = - bty
|
||||
rvec(8,ivz) = qa - qc
|
||||
|
||||
rvec(1,ipr) = alpha_f * (hp - q(ivx) * cf) + qs * vbet + aspbb
|
||||
rvec(2,ipr) = -(q(ivy) * bet3 - q(ivz) * bet2)
|
||||
rvec(3,ipr) = alpha_s * (hp - q(ivx) * cs) - qf * vbet - afpbb
|
||||
rvec(4,ipr) = 0.5d+00 * vsq + gammam2 * x / gammam1
|
||||
rvec(6,ipr) = alpha_s * (hp + q(ivx) * cs) + qf * vbet - afpbb
|
||||
rvec(1,ipr) = alf * (hp - q(ivx) * cf) + qs * vbet + aspbb
|
||||
rvec(2,ipr) = -(q(ivy) * btz - q(ivz) * bty)
|
||||
rvec(3,ipr) = als * (hp - q(ivx) * cs) - qf * vbet - afpbb
|
||||
rvec(4,ipr) = 0.5d+00 * v2 + gammam2 * x / gammam1
|
||||
rvec(6,ipr) = als * (hp + q(ivx) * cs) + qf * vbet - afpbb
|
||||
rvec(7,ipr) = - rvec(2,ipr)
|
||||
rvec(8,ipr) = alpha_f * (hp + q(ivx) * cf) - qs * vbet + aspbb
|
||||
rvec(8,ipr) = alf * (hp + q(ivx) * cf) - qs * vbet + aspbb
|
||||
|
||||
rvec(1,iby) = as_prime * bet2_star
|
||||
rvec(2,iby) = - bet3 * s * isqrtd
|
||||
rvec(3,iby) = - af_prime * bet2_star
|
||||
rvec(1,iby) = as_prime * btys
|
||||
rvec(2,iby) = - btz * sgn / sqrtd
|
||||
rvec(3,iby) = - af_prime * btys
|
||||
rvec(6,iby) = rvec(3,iby)
|
||||
rvec(7,iby) = rvec(2,iby)
|
||||
rvec(8,iby) = rvec(1,iby)
|
||||
|
||||
rvec(1,ibz) = as_prime * bet3_star
|
||||
rvec(2,ibz) = bet2 * s * isqrtd
|
||||
rvec(3,ibz) = - af_prime * bet3_star
|
||||
rvec(1,ibz) = as_prime * btzs
|
||||
rvec(2,ibz) = bty * sgn / sqrtd
|
||||
rvec(3,ibz) = - af_prime * btzs
|
||||
rvec(6,ibz) = rvec(3,ibz)
|
||||
rvec(7,ibz) = rvec(2,ibz)
|
||||
rvec(8,ibz) = rvec(1,ibz)
|
||||
|
||||
! update the varying elements of the matrix of left eigenvectors
|
||||
! === update the varying elements of the left eigenvectors matrix
|
||||
!
|
||||
norm = 0.5d+00 / twid_asq
|
||||
cff = norm * alpha_f * cf
|
||||
css = norm * alpha_s * cs
|
||||
qf = qf * norm
|
||||
qs = qs * norm
|
||||
af = norm * af_prime * q(idn)
|
||||
as = norm * as_prime * q(idn)
|
||||
afpb = norm * af_prime * bt_star
|
||||
aspb = norm * as_prime * bt_star
|
||||
norm = 2.0d+00 * twid_a2
|
||||
cff = alf * cf / norm
|
||||
css = als * cs / norm
|
||||
qf = qf / norm
|
||||
qs = qs / norm
|
||||
af = af_prime * q(idn) / norm
|
||||
as = as_prime * q(idn) / norm
|
||||
afpb = af_prime * brs / norm
|
||||
aspb = as_prime * brs / norm
|
||||
|
||||
norm = norm * gammam1
|
||||
alpha_f = alpha_f * norm
|
||||
alpha_s = alpha_s * norm
|
||||
q2_star = bet2_star / bet_starsq
|
||||
q3_star = bet3_star / bet_starsq
|
||||
vqstr = (q(ivy) * q2_star + q(ivz) * q3_star)
|
||||
norm = 2.0d+00 * norm
|
||||
norm = norm / gammam1
|
||||
alf = alf / norm
|
||||
als = als / norm
|
||||
q2s = btys / bt2s
|
||||
q3s = btzs / bt2s
|
||||
vqstr = (q(ivy) * q2s + q(ivz) * q3s)
|
||||
|
||||
! left-going fast wave
|
||||
!
|
||||
lvec(idn,1) = alpha_f * (vsq - hp) + cff * (cf + q(ivx)) - qs * vqstr - aspb
|
||||
lvec(ivx,1) = - alpha_f * q(ivx) - cff
|
||||
lvec(ivy,1) = - alpha_f * q(ivy) + qs * q2_star
|
||||
lvec(ivz,1) = - alpha_f * q(ivz) + qs * q3_star
|
||||
lvec(ipr,1) = alpha_f
|
||||
lvec(iby,1) = as * q2_star - alpha_f * q(iby)
|
||||
lvec(ibz,1) = as * q3_star - alpha_f * q(ibz)
|
||||
lvec(idn,1) = alf * (v2 - hp) + cff * (cf + q(ivx)) - qs * vqstr - aspb
|
||||
lvec(ipr,1) = alf
|
||||
lvec(ivx,1) = - alf * q(ivx) - cff
|
||||
lvec(ivy,1) = - alf * q(ivy) + qs * q2s
|
||||
lvec(ivz,1) = - alf * q(ivz) + qs * q3s
|
||||
lvec(iby,1) = as * q2s - alf * q(iby)
|
||||
lvec(ibz,1) = as * q3s - alf * q(ibz)
|
||||
|
||||
! left-going Alfvèn wave
|
||||
!
|
||||
lvec(idn,2) = 0.5d+00 * (q(ivy) * bet3 - q(ivz) * bet2)
|
||||
lvec(ivy,2) = - 0.5d+00 * bet3
|
||||
lvec(ivz,2) = 0.5d+00 * bet2
|
||||
lvec(iby,2) = - 0.5d+00 * sqrtd * bet3 * s
|
||||
lvec(ibz,2) = 0.5d+00 * sqrtd * bet2 * s
|
||||
lvec(idn,2) = 5.0d-01 * (q(ivy) * btz - q(ivz) * bty)
|
||||
lvec(ivy,2) = - 5.0d-01 * btz
|
||||
lvec(ivz,2) = 5.0d-01 * bty
|
||||
lvec(iby,2) = - 5.0d-01 * sqrtd * btz * sgn
|
||||
lvec(ibz,2) = 5.0d-01 * sqrtd * bty * sgn
|
||||
|
||||
! left-going slow wave
|
||||
!
|
||||
lvec(idn,3) = alpha_s * (vsq - hp) + css * (cs + q(ivx)) + qf * vqstr + afpb
|
||||
lvec(ivx,3) = - alpha_s * q(ivx) - css
|
||||
lvec(ivy,3) = - alpha_s * q(ivy) - qf * q2_star
|
||||
lvec(ivz,3) = - alpha_s * q(ivz) - qf * q3_star
|
||||
lvec(ipr,3) = alpha_s
|
||||
lvec(iby,3) = - af * q2_star - alpha_s * q(iby)
|
||||
lvec(ibz,3) = - af * q3_star - alpha_s * q(ibz)
|
||||
lvec(idn,3) = als * (v2 - hp) + css * (cs + q(ivx)) + qf * vqstr + afpb
|
||||
lvec(ipr,3) = als
|
||||
lvec(ivx,3) = - als * q(ivx) - css
|
||||
lvec(ivy,3) = - als * q(ivy) - qf * q2s
|
||||
lvec(ivz,3) = - als * q(ivz) - qf * q3s
|
||||
lvec(iby,3) = - af * q2s - als * q(iby)
|
||||
lvec(ibz,3) = - af * q3s - als * q(ibz)
|
||||
|
||||
! entropy wave
|
||||
!
|
||||
lvec(idn,4) = 1.0d+00 - norm * (0.5d+00 * vsq - gammam2 * x / gammam1)
|
||||
lvec(ivx,4) = norm * q(ivx)
|
||||
lvec(ivy,4) = norm * q(ivy)
|
||||
lvec(ivz,4) = norm * q(ivz)
|
||||
lvec(ipr,4) = - norm
|
||||
lvec(iby,4) = norm * q(iby)
|
||||
lvec(ibz,4) = norm * q(ibz)
|
||||
lvec(idn,4) = 1.0d+00 - (5.0d-01 * v2 - gammam2 * x / gammam1) / twid_a2
|
||||
lvec(ipr,4) = - 1.0d+00 / twid_a2
|
||||
lvec(ivx,4) = q(ivx) / twid_a2
|
||||
lvec(ivy,4) = q(ivy) / twid_a2
|
||||
lvec(ivz,4) = q(ivz) / twid_a2
|
||||
lvec(iby,4) = q(iby) / twid_a2
|
||||
lvec(ibz,4) = q(ibz) / twid_a2
|
||||
|
||||
! right-going slow wave
|
||||
!
|
||||
lvec(idn,6) = alpha_s * (vsq - hp) + css * (cs - q(ivx)) - qf * vqstr + afpb
|
||||
lvec(ivx,6) = - alpha_s * q(ivx) + css
|
||||
lvec(ivy,6) = - alpha_s * q(ivy) + qf * q2_star
|
||||
lvec(ivz,6) = - alpha_s * q(ivz) + qf * q3_star
|
||||
lvec(ipr,6) = alpha_s
|
||||
lvec(idn,6) = als * (v2 - hp) + css * (cs - q(ivx)) - qf * vqstr + afpb
|
||||
lvec(ipr,6) = als
|
||||
lvec(ivx,6) = - als * q(ivx) + css
|
||||
lvec(ivy,6) = - als * q(ivy) + qf * q2s
|
||||
lvec(ivz,6) = - als * q(ivz) + qf * q3s
|
||||
lvec(iby,6) = lvec(iby,3)
|
||||
lvec(ibz,6) = lvec(ibz,3)
|
||||
|
||||
@ -3519,11 +3516,11 @@ module equations
|
||||
|
||||
! right-going fast wave
|
||||
!
|
||||
lvec(idn,8) = alpha_f * (vsq - hp) + cff * (cf - q(ivx)) + qs * vqstr - aspb
|
||||
lvec(ivx,8) = - alpha_f * q(ivx) + cff
|
||||
lvec(ivy,8) = - alpha_f * q(ivy) - qs * q2_star
|
||||
lvec(ivz,8) = - alpha_f * q(ivz) - qs * q3_star
|
||||
lvec(ipr,8) = alpha_f
|
||||
lvec(idn,8) = alf * (v2 - hp) + cff * (cf - q(ivx)) + qs * vqstr - aspb
|
||||
lvec(ipr,8) = alf
|
||||
lvec(ivx,8) = - alf * q(ivx) + cff
|
||||
lvec(ivy,8) = - alf * q(ivy) - qs * q2s
|
||||
lvec(ivz,8) = - alf * q(ivz) - qs * q3s
|
||||
lvec(iby,8) = lvec(iby,1)
|
||||
lvec(ibz,8) = lvec(ibz,1)
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user