TestBBoxSet2: New thorn

9 files changed, 980 insertions, 0 deletions

diff --git a/CarpetExtra/TestBBoxSet2/README b/CarpetExtra/TestBBoxSet2/README
new file mode 100644
index 000000000..e6a6d5935
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/README
@@ -0,0 +1,10 @@
+Cactus Code Thorn TestBBoxSet2
+Author(s)    : Erik Schnetter <schnetter@gmail.com>
+Maintainer(s): Erik Schnetter <schnetter@gmail.com>
+Licence      : GPL
+--------------------------------------------------------------------------
+
+1. Purpose
+
+Test the new bboxset class by comparing its results to the old bboxset
+class.
diff --git a/CarpetExtra/TestBBoxSet2/configuration.ccl b/CarpetExtra/TestBBoxSet2/configuration.ccl
new file mode 100644
index 000000000..bbc806dc2
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/configuration.ccl
@@ -0,0 +1,3 @@
+# Configuration definitions for thorn TestBBoxSet2
+
+REQUIRES CarpetLib CycleClock
diff --git a/CarpetExtra/TestBBoxSet2/doc/documentation.tex b/CarpetExtra/TestBBoxSet2/doc/documentation.tex
new file mode 100644
index 000000000..77fa7f337
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header$
+%
+%  Thorn documentation in the latex file doc/documentation.tex
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include
+%  pages about variables, parameters and scheduling parsed from the
+%  relevant thorn CCL files.
+%
+%  This template contains guidelines which help to assure that your
+%  documentation will be correctly added to ThornGuides. More
+%  information is available in the Cactus UsersGuide.
+%
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date, etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma.
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations',
+%     'references', and 'image' names should conform to the following
+%     convention:
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package.
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This
+%     will allow us to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command.
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as
+%     standard sections.
+%   - For the benefit of our Perl scripts, and for future extensions,
+%     please use simple latex.
+%
+% *======================================================================*
+%
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================*
+
+% If you are using CVS use this line to give version information
+% $Header$
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Erik Schnetter \textless schnetter@gmail.com\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{TestBBoxSet2}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+%    \date{$ $Date: 2004-01-07 14:12:39 -0600 (Wed, 07 Jan 2004) $ $}
+\date{March 30 2013}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/CarpetExtra/TestBBoxSet2/interface.ccl b/CarpetExtra/TestBBoxSet2/interface.ccl
new file mode 100644
index 000000000..d3df6b03e
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/interface.ccl
@@ -0,0 +1,7 @@
+# Interface definition for thorn TestBBoxSet2
+
+IMPLEMENTS: TestBBoxSet2
+
+USES INCLUDES HEADER: cycleclock.h
+
+USES INCLUDE HEADER: bboxset.hh
diff --git a/CarpetExtra/TestBBoxSet2/param.ccl b/CarpetExtra/TestBBoxSet2/param.ccl
new file mode 100644
index 000000000..91cfad826
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/param.ccl
@@ -0,0 +1,10 @@
+# Parameter definitions for thorn TestBBoxSet2
+
+BOOLEAN verbose "Provide verbose output" STEERABLE=always
+{
+} "yes"
+
+INT components_goal "Approximate number of components to simulate" STEERABLE=always
+{
+  1:* :: ""
+} 10000
diff --git a/CarpetExtra/TestBBoxSet2/schedule.ccl b/CarpetExtra/TestBBoxSet2/schedule.ccl
new file mode 100644
index 000000000..01a173a31
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/schedule.ccl
@@ -0,0 +1,7 @@
+# Schedule definitions for thorn TestBBoxSet2
+
+SCHEDULE TestBBoxSet2_test AT paramcheck
+{
+  LANG: C
+  OPTIONS: meta
+} "Test bboxset2"
diff --git a/CarpetExtra/TestBBoxSet2/src/make.code.defn b/CarpetExtra/TestBBoxSet2/src/make.code.defn
new file mode 100644
index 000000000..d403969cf
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/src/make.code.defn
@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn TestBBoxSet2
+
+# Source files in this directory
+SRCS = test.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/CarpetExtra/TestBBoxSet2/src/test.cc b/CarpetExtra/TestBBoxSet2/src/test.cc
new file mode 100644
index 000000000..8ab99d13e
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/src/test.cc
@@ -0,0 +1,786 @@
+#include <cctk_Parameters.h>
+
+#include <cycleclock.h>
+
+#include <bboxset.hh>
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <iostream>
+#include <set>
+#include <sstream>
+#include <string>
+#include <vector>
+using namespace std;
+
+
+
+namespace {
+  
+  
+  
+  const double nano = 1.0e-9;
+  
+  // Time a region of code
+  template<typename F>
+  auto xtime(double& time, const F& f) ->
+    typename enable_if<is_void<decltype(f())>::value, void>::type
+  {
+    const ticks tb = getticks();
+    f();
+    const ticks te = getticks();
+    time += seconds_per_tick() * elapsed(te, tb);
+  }
+  
+  template<typename F>
+  auto xtime(double& time, const F& f) ->
+    typename enable_if<not is_void<decltype(f())>::value, decltype(f())>::type
+  {
+    const ticks tb = getticks();
+    const decltype(f()) ret = f();
+    const ticks te = getticks();
+    time += seconds_per_tick() * elapsed(te, tb);
+    return ret;
+  }
+  
+  
+  
+  // Convert to string
+  // (Would like to use to_string instead, but this doesn't seem to
+  // compile on Stampede)
+  template<typename T>
+  string to_str(const T& x)
+  {
+    ostringstream buf;
+    buf << x;
+    return buf.str();
+  }
+  
+  
+  
+  template<typename A1, typename A2, typename T, int D>
+  void check_equal(const string descr,
+                   const A1& x1,
+                   const A2& x2,
+                   const bboxset1::bboxset<T,D>& bset1,
+                   const bboxset2::bboxset<T,D>& bset2)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    if (verbose) {
+      cout << "  <" << D << "> " << descr << "...";
+      cout.flush();
+    }
+    
+    set<bbox<T,D>> v1, v2;
+    bset1.serialise(v1);
+    bset2.serialise(v2);
+    
+    bool eq = v1 == v2;
+    if (not eq) {
+      eq = (bset1 == bboxset1::bboxset<T,D>(v2) and
+            bset2 == bboxset2::bboxset<T,D>(v1));
+    }
+    
+    if (eq) {
+      if (verbose) cout << " [success]\n";
+    } else {
+      if (verbose) cout << " [failure]\n";
+      cout.flush();
+      cerr.flush();
+      fflush(stdout);
+      fflush(stderr);
+      CCTK_VParamWarn(CCTK_THORNSTRING,
+                      "Found difference in operation \"%s\":", descr.c_str());
+      fflush(stdout);
+      fflush(stderr);
+      cout << "\n"
+           << "input1=" << x1 << "\n"
+           << "input2=" << x2 << "\n"
+           << "bset1=" << bset1 << "\n"
+           << "bset2=" << bset2 << "\n"
+           << "bset1.serialised=" << v1 << "\n"
+           << "bset2.serialised=" << v2 << "\n";
+      cout.flush();
+      cerr.flush();
+      
+      const bboxset1::bboxset<T,D> bs(v2);
+      cout << "bset1-bset2=" << bset1-bs << "\n"
+           << "bset2-bset1=" << bs-bset1 << "\n";
+    }
+  }
+  
+  template<typename T, int D>
+  void check_equal(const string descr,
+                   const bboxset1::bboxset<T,D>& bset1,
+                   const bboxset2::bboxset<T,D>& bset2)
+  {
+    check_equal(descr, "", "", bset1, bset2);
+  }
+  
+  
+  
+#define bvect bvect_
+#define b2vect b2vect_
+#define ivect ivect_
+#define i2vect i2vect_
+#define ibbox ibbox_
+  template<int D>
+  struct full_boxes_t {
+    typedef ::vect<bool,D> bvect;
+    typedef ::vect<bvect,2> b2vect;
+    typedef ::vect<int,D> ivect;
+    typedef ::vect<ivect,2> i2vect;
+    typedef ::bbox<int,D> ibbox;
+    typedef bboxset1::bboxset<int,D> ibset1;
+    typedef bboxset2::bboxset<int,D> ibset2;
+    ibbox interior;
+    b2vect is_outer_boundary;
+    ibbox exterior;
+    ibset1 ghosts1;
+    ibset2 ghosts2;
+    ibbox communicated;
+    ibset1 outer_boundaries1;
+    ibset2 outer_boundaries2;
+    ibbox owned;
+    ibset1 boundaries1;
+    ibset2 boundaries2;
+    ibset1 buffers1;
+    ibset2 buffers2;
+    ibset1 overlaps1;
+    ibset2 overlaps2;
+    ibset1 active1;
+    ibset2 active2;
+  };
+  template<int D>
+  struct level_boxes_t {
+    typedef ::vect<bool,D> bvect;
+    typedef ::vect<bvect,2> b2vect;
+    typedef ::vect<int,D> ivect;
+    typedef ::vect<ivect,2> i2vect;
+    typedef ::bbox<int,D> ibbox;
+    typedef bboxset1::bboxset<int,D> ibset1;
+    typedef bboxset2::bboxset<int,D> ibset2;
+    ibset1 buffers1;
+    ibset2 buffers2;
+    ibset1 active1;
+    ibset2 active2;
+  };
+  
+  template<int D>
+  void test()
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    CCTK_VInfo(CCTK_THORNSTRING, "Testing bboxset2<%d>...", D);
+    
+    // Prevent warnings about shadowing declarations
+    typedef ::vect<bool,D> bvect;
+    typedef ::vect<bvect,2> b2vect;
+    typedef ::vect<int,D> ivect;
+    typedef ::vect<ivect,2> i2vect;
+    typedef ::bbox<int,D> ibbox;
+    typedef bboxset1::bboxset<int,D> ibset1;
+    typedef bboxset2::bboxset<int,D> ibset2;
+    
+    const int str = 3;          // a non-trivial stride
+    
+    double time1, time2, time1a, time2a, time_check;
+    
+    
+    
+    // Some warm-up gymnastics with emtpy sets
+    
+    ibset1 bs1;
+    ibset2 bs2;
+    xtime(time_check, [&]() { check_equal("create ibset", bs1, bs2); });
+    assert(bs1.empty());
+    assert(bs2.empty());
+    
+    ibbox b;
+    bs1 |= b;
+    bs2 |= b;
+
+    xtime(time_check, [&]() {
+        check_equal("union with empty ibbox", bs1, bs2); });
+    assert(bs1.empty());
+    assert(bs2.empty());
+    
+    bs1 &= b;
+    bs2 &= b;
+    xtime(time_check, [&]() {
+        check_equal("intersection with empty ibbox", bs1, bs2); });
+    assert(bs1.empty());
+    assert(bs2.empty());
+    
+    bs1 += b;
+    bs2 += b;
+    xtime(time_check, [&]() {
+        check_equal("symmetric union with empty ibbox", bs1, bs2); });
+    assert(bs1.empty());
+    assert(bs2.empty());
+    
+    bs1 -= b;
+    bs2 -= b;
+    xtime(time_check, [&]() {
+        check_equal("difference with empty ibbox", bs1, bs2); });
+    assert(bs1.empty());
+    assert(bs2.empty());
+    
+    
+    
+    // Some basic tests
+    
+    b = ibbox(ivect(0*str), ivect(9*str), ivect(str));
+    bs1 += b;
+    bs2 += b;
+    xtime(time_check, [&]() { check_equal("add ibbox", bs1, bs2); });
+    
+    bs1 -= b;
+    bs2 -= b;
+    xtime(time_check, [&]() { check_equal("subtract ibbox", bs1, bs2); });
+    
+    bs1 |= b;
+    bs2 |= b;
+    xtime(time_check, [&]() { check_equal("union", bs1, bs2); });
+    
+    bs1 &= b;
+    bs2 &= b;
+    xtime(time_check, [&]() { check_equal("intersection", bs1, bs2); });
+    
+    ibset1 bs1a = bs1.shift(ivect(1));
+    ibset2 bs2a = bs2.shift(ivect(1));
+    xtime(time_check, [&]() { check_equal("shift", bs1a, bs2a); });
+    
+    ibset1 bs1b = bs1.expand(ivect(1), ivect(1));
+    ibset2 bs2b = bs2.expand(ivect(1), ivect(1));
+    xtime(time_check, [&]() { check_equal("expand", bs1b, bs2b); });
+    
+    ibset1 bs1c = bs1 & bs1a;
+    ibset2 bs2c = bs2 & bs2a;
+    xtime(time_check, [&]() {
+        check_equal("non-trivial intersection", bs1c, bs2c); });
+    
+    ibset1 bs1d = bs1b - bs1;
+    ibset2 bs2d = bs2b - bs2;
+    xtime(time_check, [&]() {
+        check_equal("non-trivial difference", bs1d, bs2d); });
+    
+    
+    
+    // Many tests with random sets
+    
+    bs1 = ibset1();
+    bs2 = ibset2();
+    xtime(time_check, [&]() {
+        check_equal("many dense bboxes (init)", bs1, bs2); });
+    
+    time1 = time2 = time_check = 0.0;
+    for (int n=0; n<40; ++n) {
+      ivect lo, up;
+      for (int d=0; d<D; ++d) lo[d] = random() % 10 * str;
+      for (int d=0; d<D; ++d) up[d] = lo[d] + random() % 10 * str;
+      b = ibbox(lo, up, ivect(str));
+      xtime(time1, [&](){ bs1 |= b; });
+      xtime(time2, [&](){ bs2 |= b; });
+      xtime(time_check, [&]() {
+          check_equal("many dense bboxes (union) #" + to_str(n),
+                      bs1, bs2); });
+    }
+    cout << "  "
+         << "time1: " << time1 << " s, "
+         << "time2: " << time2 << " s, "
+         << "time_check: " << time_check << " s\n";
+    
+    time1 = time2 = time_check = 0.0;
+    for (int n=0; n<100; ++n) {
+      ivect lo, up;
+      for (int d=0; d<D; ++d) lo[d] = random() % 10 * str;
+      for (int d=0; d<D; ++d) up[d] = lo[d] + random() % 10 * str;
+      b = ibbox(lo, up, ivect(str));
+      xtime(time1, [&](){ bs1 -= b; });
+      xtime(time2, [&](){ bs2 -= b; });
+      xtime(time_check, [&]() {
+          check_equal("many dense bboxes (difference) #" + to_str(n),
+                      bs1, bs2); });
+    }
+    cout << "  "
+         << "time1: " << time1 << " s, "
+         << "time2: " << time2 << " s, "
+         << "time_check: " << time_check << " s\n";
+    
+    bs1 = ibset1();
+    bs2 = ibset2();
+    xtime(time_check, [&]() {
+        check_equal("many sparse bboxes (init)", bs1, bs2); });
+    
+    time1 = time2 = time_check = 0.0;
+    for (int n=0; n<10; ++n) {
+      ivect lo, up;
+      for (int d=0; d<D; ++d) lo[d] = random() % 1000 * str;
+      for (int d=0; d<D; ++d) up[d] = lo[d] + random() % 1000 * str;
+      b = ibbox(lo, up, ivect(str));
+      xtime(time1, [&](){ bs1 |= b; });
+      xtime(time2, [&](){ bs2 |= b; });
+      xtime(time_check, [&]() {
+          check_equal("many sparse bboxes (union) #" + to_str(n),
+                      bs1, bs2); });
+    }
+    cout << "  "
+         << "time1: " << time1 << " s, "
+         << "time2: " << time2 << " s, "
+         << "time_check: " << time_check << " s\n";
+    
+    time1 = time2 = time_check = 0.0;
+    for (int n=0; n<20; ++n) {
+      ivect lo, up;
+      for (int d=0; d<D; ++d) lo[d] = random() % 1000 * str;
+      for (int d=0; d<D; ++d) up[d] = lo[d] + random() % 1000 * str;
+      b = ibbox(lo, up, ivect(str));
+      xtime(time1, [&](){ bs1 -= b; });
+      xtime(time2, [&](){ bs2 -= b; });
+      xtime(time_check, [&]() {
+          check_equal("many sparse bboxes (difference) #" + to_str(n),
+                      bs1, bs2); });
+    }
+    cout << "  "
+         << "time1: " << time1 << " s, "
+         << "time2: " << time2 << " s, "
+         << "time_check: " << time_check << " s\n";
+    
+    
+    
+    // Test shifting, expanding, contracting
+    
+    int sec=0;
+    for (int ts=1; ts<=3; ++ts) {
+      const ivect target_stride(1<<ts);
+      for (int to=0; to<target_stride[0]; ++to) {
+        const ivect target_offset(to);
+        const ibbox target(target_offset, target_offset, target_stride);
+        
+        for (int s=1; s<=3; ++s) {
+          const ivect stride(1<<s);
+          for (int o=0; o<stride[0]; ++o) {
+            const ivect offset(o);
+            
+            bs1 = ibset1();
+            bs2 = ibset2();
+            for (int n=0; n<10; ++n) {
+              ivect lo, up;
+              for (int d=0; d<D; ++d) {
+                lo[d] = offset[d] + random() % 10 * stride[d];
+                up[d] = lo[d] + random() % 10 * stride[d];
+              }
+              b = ibbox(lo, up, stride);
+              bs1 |= b;
+              bs2 |= b;
+              xtime(time_check, [&]() {
+                  check_equal("repeated test #" +
+                              to_str(sec) + "." + to_str(n),
+                              bs1, bs2); });
+              
+              bs1a = bs1.expanded_for(target);
+              bs2a = bs2.expanded_for(target);
+              xtime(time_check, [&]() {
+                  check_equal("expanded_for #" +
+                              to_str(sec) + "." + to_str(n),
+                              bs1, target, bs1a, bs2a); });
+              bs1b = bs1.contracted_for(target);
+              bs2b = bs2.contracted_for(target);
+              // Omitting check since bboxset1 has known bugs here
+              // xtime(time_check, [&]() { check_equal("contracted_for", bs1, target, bs1b, bs2b); });
+              
+            } // n
+            ++sec;
+            
+          } // offset
+        }   // stride
+        
+      } // target_offset
+    }   // target_stride
+    
+    
+    
+    // Some real-world tests
+    
+    {
+      
+      // const auto minus1 =
+      //   (ibset1(*)(const ibbox&, const ibbox&))bboxset1::operator-;
+      // const auto minus2 =
+      //   (ibset2(*)(const ibbox&, const ibbox&))bboxset2::operator-;
+      
+      // ibset1(*const minus1)(const ibbox&, const ibbox&)(bboxset1::operator-);
+      // ibset2(*const minus2)(const ibbox&, const ibbox&)(bboxset2::operator-);
+      
+      typedef ibset1 (minus1_t)(const ibbox&, const ibbox&);
+      typedef ibset2 (minus2_t)(const ibbox&, const ibbox&);
+      // typedef ibset2 (&xor2_t)(const ibbox&, const ibbox&);
+      minus1_t& minus1(bboxset1::operator-);
+      minus2_t& minus2(bboxset2::operator-);
+      // const xor2_t xor2(bboxset2::operator^);
+      
+      const i2vect ghost_width(3);
+      const i2vect buffer_width(9);
+      const i2vect overlap_width(0);
+      const i2vect boundary_width(3);
+      
+      const int components_per_direction = lrint(pow(components_goal, 1.0/D));
+      const int components = lrint(pow(components_per_direction, D));
+      const ivect origin(0);
+      const ivect stride(1);
+      const ivect component_size(30);
+      
+      cout << "  components: " << components << "\n";
+      
+      // Domain
+      
+      time1 = time2 = time1a = time2a = time_check = 0.0;
+      
+      const ibbox domain_exterior
+        (origin,
+         origin +
+         stride * component_size * ivect(components_per_direction) - stride,
+         stride);
+      
+      const ibbox domain_active = domain_exterior.expand(-boundary_width);
+      assert(domain_active <= domain_exterior);
+      ibset1 domain_boundary1;
+      ibset2 domain_boundary2;
+      xtime(time1, [&]() {
+          domain_boundary1 = minus1(domain_exterior, domain_active); });
+      xtime(time2, [&]() {
+          domain_boundary2 = minus2(domain_exterior, domain_active); });
+      xtime(time_check, [&]() {
+          check_equal("domain_boundary",
+                      domain_boundary1, domain_boundary2); });
+      
+      cout << "  Domain:\n"
+           << "    time1:        " << time1 << " s\n"
+           << "    time2:        " << time2 << " s\n"
+           << "    time1_assert: " << time1a << " s\n"
+           << "    time2_assert: " << time2a << " s\n"
+           << "    time_check:   " << time_check << " s\n";
+      
+      // Region
+      
+      time1 = time2 = time1a = time2a = time_check = 0.0;
+      
+      vector<full_boxes_t<D>> full_boxes(components);
+      level_boxes_t<D> level_boxes;
+      for (int c=0; c<components; ++c) {
+        
+        // Interior:
+        ibbox& intr = full_boxes.AT(c).interior;
+        ivect ipos;
+        int ic = c;
+        for (int d=0; d<D; ++d) {
+          ipos[d] = ic % components_per_direction;
+          ic /= components_per_direction;
+        }
+        assert(ic==0);
+        intr = ibbox(origin + stride * component_size * ipos,
+                     origin + stride * component_size * (ipos+1) - stride,
+                     stride);
+        assert(intr <= domain_exterior);
+        for (int cc=0; cc<c; ++cc) {
+          assert(not intr.intersects(full_boxes.AT(cc).interior));
+        }
+        
+        // Outer boundary faces:
+        b2vect& is_outer_boundary = full_boxes.AT(c).is_outer_boundary;
+        is_outer_boundary[0] = intr.lower() == domain_exterior.lower(); 
+        is_outer_boundary[1] = intr.upper() == domain_exterior.upper(); 
+        
+        // Exterior:
+        ibbox& extr = full_boxes.AT(c).exterior;
+        extr = intr.expand(i2vect(not is_outer_boundary) * ghost_width);
+        assert(extr <= domain_exterior);
+        
+        // Cactus ghost zones (which include outer boundaries):
+        ibset1& ghosts1 = full_boxes.AT(c).ghosts1;
+        ibset2& ghosts2 = full_boxes.AT(c).ghosts2;
+        xtime(time1, [&]() { ghosts1 = minus1(extr, intr); });
+        xtime(time2, [&]() { ghosts2 = minus2(extr, intr); });
+        xtime(time_check, [&]() { check_equal("ghosts", ghosts1, ghosts2); });
+        xtime(time1a, [&]() { assert(ghosts1 <= domain_exterior); });
+        xtime(time2a, [&]() { assert(ghosts2 <= domain_exterior); });
+        
+        // Communicated region:
+        ibbox& comm = full_boxes.AT(c).communicated;
+        comm = extr.expand(i2vect(is_outer_boundary) * (-boundary_width));
+        assert(comm <= domain_active);
+        
+        // Outer boundary:
+        ibset1& outer_boundaries1 = full_boxes.AT(c).outer_boundaries1;
+        ibset2& outer_boundaries2 = full_boxes.AT(c).outer_boundaries2;
+        xtime(time1, [&]() { outer_boundaries1 = minus1(extr, comm); });
+        xtime(time2, [&]() { outer_boundaries2 = minus2(extr, comm); });
+        xtime(time_check, [&]() {
+            check_equal("outer_boundaries",
+                        outer_boundaries1, outer_boundaries2); });
+        xtime(time1a, [&]() { assert(outer_boundaries1 <= domain_boundary1); });
+        xtime(time2a, [&]() { assert(outer_boundaries2 <= domain_boundary2); });
+        
+        // Owned region:
+        ibbox& owned = full_boxes.AT(c).owned;
+        owned = intr.expand(i2vect(is_outer_boundary) * (-boundary_width));
+        assert(owned <= domain_active);
+        for (int cc=0; cc<c; ++cc) {
+          assert(not owned.intersects(full_boxes.AT(cc).owned));
+        }
+        
+        // Boundary (Carpet ghost zones, which do not include outer
+        // boundaries):
+        ibset1& boundaries1 = full_boxes.AT(c).boundaries1;
+        ibset2& boundaries2 = full_boxes.AT(c).boundaries2;
+        xtime(time1, [&]() { boundaries1 = minus1(comm, owned); });
+        xtime(time2, [&]() { boundaries2 = minus2(comm, owned); });
+        xtime(time_check, [&]() {
+            check_equal("boundaries", boundaries1, boundaries2); });
+        xtime(time1a, [&]() { assert(boundaries1 <= domain_active); });
+        xtime(time2a, [&]() { assert(boundaries2 <= domain_active); });
+        
+      } // for c
+      
+      cout << "  Region:\n"
+           << "    time1:        " << time1 << " s\n"
+           << "    time2:        " << time2 << " s\n"
+           << "    time1_assert: " << time1a << " s\n"
+           << "    time2_assert: " << time2a << " s\n"
+           << "    time_check:   " << time_check << " s\n";
+      
+      // Buffer zones
+      
+      time1 = time2 = time1a = time2a = time_check = 0.0;
+      
+      // All owned regions
+      ibset1 allowned1;
+      ibset2 allowned2;
+      xtime(time1, [&]() {
+          allowned1 = ibset1(full_boxes, &full_boxes_t<D>::owned); });
+      xtime(time2, [&]() {
+          allowned2 = ibset2(full_boxes, &full_boxes_t<D>::owned); });
+      xtime(time_check, [&]() {
+          check_equal("allowned", allowned1, allowned2); });
+      xtime(time1a, [&]() { assert(allowned1 <= domain_active); });
+      xtime(time2a, [&]() { assert(allowned2 <= domain_active); });
+      
+      // All not-owned regions
+      ibset1 notowned1;
+      ibset2 notowned2;
+      xtime(time1, [&]() { notowned1 = domain_active - allowned1; });
+      xtime(time2, [&]() { notowned2 = domain_active - allowned2; });
+      xtime(time_check, [&]() {
+          check_equal("notowned", notowned1, notowned2); });
+      
+      // All not-active points
+      ibset1 notactive1;
+      ibset2 notactive2;
+      xtime(time1, [&]() {
+          notactive1 = notowned1.expand(buffer_width + overlap_width); });
+      xtime(time2, [&]() {
+          notactive2 = notowned2.expand(buffer_width + overlap_width); });
+      xtime(time_check, [&]() {
+          check_equal("notactive", notactive1, notactive2); });
+      
+      // All not-active points, in stages
+      int const num_substeps =
+        any(any(ghost_width == 0)) ?
+        0 :
+        minval(minval(buffer_width / ghost_width));
+      assert(all(all(buffer_width == num_substeps * ghost_width)));
+      vector<ibset1> notactive_stepped1(num_substeps+1);
+      vector<ibset2> notactive_stepped2(num_substeps+1);
+      notactive_stepped1.AT(0) = notowned1;
+      notactive_stepped2.AT(0) = notowned2;
+      for (int substep=1; substep<=num_substeps; ++substep) {
+        xtime(time1, [&]() {
+            notactive_stepped1.AT(substep) =
+              notactive_stepped1.AT(substep-1).expand(ghost_width); });
+        xtime(time2, [&]() {
+            notactive_stepped2.AT(substep) =
+              notactive_stepped2.AT(substep-1).expand(ghost_width); });
+        xtime(time_check, [&]() {
+            check_equal("notactive_stepped[" + to_str(substep) + "]",
+                        notactive_stepped1.AT(substep),
+                        notactive_stepped2.AT(substep)); });
+      }
+      ibset1 notactive_overlaps1;
+      ibset2 notactive_overlaps2;
+      xtime(time1a, [&]() {
+          notactive_overlaps1 =
+            notactive_stepped1.AT(num_substeps).expand(overlap_width); });
+      xtime(time2a, [&]() {
+          notactive_overlaps2 =
+            notactive_stepped2.AT(num_substeps).expand(overlap_width); });
+      xtime(time_check, [&]() {
+          check_equal("notactive_overlaps",
+                      notactive_overlaps1, notactive_overlaps2); });
+      if (all(all(buffer_width == num_substeps * ghost_width))) {
+        xtime(time1a, [&]() { assert(notactive_overlaps1 == notactive1); });
+        xtime(time2a, [&]() { assert(notactive_overlaps2 == notactive2); });
+      }
+      
+      // All buffer zones
+      ibset1& allbuffers1 = level_boxes.buffers1;
+      ibset2& allbuffers2 = level_boxes.buffers2;
+      xtime(time1, [&]() {
+          allbuffers1 = allowned1 & notowned1.expand(buffer_width); });
+      xtime(time2, [&]() {
+          allbuffers2 = allowned2 & notowned2.expand(buffer_width); });
+      xtime(time_check, [&]() {
+          check_equal("allbuffers", allbuffers1, allbuffers2); });
+      
+      // All overlap zones
+      ibset1 alloverlaps1;
+      ibset2 alloverlaps2;
+      xtime(time1, [&]() {
+          alloverlaps1 = (allowned1 & notactive1) - allbuffers1; });
+      xtime(time2, [&]() {
+          alloverlaps2 = (allowned2 & notactive2) - allbuffers2; });
+      xtime(time_check, [&]() {
+          check_equal("alloverlaps", alloverlaps1, alloverlaps2); });
+      
+      // All active points
+      ibset1& allactive1 = level_boxes.active1;
+      ibset2& allactive2 = level_boxes.active2;
+      xtime(time1, [&]() { allactive1 = allowned1 - notactive1; });
+      xtime(time2, [&]() { allactive2 = allowned2 - notactive2; });
+      xtime(time_check, [&]() {
+          check_equal("allactive", allactive1, allactive2); });
+      
+      // All stepped buffer zones
+      vector<ibset1> allbuffers_stepped1(num_substeps);
+      vector<ibset2> allbuffers_stepped2(num_substeps);
+      ibset1 allbuffers_stepped_combined1;
+      ibset2 allbuffers_stepped_combined2;
+      for (int substep=0; substep<num_substeps; ++substep) {
+        xtime(time1, [&]() {
+            allbuffers_stepped1.AT(substep) =
+              allowned1 &
+              (notactive_stepped1.AT(substep+1) -
+               notactive_stepped1.AT(substep));
+          });
+        xtime(time2, [&]() {
+            allbuffers_stepped2.AT(substep) =
+              allowned2 &
+              (notactive_stepped2.AT(substep+1) -
+               notactive_stepped2.AT(substep));
+          });
+        xtime(time_check, [&]() {
+            check_equal("notactive_stepped[" + to_str(substep) + "]",
+                        notactive_stepped1.AT(substep),
+                        notactive_stepped2.AT(substep)); });
+        xtime(time1a, [&]() {
+            allbuffers_stepped_combined1 += allbuffers_stepped1.AT(substep); });
+        xtime(time2a, [&]() {
+            allbuffers_stepped_combined2 += allbuffers_stepped2.AT(substep); });
+        xtime(time_check, [&]() {
+            check_equal("allbuffers_stepped_combined",
+                        allbuffers_stepped_combined1,
+                        allbuffers_stepped_combined2); });
+      }
+      if (all(all(buffer_width == num_substeps * ghost_width))) {
+        xtime(time1a, [&]() {
+            assert(allbuffers_stepped_combined1 == allbuffers1); });
+        xtime(time2a, [&]() {
+            assert(allbuffers_stepped_combined2 == allbuffers2); });
+      }
+      
+      // Overlap zones and buffer zones must be in the active part of
+      // the domain
+      xtime(time1a, [&]() { assert(allactive1 <= domain_active); });
+      xtime(time2a, [&]() { assert(allactive2 <= domain_active); });
+      xtime(time1a, [&]() { assert(alloverlaps1 <= domain_active); });
+      xtime(time2a, [&]() { assert(alloverlaps2 <= domain_active); });
+      xtime(time1a, [&]() { assert(allbuffers1 <= domain_active); });
+      xtime(time2a, [&]() { assert(allbuffers2 <= domain_active); });
+      xtime(time1a, [&]() { assert((allactive1 & alloverlaps1).empty()); });
+      xtime(time2a, [&]() { assert((allactive2 & alloverlaps2).empty()); });
+      xtime(time1a, [&]() { assert((allactive1 & allbuffers1).empty()); });
+      xtime(time2a, [&]() { assert((allactive2 & allbuffers2).empty()); });
+      xtime(time1a, [&]() { assert((alloverlaps1 & allbuffers1).empty()); });
+      xtime(time2a, [&]() { assert((alloverlaps2 & allbuffers2).empty()); });
+      xtime(time1a, [&]() {
+          assert(allactive1 + alloverlaps1 + allbuffers1 == allowned1); });
+      xtime(time2a, [&]() {
+          assert(allactive2 + alloverlaps2 + allbuffers2 == allowned2); });
+      
+      for (int c=0; c<components; ++c) {
+        auto& box = full_boxes.AT(c);
+        
+        // Buffer zones:
+        xtime(time1, [&]() { box.buffers1 = box.owned & allbuffers1; });
+        xtime(time2, [&]() { box.buffers2 = box.owned & allbuffers2; });
+        xtime(time_check, [&]() {
+            check_equal("buffers", box.buffers1, box.buffers2); });
+        
+        // Overlap zones:
+        xtime(time1, [&]() { box.overlaps1 = box.owned & alloverlaps1; });
+        xtime(time2, [&]() { box.overlaps2 = box.owned & alloverlaps2; });
+        xtime(time_check, [&]() {
+            check_equal("overlaps", box.overlaps1, box.overlaps2); });
+        
+        // Active region:
+        xtime(time1, [&]() { box.active1 = box.owned & allactive1; });
+        xtime(time2, [&]() { box.active2 = box.owned & allactive2; });
+        xtime(time_check, [&]() {
+            check_equal("active", box.active1, box.active2); });
+        xtime(time1a, [&]() {
+            assert(box.active1 ==
+                   box.owned - (box.buffers1 + box.overlaps1)); });
+        xtime(time2a, [&]() {
+            assert(box.active2 ==
+                   box.owned - (box.buffers2 + box.overlaps2)); });
+      } // for c
+      
+#if 0
+      for (int lc = 0; lc < h.local_components(rl); ++ lc) {
+        int const c = h.get_component (rl, lc);
+        local_dboxes & local_box = local_level.AT(lc);
+        full_dboxes const& box = full_boxes.AT(c);
+        
+        local_box.buffers = box.buffers;
+        local_box.overlaps = box.overlaps;
+        
+        local_box.active = box.active;
+      } // for lc
+      
+      // The conjunction of all buffer zones must equal allbuffers
+      ibset const allbuffers1 (full_boxes, & full_dboxes::buffers);
+      ASSERT_rl (allbuffers1 == allbuffers,
+                 "Buffer zone consistency check");
+#endif
+      
+      cout << "  Buffer zones:\n"
+           << "    time1:        " << time1 << " s\n"
+           << "    time2:        " << time2 << " s\n"
+           << "    time1_assert: " << time1a << " s\n"
+           << "    time2_assert: " << time2a << " s\n"
+           << "    time_check:   " << time_check << " s\n";
+
+    }
+    
+    
+    
+    CCTK_VInfo(CCTK_THORNSTRING, "Done testing bboxset2<%d>.", D);
+  }
+  
+} // namespace
+
+
+
+extern "C"
+void TestBBoxSet2_test(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  
+  test<1>();
+  test<2>();
+  test<3>();
+}
diff --git a/CarpetExtra/TestBBoxSet2/test/bboxset2.par b/CarpetExtra/TestBBoxSet2/test/bboxset2.par
new file mode 100644
index 000000000..cb54ea8ed
--- /dev/null
+++ b/CarpetExtra/TestBBoxSet2/test/bboxset2.par
@@ -0,0 +1,6 @@
+ActiveThorns = "Carpet CarpetLib CoordBase IOUtil TestBBoxSet2"
+
+IOUtil::out_dir = $parfile
+
+TestBBoxSet2::verbose         = no
+TestBBoxSet2::components_goal = 20000