diff options
Diffstat (limited to 'src/GRHydro_ReconstructPoly.F90')
| -rw-r--r-- | src/GRHydro_ReconstructPoly.F90 | 435 |
1 files changed, 358 insertions, 77 deletions
diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90 index 06c969b..cd735e1 100644 --- a/src/GRHydro_ReconstructPoly.F90 +++ b/src/GRHydro_ReconstructPoly.F90 @@ -20,6 +20,9 @@ #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) /*@@ @@ -56,7 +59,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) &type_bitsz, trivialz, not_trivialz CCTK_REAL, dimension(:,:,:),allocatable :: & - &psi4, lbetax, lbetay, lbetaz + &psi4, lbetax, lbetay, lbetaz, dum, dump, dumm !!$ CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: & !!$ &psi4, lbetax, lbetay, lbetaz @@ -72,7 +75,11 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) allocate(psi4(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& lbetax(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& lbetay(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& - lbetaz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + 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 @@ -125,6 +132,18 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) velyminus = 0.0d0 velzplus = 0.0d0 velzminus = 0.0d0 + if(evolve_mhd.ne.0) then + Bvecxplus = 0.0d0 + Bvecxminus = 0.0d0 + Bvecyplus = 0.0d0 + Bvecyminus = 0.0d0 + Bveczplus = 0.0d0 + Bveczminus = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus = 0.0d0 + psidcminus=0.0d0 + endif + endif if (evolve_tracer .ne. 0) then tracerplus = 0.0d0 @@ -178,6 +197,20 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") end if +!!$ B-field always needs reconstruction, as it is both prim and con + 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) + if(clean_divergence.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask) + endif + endif + !$OMP PARALLEL DO PRIVATE(i, j) do k = 1, nz do j = 1, ny @@ -207,31 +240,91 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) else if (CCTK_EQUALS(recon_method,"ppm")) then if (flux_direction == 1) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& - rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& - press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& - velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& - velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& - trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& - gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & - gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& - w_lorentz(:,j,k), & - 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & - GRHydro_mppm_eigenvalue_x_right, & - GRHydro_mppm_xwind) - do i = 1, nx - if (trivial_rp(i,j,k)) then - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) - else - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) - end if + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,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) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,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) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& + rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& + press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& + velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& + velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& + gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & + gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& + w_lorentz(:,j,k), & + 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & + GRHydro_mppm_eigenvalue_x_right, & + GRHydro_mppm_xwind) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 @@ -258,31 +351,91 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) end if else if (flux_direction == 2) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& - rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& - press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& - velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& - velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& - trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& - gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & - gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& - w_lorentz(j,:,k), & - 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & - GRHydro_mppm_eigenvalue_y_right, & - GRHydro_mppm_xwind) - do i = 1, ny - if (trivial_rp(j,i,k)) then - SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) - else - SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) - end if + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,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) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,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) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& + rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& + press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& + velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& + velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& + gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & + gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& + w_lorentz(j,:,k), & + 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & + GRHydro_mppm_eigenvalue_y_right, & + GRHydro_mppm_xwind) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 @@ -309,31 +462,91 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) end if else if (flux_direction == 3) then - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& - rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& - press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& - velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& - velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& - trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& - gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & - gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& - w_lorentz(j,k,:), & - 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & - GRHydro_mppm_eigenvalue_z_right, & - GRHydro_mppm_xwind) - do i = 1, nz - if (trivial_rp(j,k,i)) then - SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) - else - SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) - end if + if(evolve_mhd.ne.0) then + if(clean_divergence.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,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) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1dM(GRHydro_eos_handle,1,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) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + call SimplePPM_1d(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& + rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& + press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& + velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& + velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& + gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & + gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& + w_lorentz(j,k,:), & + 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & + GRHydro_mppm_eigenvalue_z_right, & + GRHydro_mppm_xwind) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif !evolve_mhd if(evolve_tracer.ne.0) then !$OMP PARALLEL DO PRIVATE(j) do k = GRHydro_stencil, ny - GRHydro_stencil + 1 @@ -415,6 +628,25 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + + if(evolve_mhd.ne.0) then +!!$ Always do B-field + 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)) + if(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 + endif + do i = 1, cctk_lsh(1) if (trivial_rp(i,j,k)) then SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) @@ -460,6 +692,24 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + + 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)) + if(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 + endif + do i = 1, cctk_lsh(2) if (trivial_rp(j,i,k)) then SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) @@ -505,6 +755,24 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Variable type to reconstruct not recognized.") !$OMP END CRITICAL end if + + 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,:)) + if(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 + endif + do i = 1, cctk_lsh(3) if (trivial_rp(j,k,i)) then SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) @@ -524,6 +792,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) deallocate(trivial_rp) deallocate(psi4, lbetax, lbetay, lbetaz) + deallocate(dum, dump, dumm) !$OMP WORKSHARE where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min).or.& @@ -563,17 +832,29 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive").or.& CCTK_EQUALS(recon_method,"ppm")) then if (use_eosgeneral == 0) then - call Prim2ConservativePolytype(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call Prim2ConservativePolytypeM(CCTK_PASS_FTOF) + else + call Prim2ConservativePolytype(CCTK_PASS_FTOF) + endif else - call primitive2conservativegeneral(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call primitive2conservativegeneral(CCTK_PASS_FTOF) + else + call primitive2conservativegeneral(CCTK_PASS_FTOF) + endif end if else if (CCTK_EQUALS(recon_vars,"conservative")) then - call Con2PrimBoundsPolytype(CCTK_PASS_FTOF) + if(evolve_mhd.ne.0) then + call Con2PrimBoundsPolytype(CCTK_PASS_FTOF) + else + call Con2PrimBoundsPolytype(CCTK_PASS_FTOF) + endif else - call CCTK_WARN(0,"Variable type to reconstruct not recognized.") + call CCTK_WARN(0,"Variable type to reconstruct not recognized.") end if - + return - + end subroutine ReconstructionPolytype |