CarpetLib: Add new transport operator "accumulate"

Add a new transport operator "accumulate", which accumulates ghost
zones into grid functions. (This is in a sense the opposite of regular
synchronisation, which copies grid function values into ghost zones.)

8 files changed, 229 insertions, 27 deletions

diff --git a/Carpet/CarpetLib/src/accumulate_3d.cc b/Carpet/CarpetLib/src/accumulate_3d.cc
new file mode 100644
index 000000000..f57cc06c0
--- /dev/null
+++ b/Carpet/CarpetLib/src/accumulate_3d.cc
@@ -0,0 +1,146 @@
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+
+#include <cctk.h>
+#include <cctk_Parameters.h>
+
+#include "operator_prototypes_3d.hh"
+#include "typeprops.hh"
+
+using namespace std;
+
+
+
+namespace CarpetLib {
+
+
+  
+#define SRCIND3(i,j,k)                                  \
+  index3 (srcioff + (i), srcjoff + (j), srckoff + (k),  \
+          srciext, srcjext, srckext)
+#define DSTIND3(i,j,k)                                  \
+  index3 (dstioff + (i), dstjoff + (j), dstkoff + (k),  \
+          dstiext, dstjext, dstkext)
+  
+  
+  
+  template <typename T>
+  void
+  accumulate_3d (T const * restrict const src,
+                 ivect3 const & restrict srcext,
+                 T * restrict const dst,
+                 ivect3 const & restrict dstext,
+                 ibbox3 const & restrict srcbbox,
+                 ibbox3 const & restrict dstbbox,
+                 ibbox3 const & restrict regbbox)
+  {
+    if (any (srcbbox.stride() != regbbox.stride() or
+             dstbbox.stride() != regbbox.stride()))
+    {
+#pragma omp critical
+      {
+        cout << "accumulate_3d.cc:" << endl
+             << "srcbbox=" << srcbbox << endl
+             << "dstbbox=" << dstbbox << endl
+             << "regbbox=" << regbbox << endl;
+        CCTK_WARN (0, "Internal error: strides disagree");
+      }
+    }
+    
+    if (any (srcbbox.stride() != dstbbox.stride())) {
+      CCTK_WARN (0, "Internal error: strides disagree");
+    }
+    
+    // This could be handled, but is likely to point to an error
+    // elsewhere
+    if (regbbox.empty()) {
+#pragma omp critical
+      CCTK_WARN (0, "Internal error: region extent is empty");
+    }
+    
+    if (not regbbox.is_contained_in(srcbbox) or
+        not regbbox.is_contained_in(dstbbox))
+    {
+#pragma omp critical
+      {
+        cout << "accumulate_3d.cc:" << endl
+             << "srcbbox=" << srcbbox << endl
+             << "dstbbox=" << dstbbox << endl
+             << "regbbox=" << regbbox << endl;
+        CCTK_WARN (0, "Internal error: region extent is not contained in array extent");
+      }
+    }
+    
+    if (any (srcext != srcbbox.shape() / srcbbox.stride() or
+             dstext != dstbbox.shape() / dstbbox.stride()))
+    {
+#pragma omp critical
+      CCTK_WARN (0, "Internal error: array sizes don't agree with bounding boxes");
+    }
+    
+    
+    
+    ivect3 const regext = regbbox.shape() / regbbox.stride();
+    assert (all ((regbbox.lower() - srcbbox.lower()) % srcbbox.stride() == 0));
+    ivect3 const srcoff = (regbbox.lower() - srcbbox.lower()) / srcbbox.stride();
+    assert (all ((regbbox.lower() - dstbbox.lower()) % dstbbox.stride() == 0));
+    ivect3 const dstoff = (regbbox.lower() - dstbbox.lower()) / dstbbox.stride();
+    
+    
+    
+    ptrdiff_t const srciext = srcext[0];
+    ptrdiff_t const srcjext = srcext[1];
+    ptrdiff_t const srckext = srcext[2];
+    
+    ptrdiff_t const dstiext = dstext[0];
+    ptrdiff_t const dstjext = dstext[1];
+    ptrdiff_t const dstkext = dstext[2];
+    
+    ptrdiff_t const regiext = regext[0];
+    ptrdiff_t const regjext = regext[1];
+    ptrdiff_t const regkext = regext[2];
+    
+    ptrdiff_t const srcioff = srcoff[0];
+    ptrdiff_t const srcjoff = srcoff[1];
+    ptrdiff_t const srckoff = srcoff[2];
+    
+    ptrdiff_t const dstioff = dstoff[0];
+    ptrdiff_t const dstjoff = dstoff[1];
+    ptrdiff_t const dstkoff = dstoff[2];
+    
+    
+    
+    // Loop over region
+    for (int k=0; k<regkext; ++k) {
+      for (int j=0; j<regjext; ++j) {
+        for (int i=0; i<regiext; ++i) {
+          
+          dst [DSTIND3(i, j, k)] += src [SRCIND3(i, j, k)];
+          
+        }
+      }
+    }
+    
+  }
+  
+  
+  
+#define TYPECASE(N,T)                                   \
+  template                                              \
+  void                                                  \
+  accumulate_3d (T const * restrict const src,          \
+                 ivect3 const & restrict srcext,        \
+                 T * restrict const dst,                \
+                 ivect3 const & restrict dstext,        \
+                 ibbox3 const & restrict srcbbox,       \
+                 ibbox3 const & restrict dstbbox,       \
+                 ibbox3 const & restrict regbbox);
+#include "typecase.hh"
+#undef TYPECASE
+  
+  
+  
+} // namespace CarpetLib
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 49229008a..2c38fcd9a 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -429,27 +429,38 @@ copy_from_innerloop (gdata const * const gsrc,
     assert (0);
   }
   
