From e339cf15e8881b344ee99cac3cf9ccfcb109602d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 30 Jul 2024 16:32:16 -0300 Subject: [PATCH 1/2] SCHEMES: Add symmetry preserving optimizations to isothermal KEPES Riemann solver - Update riemann_mhd_iso_kepes() subroutine to include missing symmetry preserving optimizations. - Ensure consistency with the adiabatic version by adding and utilizing the ch variable. - Correct the initialization and usage of transformation matrices (rm and tm). - Adjust calculations to ensure symmetric properties, particularly in the handling of variables cs, ca, and cglm. - Include appropriate sign functions for consistent behavior. - Make changes to the calculation of intermediate variables to maintain symmetry and consistency. - Adjusted flux calculation to correctly apply symmetry preserving terms. Signed-off-by: Grzegorz Kowal --- sources/schemes.F90 | 44 ++++++++++++++++++++++++-------------------- 1 file changed, 24 insertions(+), 20 deletions(-) diff --git a/sources/schemes.F90 b/sources/schemes.F90 index eca68ed..6b215a7 100644 --- a/sources/schemes.F90 +++ b/sources/schemes.F90 @@ -2155,7 +2155,7 @@ module schemes subroutine riemann_mhd_iso_kepes(ql, qr, fi) use coordinates, only : nn => bcells - use equations , only : nf, csnd, csnd2, ch => cglm + use equations , only : nf, csnd, csnd2, cglm use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, & ibx, iby, ibz, ibp @@ -2167,7 +2167,7 @@ module schemes integer :: i real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa real(kind=8) :: dnl, v2l, v2r, b2l, b2r, bp2, das, daf, csq - real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3 + real(kind=8) :: ch, ca, cs, cf, ca2, cs2, cf2, x2, x3 real(kind=8) :: alf, als, f1, f2, f3, f4 real(kind=8), dimension(nf) :: v, lm, tm @@ -2181,9 +2181,7 @@ module schemes if (first) then rm(:,:) = 0.0d+00 rm(5,4) = 1.0d+00 - rm(8,4) = 1.0d+00 rm(5,5) = 1.0d+00 - rm(8,5) = - 1.0d+00 first = .false. end if @@ -2219,6 +2217,7 @@ module schemes v(iby) = qr(iby,i) - ql(iby,i) v(ibz) = qr(ibz,i) - ql(ibz,i) v(ibp) = qr(ibp,i) - ql(ibp,i) + v( : ) = v( : ) / csnd2 bp2 = bya * bya + bza * bza if (bp2 > 0.0d+00) then @@ -2251,8 +2250,9 @@ module schemes end if cf = sign(sqrt(cf2), vxa) - ca = sqrt(ca2) cs = sign(sqrt(cs2), vxa) + ca = sign(sqrt(ca2), vxa) + ch = sign( cglm, vxa) lm(1) = vxa + cf lm(2) = vxa + ca @@ -2263,10 +2263,10 @@ module schemes lm(7) = vxa - ca lm(8) = vxa - cf - tm(1) = 5.0d-01 / (dnl * csnd2) - tm(2) = 5.0d-01 / dnl**2 + tm(1) = 5.0d-01 / dnl + tm(2) = 5.0d-01 * csnd2 / dnl**2 tm(3) = tm(1) - tm(4) = 5.0d-01 + tm(4) = 5.0d-01 * csnd2 tm(5) = tm(4) tm(6) = tm(1) tm(7) = tm(2) @@ -2276,24 +2276,24 @@ module schemes daf = dnl * alf csq = csnd * sqrt(dnl) - f1 = daf - f2 = sign(das, bxa) * cs - f3 = als * csq + f1 = daf + f2 = sign(das, bxa) * cs + f3 = als * csq rm(1,1) = f1 - rm(2,1) = f1 * (vxa + cf) + rm(2,1) = f1 * lm(1) rm(3,1) = f1 * vya - f2 * x2 rm(4,1) = f1 * vza - f2 * x3 rm(6,1) = f3 * x2 rm(7,1) = f3 * x3 rm(1,8) = f1 - rm(2,8) = f1 * (vxa - cf) + rm(2,8) = f1 * lm(8) rm(3,8) = f1 * vya + f2 * x2 rm(4,8) = f1 * vza + f2 * x3 rm(6,8) = rm(6,1) rm(7,8) = rm(7,1) - f1 = dnl * sqrt(dna) + f1 = sign(dnl * sqrt(dna), vxa) rm(3,2) = f1 * x3 rm(4,2) = - f1 * x2 rm(6,2) = - dnl * x3 @@ -2304,31 +2304,35 @@ module schemes rm(6,7) = rm(6,2) rm(7,7) = rm(7,2) - f1 = das - f2 = sign(daf, bxa) * cf + f1 = das + f2 = sign(daf, bxa) * cf f3 = - alf * csq rm(1,3) = f1 - rm(2,3) = f1 * (vxa + cs) + rm(2,3) = f1 * lm(3) rm(3,3) = f1 * vya + f2 * x2 rm(4,3) = f1 * vza + f2 * x3 rm(6,3) = f3 * x2 rm(7,3) = f3 * x3 rm(1,6) = f1 - rm(2,6) = f1 * (vxa - cs) + rm(2,6) = f1 * lm(6) rm(3,6) = f1 * vya - f2 * x2 rm(4,6) = f1 * vza - f2 * x3 rm(6,6) = rm(6,3) rm(7,6) = rm(7,3) + f1 = sign(1.0d+00, vxa) + rm(8,4) = f1 + rm(8,5) = - f1 + fi(idn,i) = dnl * vxa fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta fi(imy,i) = fi(idn,i) * vya - bxa * bya fi(imz,i) = fi(idn,i) * vza - bxa * bza - fi(ibx,i) = ch * bpa + fi(ibx,i) = cglm * bpa fi(iby,i) = vxa * bya - bxa * vya fi(ibz,i) = vxa * bza - bxa * vza - fi(ibp,i) = ch * bxa + fi(ibp,i) = cglm * bxa fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm)) From a45a380d0c1c2bdf45d8fd92f44c00b586a032db Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 30 Jul 2024 16:46:29 -0300 Subject: [PATCH 2/2] SOURCES: Consider -B.div(B) terms for isothermal KEPES too. The -B.div(B) source term in the momentum equation was dropped for the isothermal MHD case. Apparently, its lack causes some numerical instabilities related to the accumulation of the divergence of B. Therefore, take it into account for the isothermal case too. Signed-off-by: Grzegorz Kowal --- sources/sources.F90 | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/sources/sources.F90 b/sources/sources.F90 index 201a456..be5da5e 100644 --- a/sources/sources.F90 +++ b/sources/sources.F90 @@ -571,21 +571,16 @@ module sources ! call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:)) -! update the momentum component increments (only in the adiabatic case for -! KEPES), i.e. +! update the momentum component increments, i.e. ! d/dt (ρv) + ∇.F = - (∇.B)B ! - if (glm_type < 3 .or. ien > 0) then - - pdata%du(imx,:,:,:) = pdata%du(imx,:,:,:) & + pdata%du(imx,:,:,:) = pdata%du(imx,:,:,:) & - db(:,:,:) * pdata%q(ibx,:,:,:) - pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) & + pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) & - db(:,:,:) * pdata%q(iby,:,:,:) - pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) & + pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) & - db(:,:,:) * pdata%q(ibz,:,:,:) - end if ! ien > 0 - ! add the HEGLM-MHD and KEPES source terms ! if (glm_type > 1) then