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.

1 files changed, 18 insertions, 2 deletions

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;
           }