+  if (transport_operator != op_accumulate) {
 #if CARPET_DIM == 3
-  call_operator<T> (& copy_3d,
-                    static_cast <T const *> (src->storage()),
-                    src->shape(),
-                    static_cast <T *> (this->storage()),
-                    this->shape(),
-                    srcbox,
-                    dstbox,
-                    box);
+    call_operator<T> (& copy_3d,
+                      static_cast <T const *> (src->storage()),
+                      src->shape(),
+                      static_cast <T *> (this->storage()),
+                      this->shape(),
+                      srcbox,
+                      dstbox,
+                      box);
 #elif CARPET_DIM == 4
-  call_operator<T> (& copy_4d,
-                    static_cast <T const *> (src->storage()),
-                    src->shape(),
-                    static_cast <T *> (this->storage()),
-                    this->shape(),
-                    srcbox,
-                    dstbox,
-                    box);
+    call_operator<T> (& copy_4d,
+                      static_cast <T const *> (src->storage()),
+                      src->shape(),
+                      static_cast <T *> (this->storage()),
+                      this->shape(),
+                      srcbox,
+                      dstbox,
+                      box);
 #else
 #  error "Value for CARPET_DIM not supported"
 #endif
+  } else {
+    call_operator<T> (& accumulate_3d,
+                      static_cast <T const *> (src->storage()),
+                      src->shape(),
+                      static_cast <T *> (this->storage()),
+                      this->shape(),
+                      srcbox,
+                      dstbox,
+                      box);
+  }
 }
 
 
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index c9f9f8ad2..3d3c06524 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -305,7 +305,7 @@ find_source_timelevel (vector <CCTK_REAL> const & times,
   CCTK_REAL const max_time = * max_element (times.begin(), times.end());
   // TODO: Use a real delta-time from somewhere instead of 1.0
   CCTK_REAL const some_time = abs (min_time) + abs (max_time) + 1.0;
-  if (op != op_copy) {
+  if (op != op_copy and op != op_accumulate) {
     if (time < min_time - eps * some_time or
         time > max_time + eps * some_time)
     {
@@ -333,7 +333,7 @@ find_source_timelevel (vector <CCTK_REAL> const & times,
     }
   }
   if (timelevel == -1) {
-    if (op == op_copy) {
+    if (op == op_copy or op == op_accumulate) {
       timelevel = 0;
     }
   }
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index ce6fe3ee3..314a04843 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -346,6 +346,26 @@ sync_all (comm_state & state,
 
 
 
+// Accumulate the boundaries of all components
+void
+ggf::
+accumulate_all (comm_state & state,
+                int const tl, int const rl, int const ml)
+{
+  assert (transport_operator == op_accumulate);
+  // Accumulate
+  static Timer timer ("accumulate_all");
+  timer.start ();
+  transfer_from_all (state,
+                     tl,rl,ml,
+                     & dh::fast_dboxes::fast_sync_sendrecv,
+                     tl,rl,ml,
+                     NULL, true);
+  timer.stop (0);
+}
+
+
+
 // Prolongate the boundaries of all components
 void
 ggf::
@@ -361,7 +381,7 @@ ref_bnd_prolongate_all (comm_state & state,
   vector<int> tl2s;
   static Timer timer ("ref_bnd_prolongate_all");
   timer.start ();
-  if (transport_operator != op_copy) {
+  if (transport_operator != op_copy and transport_operator != op_accumulate) {
     // Interpolation in time
     if (not (timelevels(ml,rl) >= prolongation_order_time+1)) {
       char * const fullname = CCTK_FullName (varindex);
@@ -518,7 +538,8 @@ transfer_from_all (comm_state & state,
                    srpvect const dh::fast_dboxes::* sendrecvs,
                    vector<int> const & tl2s, int const rl2, int const ml2,
                    CCTK_REAL const & time,
-                   mdata * const srcstorage_)
+                   mdata * const srcstorage_,
+                   bool const flip_send_recv)
 {
   assert (rl1>=0 and rl1<h.reflevels());
   assert (ml1>=0 and ml1<h.mglevels());
@@ -558,9 +579,11 @@ transfer_from_all (comm_state & state,
   // Interpolation orders
   assert (transport_operator != op_none);
   int const pos =
-    transport_operator == op_copy ? 0 : d.prolongation_orders_space.AT(rl2);
+    transport_operator == op_copy or transport_operator == op_accumulate
+    ? 0 : d.prolongation_orders_space.AT(rl2);
   int const pot =
-    transport_operator == op_copy ? 0 : prolongation_order_time;
+    transport_operator == op_copy or transport_operator == op_accumulate
+    ? 0 : prolongation_order_time;
   
   vector<const gdata*> gsrcs(tl2s.size());
   
@@ -568,8 +591,12 @@ transfer_from_all (comm_state & state,
   for (srpvect::const_iterator ipsendrecv = psendrecvs.begin();
        ipsendrecv!=psendrecvs.end(); ++ ipsendrecv)
   {
-    pseudoregion_t const & psend = (* ipsendrecv).send;
-    pseudoregion_t const & precv = (* ipsendrecv).recv;
+    typedef sendrecv_pseudoregion_t srp;
+    typedef pseudoregion_t sendrecv_pseudoregion_t::* srp_pr;
+    srp_pr const send_field = flip_send_recv ? &srp::recv : &srp::send;
+    srp_pr const recv_field = flip_send_recv ? &srp::send : &srp::recv;
+    pseudoregion_t const & psend = (* ipsendrecv).*send_field;
+    pseudoregion_t const & precv = (* ipsendrecv).*recv_field;
     ibbox const & send = psend.extent;
     ibbox const & recv = precv.extent;
     assert (all (send.stride() == h.baseextent(ml2,rl2).stride()));
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index d9412f822..56263fa80 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -124,6 +124,9 @@ public:
   // Synchronise the boundaries of a component
   void sync_all (comm_state& state, int tl, int rl, int ml);
 
+  // Accumulate from boundaries of a component
+  void accumulate_all (comm_state& state, int tl, int rl, int ml);
+
   // Prolongate the boundaries of a component
   void ref_bnd_prolongate_all (comm_state& state,
                                int tl, int rl, int ml, CCTK_REAL time);
@@ -168,14 +171,16 @@ protected:
                      srpvect const dh::fast_dboxes::* sendrecvs,
                      vector<int> const & tl2s, int rl2, int ml2,
                      CCTK_REAL const & time,
-                     mdata * srcstorage = 0);
+                     mdata * srcstorage = 0,
+                     bool flip_send_recv = false);
   
   void
   transfer_from_all (comm_state & state,
                      int tl1, int rl1, int ml1,
                      srpvect const dh::fast_dboxes::* sendrecvs,
                      int tl2, int rl2, int ml2,
-                     mdata * srcstorage = 0)
+                     mdata * srcstorage = 0,
+                     bool flip_send_recv = false)
   {
     vector <int> tl2s(1);
     tl2s.AT(0) = tl2;
@@ -185,7 +190,8 @@ protected:
                        sendrecvs,
                        tl2s, rl2, ml2,
                        time,
-                       srcstorage);
+                       srcstorage,
+                       flip_send_recv);
   }
 
 
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index 96d621172..f63ce2b62 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -25,6 +25,7 @@ SRCS =	balance.cc				\
 	th.cc					\
 	timestat.cc				\
 	vect.cc					\
+	accumulate_3d.cc			\
 	copy_3d.cc				\
 	copy_4d.cc				\
 	interpolate_3d_2tl.cc			\
diff --git a/Carpet/CarpetLib/src/operator_prototypes_3d.hh b/Carpet/CarpetLib/src/operator_prototypes_3d.hh
index 81a384d76..c4793f606 100644
--- a/Carpet/CarpetLib/src/operator_prototypes_3d.hh
+++ b/Carpet/CarpetLib/src/operator_prototypes_3d.hh
@@ -58,6 +58,16 @@ namespace CarpetLib {
            ibbox3 const & restrict dstbbox,
            ibbox3 const & restrict regbbox);
   
+  template <typename T>
+  void
+  accumulate_3d (T const * restrict const src,
+                 ivect3 const & restrict srcext,
+                 T * restrict const dst,
+                 ivect3 const & restrict dstext,
+                 ibbox3 const & restrict srcbbox,
+                 ibbox3 const & restrict dstbbox,
+                 ibbox3 const & restrict regbbox);
+  
   
   
   template <typename T>
diff --git a/Carpet/CarpetLib/src/operators.hh b/Carpet/CarpetLib/src/operators.hh
index 07e4e5616..df3828d44 100644
--- a/Carpet/CarpetLib/src/operators.hh
+++ b/Carpet/CarpetLib/src/operators.hh
@@ -12,6 +12,7 @@ enum operator_type
   op_restrict,                  // restrict only, do not prolongate
   op_copy,                      // use simple copying for prolongation
                                 // (needs only one time level)
+  op_accumulate,                // accumulate (sum) into the destination
   op_Lagrange,                  // Lagrange interpolation (standard)
   op_ENO,                       // use ENO stencils (for hydro)
   op_WENO,                      // use WENO stencils (for hydro)