diff --git a/src/schemes.F90 b/src/schemes.F90 index b939795..59482cd 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -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 !