From d822d7044d70a9fd43c8be48b0c5bdce51789ed9 Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 5 Feb 2008 22:44:56 +0000 Subject: 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 --- interface.ccl | 6 +++- schedule.ccl | 10 ++++++- src/cg.F90 | 84 ++++++++++++++++++++---------------------------------- src/copy_mask.cc | 47 ++++++++++++++++++++++++++++++ src/copy_mask.hh | 11 +++++++ src/make.code.defn | 2 +- 6 files changed, 104 insertions(+), 56 deletions(-) create mode 100644 src/copy_mask.cc create mode 100644 src/copy_mask.hh 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 + +#include +#include +#include + +#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_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 +#include + +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 = -- cgit v1.2.3