SCHEMES: Rewrite riemann_mhd_iso_roe().

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-01-26 13:10:55 -03:00
parent 5c96f3fd5a
commit 2fdb52d972

View File

@ -2040,117 +2040,230 @@ module schemes
subroutine riemann_mhd_iso_roe(ql, qr, fi) subroutine riemann_mhd_iso_roe(ql, qr, fi)
use coordinates, only : nn => bcells use coordinates, only : nn => bcells
use equations , only : nf, idn, ivx, ivz, imx, imy, imz, ibx, iby, ibz, ibp use equations , only : nf, csnd2
use equations , only : prim2cons, fluxspeed, eigensystem_roe use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
ibx, iby, ibz, ibp
use equations , only : prim2cons, fluxspeed
implicit none implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i, p integer :: i
real(kind=8) :: sdl, sdr, sds real(kind=8) :: sdl, sdr, sds, xx, yy, bty, btz, br, br2
real(kind=8) :: xx, yy real(kind=8) :: cc, ca, cs, cf, cc2, ca2, cs2, cf2, alf, als
real(kind=8) :: vqstr, sqrtd, qs, qf, norm, css, cff
real(kind=8) :: as_prime, af_prime, as, af, aspb, afpb
real(kind=8) :: f1, f2, f3
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr real(kind=8), dimension(nf ) :: qi
real(kind=8), dimension(nf ) :: qi, ci, al real(kind=8), dimension(6 ) :: lm, al, df
real(kind=8), dimension(nf,nf) :: li, ri
logical , save :: first = .true.
integer , dimension(6) , save :: ivs
real(kind=8), dimension(6,6), save :: rvec, lvec
!$omp threadprivate(first, ivs, rvec, lvec)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
if (first) then
rvec(:,:) = 0.0d+00
lvec(:,:) = 0.0d+00
ivs = [ idn, imx, imy, imz, iby, ibz ]
first = .false.
end if
call prim2cons(ql, ul) call prim2cons(ql, ul)
call prim2cons(qr, ur) call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl) call fluxspeed(ql, ul, fl)
call fluxspeed(qr, ur, fr, cr) call fluxspeed(qr, ur, fr)
do i = 1, nn do i = 1, nn
! calculate variables for the Roe intermediate state
!
sdl = sqrt(ql(idn,i)) sdl = sqrt(ql(idn,i))
sdr = sqrt(qr(idn,i)) sdr = sqrt(qr(idn,i))
sds = sdl + sdr sds = sdl + sdr
! prepare the Roe intermediate state vector (eq. 11.60 in [2])
!
qi(idn) = sdl * sdr qi(idn) = sdl * sdr
qi(ivx:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds
qi(ibx) = ql(ibx,i) qi(ibx) = ql(ibx,i)
qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby:ibz,i) / sds
qi(ibp) = ql(ibp,i) qi(ibp) = ql(ibp,i)
! prepare coefficients f1 = ql(iby,i) - qr(iby,i)
! f2 = ql(ibz,i) - qr(ibz,i)
xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2 & xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds)
+ (ql(ibz,i) - qr(ibz,i))**2) / sds**2 yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn)
yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn)
! obtain eigenvalues and eigenvectors
!
call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:))
! return upwind fluxes
!
if (ci(1) >= 0.0d+00) then
fi(:,i) = fl(:,i)
else if (ci(nf) <= 0.0d+00) then
fi(:,i) = fr(:,i)
br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz)
br = sqrt(br2)
if (br2 > 0.0d+00) then
bty = qi(iby) / br
btz = qi(ibz) / br
else else
bty = 0.0d+00
! prepare fluxes which do not change across the states btz = 0.0d+00
!
fi(ibx,i) = fl(ibx,i)
fi(ibp,i) = fl(ibp,i)
! calculate wave amplitudes α = L.ΔU
!
al(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
end do
! calculate the flux average
!
fi(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i))
fi(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i))
fi(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i))
fi(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i))
fi(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i))
fi(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i))
! correct the flux by adding the characteristic wave contribution: (α|λ|K)
!
if (qi(ivx) >= 0.0d+00) then
do p = 1, nf
xx = - 0.5d+00 * al(p) * abs(ci(p))
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
end do
else
do p = nf, 1, -1
xx = - 0.5d+00 * al(p) * abs(ci(p))
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
end do
end if end if
end if ! sl < 0 < sr cc2 = csnd2 + xx
ca2 = qi(ibx) * qi(ibx)
f1 = (ca2 + br2) / qi(idn) + cc2
ca2 = ca2 / qi(idn)
f2 = f1 * f1
f3 = 4.0d+00 * cc2 * ca2
if (f2 > f3) then
f2 = sqrt(f2 - f3)
cf2 = 5.0d-01 * (f1 + f2)
cs2 = 5.0d-01 * (f1 - f2)
else
cf2 = 5.0d-01 * f1
cs2 = cf2
end if
if (cf2 > cs2) then
f2 = cf2 - cs2
if (cc2 > cs2) then
alf = sqrt((cc2 - cs2) / f2)
else
alf = 0.0d+00
end if
if (cf2 > cc2) then
als = sqrt((cf2 - cc2) / f2)
else
als = 0.0d+00
end if
else
alf = 1.0d+00
als = 1.0d+00
end if
cc = sqrt(cc2)
ca = sqrt(ca2)
cf = sign(sqrt(cf2), qi(ivx))
cs = sign(sqrt(cs2), qi(ivx))
lm(1) = qi(ivx) + cf
lm(2) = qi(ivx) + ca
lm(3) = qi(ivx) + cs
lm(4) = qi(ivx) - cs
lm(5) = qi(ivx) - ca
lm(6) = qi(ivx) - cf
sqrtd = sqrt(qi(idn))
f1 = sqrt(yy)
qf = cf * sign(alf, qi(ibx)) / f1
qs = cs * sign(als, qi(ibx)) / f1
f1 = cc / (f1 * sqrtd)
af_prime = alf * f1
as_prime = als * f1
rvec(1,1) = alf
rvec(1,2) = alf * lm(1)
rvec(1,3) = alf * qi(ivy) - qs * bty
rvec(1,4) = alf * qi(ivz) - qs * btz
rvec(1,5) = as_prime * bty
rvec(1,6) = as_prime * btz
rvec(6,1) = alf
rvec(6,2) = alf * lm(6)
rvec(6,3) = alf * qi(ivy) + qs * bty
rvec(6,4) = alf * qi(ivz) + qs * btz
rvec(6,5) = rvec(1,5)
rvec(6,6) = rvec(1,6)
rvec(3,1) = als
rvec(3,2) = als * lm(3)
rvec(3,3) = als * qi(ivy) + qf * bty
rvec(3,4) = als * qi(ivz) + qf * btz
rvec(3,5) = - af_prime * bty
rvec(3,6) = - af_prime * btz
rvec(4,1) = als
rvec(4,2) = als * lm(4)
rvec(4,3) = als * qi(ivy) - qf * bty
rvec(4,4) = als * qi(ivz) - qf * btz
rvec(4,5) = rvec(3,5)
rvec(4,6) = rvec(3,6)
f1 = sign(1.0d+00 / sqrtd, qi(ibx))
rvec(2,3) = btz
rvec(2,4) = - bty
rvec(2,5) = - btz * f1
rvec(2,6) = bty * f1
rvec(5,3) = - rvec(2,3)
rvec(5,4) = - rvec(2,4)
rvec(5,5) = rvec(2,5)
rvec(5,6) = rvec(2,6)
norm = 2.0d+00 * cc2
cff = alf * cf / norm
css = als * cs / norm
norm = norm / yy
qf = qf / norm
qs = qs / norm
f1 = qi(idn) / norm
af = af_prime * f1
as = as_prime * f1
f1 = br / norm
afpb = af_prime * f1
aspb = as_prime * f1
vqstr = (qi(ivy) * bty + qi(ivz) * btz)
f1 = qs * vqstr
lvec(1,1) = - cff * lm(6) + f1 - aspb
lvec(1,2) = cff
lvec(1,3) = - qs * bty
lvec(1,4) = - qs * btz
lvec(1,5) = as * bty
lvec(1,6) = as * btz
lvec(6,1) = cff * lm(1) - f1 - aspb
lvec(6,2) = - lvec(1,2)
lvec(6,3) = - lvec(1,3)
lvec(6,4) = - lvec(1,4)
lvec(6,5) = lvec(1,5)
lvec(6,6) = lvec(1,6)
f1 = qf * vqstr
lvec(3,1) = - css * lm(4) - f1 + afpb
lvec(3,2) = css
lvec(3,3) = qf * bty
lvec(3,4) = qf * btz
lvec(3,5) = - af * bty
lvec(3,6) = - af * btz
lvec(4,1) = css * lm(3) + f1 + afpb
lvec(4,2) = - lvec(3,2)
lvec(4,3) = - lvec(3,3)
lvec(4,4) = - lvec(3,4)
lvec(4,5) = lvec(3,5)
lvec(4,6) = lvec(3,6)
f1 = sign(5.0d-01 * sqrtd, qi(ibx))
lvec(2,1) = - 5.0d-01 * (qi(ivy) * btz - qi(ivz) * bty)
lvec(2,3) = 5.0d-01 * btz
lvec(2,4) = - 5.0d-01 * bty
lvec(2,5) = - f1 * btz
lvec(2,6) = f1 * bty
lvec(5,1) = - lvec(2,1)
lvec(5,3) = - lvec(2,3)
lvec(5,4) = - lvec(2,4)
lvec(5,5) = lvec(2,5)
lvec(5,6) = lvec(2,6)
al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i))
df = matmul(al, rvec)
fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:))
fi(ibx,i) = fl(ibx,i)
fi(ibp,i) = fl(ibp,i)
end do end do