CarpetMask: Keep track of the volume that is masked out

Keep track of the volume that is masked out by CarpetMask, and take
this volume into account when checking in CarpetReduce that the
integral over the simulation domain equals the domain volume.

3 files changed, 40 insertions, 4 deletions

diff --git a/Carpet/CarpetMask/interface.ccl b/Carpet/CarpetMask/interface.ccl
index 7297a1c46..9cd31f3bd 100644
--- a/Carpet/CarpetMask/interface.ccl
+++ b/Carpet/CarpetMask/interface.ccl
@@ -3,3 +3,5 @@
 IMPLEMENTS: CarpetMask
 
 INHERITS: grid SphericalSurface
+
+USES INCLUDE HEADER: CarpetReduce_bits.h
diff --git a/Carpet/CarpetMask/src/mask_excluded.cc b/Carpet/CarpetMask/src/mask_excluded.cc
index 1b1c0beb5..aec5ffa6c 100644
--- a/Carpet/CarpetMask/src/mask_excluded.cc
+++ b/Carpet/CarpetMask/src/mask_excluded.cc
@@ -8,6 +8,8 @@
 
 #include "mask_excluded.hh"
 
+#include <CarpetReduce_bits.h>
+
 
 
 namespace CarpetMask {
@@ -26,15 +28,27 @@ namespace CarpetMask {
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
     
-    CCTK_INT * const iweight =
+    CCTK_INT * restrict const iweight =
       static_cast <CCTK_INT *>
       (CCTK_VarDataPtr (cctkGH, 0, "CarpetReduce::iweight"));
     
-    if (not iweight) {
+    CCTK_REAL * restrict const excised_cells =
+      static_cast <CCTK_REAL *>
+      (CCTK_VarDataPtr (cctkGH, 0, "CarpetReduce::excised_cells"));
+
+    if (not iweight or not excised_cells) {
       CCTK_WARN (CCTK_WARN_ABORT,
                  "CarpetReduce is not active, or CarpetReduce::iweight does not have storage");
     }
     
+    // Volume of a grid cell on this level, in terms of coarse grid
+    // cells
+    CCTK_REAL const cell_volume =
+      1.0 / (cctk_levfac[0] * cctk_levfac[1] * cctk_levfac[2]);
+    
+    unsigned const bits = BMSK(cctk_dim);
+    CCTK_REAL const factor = 1.0 / bits;
+    
     for (int n = 0; n < 10; ++ n) {
       
       CCTK_REAL const r0 = excluded_radius[n];
@@ -60,6 +74,8 @@ namespace CarpetMask {
               dx2 + dy2 + dz2 >= r2 :
               dx2 + dy2 + dz2 <= r2)
           {
+            // Tally up the weight we are removing
+            * excised_cells += cell_volume * factor * BCNT(iweight[ind]);
             iweight[ind] = 0;
           }
           
diff --git a/Carpet/CarpetMask/src/mask_surface.cc b/Carpet/CarpetMask/src/mask_surface.cc
index 02d80f9f8..f2a592b3c 100644
--- a/Carpet/CarpetMask/src/mask_surface.cc
+++ b/Carpet/CarpetMask/src/mask_surface.cc
@@ -11,6 +11,8 @@
 
 #include "mask_surface.hh"
 
+#include <CarpetReduce_bits.h>
+
 
 
 namespace CarpetMask {
@@ -67,15 +69,27 @@ namespace CarpetMask {
       last_output.resize (num_excluded, -1);
     }
     
-    CCTK_INT * const iweight =
+    CCTK_INT * restrict const iweight =
       static_cast <CCTK_INT *>
       (CCTK_VarDataPtr (cctkGH, 0, "CarpetReduce::iweight"));
     
-    if (not iweight) {
+    CCTK_REAL * restrict const excised_cells =
+      static_cast <CCTK_REAL *>
+      (CCTK_VarDataPtr (cctkGH, 0, "CarpetReduce::excised_cells"));
+    
+    if (not iweight or not excised_cells) {
       CCTK_WARN (CCTK_WARN_ABORT,
                  "CarpetReduce is not active, or CarpetReduce::iweight does not have storage");
     }
     
+    // Volume of a grid cell on this level, in terms of coarse grid
+    // cells
+    CCTK_REAL const cell_volume =
+      1.0 / (cctk_levfac[0] * cctk_levfac[1] * cctk_levfac[2]);
+    
+    unsigned const bits = BMSK(cctk_dim);
+    CCTK_REAL const factor = 1.0 / bits;
+    
     for (int n = 0; n < num_excluded; ++ n) {
       
       int const sn = excluded_surface[n];
@@ -127,6 +141,8 @@ namespace CarpetMask {
               sqrt (pow (dx, 2) + pow (dy, 2) + pow (dz, 2));
             if (rho < 1.0e-12) {
               // Always excise the surface origin
+              // Tally up the weight we are removing
+              * excised_cells += cell_volume * factor * BCNT(iweight[ind]);
               iweight[ind] = 0;
             } else {
               CCTK_REAL theta =
@@ -171,6 +187,8 @@ namespace CarpetMask {
               CCTK_REAL const dr =
                 sf_radius[a + maxntheta * (b + maxnphi * sn)];
               if (rho <= dr * shrink_factor) {
+                // Tally up the weight we are removing
+                * excised_cells += cell_volume * factor * BCNT(iweight[ind]);
                 iweight[ind] = 0;
               }
             }