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

SCHEMES: Adopt update_flux_srhd_adi() to use separate reconstruction.

Browse Source 
Implement reconstruction based on 3-vector and 4-vector as well.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2016-08-14 22:52:40 -03:00
parent 5c8d60e47e
commit 7a2c83d0ea
1 changed files with 92 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

126
src/schemes.F90
View File

@ -49,6 +49,10 @@ module schemes
integer , save :: imi, imf, ims, imr
#endif /* PROFILE */
! 4-vector reconstruction flag
!
logical , save :: states_4vec = .false.
! pointer to the flux update procedure
!
procedure(update_flux_hd_iso), pointer, save :: update_flux => null()
@ -377,6 +381,7 @@ module schemes
! set pointers to subroutines
!
states => states_srhd_adi_4vec
states_4vec = .true.
! in the case of state variables, revert to primitive
!
@ -449,6 +454,7 @@ module schemes
! set pointers to subroutines
!
states => states_srmhd_adi_4vec
states_4vec = .true.
! in the case of state variables, revert to primitive
!
@ -4668,13 +4674,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
!
!-------------------------------------------------------------------------------
!
@ -4688,6 +4700,64 @@ 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)
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
@ -4695,19 +4765,15 @@ 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(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(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
!
@ -4727,19 +4793,15 @@ 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(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(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
!
@ -4760,19 +4822,15 @@ 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(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(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
!