Move Evec calculation in seperate loops

this avoids accessing unitinialized values in the ghost zones when
computing the hydro varaible updates. This also saves time in running those loops since they are smaller.

From: Philipp Moesta <pmoesta@tapir.caltech.edu>

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

1 files changed, 40 insertions, 24 deletions

diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90
index 3b08360..eee7d14 100644
--- a/src/GRHydro_CalcUpdate.F90
+++ b/src/GRHydro_CalcUpdate.F90
@@ -47,10 +47,46 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
 
     if (use_weighted_fluxes == 0) then
 
+      if(evolve_mhd.ne.0 .and. transport_constraints.ne.0) then
+        !$OMP PARALLEL DO PRIVATE(i,j,k,alp_l,alp_r,alp_tmp)
+        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
+
+              alp_l = 0.5d0 * (alp(i,j,k) + &
+                   alp(i-xoffset,j-yoffset,k-zoffset))
+              alp_r = 0.5d0 * (alp(i,j,k) + &
+                   alp(i+xoffset,j+yoffset,k+zoffset))
+             
+              ! we have to first compute all components of v\crossB = E and
+              ! combine them in the last substep into Bconshs
+              ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k)
+              if(flux_direction.eq.1) then
+                 alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1))
+                 Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i  ,j  ,k+1))
+                 alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset))
+                 Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i  ,j+1,k  ))
+              elseif(flux_direction.eq.2) then
+                 alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1))
+                 Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i  ,j  ,k+1))
+                 alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset))
+                 Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j  ,k  ))
+              elseif(flux_direction.eq.3) then
+                 alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset))
+                 Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i  ,j+1,k  ))
+                 alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset))
+                 Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j  ,k  ))
+              end if
+            enddo
+          enddo
+        enddo
+        !$OMP END PARALLEL DO
+      endif
+ 
       !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,alp_tmp,Bcons_l,Bcons_r)
-      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
+      do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined
+        do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil
+          do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil
 
             alp_l = 0.5d0 * (alp(i,j,k) + &
                  alp(i-xoffset,j-yoffset,k-zoffset))
@@ -74,27 +110,7 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
                  alp_r * tauflux(i,j,k)) * idx
             
             if(evolve_mhd.ne.0) then
-               if(transport_constraints.ne.0) then
-                  ! we have to first compute all components of v\crossB = E and
-                  ! combine them in the last substep into Bconshs
-                  ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k)
-                  if(flux_direction.eq.1) then
-                     alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1))
-                     Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i  ,j  ,k+1))
-                     alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset))
-                     Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i  ,j+1,k  ))
-                  elseif(flux_direction.eq.2) then
-                     alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1))
-                     Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i  ,j  ,k+1))
-                     alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset))
-                     Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j  ,k  ))
-                  elseif(flux_direction.eq.3) then
-                     alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset))
-                     Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i  ,j+1,k  ))
-                     alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset))
-                     Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j  ,k  ))
-                  end if
-               else
+               if(transport_constraints.ne.1) then
                   Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + &
                        (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - &
                        alp_r * Bconsxflux(i,j,k)) * idx