expand region where fluxes are computed for constraint transport

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@307 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

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.