Change the scheme for reductions. Before, reductions was done over the full

grid functions using CCTK_Reduce. Now, reductions are done alongside the
computations into local arrays and then these are reduced using
CCTK_ReduceLocArrayToArray1D. This is much more efficient. In addition
OpenMP directives are added to allow for mixed MPI/OpenMP usage.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/NoExcision/trunk@15 4ec1db94-0e4f-0410-ada3-8bed251432c9

6 files changed, 104 insertions, 56 deletions

diff --git a/interface.ccl b/interface.ccl
index 968f4cf..4f64e13 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -18,7 +18,6 @@ CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER_TO_CONST IN GH, \
                                             CCTK_INT IN table_handle, \
                                             CCTK_STRING IN var_name, \
                                             CCTK_STRING IN bc_name)
-
 USES FUNCTION Boundary_SelectGroupForBC
 
 CCTK_INT smask type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
@@ -26,6 +25,11 @@ CCTK_INT smask type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation=
   nes_mask
 } "mask for smoothing"
 
+CCTK_REAL reduction_mask type=GF tags='Prolongation="None" InterpNumTimelevels=1 checkpoint="no"'
+{
+  red_mask
+} "Copy of the weight grid function from CarpetReduce"
+
 CCTK_REAL cg_res_metric type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
 {
   resgxx, resgxy, resgxz, resgyy, resgyz, resgzz
diff --git a/schedule.ccl b/schedule.ccl
index 824d112..28e4603 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,6 +1,8 @@
 # Schedule definitions for thorn NoExcision
 # $Header$
 
+STORAGE: reduction_mask
+
 SCHEDULE NoExcision_Overwrite IN ADMBase_PostInitial
 {
   LANG: Fortran
@@ -39,7 +41,13 @@ if (smooth_regions) {
 
     } "Conjugate Gradient smoothing"
 
-    SCHEDULE NoExcision_CGInit_1 IN NoExcision_CGSmoothing
+    SCHEDULE CopyMask IN NoExcision_CGSmoothing
+    {
+      LANG: C
+      OPTIONS: global loop-local
+    } "Copy the weight function from CarpetReduce"
+
+    SCHEDULE NoExcision_CGInit_1 IN NoExcision_CGSmoothing AFTER CopyMask
     {
       LANG: Fortran
     } "Initialise the conjugate gradient method 1"
diff --git a/src/cg.F90 b/src/cg.F90
index ce95c93..316fe68 100644
--- a/src/cg.F90
+++ b/src/cg.F90
@@ -92,7 +92,8 @@ subroutine NoExcision_CGInit_1 (CCTK_ARGUMENTS)
                     res, resx, resy, resz, &
                     redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
                     redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
-                    red, redx, redy, redz, nes_mask )
+                    red, redx, redy, redz, nes_mask, red_mask, &
+                    .true., smoothing_order )
                      
     call CCTK_ReductionArrayHandle ( sum_handle, 'sum' )
     if ( sum_handle .lt. 0 ) then
@@ -129,27 +130,17 @@ subroutine NoExcision_CGInit_2 (CCTK_ARGUMENTS)
   character(len=56) :: conv_message
 
   if ( loop_control == 1 ) then
-    do i = 1, 16
-      if ( cont(i) ) then
-        call CCTK_VarIndex ( res_var, red_names(i) )
-        if ( res_var < 0 ) then
-          call CCTK_WARN ( 0, 'Could not get index to grid function red' )
-        end if
-
-        ! delta_new = r^T r.
-
-        call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
-                           CCTK_VARIABLE_REAL, delta_new(i), 1, res_var )
+    ! delta_new = r^T r.
+    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, sum_handle, lsumred, &
+                                        delta_new, 16, CCTK_VARIABLE_REAL)
+    if ( ierr < 0 ) call CCTK_WARN ( 0, 'Could not perform reduction of local 1D array' )
 
-        call CCTK_Reduce ( ierr, cctkGH, -1, infnorm_handle, 1, &
-                           CCTK_VARIABLE_REAL, infnormresid(i), 1, res_var )
+    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, infnorm_handle, &
+                                        linfred, infnormresid, 16, &
+                                        CCTK_VARIABLE_REAL)
+    if ( ierr < 0 ) call CCTK_WARN ( 0, 'Could not perform reduction of local 1D array' )
 
-        ! Since we are reducing r*r, we have to take the square root.
-
-        infnormresid(i) = sqrt(infnormresid(i))
-
-      end if
-    end do
+    where ( cont ) infnormresid = sqrt(infnormresid)
 
     ! Check if some variables has already converged. This happens when the
     ! variable is identically zero.
@@ -209,7 +200,8 @@ subroutine NoExcision_CG_1 (CCTK_ARGUMENTS)
                   q, qx, qy, qz, &
                   redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
                   redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
-                  red, redx, redy, redz, nes_mask )
+                  red, redx, redy, redz, nes_mask, red_mask,&
+                  .false., smoothing_order )
 
   sym_selector = 2
 
@@ -231,22 +223,12 @@ subroutine NoExcision_CG_2 (CCTK_ARGUMENTS)
   integer :: i
 
   if ( loop_control == 1 ) then
