diff options
-rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 2 | ||||
-rw-r--r-- | src/GRHydro_ENOReconstruct_drv.F90 | 18 | ||||
-rw-r--r-- | src/GRHydro_HLLEM.F90 | 8 | ||||
-rw-r--r-- | src/GRHydro_PPMMReconstruct_drv.F90 | 18 | ||||
-rw-r--r-- | src/GRHydro_Prim2ConM.F90 | 16 | ||||
-rw-r--r-- | src/GRHydro_TVDReconstruct.F90 | 8 |
6 files changed, 45 insertions, 25 deletions
diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90 index e53f4b8..7d9408e 100644 --- a/src/GRHydro_CalcUpdate.F90 +++ b/src/GRHydro_CalcUpdate.F90 @@ -48,7 +48,7 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) if (use_weighted_fluxes == 0) then !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,Bvec_l,Bvec_r) - do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! mean we need to compute fluxes for one more point the GRHydro_stencil would indicate + do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90 index 58592f1..207086c 100644 --- a/src/GRHydro_ENOReconstruct_drv.F90 +++ b/src/GRHydro_ENOReconstruct_drv.F90 @@ -179,9 +179,11 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS) endif if (flux_direction == 1) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces !$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 + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + transport_constraints if (CCTK_EQUALS(recon_vars,"primitive")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),& @@ -259,9 +261,11 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS) end do !$OMP END PARALLEL DO else if (flux_direction == 2) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces !$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 + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + transport_constraints if (CCTK_EQUALS(recon_vars,"primitive")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),& @@ -339,9 +343,11 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS) end do !$OMP END PARALLEL DO else if (flux_direction == 3) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces !$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 + do k = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + transport_constraints if (CCTK_EQUALS(recon_vars,"primitive")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),& diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index 7e9c541..6f1791e 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -74,9 +74,11 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) call CCTK_WARN(0, "Flux direction not x,y,z") end if - do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil - do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil - do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + transport_constraints*(1-xoffset) f1 = 0.d0 lamminus = 0.d0 diff --git a/src/GRHydro_PPMMReconstruct_drv.F90 b/src/GRHydro_PPMMReconstruct_drv.F90 index 52c8162..852f447 100644 --- a/src/GRHydro_PPMMReconstruct_drv.F90 +++ b/src/GRHydro_PPMMReconstruct_drv.F90 @@ -181,9 +181,11 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) !!$ PPM starts: if (flux_direction == 1) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints 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),& @@ -252,9 +254,11 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) end if else if (flux_direction == 2) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints 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),& @@ -322,9 +326,11 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) end if else if (flux_direction == 3) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints 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,:),& diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90 index 87ae5e7..59acf0d 100644 --- a/src/GRHydro_Prim2ConM.F90 +++ b/src/GRHydro_Prim2ConM.F90 @@ -57,9 +57,11 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) @@ -275,11 +277,13 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr) - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) diff --git a/src/GRHydro_TVDReconstruct.F90 b/src/GRHydro_TVDReconstruct.F90 index ebe1a43..9abc8f2 100644 --- a/src/GRHydro_TVDReconstruct.F90 +++ b/src/GRHydro_TVDReconstruct.F90 @@ -52,10 +52,12 @@ subroutine tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & !!$ Initially all Riemann problems are NON-trivial trivial_rp = .false. + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces !$OMP PARALLEL DO PRIVATE(i,j,k,dupw,dloc,delta,ratio,hdelta) - do k = GRHydro_stencil, nz-GRHydro_stencil+1 - do j = GRHydro_stencil, ny-GRHydro_stencil+1 - do i = GRHydro_stencil, nx-GRHydro_stencil+1 + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints*(1-xoffset) if (GRHydro_enable_internal_excision /= 0 .and. & (hydro_excision_mask(i,j,k) .ne. 0)) then trivial_rp(i-xoffset, j-yoffset, k-zoffset) = .true. |