From 8843b6f45f5991e3047378a8b19cacaf590d8464 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 15 Dec 2021 16:45:33 -0300 Subject: [PATCH] EQUATIONS: Rewrite a bit esystem_roe_mhd_adi(). Signed-off-by: Grzegorz Kowal --- sources/equations.F90 | 299 +++++++++++++++++++++--------------------- 1 file changed, 148 insertions(+), 151 deletions(-) diff --git a/sources/equations.F90 b/sources/equations.F90 index cf94e41..14abcdf 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -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,34 +3306,33 @@ 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 ! - c(1) = q(ivx) - cf - c(2) = q(ivx) - ca - c(3) = q(ivx) - cs - c(4) = q(ivx) - c(5) = q(ivx) - c(6) = q(ivx) + cs - c(7) = q(ivx) + ca - c(8) = q(ivx) + cf - c(9) = c(8) + c(1) = q(ivx) - cf + c(2) = q(ivx) - ca + c(3) = q(ivx) - cs + c(4) = q(ivx) + c(5) = q(ivx) + c(6) = q(ivx) + cs + c(7) = q(ivx) + ca + c(8) = q(ivx) + cf + c(9) = c(8) ! eigenvectors only for the case of waves propagating in both direction ! @@ -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)