-    do i = 1, 16
-      if ( cont(i) ) then
-        call CCTK_VarIndex ( res_var, red_names(i) )
-        if ( res_var < 0 ) then
-          call CCTK_WARN ( 0, 'Could not get index to grid function red' )
-        end if
-
-        ! alpha = delta_new / ( d^T A d ).
-        call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
-                           CCTK_VARIABLE_REAL, alpha(i), 1, res_var )
-      
-        alpha(i) = delta_new(i) / alpha(i)
-
-      end if
-    end do
+    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, sum_handle, lsumred, &
+                                        alpha, 16, CCTK_VARIABLE_REAL)
+    if ( ierr < 0 ) call CCTK_WARN ( 0, 'Could not perform reduction of local 1D array' )
 
+    ! alpha = delta_new / ( d^T A d ).
+    where ( cont ) alpha = delta_new / alpha
   end if
 
 end subroutine NoExcision_CG_2
@@ -315,7 +297,8 @@ subroutine NoExcision_CG_3 (CCTK_ARGUMENTS)
                   res, resx, resy, resz, &
                   redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
                   redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
-                  red, redx, redy, redz, nes_mask )
+                  red, redx, redy, redz, nes_mask, red_mask, &
+                  .true., smoothing_order )
 
   sym_selector = 3
 
@@ -337,25 +320,20 @@ subroutine NoExcision_CG_4 (CCTK_ARGUMENTS)
   integer :: i
 
   if ( loop_control == 1 ) then
-    do i = 1, 16
-      if ( cont(i) ) then
-        call CCTK_VarIndex ( res_var, red_names(i) )
-        if ( res_var < 0 ) then
-          call CCTK_WARN ( 0, 'Could not get index to grid function red' )
-        end if
-
-        ! delta_new = r^T r.
-        call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
-                           CCTK_VARIABLE_REAL, delta_new(i), 1, res_var )
-
-        call CCTK_Reduce ( ierr, cctkGH, -1, infnorm_handle, 1, &
-                           CCTK_VARIABLE_REAL, infnormresid(i), 1, res_var )
-        infnormresid(i) = sqrt(infnormresid(i))
-      end if
-    end do
+    ! delta_new = r^T r.
+    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, sum_handle, lsumred, &
+                                        delta_new, 16, CCTK_VARIABLE_REAL)
+    if ( ierr < 0 ) call CCTK_WARN ( 0, 'Could not perform reduction of local 1D array' )
 
+    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, infnorm_handle, &
+                                        linfred, infnormresid, 16, &
+                                        CCTK_VARIABLE_REAL)
+    if ( ierr < 0 ) call CCTK_WARN ( 0, 'Could not perform reduction of local 1D array' )
+
+    where ( cont ) infnormresid = sqrt(infnormresid)
   end if
 
+
 end subroutine NoExcision_CG_4
 
 
diff --git a/src/copy_mask.cc b/src/copy_mask.cc
new file mode 100644
index 0000000..4dd0538
--- /dev/null
+++ b/src/copy_mask.cc
@@ -0,0 +1,47 @@
+#include <cmath>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include "copy_mask.hh"
+
+namespace NoExcision {
+
+#ifdef HAVE_CARPET
+  using namespace Carpet;
+#endif
+
+  /**
+   * Modify the mask according to the CarpetReduce mask in order to be able to do the local reductions.
+   */
+
+  extern "C" 
+    void CopyMask (CCTK_ARGUMENTS)
+    {
+      DECLARE_CCTK_ARGUMENTS;
+      DECLARE_CCTK_PARAMETERS;
+
+      CCTK_REAL * const weight =
+        static_cast <CCTK_REAL *>
+        (CCTK_VarDataPtr (cctkGH, 0, "CarpetReduce::weight"));
+
+      if (not weight) {
+        CCTK_WARN (CCTK_WARN_ABORT,
+                   "CarpetReduce is not active, or CarpetReduce::mask does not have storage");
+      }
+
+      for (int k = 0; k < cctk_lsh[2]; ++ k) {
+        for (int j = 0; j < cctk_lsh[1]; ++ j) {
+          for (int i = 0; i < cctk_lsh[0]; ++ i) {
+            int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+
+            red_mask[ind] = weight[ind];
+          }
+        }
+      }
+
+    }
+
+
+} // namespace NoExcision
diff --git a/src/copy_mask.hh b/src/copy_mask.hh
new file mode 100644
index 0000000..bb016f9
--- /dev/null
+++ b/src/copy_mask.hh
@@ -0,0 +1,11 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+
+namespace CarpetMask {
+
+  extern "C" {
+    void
+    CopyMask (CCTK_ARGUMENTS);
+  }
+
+} // namespace CarpetMask
diff --git a/src/make.code.defn b/src/make.code.defn
index a174c46..f72509a 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = NoExcision_mod.F90 overwrite.F90 smooth.F90 reduce.F90 reduce.c getlevelinfo.cc cg.F90
+SRCS = NoExcision_mod.F90 overwrite.F90 smooth.F90 reduce.F90 reduce.c getlevelinfo.cc copy_mask.cc cg.F90
 
 # Subdirectories containing source files
 SUBDIRS =