diff options
Diffstat (limited to 'src/GRHydro_ReconstructPoly.F90')
-rw-r--r-- | src/GRHydro_ReconstructPoly.F90 | 530 |
1 files changed, 530 insertions, 0 deletions
diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90 new file mode 100644 index 0000000..0467c27 --- /dev/null +++ b/src/GRHydro_ReconstructPoly.F90 @@ -0,0 +1,530 @@ + /*@@ + @file GRHydro_ReconstructPoly.F90 + @date Sat Jan 26 02:13:25 2002 + @author + @desc + Wrapper routine to perform the reconstruction for polytropes. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" + +#include "SpaceMask.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) + + + /*@@ + @routine ReconstructionPolytype + @date Tue Mar 19 11:36:55 2002 + @author Ian Hawke + @desc + If using a polytropic type EOS, we do not have to reconstruct all the + variables. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine ReconstructionPolytype(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer :: nx, ny, nz, i, j, k, itracer + + logical, dimension(:,:,:), allocatable :: trivial_rp +!!$ logical, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: trivial_rp + + CCTK_INT :: type_bitsx, trivialx, not_trivialx, & + &type_bitsy, trivialy, not_trivialy, & + &type_bitsz, trivialz, not_trivialz + + CCTK_REAL, dimension(:,:,:),allocatable :: & + &psi4, lbetax, lbetay, lbetaz +!!$ CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: & +!!$ &psi4, lbetax, lbetay, lbetaz + + CCTK_INT :: ierr + + CCTK_REAL :: local_min_tracer + + allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + if (ierr .ne. 0) then + call CCTK_WARN(0, "Allocation problems with trivial_rp") + end if + + allocate(psi4(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + lbetax(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + lbetay(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + lbetaz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + if (ierr .ne. 0) then + call CCTK_WARN(0, "Allocation problems with lbeta") + end if + + call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialx, "Hydro_RiemannProblemX", & + &"not_trivial") + call SpaceMask_GetTypeBits(type_bitsy, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(trivialy, "Hydro_RiemannProblemY", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialy, "Hydro_RiemannProblemY", & + &"not_trivial") + call SpaceMask_GetTypeBits(type_bitsz, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(trivialz, "Hydro_RiemannProblemZ", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialz, "Hydro_RiemannProblemZ", & + &"not_trivial") + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + trivial_rp = .false. + +!!$ Currently only option is reconstruction on primitive variables. +!!$ Should change this. + + if (conformal_state .eq. 0) then + psi4 = 1.d0 + else + psi4 = psi**4 + end if + if (shift_state .ne. 0) then + lbetax = betax + lbetay = betay + lbetaz = betaz + else + lbetax = 0.d0 + lbetay = 0.d0 + lbetaz = 0.d0 + end if + +!!$ Initialize variables that store reconstructed quantities + + rhoplus = 0.0d0 + rhominus = 0.0d0 + epsplus = 0.0d0 + epsminus = 0.0d0 + velxplus = 0.0d0 + velxminus = 0.0d0 + velyplus = 0.0d0 + velyminus = 0.0d0 + velzplus = 0.0d0 + velzminus = 0.0d0 + + if (evolve_tracer .ne. 0) then + tracerplus = 0.0d0 + tracerminus = 0.0d0 + endif + + if (CCTK_EQUALS(recon_method,"tvd")) then + + if (evolve_tracer .ne. 0) then + do itracer=1,number_of_tracers + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), & + tracerminus(:,:,:,itracer), & + trivial_rp, hydro_excision_mask) + enddo + end if + + + if (CCTK_EQUALS(recon_vars,"primitive")) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) + + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + dens, densplus, densminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + scon(:,:,:,1), sxplus, sxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + scon(:,:,:,2), syplus, syminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask) + + else + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + end if + + !$OMP PARALLEL DO PRIVATE(i, j) + do k = 1, nz + do j = 1, ny + do i = 1, nx + if (trivial_rp(i,j,k)) then + if (flux_direction == 1) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else if (flux_direction == 2) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsy, trivialy) + else if (flux_direction == 3) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsz, trivialz) + end if + else + if (flux_direction == 1) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + else if (flux_direction == 2) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsy, not_trivialy) + else if (flux_direction == 3) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsz, not_trivialz) + end if + end if + end do + end do + end do + !$OMP END PARALLEL DO + + else if (CCTK_EQUALS(recon_method,"ppm")) then + + if (flux_direction == 1) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& + press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& + velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& + velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& + gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & + gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& + w_lorentz(:,j,k), & + 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & + GRHydro_mppm_eigenvalue_x_right, & + GRHydro_mppm_xwind) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + if(evolve_tracer.ne.0) then + !$OMP PARALLEL DO PRIVATE(j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_tracer_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & + velx(:,j,k),vely(:,j,k),velz(:,j,k), & + tracer(:,j,k,:),tracerminus(:,j,k,:),tracerplus(:,j,k,:), & + press(:,j,k)) + end do + end do + !$OMP END PARALLEL DO + end if + + else if (flux_direction == 2) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& + press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& + velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& + velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& + gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & + gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& + w_lorentz(j,:,k), & + 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & + GRHydro_mppm_eigenvalue_y_right, & + GRHydro_mppm_xwind) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do + end do + end do + !$OMP END PARALLEL DO + if(evolve_tracer.ne.0) then + !$OMP PARALLEL DO PRIVATE(j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_tracer_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & + vely(j,:,k),velz(j,:,k),velx(j,:,k), & + tracer(j,:,k,:),tracerminus(j,:,k,:),tracerplus(j,:,k,:), & + press(j,:,k)) + end do + end do + !$OMP END PARALLEL DO + end if + + else if (flux_direction == 3) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& + press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& + velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& + velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& + gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & + gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& + w_lorentz(j,k,:), & + 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & + GRHydro_mppm_eigenvalue_z_right, & + GRHydro_mppm_xwind) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do + end do + end do + !$OMP END PARALLEL DO + if(evolve_tracer.ne.0) then + !$OMP PARALLEL DO PRIVATE(j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + + call SimplePPM_tracer_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & + velz(j,k,:),velx(j,k,:),vely(j,k,:), & + tracer(j,k,:,:),tracerminus(j,k,:,:),tracerplus(j,k,:,:), & + press(j,k,:)) + end do + end do + !$OMP END PARALLEL DO + end if + + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + else if (CCTK_EQUALS(recon_method,"eno")) then + + if (evolve_tracer .ne. 0) then + do itracer=1,number_of_tracers + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), & + tracerminus(:,:,:,itracer), trivial_rp, & + hydro_excision_mask) + enddo + end if + + if (flux_direction == 1) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + if (CCTK_EQUALS(recon_vars,"primitive")) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + velx(:,j,k),velxminus(:,j,k),velxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + vely(:,j,k),velyminus(:,j,k),velyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + dens(:,j,k),densminus(:,j,k),densplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + sx(:,j,k),sxminus(:,j,k),sxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + sy(:,j,k),syminus(:,j,k),syplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + sz(:,j,k),szminus(:,j,k),szplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + do i = 1, cctk_lsh(1) + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else if (flux_direction == 2) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + if (CCTK_EQUALS(recon_vars,"primitive")) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + dens(j,:,k),densminus(j,:,k),densplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + sy(j,:,k),syminus(j,:,k),syplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + sz(j,:,k),szminus(j,:,k),szplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + do i = 1, cctk_lsh(2) + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else if (flux_direction == 3) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + if (CCTK_EQUALS(recon_vars,"primitive")) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + velx(j,k,:),velxminus(j,k,:),velxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + vely(j,k,:),velyminus(j,k,:),velyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + dens(j,k,:),densminus(j,k,:),densplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + sx(j,k,:),sxminus(j,k,:),sxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + sy(j,k,:),syminus(j,k,:),syplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + sz(j,k,:),szminus(j,k,:),szplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + do i = 1, cctk_lsh(3) + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + else + call CCTK_WARN(0, "Reconstruction method not recognized!") + end if + + deallocate(trivial_rp) + deallocate(psi4, lbetax, lbetay, lbetaz) + + !$OMP WORKSHARE + where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min).or.& + (epsplus < 0.d0).or.(epsminus < 0.d0) ) + rhoplus = rho + rhominus = rho + velxplus = vel(:,:,:,1) + velxminus = vel(:,:,:,1) + velyplus = vel(:,:,:,2) + velyminus = vel(:,:,:,2) + velzplus = vel(:,:,:,3) + velzminus = vel(:,:,:,3) + epsplus = eps + epsminus = eps + end where + !$OMP END WORKSHARE + + if (evolve_tracer .ne. 0) then + if (use_min_tracer .ne. 0) then + local_min_tracer = min_tracer + else + local_min_tracer = 0d0 + end if + + !$OMP WORKSHARE + where( (tracerplus .le. local_min_tracer).or.& + (tracerminus .le. local_min_tracer) ) + tracerplus = tracer + tracerminus = tracer + end where + !$OMP END WORKSHARE + ! Call the conserved tracer routine in any case because (accord. to + ! Christian Ott) this is the only way this works + call Prim2ConservativeTracer(CCTK_PASS_FTOF) + endif + + if (CCTK_EQUALS(recon_vars,"primitive").or.& + CCTK_EQUALS(recon_method,"ppm")) then + if (use_eosgeneral == 0) then + call Prim2ConservativePolytype(CCTK_PASS_FTOF) + else + call primitive2conservativegeneral(CCTK_PASS_FTOF) + end if + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call Con2PrimBoundsPolytype(CCTK_PASS_FTOF) + else + call CCTK_WARN(0,"Variable type to reconstruct not recognized.") + end if + + return + +end subroutine ReconstructionPolytype + |