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

SCHEMES: Adopt update_flux_srmhd_adi() to use separate reconstruction.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2016-08-14 23:01:16 -03:00
parent 7a2c83d0ea
commit 5f25de32f9
1 changed files with 103 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

149
src/schemes.F90
View File

@ -5414,13 +5414,19 @@ module schemes
! local variables
!
integer :: i, j, k
integer :: i, j, k, l, p
real(kind=8) :: vm
! local temporary arrays
!
real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx
real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy
real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz
real(kind=8), dimension(nv,im,jm,km) :: qq
real(kind=8), dimension(nv,im,jm,km,2,NDIMS) :: qs
real(kind=8), dimension(nv,im,2) :: qx
real(kind=8), dimension(nv,jm,2) :: qy
real(kind=8), dimension(nv,km,2) :: qz
real(kind=8), dimension(nv,im) :: fx
real(kind=8), dimension(nv,jm) :: fy
real(kind=8), dimension(nv,km) :: fz
!
!-------------------------------------------------------------------------------
!
@ -5434,6 +5440,69 @@ module schemes
!
f(1:NDIMS,1:nv,1:im,1:jm,1:km) = 0.0d+00
! apply reconstruction on variables using 4-vector or 3-vector
!
if (states_4vec) then
! convert velocities to four-velocities for physical reconstruction
!
do k = 1, km
do j = 1, jm
do i = 1, im
vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))
qq(idn,i,j,k) = q(idn,i,j,k) / vm
qq(ivx,i,j,k) = q(ivx,i,j,k) / vm
qq(ivy,i,j,k) = q(ivy,i,j,k) / vm
qq(ivz,i,j,k) = q(ivz,i,j,k) / vm
qq(ipr,i,j,k) = q(ipr,i,j,k)
qq(ibx,i,j,k) = q(ibx,i,j,k)
qq(iby,i,j,k) = q(iby,i,j,k)
qq(ibz,i,j,k) = q(ibz,i,j,k)
qq(ibp,i,j,k) = q(ibp,i,j,k)
qq(ipr,i,j,k) = q(ipr,i,j,k)
end do ! i = 1, im
end do ! j = 1, jm
end do ! k = 1, km
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), qq(1:nv,1:im,1:jm,1:km) &
, qs(1:nv,1:im,1:jm,1:km,1:2,1:NDIMS))
! convert state four-velocities back to velocities
!
do k = 1, km
do j = 1, jm
do i = 1, im
do l = 1, 2
do p = 1, NDIMS
vm = sqrt(1.0d+00 + sum(qs(ivx:ivz,i,j,k,l,p)**2))
qs(idn,i,j,k,l,p) = qs(idn,i,j,k,l,p) / vm
qs(ivx,i,j,k,l,p) = qs(ivx,i,j,k,l,p) / vm
qs(ivy,i,j,k,l,p) = qs(ivy,i,j,k,l,p) / vm
qs(ivz,i,j,k,l,p) = qs(ivz,i,j,k,l,p) / vm
end do ! p = 1, ndims
end do ! l = 1, 2
end do ! i = 1, im
end do ! j = 1, jm
end do ! k = 1, km
else
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), q (1:nv,1:im,1:jm,1:km) &
, qs(1:nv,1:im,1:jm,1:km,1:2,1:NDIMS))
end if
! calculate the flux along the X-direction
!
do k = kbl, keu
@ -5441,23 +5510,19 @@ module schemes
! copy directional variable vectors to pass to the one dimensional solver
!
qx(idn,1:im) = q(idn,1:im,j,k)
qx(ivx,1:im) = q(ivx,1:im,j,k)
qx(ivy,1:im) = q(ivy,1:im,j,k)
qx(ivz,1:im) = q(ivz,1:im,j,k)
qx(ibx,1:im) = q(ibx,1:im,j,k)
qx(iby,1:im) = q(iby,1:im,j,k)
qx(ibz,1:im) = q(ibz,1:im,j,k)
qx(ibp,1:im) = q(ibp,1:im,j,k)
qx(ipr,1:im) = q(ipr,1:im,j,k)
! reconstruct Riemann states
!
call states(im, dx(1), qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im))
qx(idn,1:im,1:2) = qs(idn,1:im,j,k,1:2,1)
qx(ivx,1:im,1:2) = qs(ivx,1:im,j,k,1:2,1)
qx(ivy,1:im,1:2) = qs(ivy,1:im,j,k,1:2,1)
qx(ivz,1:im,1:2) = qs(ivz,1:im,j,k,1:2,1)
qx(ibx,1:im,1:2) = qs(ibx,1:im,j,k,1:2,1)
qx(iby,1:im,1:2) = qs(iby,1:im,j,k,1:2,1)
qx(ibz,1:im,1:2) = qs(ibz,1:im,j,k,1:2,1)
qx(ibp,1:im,1:2) = qs(ibp,1:im,j,k,1:2,1)
qx(ipr,1:im,1:2) = qs(ipr,1:im,j,k,1:2,1)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im))
call riemann(im, qx(1:nv,1:im,1), qx(1:nv,1:im,2), fx(1:nv,1:im))
! update the array of fluxes
!
@ -5481,23 +5546,19 @@ module schemes
! copy directional variable vectors to pass to the one dimensional solver
!
qy(idn,1:jm) = q(idn,i,1:jm,k)
qy(ivx,1:jm) = q(ivy,i,1:jm,k)
qy(ivy,1:jm) = q(ivz,i,1:jm,k)
qy(ivz,1:jm) = q(ivx,i,1:jm,k)
qy(ibx,1:jm) = q(iby,i,1:jm,k)
qy(iby,1:jm) = q(ibz,i,1:jm,k)
qy(ibz,1:jm) = q(ibx,i,1:jm,k)
qy(ibp,1:jm) = q(ibp,i,1:jm,k)
qy(ipr,1:jm) = q(ipr,i,1:jm,k)
! reconstruct Riemann states
!
call states(jm, dx(2), qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm))
qy(idn,1:jm,1:2) = qs(idn,i,1:jm,k,1:2,2)
qy(ivx,1:jm,1:2) = qs(ivy,i,1:jm,k,1:2,2)
qy(ivy,1:jm,1:2) = qs(ivz,i,1:jm,k,1:2,2)
qy(ivz,1:jm,1:2) = qs(ivx,i,1:jm,k,1:2,2)
qy(ibx,1:jm,1:2) = qs(iby,i,1:jm,k,1:2,2)
qy(iby,1:jm,1:2) = qs(ibz,i,1:jm,k,1:2,2)
qy(ibz,1:jm,1:2) = qs(ibx,i,1:jm,k,1:2,2)
qy(ibp,1:jm,1:2) = qs(ibp,i,1:jm,k,1:2,2)
qy(ipr,1:jm,1:2) = qs(ipr,i,1:jm,k,1:2,2)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm))
call riemann(jm, qy(1:nv,1:jm,1), qy(1:nv,1:jm,2), fy(1:nv,1:jm))
! update the array of fluxes
!
@ -5522,23 +5583,19 @@ module schemes
! copy directional variable vectors to pass to the one dimensional solver
!
qz(idn,1:km) = q(idn,i,j,1:km)
qz(ivx,1:km) = q(ivz,i,j,1:km)
qz(ivy,1:km) = q(ivx,i,j,1:km)
qz(ivz,1:km) = q(ivy,i,j,1:km)
qz(ibx,1:km) = q(ibz,i,j,1:km)
qz(iby,1:km) = q(ibx,i,j,1:km)
qz(ibz,1:km) = q(iby,i,j,1:km)
qz(ibp,1:km) = q(ibp,i,j,1:km)
qz(ipr,1:km) = q(ipr,i,j,1:km)
! reconstruct Riemann states
!
call states(km, dx(3), qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km))
qz(idn,1:km,1:2) = qs(idn,i,j,1:km,1:2,3)
qz(ivx,1:km,1:2) = qs(ivz,i,j,1:km,1:2,3)
qz(ivy,1:km,1:2) = qs(ivx,i,j,1:km,1:2,3)
qz(ivz,1:km,1:2) = qs(ivy,i,j,1:km,1:2,3)
qz(ibx,1:km,1:2) = qs(ibz,i,j,1:km,1:2,3)
qz(iby,1:km,1:2) = qs(ibx,i,j,1:km,1:2,3)
qz(ibz,1:km,1:2) = qs(iby,i,j,1:km,1:2,3)
qz(ibp,1:km,1:2) = qs(ibp,i,j,1:km,1:2,3)
qz(ipr,1:km,1:2) = qs(ipr,i,j,1:km,1:2,3)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km))
call riemann(km, qz(1:nv,1:km,1), qz(1:nv,1:km,2), fz(1:nv,1:km))
! update the array of fluxes
!