[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 571)

rhaas at tapir.caltech.edu rhaas at tapir.caltech.edu
Tue Aug 13 09:56:21 CDT 2013


User: rhaas
Date: 2013/08/13 09:56 AM

Modified:
 /trunk/src/
  GRHydro_CalcUpdate.F90

Log:
 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 at tapir.caltech.edu>

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_CalcUpdate.F90
Delta lines: +40 -24
===================================================================
--- trunk/src/GRHydro_CalcUpdate.F90	2013-08-13 14:56:16 UTC (rev 570)
+++ trunk/src/GRHydro_CalcUpdate.F90	2013-08-13 14:56:21 UTC (rev 571)
@@ -47,10 +47,46 @@
 
     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 @@
                  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



More information about the Commits mailing list