CarpetReduce: Correct errors in calculating the weight function

10 files changed, 336 insertions, 100 deletions

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 <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <assert.h>
+
+
+
+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 <loopcontrol.h>
 
+#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<ibset,dim> const & boundaries =
+        vect<vect<ibset,2>,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<dim; ++d) {
-          assert ((notactive & boundaries[d]).empty());
+          for (int f=0; f<2; ++f) {
+            assert ((notactive & boundaries[d][f]).empty());
+          }
         }
         
         LOOP_OVER_BSET (cctkGH, notactive, box, imin, imax) {
@@ -106,34 +108,41 @@ namespace CarpetMask {
                    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
           {
             int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-            weight[ind] = 0.0;
+            iweight[ind] = 0;
           } LC_ENDLOOP3(CarpetMaskSetup_buffers);
           
         } END_LOOP_OVER_BSET;
         
         for (int d=0; d<dim; ++d) {
-          LOOP_OVER_BSET (cctkGH, boundaries[d], box, imin, imax) {
-            
-            if (verbose) {
-              ostringstream buf;
-              buf << "Setting prolongation boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione;
-              CCTK_INFO (buf.str().c_str());
-            }
-            
-            // Set weight on the boundary to 1/2
-            assert (dim == 3);
+          for (int f=0; f<2; ++f) {
+            LOOP_OVER_BSET (cctkGH, boundaries[d][f], box, imin, imax) {
+              
+              if (verbose) {
+                ostringstream buf;
+                buf << "Setting prolongation boundary on level " << reflevel << " direction " << d << " face " << f << " to weight 1/2: " << imin << ":" << imax-ione;
+                CCTK_INFO (buf.str().c_str());
+              }
+              
+              // Set weight on the boundary to 1/2
+              int bmask = 0;
+              for (unsigned b=0; b<BMSK(dim); ++b) {
+                if (BGET(b,d)==f) {
+                  bmask = BSET(bmask, b);
+                }
+              }
+              assert (dim == 3);
 #pragma omp parallel
-            LC_LOOP3(CarpetMaskSetup_prolongation,
-                     i,j,k,
-                     imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
-                     cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
-            {
-              int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-              weight[ind] *= 0.5;
-            } LC_ENDLOOP3(CarpetMaskSetup_prolongation);
-            
-            
-          } END_LOOP_OVER_BSET;
+              LC_LOOP3(CarpetMaskSetup_prolongation,
+                       i,j,k,
+                       imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+                       cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+              {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                iweight[ind] &= ~bmask;
+              } LC_ENDLOOP3(CarpetMaskSetup_prolongation);
+              
+            } END_LOOP_OVER_BSET;
+          } // for f
         } // for d
         
       } END_LOCAL_COMPONENT_LOOP;
@@ -147,7 +156,7 @@ namespace CarpetMask {
           
           ibset const & refined =
             dd.local_boxes.AT(mglevel).AT(reflevel).AT(local_component).restricted_region;
-          vect<ibset,dim> const & boundaries =
+          vect<vect<ibset,2>,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<int> mask (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]);
+          vector<int> imask (prod(ivect::ref(cctk_lsh)));
+          vector<int> 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<dim; ++d) {
-            LOOP_OVER_BSET (cctkGH, boundaries[d], box, imin, imax) {
-              
-              if (verbose) {
-                ostringstream buf;
-                buf << "Setting restriction boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione;
-                CCTK_INFO (buf.str().c_str());
-              }
-              
-              // Set weight on the boundary to 1/2
-              assert (dim == 3);
-#pragma omp parallel
-              LC_LOOP3(CarpetMaskSetup_restriction_boundary_partial,
-                       i,j,k,
-                       imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
-                       cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
-              {
-                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                if (mask[ind] == 0) {
-                  mask[ind] = 1;
+            for (int f=0; f<2; ++f) {
+              LOOP_OVER_BSET (cctkGH, boundaries[d][f], box, imin, imax) {
+                
+                if (verbose) {
+                  ostringstream buf;
+                  buf << "Setting restriction boundary on level " << reflevel << " direction " << d << " face " << f << " to weight 1/2: " << imin << ":" << imax-ione;
+                  CCTK_INFO (buf.str().c_str());
                 }
-                mask[ind] *= 2;
-              } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_partial);
-              
-            } END_LOOP_OVER_BSET;
+                
+                // Set weight on the boundary to 1/2
+                int bmask = 0;
+                for (unsigned b=0; b<BMSK(dim); ++b) {
+                  if (BGET(b,d)==f) {
+                    bmask = BSET(bmask, b);
+                  }
+                }
+                assert (dim == 3);
+#pragma omp parallel
+                LC_LOOP3(CarpetMaskSetup_restriction_boundary_partial,
+                         i,j,k,
+                         imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+                         cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+                {
+                  int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                  imask[ind] |= bmask;
+                  if (mask[ind] == 0) {
+                    mask[ind] = 1;
+                  }
+                  mask[ind] *= 2;
+                } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_partial);
+                
+              } END_LOOP_OVER_BSET;
+            } // for f
           } // for d
           
           assert (dim == 3);
@@ -222,10 +241,7 @@ namespace CarpetMask {
           {
             int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
             if (mask[ind] > 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 <loopcontrol.h>
 
+#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<BMSK(cctk_dim); ++b) {
+                if (BGET(b,d)==f) {
+                  bmask = BSET(bmask, b);
+                }
               }
 #pragma omp parallel
               LC_LOOP3(CoordBase_SetupMask_boundary2,
@@ -167,10 +179,8 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
                        bmin[0],bmin[1],bmin[2], bmax[0],bmax[1],bmax[2],
                        cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
               {
-                
                 int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                weight[ind] *= 0.5;
-                
+                iweight[ind] &= ~bmask;
               } LC_ENDLOOP3(CoordBase_SetupMask_boundary2);
               
             } /* if the domain is not empty */
diff --git a/Carpet/CarpetReduce/src/mask_init.c b/Carpet/CarpetReduce/src/mask_init.c
index 1ec321773..806a4b3d0 100644
--- a/Carpet/CarpetReduce/src/mask_init.c
+++ b/Carpet/CarpetReduce/src/mask_init.c
@@ -4,6 +4,8 @@
 
 #include <loopcontrol.h>
 
+#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 <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <loopcontrol.h>
+
+#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 <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <assert.h>
+#include <math.h>
+
+
+
+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<cctk_dim; ++d) {
+      domain_volume *= (physical_max[d] - physical_min[d]) / spacing[d];
+    }
+    
+    
+    
+    int const there_is_a_problem =
+      fabs(sum_weight - domain_volume) > 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);
+}