From 673d981e37b81792b1c6b304cc3078daebb6dbc4 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 1 Oct 2010 15:20:07 -0500 Subject: CarpetReduce: Correct errors in calculating the weight function --- Carpet/CarpetReduce/interface.ccl | 30 +++++++- Carpet/CarpetReduce/schedule.ccl | 67 +++++++++-------- Carpet/CarpetReduce/src/bits.h | 25 +++++++ Carpet/CarpetReduce/src/make.code.defn | 2 +- Carpet/CarpetReduce/src/mask_allocate.c | 24 ++++++ Carpet/CarpetReduce/src/mask_carpet.cc | 128 ++++++++++++++++++-------------- Carpet/CarpetReduce/src/mask_coords.c | 22 ++++-- Carpet/CarpetReduce/src/mask_init.c | 9 ++- Carpet/CarpetReduce/src/mask_set.c | 34 +++++++++ Carpet/CarpetReduce/src/mask_test.c | 95 ++++++++++++++++++++++++ 10 files changed, 336 insertions(+), 100 deletions(-) create mode 100644 Carpet/CarpetReduce/src/bits.h create mode 100644 Carpet/CarpetReduce/src/mask_allocate.c create mode 100644 Carpet/CarpetReduce/src/mask_set.c create mode 100644 Carpet/CarpetReduce/src/mask_test.c (limited to 'Carpet/CarpetReduce') diff --git a/Carpet/CarpetReduce/interface.ccl b/Carpet/CarpetReduce/interface.ccl index 73a0ed7ef..641ede55a 100644 --- a/Carpet/CarpetReduce/interface.ccl +++ b/Carpet/CarpetReduce/interface.ccl @@ -14,12 +14,14 @@ uses include header: typeprops.hh uses include header: loopcontrol.h -CCTK_INT FUNCTION \ + + +CCTK_INT FUNCTION \ SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH) REQUIRES FUNCTION SymmetryTableHandleForGrid -CCTK_INT FUNCTION \ - MultiPatch_GetMap \ +CCTK_INT FUNCTION \ + MultiPatch_GetMap \ (CCTK_POINTER_TO_CONST IN cctkGH) USES FUNCTION MultiPatch_GetMap @@ -42,6 +44,28 @@ CCTK_INT FUNCTION \ CCTK_INT OUT ARRAY shiftout) REQUIRES FUNCTION GetBoundarySpecification +CCTK_INT FUNCTION \ + GetDomainSpecification \ + (CCTK_INT IN size, \ + CCTK_REAL OUT ARRAY physical_min, \ + CCTK_REAL OUT ARRAY physical_max, \ + CCTK_REAL OUT ARRAY interior_min, \ + CCTK_REAL OUT ARRAY interior_max, \ + CCTK_REAL OUT ARRAY exterior_min, \ + CCTK_REAL OUT ARRAY exterior_max, \ + CCTK_REAL OUT ARRAY spacing) +USES FUNCTION GetDomainSpecification + +# Get current refinement level +CCTK_INT FUNCTION \ + GetRefinementLevel \ + (CCTK_POINTER_TO_CONST IN cctkGH) +REQUIRES FUNCTION GetRefinementLevel + +INT iweight TYPE=gf TAGS='prolongation="none" InterpNumTimelevels=1 checkpoint="no"' "Integer weight mask, using 2^D bits" + REAL weight TYPE=gf TAGS='prolongation="none" InterpNumTimelevels=1 checkpoint="no"' "Weight function" + +REAL one TYPE=gf TAGS='prolongation="none" InterpNumTimelevels=1 checkpoint="no"' "Constant one" diff --git a/Carpet/CarpetReduce/schedule.ccl b/Carpet/CarpetReduce/schedule.ccl index 7610f1093..a0b6614a2 100644 --- a/Carpet/CarpetReduce/schedule.ccl +++ b/Carpet/CarpetReduce/schedule.ccl @@ -7,7 +7,7 @@ schedule CarpetReduceStartup at STARTUP -# This might move to MaskBase +# Should this move to a new thorn MaskBase? STORAGE: weight @@ -28,53 +28,56 @@ SCHEDULE GROUP MaskBase_SetupMask AT post_recover_variables { } "Set up the weight function" -SCHEDULE MaskBase_InitMask IN MaskBase_SetupMask + + +SCHEDULE GROUP MaskBase_SetupMaskAll IN MaskBase_SetupMask +{ + STORAGE: iweight + STORAGE: one +} "Set up the weight function" + + + +SCHEDULE MaskBase_AllocateMask IN MaskBase_SetupMaskAll +{ + LANG: C + OPTIONS: global +} "Allocate the weight function" + +SCHEDULE MaskBase_InitMask IN MaskBase_SetupMaskAll AFTER MaskBase_AllocateMask { LANG: C OPTIONS: global loop-local } "Initialise the weight function" -SCHEDULE GROUP SetupMask IN MaskBase_SetupMask AFTER MaskBase_InitMask +SCHEDULE GROUP SetupMask IN MaskBase_SetupMaskAll AFTER MaskBase_InitMask { } "Set up the weight function (schedule other routines in here)" -# This might move to CoordBase +SCHEDULE MaskBase_SetMask IN MaskBase_SetupMaskAll AFTER SetupMask +{ + LANG: C + OPTIONS: global loop-local +} "Set the weight function" + +SCHEDULE MaskBase_TestMask IN MaskBase_SetupMaskAll AFTER MaskBase_SetMask +{ + LANG: C + OPTIONS: global +} "Test the weight function" + + + +# Should this move to CoordBase? SCHEDULE CoordBase_SetupMask IN SetupMask { LANG: C OPTIONS: global loop-local } "Set up the outer boundaries of the weight function" -# This might move to CarpetMask +# Should this move to CarpetMask? SCHEDULE CarpetMaskSetup IN SetupMask { LANG: C OPTIONS: global loop-singlemap } "Set up the weight function for the restriction regions" - - - -#SCHEDULE GROUP MaskBase_SetupMask_LevelMode AT basegrid AFTER (SpatialCoordinates SphericalSurface_Setup) -#{ -#} "Set up the weight function" -# -#SCHEDULE MaskBase_InitMask IN MaskBase_SetupMask_LevelMode -#{ -# LANG: C -#} "Initialise the weight function" -# -#SCHEDULE GROUP SetupMask_LevelMode IN MaskBase_SetupMask_LevelMode AFTER MaskBase_InitMask -#{ -#} "Set up the weight function (schedule other routines in here)" -# -## This might move to CoordBase -#SCHEDULE CoordBase_SetupMask IN SetupMask_LevelMode -#{ -# LANG: C -#} "Set up the outer boundaries of the weight function" -# -## This might move to CarpetMask -#SCHEDULE CarpetMaskSetup IN SetupMask_LevelMode -#{ -# LANG: C -#} "Set up the weight function for the restriction regions" diff --git a/Carpet/CarpetReduce/src/bits.h b/Carpet/CarpetReduce/src/bits.h new file mode 100644 index 000000000..65bd450da --- /dev/null +++ b/Carpet/CarpetReduce/src/bits.h @@ -0,0 +1,25 @@ +#ifndef BITS_H +#define BITS_H + +/* a mask with exactly bit n set */ +#define BMSK(n) (1U << (n)) + +/* whether bit n in value v is set */ +#define BGET(v,n) (!!((v) & BMSK(n))) + +/* set, clear, or invert bit n in value v */ +#define BSET(v,n) ((v) | BMSK(n)) +#define BCLR(v,n) ((v) & ~BMSK(n)) +#define BINV(v,n) ((v) ^ BMSK(n)) + +/* set bit n in value v to b */ +#define BCPY(v,b,n) ((b) ? BSET(v,n) : BCLR(v,n)) + +/* count the number of set bits */ +#define BCNT(x) (((BX_(x)+(BX_(x)>>4)) & 0x0F0F0F0F) % 0xFF) +#define BX_(x) ((x) \ + - (((x)>>1)&0x77777777) \ + - (((x)>>2)&0x33333333) \ + - (((x)>>3)&0x11111111)) + +#endif /* #ifndef BITS_H */ diff --git a/Carpet/CarpetReduce/src/make.code.defn b/Carpet/CarpetReduce/src/make.code.defn index d68722520..35ae77a52 100644 --- a/Carpet/CarpetReduce/src/make.code.defn +++ b/Carpet/CarpetReduce/src/make.code.defn @@ -1,7 +1,7 @@ # Main make.code.defn file for thorn CarpetReduce # Source files in this directory -SRCS = mask_carpet.cc mask_coords.c mask_init.c reduce.cc +SRCS = mask_allocate.c mask_init.c mask_set.c mask_test.c mask_carpet.cc mask_coords.c reduce.cc # Subdirectories containing source files SUBDIRS = diff --git a/Carpet/CarpetReduce/src/mask_allocate.c b/Carpet/CarpetReduce/src/mask_allocate.c new file mode 100644 index 000000000..24029942d --- /dev/null +++ b/Carpet/CarpetReduce/src/mask_allocate.c @@ -0,0 +1,24 @@ +#include +#include +#include + +#include + + + +void +MaskBase_AllocateMask (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + /* Allocate helpers for the weight function */ + if (verbose) { + CCTK_INFO ("Allocating weight function helpers"); + } + + int const ierr1 = CCTK_EnableGroupStorage (cctkGH, "CarpetReduce::iweight"); + int const ierr2 = CCTK_EnableGroupStorage (cctkGH, "CarpetReduce::one"); + assert (!ierr1); + assert (!ierr2); +} diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc index 3e41f1846..6db04bc9a 100644 --- a/Carpet/CarpetReduce/src/mask_carpet.cc +++ b/Carpet/CarpetReduce/src/mask_carpet.cc @@ -11,6 +11,7 @@ #include +#include "bits.h" #include "mask_carpet.hh" @@ -62,7 +63,6 @@ namespace CarpetMask { ivect const ione = ivect(1); - gh const & hh = *vhh.AT(Carpet::map); dh const & dd = *vdd.AT(Carpet::map); ivect const reffact = @@ -80,13 +80,15 @@ namespace CarpetMask { ibset const & active = dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).active; - vect const & boundaries = + vect,dim> const & boundaries = dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).prolongation_boundaries; ibset const notactive = ext - active; for (int d=0; d const & boundaries = + vect,dim> const & boundaries = dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).restriction_boundaries; LOOP_OVER_BSET (cctkGH, refined, box, imin, imax) { @@ -167,13 +176,13 @@ namespace CarpetMask { cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - weight[ind] = 0; + iweight[ind] = 0; } LC_ENDLOOP3(CarpetMaskSetup_restriction); } END_LOOP_OVER_BSET; - assert (dim == 3); - vector mask (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]); + vector imask (prod(ivect::ref(cctk_lsh))); + vector mask (prod(ivect::ref(cctk_lsh))); assert (dim == 3); #pragma omp parallel @@ -183,34 +192,44 @@ namespace CarpetMask { cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + imask[ind] = 0; mask[ind] = 0; } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_init); for (int d=0; d 0) { -#if 0 - weight[ind] *= 1.0 - 1.0 / mask[ind]; -#endif - weight[ind] -= 1.0 / mask[ind]; + iweight[ind] &= imask[ind]; } } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_apply); diff --git a/Carpet/CarpetReduce/src/mask_coords.c b/Carpet/CarpetReduce/src/mask_coords.c index 6749df74d..b018dbce4 100644 --- a/Carpet/CarpetReduce/src/mask_coords.c +++ b/Carpet/CarpetReduce/src/mask_coords.c @@ -6,6 +6,8 @@ #include +#include "bits.h" + void @@ -30,6 +32,10 @@ CoordBase_SetupMask (CCTK_ARGUMENTS) + int const reflevel = GetRefinementLevel(cctkGH); + + + if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) { int const m = MultiPatch_GetMap (cctkGH); assert (m >= 0); @@ -113,7 +119,7 @@ CoordBase_SetupMask (CCTK_ARGUMENTS) /* Loop over the boundary */ if (verbose) { CCTK_VInfo (CCTK_THORNSTRING, - "Setting boundary points in direction %d face %d to weight 0", d, f); + "Setting boundary points in direction %d face %d to weight 0 on level %d", d, f, reflevel); } #pragma omp parallel LC_LOOP3(CoordBase_SetupMask_boundary, @@ -123,7 +129,7 @@ CoordBase_SetupMask (CCTK_ARGUMENTS) { int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - weight[ind] = 0.0; + iweight[ind] = 0; } LC_ENDLOOP3(CoordBase_SetupMask); @@ -159,7 +165,13 @@ CoordBase_SetupMask (CCTK_ARGUMENTS) /* Loop over the points next to boundary */ if (verbose) { CCTK_VInfo (CCTK_THORNSTRING, - "Setting non-staggered boundary points in direction %d face %d to weight 1/2", d, f); + "Setting non-staggered boundary points in direction %d face %d to weight 1/2 on level %d", d, f, reflevel); + } + int bmask = 0; + for (unsigned b=0; b +#include "bits.h" + void @@ -14,9 +16,12 @@ MaskBase_InitMask (CCTK_ARGUMENTS) /* Initialise the weight to 1 everywhere */ if (verbose) { - CCTK_INFO ("Initialising weight to 1"); + int const reflevel = GetRefinementLevel(cctkGH); + CCTK_VInfo (CCTK_THORNSTRING, + "Initialising weight to 1 on level %d", reflevel); } + int const allbits = BMSK(BMSK(cctk_dim)) - 1; #pragma omp parallel LC_LOOP3(MaskBase_InitMask_interior, i,j,k, @@ -24,6 +29,6 @@ MaskBase_InitMask (CCTK_ARGUMENTS) cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - weight[ind] = 1.0; + iweight[ind] = allbits; } LC_ENDLOOP3(MaskBase_InitMask_interior); } diff --git a/Carpet/CarpetReduce/src/mask_set.c b/Carpet/CarpetReduce/src/mask_set.c new file mode 100644 index 000000000..2444bda3d --- /dev/null +++ b/Carpet/CarpetReduce/src/mask_set.c @@ -0,0 +1,34 @@ +#include +#include +#include + +#include + +#include "bits.h" + + + +void +MaskBase_SetMask (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (verbose) { + int const reflevel = GetRefinementLevel(cctkGH); + CCTK_VInfo (CCTK_THORNSTRING, + "Finalise the weight on level %d", reflevel); + } + + CCTK_REAL const factor = 1.0 / BMSK(cctk_dim); +#pragma omp parallel + LC_LOOP3(MaskBase_InitMask_interior, + i,j,k, + 0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + weight[ind] = factor * BCNT(iweight[ind]); + one[ind] = 1.0; + } LC_ENDLOOP3(MaskBase_InitMask_interior); +} diff --git a/Carpet/CarpetReduce/src/mask_test.c b/Carpet/CarpetReduce/src/mask_test.c new file mode 100644 index 000000000..cbf57b3b3 --- /dev/null +++ b/Carpet/CarpetReduce/src/mask_test.c @@ -0,0 +1,95 @@ +#include +#include +#include + +#include +#include + + + +void +MaskBase_TestMask (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (verbose) { + CCTK_INFO ("Testing weight"); + } + + + + int const sum = CCTK_ReductionHandle ("sum"); + assert (sum >= 0); + + int const proc = 0; + + int const weight_var = CCTK_VarIndex ("CarpetReduce::weight"); + int const one_var = CCTK_VarIndex ("CarpetReduce::one"); + assert (weight_var >= 0); + + CCTK_REAL sum_weight; + + { + int const ierr = CCTK_Reduce (cctkGH, + proc, + sum, + 1, CCTK_VARIABLE_REAL, &sum_weight, + 1, one_var); + assert (ierr >= 0); + } + + + + if (CCTK_MyProc(cctkGH) == proc) { + + CCTK_REAL physical_min[cctk_dim]; + CCTK_REAL physical_max[cctk_dim]; + CCTK_REAL interior_min[cctk_dim]; + CCTK_REAL interior_max[cctk_dim]; + CCTK_REAL exterior_min[cctk_dim]; + CCTK_REAL exterior_max[cctk_dim]; + CCTK_REAL spacing [cctk_dim]; + int const ierr = GetDomainSpecification (cctk_dim, + physical_min, + physical_max, + interior_min, + interior_max, + exterior_min, + exterior_max, + spacing); + assert (!ierr); + + CCTK_REAL domain_volume = 1.0; + for (int d=0; d 1.0e-12 * (sum_weight + domain_volume); + + + + if (verbose || there_is_a_problem) { + CCTK_VInfo (CCTK_THORNSTRING, + "Simulation domain volume: %.15g", (double)domain_volume); + CCTK_VInfo (CCTK_THORNSTRING, + "Reduction weight sum: %.15g", (double)sum_weight); + } + + if (there_is_a_problem) { + CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Simulation domain volume and reduction weight sum differ"); + } + + } + + + + int const iret1 = CCTK_DisableGroupStorage (cctkGH, "CarpetReduce::iweight"); + int const iret2 = CCTK_DisableGroupStorage (cctkGH, "CarpetReduce::one"); + assert (iret1==1); + assert (iret2==1); +} -- cgit v1.2.3