diff options
| -rw-r--r-- | src/GRHydro_ENOReconstruct_drv.F90 | 403 | ||||
| -rw-r--r-- | src/GRHydro_PPMReconstruct_drv.F90 | 482 | ||||
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 786 | ||||
| -rw-r--r-- | src/GRHydro_TVDReconstruct_drv.F90 | 238 | ||||
| -rw-r--r-- | src/make.code.defn | 7 |
5 files changed, 1132 insertions, 784 deletions
diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90 new file mode 100644 index 0000000..656fd6e --- /dev/null +++ b/src/GRHydro_ENOReconstruct_drv.F90 @@ -0,0 +1,403 @@ + /*@@ + @file GRHydro_ENOReconstruct_drv.F90 + @date Tue Jul 19 13:22:03 EDT 2011 + @author Bruno C. Mundim, Joshua Faber + @desc + Driver routine to perform the ENO reconstruction. + @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) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) + + + /*@@ + @routine GRHydro_ENOReconstruct_drv + @date Tue Jul 19 13:24:34 EDT 2011 + @author Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber + @desc + A driver routine to do ENO reconstruction. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_ENOReconstruct_drv(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, dum, dump, dumm +!!$ 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 + +!!$ The dum variables are used as dummies if MHD on but divergence cleaning off + 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)),& + dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dumm(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) + + !$OMP PARALLEL + !$OMP WORKSHARE + trivial_rp = .false. + !$OMP END WORKSHARE NOWAIT + +!!$ Currently only option is reconstruction on primitive variables. +!!$ Should change this. + + !$OMP WORKSHARE + psi4 = 1.d0 + !$OMP END WORKSHARE NOWAIT + if (shift_state .ne. 0) then + !$OMP WORKSHARE + lbetax = betax + lbetay = betay + lbetaz = betaz + !$OMP END WORKSHARE NOWAIT + else + !$OMP WORKSHARE + lbetax = 0.d0 + lbetay = 0.d0 + lbetaz = 0.d0 + !$OMP END WORKSHARE NOWAIT + end if + !$OMP END PARALLEL + +!!$ ENO starts: + 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 (evolve_Y_e .ne. 0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Y_e(:,:,:), Y_e_plus(:,:,:), & + Y_e_minus(:,:,:), & + trivial_rp, hydro_excision_mask) + endif + + if (flux_direction == 1) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + 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)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + 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)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& + psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + + 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, k) + 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)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + 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)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& + psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + + 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, k) + 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,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + 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,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + if(evolve_mhd.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then + call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& + psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + + 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 +!!$ ENO ends. + + deallocate(trivial_rp) + deallocate(psi4, lbetax, lbetay, lbetaz) + deallocate(dum,dump,dumm) + +end subroutine GRHydro_ENOReconstruct_drv + diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90 new file mode 100644 index 0000000..ddd68cd --- /dev/null +++ b/src/GRHydro_PPMReconstruct_drv.F90 @@ -0,0 +1,482 @@ + /*@@ + @file GRHydro_PPMReconstruct_drv.F90 + @date Tue Jul 19 13:22:03 EDT 2011 + @author Bruno C. Mundim, Joshua Faber + @desc + Driver routine to perform the PPM reconstructions. + @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) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) + + + /*@@ + @routine GRHydro_PPMReconstruct_drv + @date Tue Jul 19 13:24:34 EDT 2011 + @author Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber + @desc + A driver routine to do PPM reconstructions. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_PPMReconstruct_drv(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, dum, dump, dumm +!!$ 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 + +!!$ The dum variables are used as dummies if MHD on but divergence cleaning off + 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)),& + dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dumm(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) + + !$OMP PARALLEL + !$OMP WORKSHARE + trivial_rp = .false. + !$OMP END WORKSHARE NOWAIT + +!!$ Currently only option is reconstruction on primitive variables. +!!$ Should change this. + + !$OMP WORKSHARE + psi4 = 1.d0 + !$OMP END WORKSHARE NOWAIT + if (shift_state .ne. 0) then + !$OMP WORKSHARE + lbetax = betax + lbetay = betay + lbetaz = betaz + !$OMP END WORKSHARE NOWAIT + else + !$OMP WORKSHARE + lbetax = 0.d0 + lbetay = 0.d0 + lbetaz = 0.d0 + !$OMP END WORKSHARE NOWAIT + end if + !$OMP END PARALLEL + + +!!$ PPM starts: + if (flux_direction == 1) then + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& + Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),psidc(:,j,k),eps(:,j,k),press(:,j,k),& + rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& + Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),psidcminus(:,j,k),epsminus(:,j,k),& + rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& + Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),& + 1,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 + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& + Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),dum(:,j,k),eps(:,j,k),press(:,j,k),& + rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& + Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),dumm(:,j,k),epsminus(:,j,k),& + rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& + Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),& + 0,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 + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,0,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 + endif !evolve_mhd + if(evolve_tracer.ne.0) then + !$OMP PARALLEL DO PRIVATE(j, k) + 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 + if(evolve_Y_e.ne.0) then + !$OMP PARALLEL DO PRIVATE(j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & + velx(:,j,k),vely(:,j,k),velz(:,j,k), & + Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), & + press(:,j,k)) + end do + end do + !$OMP END PARALLEL DO + end if + + else if (flux_direction == 2) then + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& + Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),& + rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& + Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),psidcminus(j,:,k),epsminus(j,:,k),& + rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& + Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),& + 1,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 + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& + Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),& + rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& + Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),dumm(j,:,k),epsminus(j,:,k),& + rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& + Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),& + 0,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 + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,0,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 + endif !evolve_mhd + if(evolve_tracer.ne.0) then + !$OMP PARALLEL DO PRIVATE(j, k) + 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 + if(evolve_Y_e.ne.0) then + !$OMP PARALLEL DO PRIVATE(j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & + velx(j,:,k),vely(j,:,k),velz(j,:,k), & + Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & + press(j,:,k)) + end do + end do + !$OMP END PARALLEL DO + end if + + else if (flux_direction == 3) then + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& + Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(j,k,:),eps(j,k,:),press(j,k,:),& + rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& + Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(j,k,:),epsminus(j,k,:),& + rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& + Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),& + 1,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 + else !clean_divergence + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& + Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),dum(j,k,:),eps(j,k,:),press(j,k,:),& + rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& + Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(j,k,:),epsminus(j,k,:),& + rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& + Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),& + 0,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 + endif !clean_divergence + else !evolve_mhd + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,0,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 + endif !evolve_mhd + if(evolve_tracer.ne.0) then + !$OMP PARALLEL DO PRIVATE(j, k) + 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 + if(evolve_Y_e.ne.0) then + !$OMP PARALLEL DO PRIVATE(j, k) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & + velx(j,k,:),vely(j,k,:),velz(j,k,:), & + Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(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 +!!$ PPM ends. + + deallocate(trivial_rp) + deallocate(psi4, lbetax, lbetay, lbetaz) + deallocate(dum,dump,dumm) + +end subroutine GRHydro_PPMReconstruct_drv + diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index 5b447c0..a8ac2ec 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -14,20 +14,6 @@ #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) -#define Bvecx(i,j,k) Bvec(i,j,k,1) -#define Bvecy(i,j,k) Bvec(i,j,k,2) -#define Bvecz(i,j,k) Bvec(i,j,k,3) -#define Bconsx(i,j,k) Bcons(i,j,k,1) -#define Bconsy(i,j,k) Bcons(i,j,k,2) -#define Bconsz(i,j,k) Bcons(i,j,k,3) - - /*@@ @routine Reconstruction @date Sat Jan 26 02:13:47 2002 @@ -52,89 +38,11 @@ subroutine Reconstruction(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, dum, dump, dumm -!!$ 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 - -!!$ The dum variables are used as dummies if MHD on but divergence cleaning off - 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)),& - dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& - dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& - dumm(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) - - !$OMP PARALLEL - !$OMP WORKSHARE - trivial_rp = .false. - !$OMP END WORKSHARE NOWAIT - -!!$ Currently only option is reconstruction on primitive variables. -!!$ Should change this. - - !$OMP WORKSHARE - psi4 = 1.d0 - !$OMP END WORKSHARE NOWAIT - if (shift_state .ne. 0) then - !$OMP WORKSHARE - lbetax = betax - lbetay = betay - lbetaz = betaz - !$OMP END WORKSHARE NOWAIT - else - !$OMP WORKSHARE - lbetax = 0.d0 - lbetay = 0.d0 - lbetaz = 0.d0 - !$OMP END WORKSHARE NOWAIT - end if - !!$ Initialize variables that store reconstructed quantities + !$OMP PARALLEL !$OMP WORKSHARE rhoplus = 0.0d0 rhominus = 0.0d0 @@ -183,706 +91,20 @@ subroutine Reconstruction(CCTK_ARGUMENTS) 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 (evolve_Y_e .ne. 0) then - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Y_e(:,:,:), Y_e_plus(:,:,:), & - Y_e_minus(:,:,:), & - trivial_rp, hydro_excision_mask) - endif - - 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) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) - if(evolve_mhd.ne.0) then - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask) - endif - - 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) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - tau, tauplus, tauminus, trivial_rp, hydro_excision_mask) - if(evolve_mhd.ne.0) then - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bcons(:,:,:,1), Bconsxplus, Bconsxminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bcons(:,:,:,2), Bconsyplus, Bconsyminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bcons(:,:,:,3), Bconszplus, Bconszminus, trivial_rp, hydro_excision_mask) - endif - - else - call CCTK_WARN(0, "Variable type to reconstruct not recognized.") - end if - - if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask) - endif - - !$OMP PARALLEL DO PRIVATE(i, j, k) - 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 + call GRHydro_TVDReconstruct_drv(CCTK_PASS_FTOF) else if (CCTK_EQUALS(recon_method,"ppm")) then - if (flux_direction == 1) then - if(evolve_mhd.ne.0) then - if(clean_divergence.ne.0) then - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& - rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& - Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),psidc(:,j,k),eps(:,j,k),press(:,j,k),& - rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& - Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),psidcminus(:,j,k),epsminus(:,j,k),& - rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& - Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),& - 1,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 - else !clean_divergence - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& - rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& - Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),dum(:,j,k),eps(:,j,k),press(:,j,k),& - rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& - Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),dumm(:,j,k),epsminus(:,j,k),& - rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& - Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),& - 0,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 - endif !clean_divergence - else !evolve_mhd - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,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 - endif !evolve_mhd - if(evolve_tracer.ne.0) then - !$OMP PARALLEL DO PRIVATE(j, k) - 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 - if(evolve_Y_e.ne.0) then - !$OMP PARALLEL DO PRIVATE(j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & - velx(:,j,k),vely(:,j,k),velz(:,j,k), & - Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), & - press(:,j,k)) - end do - end do - !$OMP END PARALLEL DO - end if - - else if (flux_direction == 2) then - if(evolve_mhd.ne.0) then - if(clean_divergence.ne.0) then - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& - rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& - Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),& - rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& - Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),psidcminus(j,:,k),epsminus(j,:,k),& - rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& - Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),& - 1,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 - else !clean_divergence - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& - rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& - Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),& - rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& - Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),dumm(j,:,k),epsminus(j,:,k),& - rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& - Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),& - 0,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 - endif !clean_divergence - else !evolve_mhd - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,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 - endif !evolve_mhd - if(evolve_tracer.ne.0) then - !$OMP PARALLEL DO PRIVATE(j, k) - 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 - if(evolve_Y_e.ne.0) then - !$OMP PARALLEL DO PRIVATE(j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & - velx(j,:,k),vely(j,:,k),velz(j,:,k), & - Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & - press(j,:,k)) - end do - end do - !$OMP END PARALLEL DO - end if + call GRHydro_PPMReconstruct_drv(CCTK_PASS_FTOF) - else if (flux_direction == 3) then - if(evolve_mhd.ne.0) then - if(clean_divergence.ne.0) then - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& - rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& - Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(j,k,:),eps(j,k,:),press(j,k,:),& - rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& - Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(j,k,:),epsminus(j,k,:),& - rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& - Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),& - 1,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 - else !clean_divergence - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& - rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& - Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),dum(j,k,:),eps(j,k,:),press(j,k,:),& - rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& - Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(j,k,:),epsminus(j,k,:),& - rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& - Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),& - 0,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 - endif !clean_divergence - else !evolve_mhd - !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,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 - endif !evolve_mhd - if(evolve_tracer.ne.0) then - !$OMP PARALLEL DO PRIVATE(j, k) - 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 - if(evolve_Y_e.ne.0) then - !$OMP PARALLEL DO PRIVATE(j, k) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & - velx(j,k,:),vely(j,k,:),velz(j,k,:), & - Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(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 (evolve_Y_e .ne. 0) then - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Y_e(:,:,:), Y_e_plus(:,:,:), & - Y_e_minus(:,:,:), & - trivial_rp, hydro_excision_mask) - endif - - if (flux_direction == 1) then - !$OMP PARALLEL DO PRIVATE(i, j, k) - 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)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - if(evolve_mhd.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - endif - 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)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - if(evolve_mhd.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - endif - else - !$OMP CRITICAL - call CCTK_WARN(0, "Variable type to reconstruct not recognized.") - !$OMP END CRITICAL - end if - - if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& - psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) - endif - - 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, k) - 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)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - if(evolve_mhd.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - endif - 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)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - if(evolve_mhd.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - endif - else - !$OMP CRITICAL - call CCTK_WARN(0, "Variable type to reconstruct not recognized.") - !$OMP END CRITICAL - end if + call GRHydro_ENOReconstruct_drv(CCTK_PASS_FTOF) - if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& - psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) - endif - - 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, k) - 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,:)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - if(evolve_mhd.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - endif - 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,:)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - if(evolve_mhd.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - endif - else - !$OMP CRITICAL - call CCTK_WARN(0, "Variable type to reconstruct not recognized.") - !$OMP END CRITICAL - end if - - if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then - call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& - psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) - endif - - 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) - deallocate(dum,dump,dumm) - !$OMP PARALLEL WORKSHARE where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min) ) rhoplus = rho diff --git a/src/GRHydro_TVDReconstruct_drv.F90 b/src/GRHydro_TVDReconstruct_drv.F90 new file mode 100644 index 0000000..6fc2356 --- /dev/null +++ b/src/GRHydro_TVDReconstruct_drv.F90 @@ -0,0 +1,238 @@ + /*@@ + @file GRHydro_TVDReconstruct_drv.F90 + @date Tue Jul 19 13:22:03 EDT 2011 + @author Bruno C. Mundim, Joshua Faber + @desc + Driver routine to perform the TVD reconstruction. + @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) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) + + + /*@@ + @routine GRHydro_TVDReconstruct_drv + @date Tue Jul 19 13:24:34 EDT 2011 + @author Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber + @desc + A driver routine to do TVD reconstruction. Currently just does + TVD on the primitive variables. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_TVDReconstruct_drv(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, dum, dump, dumm +!!$ 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 + +!!$ The dum variables are used as dummies if MHD on but divergence cleaning off + 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)),& + dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dumm(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) + + !$OMP PARALLEL + !$OMP WORKSHARE + trivial_rp = .false. + !$OMP END WORKSHARE NOWAIT + +!!$ Currently only option is reconstruction on primitive variables. +!!$ Should change this. + + !$OMP WORKSHARE + psi4 = 1.d0 + !$OMP END WORKSHARE NOWAIT + if (shift_state .ne. 0) then + !$OMP WORKSHARE + lbetax = betax + lbetay = betay + lbetaz = betaz + !$OMP END WORKSHARE NOWAIT + else + !$OMP WORKSHARE + lbetax = 0.d0 + lbetay = 0.d0 + lbetaz = 0.d0 + !$OMP END WORKSHARE NOWAIT + end if + !$OMP END PARALLEL + + +!!$ TVD starts: + 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 (evolve_Y_e .ne. 0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Y_e(:,:,:), Y_e_plus(:,:,:), & + Y_e_minus(:,:,:), & + trivial_rp, hydro_excision_mask) + endif + + 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) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) + if(evolve_mhd.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask) + endif + + 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) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + tau, tauplus, tauminus, trivial_rp, hydro_excision_mask) + if(evolve_mhd.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,1), Bconsxplus, Bconsxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,2), Bconsyplus, Bconsyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Bcons(:,:,:,3), Bconszplus, Bconszminus, trivial_rp, hydro_excision_mask) + endif + + else + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + end if + + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask) + endif + + !$OMP PARALLEL DO PRIVATE(i, j, k) + 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 +!!$ TVD ends. + + deallocate(trivial_rp) + deallocate(psi4, lbetax, lbetay, lbetaz) + deallocate(dum,dump,dumm) + +end subroutine GRHydro_TVDReconstruct_drv + + diff --git a/src/make.code.defn b/src/make.code.defn index 27c2ecb..93c6b5d 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -7,9 +7,10 @@ SRCS = Utils.F90 \ GRHydro_Boundaries.F90 \ GRHydro_CalcUpdate.F90 \ GRHydro_Con2Prim.F90 \ - GRHydro_DivergenceClean.F90 \ + GRHydro_DivergenceClean.F90 \ GRHydro_Eigenproblem.F90 \ GRHydro_Eigenproblem_Marquina.F90 \ + GRHydro_ENOReconstruct_drv.F90 \ GRHydro_ENOReconstruct.F90 \ GRHydro_ENOScalars.F90 \ GRHydro_EOS.c \ @@ -19,6 +20,7 @@ SRCS = Utils.F90 \ GRHydro_Loop.F90 \ GRHydro_Minima.F90 \ GRHydro_Minima.cc \ + GRHydro_PPMReconstruct_drv.F90 \ GRHydro_PPM.F90 \ GRHydro_ParamCheck.F90 \ GRHydro_Particle.F90 \ @@ -36,6 +38,7 @@ SRCS = Utils.F90 \ GRHydro_SlopeLimiter.F90 \ GRHydro_Source.F90 \ GRHydro_Startup.F90 \ + GRHydro_TVDReconstruct_drv.F90 \ GRHydro_TVDReconstruct.F90 \ GRHydro_Marquina.F90 \ GRHydro_UpdateMask.F90 \ @@ -49,7 +52,7 @@ SRCS = Utils.F90 \ GRHydro_EigenproblemM.F90 \ GRHydro_FluxM.F90 \ GRHydro_HLLEM.F90 \ - GRHydro_PPMM.F90 \ + GRHydro_PPMM.F90 \ GRHydro_Prim2ConM.F90 \ GRHydro_RiemannSolveM.F90 \ GRHydro_SourceM.F90 \ |