CarpetLib: Implement padding for grid variables

Ignore-this: 1a389f0dd3f40a0c0edb3fdabd6e7d40

Padding grid variables means that e.g. a component of size 32x32x32 is
allocated as 33x33x33 instead, but only 32x32x32 of this storage is
used.  This can improve cache performance considerably.  This requires
corresponding changes to the cGH entries.

5 files changed, 44 insertions, 2 deletions

diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl
index 31a43a135..0dd901181 100644
--- a/Carpet/CarpetLib/param.ccl
+++ b/Carpet/CarpetLib/param.ccl
@@ -120,6 +120,16 @@ BOOLEAN combine_recompose "Recompose all grid functions of one refinement levels
 
 
 
+# Memory allocation parameters
+
+INT avoid_arraysize_bytes "Avoid array sizes that are multiples of this" STEERABLE=recover
+{
+  0   :: "don't avoid anything"
+  2:* :: ""
+} 0
+
+
+
 # Communication experiment parameters
 
 INT message_size_multiplier "Enlarge size of transmitted messages by this factor" STEERABLE=always
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index da32ee4ea..9482ad518 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -336,15 +336,20 @@ void data<T>::allocate (const ibbox& extent_,
   assert (all(extent_.lower() > numeric_limits<int>::min() / 2));
   assert (all(extent_.upper() < numeric_limits<int>::max() / 2));
   assert (all(extent_.upper() > numeric_limits<int>::min() / 2));
+  
   // data
   _extent = extent_;
-  _shape = max(ivect(0), _extent.shape() / _extent.stride());
+  //_shape = max (ivect(0), _extent.shape() / _extent.stride());
+  _shape = allocated_memory_shape (_extent);
+  assert (all (_shape >= max (ivect(0), _extent.shape() / _extent.stride())));
+  
   _size = 1;
   for (int d=0; d<dim; ++d) {
     _stride[d] = _size;
     assert (_shape[d]==0 or _size <= numeric_limits<int>::max() / _shape[d]);
     _size *= _shape[d];
   }
+  
   _proc = proc_;
   if (dist::rank() == _proc) {
     if (vectorindex == 0) {
@@ -380,7 +385,9 @@ size_t data<T>::allocsize (const ibbox & extent_, const int proc_) const
   if (dist::rank() != proc_) return 0;
   if (vectorindex != 0) return 0;
   assert (not vectorleader);
-  return vectorlength * extent_.size() * sizeof (T);
+  ivect const shape_ = allocated_memory_shape (extent_);
+  //return vectorlength * extent_.size() * sizeof (T);
+  return vectorlength * prod(shape_) * sizeof (T);
 }
 
 
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 36b3c310f..e107ab306 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,4 +1,5 @@
 #include <cassert>
+#include <cstddef>
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index e1888835c..4fa5e889f 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -64,6 +64,29 @@ gdata::~gdata ()
 
 
 
+// Storage management
+
+ivect
+gdata::
+allocated_memory_shape (ibbox const& extent)
+{
+  DECLARE_CCTK_PARAMETERS;
+  ivect shape = max (ivect(0), extent.shape() / extent.stride());
+  // Enlarge shape to avoid multiples of cache line colours
+  if (avoid_arraysize_bytes > 0) {
+    for (int d=0; d<dim; ++d) {
+      if (shape[d] > 0 and
+          shape[d] * sizeof(CCTK_REAL) % avoid_arraysize_bytes == 0)
+      {
+        ++shape[d];
+      }
+    }
+  }
+  return shape;
+}
+
+
+
 // Data manipulators
 
 void
diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh
index 4b62cf564..5836f8f4e 100644
--- a/Carpet/CarpetLib/src/gdata.hh
+++ b/Carpet/CarpetLib/src/gdata.hh
@@ -78,6 +78,7 @@ public:
 			 void* const memptr = NULL, size_t const memsize = 0) = 0;
   virtual void free () = 0;
   virtual size_t allocsize (const ibbox& extent, const int proc) const = 0;
+  static ivect allocated_memory_shape (ibbox const& extent);
   
   // Accessors
   bool has_storage () const {