RIT GRMHD dev:

1) Further splitting of PPM Reconstruction routine into magnetic
and non-magnetic part.

2) Merge the divergence cleaning loop into the non-divergence one,
since profiling indicated no substantial difference between the
two cases. Indeed it was a little bit better after the merger
(0.14%), but that it is not significant. We decided then to apply 
Occam's razor and choose the simplest form. Branch prediction 
seems to work fine in this case.

3) Move reconstruction initialization statements to their repective 
drivers.

4) Get rid of WORKSHARE in the reconstruction routines. Profiling 
showed a 4.16% performance improvement for the hydro ppm reconstruction 
routine when using 1 processor and 2 threads only. Expect a 
better improvement for a larger number of threads.

5) Introduce a parameter to control characteristic speeds
for psidc in HLLE.



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

1 files changed, 56 insertions, 27 deletions

diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90
index 656fd6e..65ff19c 100644
--- a/src/GRHydro_ENOReconstruct_drv.F90
+++ b/src/GRHydro_ENOReconstruct_drv.F90
@@ -67,8 +67,6 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS)
 
   CCTK_INT :: ierr
 
-  CCTK_REAL :: local_min_tracer
-
   allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr)
   if (ierr .ne. 0) then
     call CCTK_WARN(0, "Allocation problems with trivial_rp")
@@ -107,31 +105,62 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS)
   ny = cctk_lsh(2)
   nz = cctk_lsh(3)
 
-  !$OMP PARALLEL
-  !$OMP WORKSHARE
-  trivial_rp = .false.
-  !$OMP END WORKSHARE NOWAIT
-
-!!$  Currently only option is reconstruction on primitive variables.
-!!$  Should change this.
-
-  !$OMP WORKSHARE
-  psi4 = 1.d0
-  !$OMP END WORKSHARE NOWAIT
-  if (shift_state .ne. 0) then
-    !$OMP WORKSHARE
-    lbetax = betax
-    lbetay = betay
-    lbetaz = betaz
-    !$OMP END WORKSHARE NOWAIT
-  else
-    !$OMP WORKSHARE
-    lbetax = 0.d0
-    lbetay = 0.d0
-    lbetaz = 0.d0
-    !$OMP END WORKSHARE NOWAIT
-  end if
-  !$OMP END PARALLEL
+!!$ Initialize variables that store reconstructed quantities:
+
+  !$OMP PARALLEL DO PRIVATE(i,j,k)
+  do k=1,cctk_lsh(3)
+     do j=1,cctk_lsh(2)
+        do i=1,cctk_lsh(1)
+           trivial_rp(i,j,k) = .false.
+           psi4(i,j,k) = 1.0d0
+           rhoplus(i,j,k) = 0.0d0
+           rhominus(i,j,k)= 0.0d0
+           epsplus(i,j,k) = 0.0d0
+           epsminus(i,j,k) = 0.0d0
+           velxplus(i,j,k) = 0.0d0
+           velxminus(i,j,k) = 0.0d0
+           velyplus(i,j,k) = 0.0d0
+           velyminus(i,j,k) = 0.0d0
+           velzplus(i,j,k) = 0.0d0
+           velzminus(i,j,k) = 0.0d0
+
+           if(evolve_mhd.ne.0) then
+             Bvecxplus(i,j,k) = 0.0d0
+             Bvecxminus(i,j,k) = 0.0d0
+             Bvecyplus(i,j,k) = 0.0d0
+             Bvecyminus(i,j,k) = 0.0d0
+             Bveczplus(i,j,k) = 0.0d0
+             Bveczminus(i,j,k) = 0.0d0
+             if(clean_divergence.ne.0) then
+                psidcplus(i,j,k) = 0.0d0
+                psidcminus(i,j,k) = 0.0d0
+             endif
+           endif
+
+           if (evolve_tracer .ne. 0) then
+              tracerplus(i,j,k,:) = 0.0d0
+              tracerminus(i,j,k,:) = 0.0d0
+           endif
+
+           if (evolve_Y_e .ne. 0) then
+              Y_e_plus(i,j,k) = 0.0d0
+              Y_e_minus(i,j,k) = 0.0d0
+           endif
+
+           if (shift_state .ne. 0) then
+             lbetax(i,j,k) = betax(i,j,k)
+             lbetay(i,j,k) = betay(i,j,k)
+             lbetaz(i,j,k) = betaz(i,j,k)
+           else
+             lbetax(i,j,k) = 0.d0
+             lbetay(i,j,k) = 0.d0
+             lbetaz(i,j,k) = 0.d0
+           end if
+
+        enddo
+     enddo
+  enddo
+  !$OMP END PARALLEL DO
 
 !!$ ENO starts:
      if (evolve_tracer .ne. 0) then