gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

SCHEMES: Pass only flux variables to Riemann solvers.

Browse Source 
There is no need to interpolate or pass all variables to the Riemann
solvers. Just pass those related to flux calculation.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2019-10-03 21:48:04 -03:00
parent e26d9b23e4
commit 92b905cecd
1 changed files with 173 additions and 179 deletions

352
sources/schemes.F90

@ -685,7 +685,7 @@ module schemes
! include external procedures
!
use equations , only : nv, idn, ipr
use equations , only : nf, idn, ipr
use interpolations , only : interfaces
! local variables are not implicit by default
@ -705,9 +705,9 @@ module schemes
!
!-------------------------------------------------------------------------------
!
! iterate over all variables
! iterate over all flux variables
!
do p = 1, nv
do p = 1, nf
! determine if the variable is positive
!
@ -718,7 +718,7 @@ module schemes
call interfaces(positive, dx(1:NDIMS), q (p,:,:,:), &
qi(p,:,:,:,1:2,1:NDIMS))
end do ! p = 1, nv
end do ! p = 1, nf
!-------------------------------------------------------------------------------
!
@ -750,7 +750,7 @@ module schemes
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
! local variables are not implicit by default
@ -770,12 +770,12 @@ module schemes
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nv,nn,nn,nn,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,nn, 1,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,2) :: qi
real(kind=8), dimension(nv,nn) :: fi
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
@ -911,7 +911,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : ivx
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
@ -932,8 +932,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -1061,7 +1061,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz
use equations , only : prim2cons, fluxspeed, eigensystem_roe
use interpolations, only : order
@ -1079,13 +1079,12 @@ module schemes
!
integer :: n, p, i
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: xx
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: qi, ci, al
real(kind=8), dimension(nv,nv) :: li, ri
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!
!-------------------------------------------------------------------------------
!
@ -1141,7 +1140,7 @@ module schemes
f(:,i) = fl(:,i)
else if (ci(nv) <= 0.0d+00) then
else if (ci(nf) <= 0.0d+00) then
f(:,i) = fr(:,i)
@ -1149,26 +1148,24 @@ module schemes
! calculate wave amplitudes α = L.ΔU
!
al(1:nv) = 0.0d+00
do p = 1, nv
al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i))
al(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
end do
! calculate the flux average
!
f(1:nv,i) = 0.5d+00 * (fl(1:nv,i) + fr(1:nv,i))
f(:,i) = 0.5d+00 * (fl(:,i) + fr(:,i))
! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
if (qi(ivx) >= 0.0d+00) then
do p = 1, nv
xx = - 0.5d+00 * al(p) * abs(ci(p))
f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv)
do p = 1, nf
f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
end do
else
do p = nv, 1, - 1
xx = - 0.5d+00 * al(p) * abs(ci(p))
f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv)
do p = nf, 1, - 1
f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
end do
end if
@ -1233,7 +1230,7 @@ module schemes
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
! local variables are not implicit by default
@ -1253,12 +1250,12 @@ module schemes
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nv,nn,nn,nn,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,nn, 1,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,2) :: qi
real(kind=8), dimension(nv,nn) :: fi
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
@ -1400,7 +1397,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : ivx
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
@ -1421,8 +1418,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -1465,11 +1462,11 @@ module schemes
!
if (sl >= 0.0d+00) then
f(1:nv,i) = fl(1:nv,i)
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(1:nv,i) = fr(1:nv,i)
f(:,i) = fr(:,i)
else ! sl < 0 < sr
@ -1479,12 +1476,12 @@ module schemes
! calculate vectors of the left and right-going waves
!
wl(1:nv) = sl * ul(1:nv,i) - fl(1:nv,i)
wr(1:nv) = sr * ur(1:nv,i) - fr(1:nv,i)
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
f(1:nv,i) = (sl * wr(1:nv) - sr * wl(1:nv)) / srml
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
@ -1547,7 +1544,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
@ -1570,8 +1567,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr, ui
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, ui
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -1614,18 +1611,18 @@ module schemes
!
if (sl >= 0.0d+00) then
f(1:nv,i) = fl(1:nv,i)
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(1:nv,i) = fr(1:nv,i)
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate vectors of the left and right-going waves
!
wl(1:nv) = sl * ul(1:nv,i) - fl(1:nv,i)
wr(1:nv) = sr * ur(1:nv,i) - fr(1:nv,i)
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! the speed of contact discontinuity
!
@ -1654,7 +1651,7 @@ module schemes
! the left intermediate flux
!
f(1:nv,i) = sl * ui(1:nv) - wl(1:nv)
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
@ -1672,7 +1669,7 @@ module schemes
! the right intermediate flux
!
f(1:nv,i) = sr * ui(1:nv) - wr(1:nv)
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
@ -1751,7 +1748,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ipr, ien
use equations , only : prim2cons, fluxspeed, eigensystem_roe
use interpolations, only : order
@ -1769,13 +1766,12 @@ module schemes
!
integer :: n, p, i
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: xx
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: qi, ci, al
real(kind=8), dimension(nv,nv) :: li, ri
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!
!-------------------------------------------------------------------------------
!
@ -1833,7 +1829,7 @@ module schemes
f(:,i) = fl(:,i)
else if (ci(nv) <= 0.0d+00) then
else if (ci(nf) <= 0.0d+00) then
f(:,i) = fr(:,i)
@ -1841,26 +1837,24 @@ module schemes
! calculate wave amplitudes α = L.ΔU
!
al(1:nv) = 0.0d+00
do p = 1, nv
al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i))
al(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
end do
! calculate the flux average
!
f(1:nv,i) = 0.5d+00 * (fl(1:nv,i) + fr(1:nv,i))
f(:,i) = 0.5d+00 * (fl(:,i) + fr(:,i))
! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
if (qi(ivx) >= 0.0d+00) then
do p = 1, nv
xx = - 0.5d+00 * al(p) * abs(ci(p))
f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv)
do p = 1, nf
f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
end do
else
do p = nv, 1, - 1
xx = - 0.5d+00 * al(p) * abs(ci(p))
f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv)
do p = nf, 1, - 1
f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
end do
end if
@ -1925,7 +1919,7 @@ module schemes
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
use equations , only : ibx, iby, ibz, ibp
@ -1946,12 +1940,12 @@ module schemes
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nv,nn,nn,nn,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,nn, 1,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,2) :: qi
real(kind=8), dimension(nv,nn) :: fi
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
@ -2111,7 +2105,7 @@ module schemes
! include external variables and procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : ivx, ibx, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
@ -2134,8 +2128,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -2276,7 +2270,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
@ -2299,8 +2293,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -2686,7 +2680,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
@ -2709,8 +2703,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -3105,7 +3099,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp
use equations , only : imx, imy, imz
use equations , only : cmax
@ -3130,9 +3124,9 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: qi, ci, al
real(kind=8), dimension(nv,nv) :: li, ri
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!
!-------------------------------------------------------------------------------
!
@ -3214,7 +3208,7 @@ module schemes
f(:,i) = fl(:,i)
else if (ci(nv) <= 0.0d+00) then
else if (ci(nf) <= 0.0d+00) then
f(:,i) = fr(:,i)
@ -3227,9 +3221,9 @@ module schemes
! calculate wave amplitudes α = L.ΔU
!
al(1:nv) = 0.0d+00
do p = 1, nv
al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i))
al(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
end do
! calculate the flux average
@ -3244,7 +3238,7 @@ module schemes
! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
if (qi(ivx) >= 0.0d+00) then
do p = 1, nv
do p = 1, nf
xx = - 0.5d+00 * al(p) * abs(ci(p))
f(idn,i) = f(idn,i) + xx * ri(p,idn)
@ -3255,7 +3249,7 @@ module schemes
f(ibz,i) = f(ibz,i) + xx * ri(p,ibz)
end do
else
do p = nv, 1, -1
do p = nf, 1, -1
xx = - 0.5d+00 * al(p) * abs(ci(p))
f(idn,i) = f(idn,i) + xx * ri(p,idn)
@ -3328,7 +3322,7 @@ module schemes
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
use equations , only : ibx, iby, ibz, ibp
@ -3349,12 +3343,12 @@ module schemes
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nv,nn,nn,nn,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,nn, 1,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,2) :: qi
real(kind=8), dimension(nv,nn) :: fi
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
@ -3520,7 +3514,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : ivx, ibx, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
@ -3543,8 +3537,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -3691,7 +3685,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr
use equations , only : imx, imy, imz, ien
use equations , only : cmax
@ -3715,8 +3709,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr, ui
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, ui
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -4011,7 +4005,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr
use equations , only : imx, imy, imz, ien
use equations , only : cmax
@ -4037,8 +4031,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -4674,7 +4668,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
use equations , only : imx, imy, imz, ien
use equations , only : cmax
@ -4700,9 +4694,9 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: qi, ci, al
real(kind=8), dimension(nv,nv) :: li, ri
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!
!-------------------------------------------------------------------------------
!
@ -4791,7 +4785,7 @@ module schemes
f(:,i) = fl(:,i)
else if (ci(nv) <= 0.0d+00) then
else if (ci(nf) <= 0.0d+00) then
f(:,i) = fr(:,i)
@ -4804,9 +4798,9 @@ module schemes
! calculate wave amplitudes α = L.ΔU
!
al(1:nv) = 0.0d+00
do p = 1, nv
al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i))
al(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
end do
! calculate the flux average
@ -4822,7 +4816,7 @@ module schemes
! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
if (qi(ivx) >= 0.0d+00) then
do p = 1, nv
do p = 1, nf
xx = - 0.5d+00 * al(p) * abs(ci(p))
f(idn,i) = f(idn,i) + xx * ri(p,idn)
@ -4834,7 +4828,7 @@ module schemes
f(ibz,i) = f(ibz,i) + xx * ri(p,ibz)
end do
else
do p = nv, 1, -1
do p = nf, 1, -1
xx = - 0.5d+00 * al(p) * abs(ci(p))
f(idn,i) = f(idn,i) + xx * ri(p,idn)
@ -4908,7 +4902,7 @@ module schemes
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
! local variables are not implicit by default
@ -4929,14 +4923,14 @@ module schemes
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nv,nn,nn,nn) :: qq
real(kind=8), dimension(nv,nn,nn,nn,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn,nn) :: qq
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,nn, 1) :: qq
real(kind=8), dimension(nv,nn,nn, 1,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn, 1) :: qq
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,2) :: qi
real(kind=8), dimension(nv,nn) :: fi
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
@ -5138,7 +5132,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : ivx
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
@ -5159,8 +5153,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -5203,11 +5197,11 @@ module schemes
!
if (sl >= 0.0d+00) then
f(1:nv,i) = fl(1:nv,i)
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(1:nv,i) = fr(1:nv,i)
f(:,i) = fr(:,i)
else ! sl < 0 < sr
@ -5217,12 +5211,12 @@ module schemes
! calculate vectors of the left and right-going waves
!
wl(1:nv) = sl * ul(1:nv,i) - fl(1:nv,i)
wr(1:nv) = sr * ur(1:nv,i) - fr(1:nv,i)
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
f(1:nv,i) = (sl * wr(1:nv) - sr * wl(1:nv)) / srml
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
@ -5287,7 +5281,7 @@ module schemes
! include external procedures
!
use algebra , only : quadratic
use equations , only : nv
use equations , only : nf
use equations , only : ivx, idn, imx, imy, imz, ien
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
@ -5309,8 +5303,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: uh, us, fh, wl, wr
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: uh, us, fh, wl, wr
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
real(kind=8), dimension(3) :: a
real(kind=8), dimension(2) :: x
@ -5355,11 +5349,11 @@ module schemes
!
if (sl >= 0.0d+00) then
f(1:nv,i) = fl(1:nv,i)
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(1:nv,i) = fr(1:nv,i)
f(:,i) = fr(:,i)
else ! sl < 0 < sr
@ -5369,13 +5363,13 @@ module schemes
! calculate vectors of the left and right-going waves
!
wl(1:nv) = sl * ul(1:nv,i) - fl(1:nv,i)
wr(1:nv) = sr * ur(1:nv,i) - fr(1:nv,i)
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
uh(1:nv) = ( wr(1:nv) - wl(1:nv)) / srml
fh(1:nv) = (sl * wr(1:nv) - sr * wl(1:nv)) / srml
uh(:) = ( wr(:) - wl(:)) / srml
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
! correct the energy waves
!
@ -5395,7 +5389,7 @@ module schemes
! if Δ < 0, just use the HLL flux
!
if (nr < 1) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! get the contact dicontinuity speed
@ -5405,7 +5399,7 @@ module schemes
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
if ((sm <= sl) .or. (sm >= sr)) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! calculate total pressure (eq. 17 in [1])
@ -5415,7 +5409,7 @@ module schemes
! if the pressure is negative, use the HLL flux
!
if (pr <= 0.0d+00) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! depending in the sign of the contact dicontinuity speed, calculate the proper
@ -5435,7 +5429,7 @@ module schemes
! calculate the flux (eq. 14 in [1])
!
f(1:nv,i) = fl(1:nv,i) + sl * (us(1:nv) - ul(1:nv,i))
f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
else if (sm < 0.0d+00) then
@ -5451,7 +5445,7 @@ module schemes
! calculate the flux (eq. 14 in [1])
!
f(1:nv,i) = fr(1:nv,i) + sr * (us(1:nv) - ur(1:nv,i))
f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
else
@ -5533,7 +5527,7 @@ module schemes
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nv
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
use equations , only : ibx, iby, ibz, ibp
@ -5555,14 +5549,14 @@ module schemes
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nv,nn,nn,nn) :: qq
real(kind=8), dimension(nv,nn,nn,nn,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn,nn) :: qq
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,nn, 1) :: qq
real(kind=8), dimension(nv,nn,nn, 1,2,NDIMS) :: qs
real(kind=8), dimension(nf,nn,nn, 1) :: qq
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,2) :: qi
real(kind=8), dimension(nv,nn) :: fi
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
@ -5793,7 +5787,7 @@ module schemes
! include external procedures
!
use equations , only : nv
use equations , only : nf
use equations , only : ivx, ibx, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
@ -5816,8 +5810,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: wl, wr
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
!
!-------------------------------------------------------------------------------
@ -5876,11 +5870,11 @@ module schemes
!
if (sl >= 0.0d+00) then
f(1:nv,i) = fl(1:nv,i)
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(1:nv,i) = fr(1:nv,i)
f(:,i) = fr(:,i)
else ! sl < 0 < sr
@ -5890,12 +5884,12 @@ module schemes
! calculate vectors of the left and right-going waves
!
wl(1:nv) = sl * ul(1:nv,i) - fl(1:nv,i)
wr(1:nv) = sr * ur(1:nv,i) - fr(1:nv,i)
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
f(1:nv,i) = (sl * wr(1:nv) - sr * wl(1:nv)) / srml
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
@ -5961,7 +5955,7 @@ module schemes
! include external procedures
!
use algebra , only : quadratic
use equations , only : nv
use equations , only : nf
use equations , only : ivx, idn, imx, imy, imz, ien, ibx, iby, ibz, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
@ -5986,8 +5980,8 @@ module schemes
! local arrays to store the states
!
real(kind=8), dimension(nv,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nv) :: uh, us, fh, wl, wr
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: uh, us, fh, wl, wr
real(kind=8), dimension(size(q,2)) :: clm, clp, crm, crp
real(kind=8), dimension(3) :: a, vs
real(kind=8), dimension(2) :: x
@ -6048,11 +6042,11 @@ module schemes
!
if (sl >= 0.0d+00) then
f(1:nv,i) = fl(1:nv,i)
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(1:nv,i) = fr(1:nv,i)
f(:,i) = fr(:,i)
else ! sl < 0 < sr
@ -6062,13 +6056,13 @@ module schemes
! calculate vectors of the left and right-going waves
!
wl(1:nv) = sl * ul(1:nv,i) - fl(1:nv,i)
wr(1:nv) = sr * ur(1:nv,i) - fr(1:nv,i)
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
uh(1:nv) = ( wr(1:nv) - wl(1:nv)) / srml
fh(1:nv) = (sl * wr(1:nv) - sr * wl(1:nv)) / srml
uh(:) = ( wr(:) - wl(:)) / srml
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
! correct the energy waves
!
@ -6100,7 +6094,7 @@ module schemes
! if Δ < 0, just use the HLL flux
!
if (nr < 1) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! get the contact dicontinuity speed
@ -6110,7 +6104,7 @@ module schemes
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
if ((sm <= sl) .or. (sm >= sr)) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! substitute magnetic field components from the HLL state (eq. 37 in [1])
@ -6142,7 +6136,7 @@ module schemes
! if the pressure is negative, use the HLL flux
!
if (pt <= 0.0d+00) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! depending in the sign of the contact dicontinuity speed, calculate the proper
@ -6169,7 +6163,7 @@ module schemes
! calculate the flux (eq. 14 in [1])
!
f(1:nv,i) = fl(1:nv,i) + sl * (us(1:nv) - ul(1:nv,i))
f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
else if (sm < 0.0d+00) then
@ -6192,13 +6186,13 @@ module schemes
! calculate the flux (eq. 14 in [1])
!
f(1:nv,i) = fr(1:nv,i) + sr * (us(1:nv) - ur(1:nv,i))
f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
else
! intermediate flux is constant across the contact discontinuity so use HLL flux
!
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
end if ! sm == 0
@ -6223,7 +6217,7 @@ module schemes
! if Δ < 0, just use the HLL flux
!
if (nr < 1) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! get the contact dicontinuity speed
@ -6233,7 +6227,7 @@ module schemes
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
if ((sm <= sl) .or. (sm >= sr)) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! calculate total pressure (eq. 48 in [1])
@ -6243,7 +6237,7 @@ module schemes
! if the pressure is negative, use the HLL flux
!
if (pt <= 0.0d+00) then
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
else
! depending in the sign of the contact dicontinuity speed, calculate the proper
@ -6267,7 +6261,7 @@ module schemes
! calculate the flux (eq. 27 in [1])
!
f(1:nv,i) = fl(1:nv,i) + sl * (us(1:nv) - ul(1:nv,i))
f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
else if (sm < 0.0d+00) then
@ -6287,13 +6281,13 @@ module schemes
! calculate the flux (eq. 27 in [1])
!
f(1:nv,i) = fr(1:nv,i) + sr * (us(1:nv) - ur(1:nv,i))
f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
else
! intermediate flux is constant across the contact discontinuity so use HLL flux
!
f(1:nv,i) = fh(1:nv)
f(:,i) = fh(:)
end if ! sm == 0