diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/GRHydro_ENOReconstruct_drv.F90 | 83 | ||||
| -rw-r--r-- | src/GRHydro_HLLEM.F90 | 23 | ||||
| -rw-r--r-- | src/GRHydro_PPMMReconstruct_drv.F90 | 386 | ||||
| -rw-r--r-- | src/GRHydro_PPMReconstruct_drv.F90 | 262 | ||||
| -rw-r--r-- | src/GRHydro_ParamCheck.F90 | 13 | ||||
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 121 | ||||
| -rw-r--r-- | src/GRHydro_TVDReconstruct_drv.F90 | 90 | ||||
| -rw-r--r-- | src/make.code.defn | 1 |
8 files changed, 614 insertions, 365 deletions
diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90 index 656fd6e..65ff19c 100644 --- a/src/GRHydro_ENOReconstruct_drv.F90 +++ b/src/GRHydro_ENOReconstruct_drv.F90 @@ -67,8 +67,6 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS) 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") @@ -107,31 +105,62 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS) 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 +!!$ Initialize variables that store reconstructed quantities: + + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + trivial_rp(i,j,k) = .false. + psi4(i,j,k) = 1.0d0 + rhoplus(i,j,k) = 0.0d0 + rhominus(i,j,k)= 0.0d0 + epsplus(i,j,k) = 0.0d0 + epsminus(i,j,k) = 0.0d0 + velxplus(i,j,k) = 0.0d0 + velxminus(i,j,k) = 0.0d0 + velyplus(i,j,k) = 0.0d0 + velyminus(i,j,k) = 0.0d0 + velzplus(i,j,k) = 0.0d0 + velzminus(i,j,k) = 0.0d0 + + if(evolve_mhd.ne.0) then + Bvecxplus(i,j,k) = 0.0d0 + Bvecxminus(i,j,k) = 0.0d0 + Bvecyplus(i,j,k) = 0.0d0 + Bvecyminus(i,j,k) = 0.0d0 + Bveczplus(i,j,k) = 0.0d0 + Bveczminus(i,j,k) = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus(i,j,k) = 0.0d0 + psidcminus(i,j,k) = 0.0d0 + endif + endif + + if (evolve_tracer .ne. 0) then + tracerplus(i,j,k,:) = 0.0d0 + tracerminus(i,j,k,:) = 0.0d0 + endif + + if (evolve_Y_e .ne. 0) then + Y_e_plus(i,j,k) = 0.0d0 + Y_e_minus(i,j,k) = 0.0d0 + endif + + if (shift_state .ne. 0) then + lbetax(i,j,k) = betax(i,j,k) + lbetay(i,j,k) = betay(i,j,k) + lbetaz(i,j,k) = betaz(i,j,k) + else + lbetax(i,j,k) = 0.d0 + lbetay(i,j,k) = 0.d0 + lbetaz(i,j,k) = 0.d0 + end if + + enddo + enddo + enddo + !$OMP END PARALLEL DO !!$ ENO starts: if (evolve_tracer .ne. 0) then diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index e841450..b85efb6 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -426,19 +426,24 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) end do if(clean_divergence.ne.0) then - psidcdiff = psidcm - psidcp - -!!$ psidcf = (charmax * psidcfp - charmin * psidcfm + & -!!$ charmax * charmin * psidcdiff) / charpm - -!!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic? - + psidcdiff = psidcm - psidcp + select case(whichpsidcspeed) + case(0) + psidcf = (charmax * psidcfp - charmin * psidcfm + & + charmax * charmin * psidcdiff) / charpm + case(1) + !!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic? psidcf = (1.d0 * psidcfp - (-1.d0) * psidcfm + & 1.d0 * (-1.d0) * psidcdiff) / 2.d0 - - + case(2) + charmax = setcharmax + charmin = setcharmin + psidcf = (charmax * psidcfp - charmin * psidcfm + & + charmax * charmin * psidcdiff) / charpm + end select endif + end if !!! The end of the SpaceMask check for a trivial RP. densflux(i, j, k) = f1(1) diff --git a/src/GRHydro_PPMMReconstruct_drv.F90 b/src/GRHydro_PPMMReconstruct_drv.F90 new file mode 100644 index 0000000..9486a77 --- /dev/null +++ b/src/GRHydro_PPMMReconstruct_drv.F90 @@ -0,0 +1,386 @@ + /*@@ + @file GRHydro_PPMMReconstruct_drv.F90 + @date Wed Jul 27 15:17:03 EDT 2011 + @author Bruno C. Mundim, Joshua Faber + @desc + Driver routine to perform the magnetic version of PPM 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_PPMMReconstruct_drv + @date Wed Jul 27 15:17:55 EDT 2011 + @author Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber + @desc + A driver routine to do magnetic version of PPM reconstructions. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_PPMMReconstruct_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 + + 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) + +!!$ Initialize variables that store reconstructed quantities: + + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + trivial_rp(i,j,k) = .false. + psi4(i,j,k) = 1.0d0 + rhoplus(i,j,k) = 0.0d0 + rhominus(i,j,k)= 0.0d0 + epsplus(i,j,k) = 0.0d0 + epsminus(i,j,k) = 0.0d0 + velxplus(i,j,k) = 0.0d0 + velxminus(i,j,k) = 0.0d0 + velyplus(i,j,k) = 0.0d0 + velyminus(i,j,k) = 0.0d0 + velzplus(i,j,k) = 0.0d0 + velzminus(i,j,k) = 0.0d0 + + Bvecxplus(i,j,k) = 0.0d0 + Bvecxminus(i,j,k) = 0.0d0 + Bvecyplus(i,j,k) = 0.0d0 + Bvecyminus(i,j,k) = 0.0d0 + Bveczplus(i,j,k) = 0.0d0 + Bveczminus(i,j,k) = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus(i,j,k) = 0.0d0 + psidcminus(i,j,k) = 0.0d0 + endif + + if (evolve_tracer .ne. 0) then + tracerplus(i,j,k,:) = 0.0d0 + tracerminus(i,j,k,:) = 0.0d0 + endif + + if (evolve_Y_e .ne. 0) then + Y_e_plus(i,j,k) = 0.0d0 + Y_e_minus(i,j,k) = 0.0d0 + endif + + if (shift_state .ne. 0) then + lbetax(i,j,k) = betax(i,j,k) + lbetay(i,j,k) = betay(i,j,k) + lbetaz(i,j,k) = betaz(i,j,k) + else + lbetax(i,j,k) = 0.d0 + lbetay(i,j,k) = 0.d0 + lbetaz(i,j,k) = 0.d0 + end if + + enddo + enddo + enddo + !$OMP END PARALLEL DO + +!!$ PPM starts: + if (flux_direction == 1) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + if(clean_divergence.ne.0) then + 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) + else !clean_divergence + 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) + endif !clean_divergence + 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, 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 + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + if(clean_divergence.ne.0) then + 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) + else !clean_divergence + 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) + endif !clean_divergence + 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, 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 + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + if(clean_divergence.ne.0) then + 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) + else !clean_divergence + 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) + endif !clean_divergence + 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, 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_PPMMReconstruct_drv + diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90 index ddd68cd..854d77e 100644 --- a/src/GRHydro_PPMReconstruct_drv.F90 +++ b/src/GRHydro_PPMReconstruct_drv.F90 @@ -20,12 +20,6 @@ #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) /*@@ @@ -61,30 +55,25 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) &type_bitsz, trivialz, not_trivialz CCTK_REAL, dimension(:,:,:),allocatable :: & - &psi4, lbetax, lbetay, lbetaz, dum, dump, dumm + &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 -!!$ 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) + STAT=ierr) if (ierr .ne. 0) then - call CCTK_WARN(0, "Allocation problems with lbeta") + call CCTK_WARN(0, "Allocation problems with psi4 or lbeta") end if call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX") @@ -107,94 +96,53 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) ny = cctk_lsh(2) nz = cctk_lsh(3) - !$OMP PARALLEL - !$OMP WORKSHARE - trivial_rp = .false. - !$OMP END WORKSHARE NOWAIT +!!$ Initialize variables that store reconstructed quantities: -!!$ Currently only option is reconstruction on primitive variables. -!!$ Should change this. + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + trivial_rp(i,j,k) = .false. + psi4(i,j,k) = 1.0d0 + rhoplus(i,j,k) = 0.0d0 + rhominus(i,j,k)= 0.0d0 + epsplus(i,j,k) = 0.0d0 + epsminus(i,j,k) = 0.0d0 + velxplus(i,j,k) = 0.0d0 + velxminus(i,j,k) = 0.0d0 + velyplus(i,j,k) = 0.0d0 + velyminus(i,j,k) = 0.0d0 + velzplus(i,j,k) = 0.0d0 + velzminus(i,j,k) = 0.0d0 - !$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 + if (evolve_tracer .ne. 0) then + tracerplus(i,j,k,:) = 0.0d0 + tracerminus(i,j,k,:) = 0.0d0 + endif + + if (evolve_Y_e .ne. 0) then + Y_e_plus(i,j,k) = 0.0d0 + Y_e_minus(i,j,k) = 0.0d0 + endif + + if (shift_state .ne. 0) then + lbetax(i,j,k) = betax(i,j,k) + lbetay(i,j,k) = betay(i,j,k) + lbetaz(i,j,k) = betaz(i,j,k) + else + lbetax(i,j,k) = 0.d0 + lbetay(i,j,k) = 0.d0 + lbetaz(i,j,k) = 0.d0 + end if + + enddo + enddo + enddo + !$OMP END PARALLEL DO !!$ 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 @@ -220,7 +168,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) 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 @@ -247,65 +194,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) 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 @@ -331,7 +219,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) 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 @@ -358,65 +245,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) 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 @@ -442,7 +270,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) 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 @@ -476,7 +303,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) deallocate(trivial_rp) deallocate(psi4, lbetax, lbetay, lbetaz) - deallocate(dum,dump,dumm) end subroutine GRHydro_PPMReconstruct_drv diff --git a/src/GRHydro_ParamCheck.F90 b/src/GRHydro_ParamCheck.F90 index f6177f7..2c29d13 100644 --- a/src/GRHydro_ParamCheck.F90 +++ b/src/GRHydro_ParamCheck.F90 @@ -112,5 +112,18 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS) call CCTK_PARAMWARN("Marquina solver is not implemented yet for MHD") end if + + if(clean_divergence.ne.0) then + if (CCTK_EQUALS(psidcspeed,"char speed")) then + whichpsidcspeed = 0 + else if (CCTK_EQUALS(psidcspeed,"light speed")) then + whichpsidcspeed = 1 + else if (CCTK_EQUALS(psidcspeed,"set speed")) then + whichpsidcspeed = 2 + else + call CCTK_PARAMWARN("Unknown type of speed setting for psidc (psidcspeed)") + end if + end if + end subroutine GRHydro_ParamCheck diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index a8ac2ec..ca93a4c 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -38,112 +38,71 @@ subroutine Reconstruction(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + CCTK_INT :: i,j,k CCTK_REAL :: local_min_tracer -!!$ Initialize variables that store reconstructed quantities - - !$OMP PARALLEL - !$OMP WORKSHARE - 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 - !$OMP END WORKSHARE NOWAIT - - if(evolve_mhd.ne.0) then - !$OMP WORKSHARE - Bvecxplus = 0.0d0 - Bvecxminus = 0.0d0 - Bvecyplus = 0.0d0 - Bvecyminus = 0.0d0 - Bveczplus = 0.0d0 - Bveczminus = 0.0d0 - !$OMP END WORKSHARE NOWAIT - if(clean_divergence.ne.0) then - !$OMP WORKSHARE - psidcplus = 0.0d0 - psidcminus = 0.0d0 - !$OMP END WORKSHARE NOWAIT - endif - endif - - if (evolve_tracer .ne. 0) then - !$OMP WORKSHARE - tracerplus = 0.0d0 - tracerminus = 0.0d0 - !$OMP END WORKSHARE NOWAIT - endif - - if (evolve_Y_e .ne. 0) then - !$OMP WORKSHARE - Y_e_plus = 0.0d0 - Y_e_minus = 0.0d0 - !$OMP END WORKSHARE NOWAIT - endif - !$OMP END PARALLEL - if (CCTK_EQUALS(recon_method,"tvd")) then - + ! this handles MHD and non-MHD call GRHydro_TVDReconstruct_drv(CCTK_PASS_FTOF) else if (CCTK_EQUALS(recon_method,"ppm")) then - call GRHydro_PPMReconstruct_drv(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call GRHydro_PPMMReconstruct_drv(CCTK_PASS_FTOF) + else + call GRHydro_PPMReconstruct_drv(CCTK_PASS_FTOF) + end if else if (CCTK_EQUALS(recon_method,"eno")) then - + ! this handles MHD and non-MHD call GRHydro_ENOReconstruct_drv(CCTK_PASS_FTOF) else call CCTK_WARN(0, "Reconstruction method not recognized!") end if - !$OMP PARALLEL WORKSHARE - where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min) ) - 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 PARALLEL WORKSHARE - if (evolve_tracer .ne. 0) then if (use_min_tracer .ne. 0) then local_min_tracer = min_tracer else local_min_tracer = 0.0d0 end if + end if - !$OMP PARALLEL WORKSHARE - where( (tracerplus .le. local_min_tracer).or.& - (tracerminus .le. local_min_tracer) ) - tracerplus = tracer - tracerminus = tracer - end where - !$OMP END PARALLEL 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 + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + if(rhoplus(i,j,k).lt.GRHydro_rho_min .or. & + rhominus(i,j,k).lt.GRHydro_rho_min) then + rhoplus(i,j,k) = rho(i,j,k) + rhominus(i,j,k) = rho(i,j,k) + epsplus(i,j,k) = eps(i,j,k) + epsminus(i,j,k) = eps(i,j,k) + velxplus(i,j,k) = vel(i,j,k,1) + velyplus(i,j,k) = vel(i,j,k,2) + velzplus(i,j,k) = vel(i,j,k,3) + velxminus(i,j,k) = vel(i,j,k,1) + velyminus(i,j,k) = vel(i,j,k,2) + velzminus(i,j,k) = vel(i,j,k,3) + end if + if(evolve_tracer.ne.0) then + where(tracerplus(i,j,k,:).le.local_min_tracer .or. & + tracerminus(i,j,k,:).le.local_min_tracer) + tracerplus(i,j,k,:) = tracer(i,j,k,:) + tracerminus(i,j,k,:) = tracer(i,j,k,:) + end where + end if + enddo + enddo + enddo + !$OMP END PARALLEL DO if (CCTK_EQUALS(recon_vars,"primitive").or.& CCTK_EQUALS(recon_method,"ppm")) then if(evolve_mhd.ne.0) then call primitive2conservativeM(CCTK_PASS_FTOF) - else call primitive2conservative(CCTK_PASS_FTOF) endif @@ -157,6 +116,12 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call CCTK_WARN(0,"Variable type to reconstruct not recognized.") end if + if(evolve_tracer.ne.0) then + ! 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) + end if + return end subroutine Reconstruction diff --git a/src/GRHydro_TVDReconstruct_drv.F90 b/src/GRHydro_TVDReconstruct_drv.F90 index 6fc2356..f5bf492 100644 --- a/src/GRHydro_TVDReconstruct_drv.F90 +++ b/src/GRHydro_TVDReconstruct_drv.F90 @@ -62,27 +62,22 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS) &type_bitsz, trivialz, not_trivialz CCTK_REAL, dimension(:,:,:),allocatable :: & - &psi4, lbetax, lbetay, lbetaz, dum, dump, dumm + &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 -!!$ 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) + STAT=ierr) if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with lbeta") @@ -108,32 +103,62 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS) 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 +!!$ Initialize variables that store reconstructed quantities: + + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + trivial_rp(i,j,k) = .false. + psi4(i,j,k) = 1.0d0 + rhoplus(i,j,k) = 0.0d0 + rhominus(i,j,k)= 0.0d0 + epsplus(i,j,k) = 0.0d0 + epsminus(i,j,k) = 0.0d0 + velxplus(i,j,k) = 0.0d0 + velxminus(i,j,k) = 0.0d0 + velyplus(i,j,k) = 0.0d0 + velyminus(i,j,k) = 0.0d0 + velzplus(i,j,k) = 0.0d0 + velzminus(i,j,k) = 0.0d0 + + if(evolve_mhd.ne.0) then + Bvecxplus(i,j,k) = 0.0d0 + Bvecxminus(i,j,k) = 0.0d0 + Bvecyplus(i,j,k) = 0.0d0 + Bvecyminus(i,j,k) = 0.0d0 + Bveczplus(i,j,k) = 0.0d0 + Bveczminus(i,j,k) = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus(i,j,k) = 0.0d0 + psidcminus(i,j,k) = 0.0d0 + endif + endif + + if (evolve_tracer .ne. 0) then + tracerplus(i,j,k,:) = 0.0d0 + tracerminus(i,j,k,:) = 0.0d0 + endif + + if (evolve_Y_e .ne. 0) then + Y_e_plus(i,j,k) = 0.0d0 + Y_e_minus(i,j,k) = 0.0d0 + endif + + if (shift_state .ne. 0) then + lbetax(i,j,k) = betax(i,j,k) + lbetay(i,j,k) = betay(i,j,k) + lbetaz(i,j,k) = betaz(i,j,k) + else + lbetax(i,j,k) = 0.d0 + lbetay(i,j,k) = 0.d0 + lbetaz(i,j,k) = 0.d0 + end if + enddo + enddo + enddo + !$OMP END PARALLEL DO !!$ TVD starts: if (evolve_tracer .ne. 0) then @@ -231,7 +256,6 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS) 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 93c6b5d..a4e0873 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -52,6 +52,7 @@ SRCS = Utils.F90 \ GRHydro_EigenproblemM.F90 \ GRHydro_FluxM.F90 \ GRHydro_HLLEM.F90 \ + GRHydro_PPMMReconstruct_drv.F90 \ GRHydro_PPMM.F90 \ GRHydro_Prim2ConM.F90 \ GRHydro_RiemannSolveM.F90 \ |