CarpetLib: Add ENO timer interpolator

darcs-hash:20070219202316-dae7b-441bdad39e75963a86bedbcaea54bbf7ced0f8a7.gz

4 files changed, 383 insertions, 34 deletions

diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 20768866d..64276b2a3 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -860,47 +860,107 @@ void data <T>
                     CCTK_REAL const time,
                     int const order_time)
 {
-  switch (order_time) {
-    
-  case 0:
-    // We could handle this, but this points to an inefficiency
-    assert (0);
+  static Timer total ("time_interpolate");
+  total.start ();
+
+  switch (transport_operator) {
     
-  case 1:
-    assert (times.size() >= 2);
-    interpolate_3d_2tl (static_cast <T const *> (srcs.AT(0)->storage()),
-                        times.AT(0),
-                        static_cast <T const *> (srcs.AT(1)->storage()),
-                        times.AT(1),
-                        srcs.AT(0)->shape(),
-                        static_cast <T *> (this->storage()),
-                        time,
-                        this->shape(),
-                        srcs.AT(0)->extent(),
-                        this->extent(),
-                        box);
+  case op_copy:
+  case op_Lagrange: {
+    static Timer timer ("time_interpolate_Lagrange");
+    timer.start ();
+    switch (order_time) {
+      
+    case 1:
+      assert (times.size() >= 2);
+      interpolate_3d_2tl (static_cast <T const *> (srcs.AT(0)->storage()),
+                          times.AT(0),
+                          static_cast <T const *> (srcs.AT(1)->storage()),
+                          times.AT(1),
+                          srcs.AT(0)->shape(),
+                          static_cast <T *> (this->storage()),
+                          time,
+                          this->shape(),
+                          srcs.AT(0)->extent(),
+                          this->extent(),
+                          box);
+      break;
+      
+    case 2:
+      assert (times.size() >= 3);
+      interpolate_3d_3tl (static_cast <T const *> (srcs.AT(0)->storage()),
+                          times.AT(0),
+                          static_cast <T const *> (srcs.AT(1)->storage()),
+                          times.AT(1),
+                          static_cast <T const *> (srcs.AT(2)->storage()),
+                          times.AT(2),
+                          srcs.AT(0)->shape(),
+                          static_cast <T *> (this->storage()),
+                          time,
+                          this->shape(),
+                          srcs.AT(0)->extent(),
+                          this->extent(),
+                          box);
+      break;
+      
+    default:
+      assert (0);
+    }
+    timer.stop (0);
     break;
+  }
     
-  case 2:
-    assert (times.size() >= 3);
-    interpolate_3d_3tl (static_cast <T const *> (srcs.AT(0)->storage()),
-                        times.AT(0),
-                        static_cast <T const *> (srcs.AT(1)->storage()),
-                        times.AT(1),
-                        static_cast <T const *> (srcs.AT(2)->storage()),
-                        times.AT(2),
-                        srcs.AT(0)->shape(),
-                        static_cast <T *> (this->storage()),
-                        time,
-                        this->shape(),
-                        srcs.AT(0)->extent(),
-                        this->extent(),
-                        box);
+  case op_ENO:
+  case op_WENO: {
+    // ENO and WENO timer interpolation is the same for order_time <= 2
+    static Timer timer ("time_interpolate_ENO");
+    timer.start ();
+    switch (order_time) {
+      
+    case 1:
+      assert (times.size() >= 2);
+      interpolate_3d_2tl (static_cast <T const *> (srcs.AT(0)->storage()),
+                          times.AT(0),
+                          static_cast <T const *> (srcs.AT(1)->storage()),
+                          times.AT(1),
+                          srcs.AT(0)->shape(),
+                          static_cast <T *> (this->storage()),
+                          time,
+                          this->shape(),
+                          srcs.AT(0)->extent(),
+                          this->extent(),
+                          box);
+      break;
+      
+    case 2:
+      assert (times.size() >= 3);
+      interpolate_eno_3d_3tl (static_cast <T const *> (srcs.AT(0)->storage()),
+                              times.AT(0),
+                              static_cast <T const *> (srcs.AT(1)->storage()),
+                              times.AT(1),
+                              static_cast <T const *> (srcs.AT(2)->storage()),
+                              times.AT(2),
+                              srcs.AT(0)->shape(),
+                              static_cast <T *> (this->storage()),
+                              time,
+                              this->shape(),
+                              srcs.AT(0)->extent(),
+                              this->extent(),
+                              box);
+      break;
+      
+    default:
+      assert (0);
+    }
+    timer.stop (0);
     break;
+  }
     
   default:
     assert (0);
-  }
+  } // switch (transport_operator)
+  
+  total.stop (0);
 }
 
 template <>
diff --git a/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc b/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
new file mode 100644
index 000000000..f0b19ed35
--- /dev/null
+++ b/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
@@ -0,0 +1,272 @@
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+
+#include <cctk.h>
+#include <cctk_Parameters.h>
+
+#include "operator_prototypes.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>
+  static inline
+  T
+  min3 (T const & x, T const & y, T const & z)
+  {
+    return min (x, min (y, z));
+  }
+  
+  template <typename T>
+  static inline
+  T
+  max3 (T const & x, T const & y, T const & z)
+  {
+    return max (x, max (y, z));
+  }
+  
+  
+  
+  template <typename T>
+  void
+  interpolate_eno_3d_3tl (T const * restrict const src1,
+                          CCTK_REAL const t1,
+                          T const * restrict const src2,
+                          CCTK_REAL const t2,
+                          T const * restrict const src3,
+                          CCTK_REAL const t3,
+                          ivect3 const & restrict srcext,
+                          T * restrict const dst,
+                          CCTK_REAL const t,
+                          ivect3 const & restrict dstext,
+                          ibbox3 const & restrict srcbbox,
+                          ibbox3 const & restrict dstbbox,
+                          ibbox3 const & restrict regbbox)
+  {
+    typedef typename typeprops<T>::real RT;
+    
+    
+    
+    if (any (srcbbox.stride() != regbbox.stride() or
+             dstbbox.stride() != regbbox.stride()))
+    {
+      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()) {
+      CCTK_WARN (0, "Internal error: region extent is empty");
+    }
+    
+    if (not regbbox.is_contained_in(srcbbox) or
+        not regbbox.is_contained_in(dstbbox))
+    {
+      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()))
+    {
+      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();
+    
+    
+    
+    size_t const srciext = srcext[0];
+    size_t const srcjext = srcext[1];
+    size_t const srckext = srcext[2];
+    
+    size_t const dstiext = dstext[0];
+    size_t const dstjext = dstext[1];
+    size_t const dstkext = dstext[2];
+    
+    size_t const regiext = regext[0];
+    size_t const regjext = regext[1];
+    size_t const regkext = regext[2];
+    
+    size_t const srcioff = srcoff[0];
+    size_t const srcjoff = srcoff[1];
+    size_t const srckoff = srcoff[2];
+    
+    size_t const dstioff = dstoff[0];
+    size_t const dstjoff = dstoff[1];
+    size_t const dstkoff = dstoff[2];
+    
+    
+    
+    // Quadratic (second order) interpolation
+    
+    RT const eps = 1.0e-10;
+    
+    if (abs (t1 - t2) < eps or abs (t1 - t3) < eps or abs (t2 - t3) < eps) {
+      CCTK_WARN (0, "Internal error: arrays have same time");
+    }
+    if (t < min3 (t1, t2, t3) - eps or t > max3 (t1, t2, t3) + eps) {
+      CCTK_WARN (0, "Internal error: extrapolation in time");
+    }
+    
+    RT const s1fac3 = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3));
+    RT const s2fac3 = (t - t1) * (t - t3) / ((t2 - t1) * (t2 - t3));
+    RT const s3fac3 = (t - t1) * (t - t2) / ((t3 - t1) * (t3 - t2));
+    
+    RT const s1fac2_12 = (t - t2) / (t1 - t2);
+    RT const s2fac2_12 = (t - t1) / (t2 - t1);
+    
+    RT const s2fac2_23 = (t - t3) / (t2 - t3);
+    RT const s3fac2_23 = (t - t2) / (t3 - t2);
+    
+    bool const use_12 =
+      t >= min (t1, t2) - eps and t <= max (t1, t2) + eps;
+    bool const use_23 =
+      t >= min (t2, t3) - eps and t <= max (t2, t3) + eps;
+    assert (use_12 or use_23);
+    
+    
+    
+    // Loop over region
+    for (size_t k=0; k<regkext; ++k) {
+      for (size_t j=0; j<regjext; ++j) {
+        for (size_t i=0; i<regiext; ++i) {
+          
+          T const s1 = src1 [SRCIND3(i, j, k)];
+          T const s2 = src2 [SRCIND3(i, j, k)];
+          T const s3 = src3 [SRCIND3(i, j, k)];
+          
+          T d = s1fac3 * s1 + s2fac3 * s2 + s3fac3 * s3;
+          
+          // if ((d - max3 (s1, s2, s3)) * (d - min3 (s1, s2, s3)) < 0)
+          if (d > max3 (s1, s2, s3) or d < min3 (s1, s2, s3)) {
+            if (use_12) {
+              d = s1fac2_12 * s1 + s2fac2_12 * s2;
+            } else {
+              d = s2fac2_23 * s2 + s3fac2_23 * s3;
+            }
+          }
+          
+          dst [DSTIND3(i, j, k)] = d;
+          
+        }
+      }
+    }
+    
+  }
+  
+  
+  
+#ifdef HAVE_CCTK_COMPLEX8
+  template <>
+  void
+  interpolate_eno_3d_3tl (CCTK_COMPLEX8 const * restrict const src1,
+                          CCTK_REAL const t1,
+                          CCTK_COMPLEX8 const * restrict const src2,
+                          CCTK_REAL const t2,
+                          CCTK_COMPLEX8 const * restrict const src3,
+                          CCTK_REAL const t3,
+                          ivect3 const & restrict srcext,
+                          CCTK_COMPLEX8 * restrict const dst,
+                          CCTK_REAL const t,
+                          ivect3 const & restrict dstext,
+                          ibbox3 const & restrict srcbbox,
+                          ibbox3 const & restrict dstbbox,
+                          ibbox3 const & restrict regbbox)
+  {
+    CCTK_WARN (CCTK_WARN_ABORT, "ENO for complex numbers is not supported");
+  }
+#endif
+  
+#ifdef HAVE_CCTK_COMPLEX16
+  template <>
+  void
+  interpolate_eno_3d_3tl (CCTK_COMPLEX16 const * restrict const src1,
+                          CCTK_REAL const t1,
+                          CCTK_COMPLEX16 const * restrict const src2,
+                          CCTK_REAL const t2,
+                          CCTK_COMPLEX16 const * restrict const src3,
+                          CCTK_REAL const t3,
+                          ivect3 const & restrict srcext,
+                          CCTK_COMPLEX16 * restrict const dst,
+                          CCTK_REAL const t,
+                          ivect3 const & restrict dstext,
+                          ibbox3 const & restrict srcbbox,
+                          ibbox3 const & restrict dstbbox,
+                          ibbox3 const & restrict regbbox)
+  {
+    CCTK_WARN (CCTK_WARN_ABORT, "ENO for complex numbers is not supported");
+  }
+#endif
+  
+#ifdef HAVE_CCTK_COMPLEX32
+  template <>
+  void
+  interpolate_eno_3d_3tl (CCTK_COMPLEX32 const * restrict const src1,
+                          CCTK_REAL const t1,
+                          CCTK_COMPLEX32 const * restrict const src2,
+                          CCTK_REAL const t2,
+                          CCTK_COMPLEX32 const * restrict const src3,
+                          CCTK_REAL const t3,
+                          ivect3 const & restrict srcext,
+                          CCTK_COMPLEX32 * restrict const dst,
+                          CCTK_REAL const t,
+                          ivect3 const & restrict dstext,
+                          ibbox3 const & restrict srcbbox,
+                          ibbox3 const & restrict dstbbox,
+                          ibbox3 const & restrict regbbox)
+  {
+    CCTK_WARN (CCTK_WARN_ABORT, "ENO for complex numbers is not supported");
+  }
+#endif
+  
+  
+  
+#define INSTANTIATE(T)                                          \
+  template                                                      \
+  void                                                          \
+  interpolate_eno_3d_3tl (T const * restrict const src1,        \
+                          CCTK_REAL const t1,                   \
+                          T const * restrict const src2,        \
+                          CCTK_REAL const t2,                   \
+                          T const * restrict const src3,        \
+                          CCTK_REAL const t3,                   \
+                          ivect3 const & restrict srcext,       \
+                          T * restrict const dst,               \
+                          CCTK_REAL const t,                    \
+                          ivect3 const & restrict dstext,       \
+                          ibbox3 const & restrict srcbbox,      \
+                          ibbox3 const & restrict dstbbox,      \
+                          ibbox3 const & restrict regbbox);
+#include "instantiate"
+#undef INSTANTIATE
+  
+  
+  
+} // namespace CarpetLib
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index 3b47c6016..8c5c4e4a9 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -20,6 +20,7 @@ SRCS =	bbox.cc					\
 	copy_3d.cc				\
 	interpolate_3d_2tl.cc			\
 	interpolate_3d_3tl.cc			\
+	interpolate_eno_3d_3tl.cc		\
 	restrict_3d_cc_rf2.cc			\
 	restrict_3d_rf2.cc			\
 	prolongate_3d_cc_rf2.cc			\
diff --git a/Carpet/CarpetLib/src/operator_prototypes.hh b/Carpet/CarpetLib/src/operator_prototypes.hh
index c2a6e8fbc..40749af1e 100644
--- a/Carpet/CarpetLib/src/operator_prototypes.hh
+++ b/Carpet/CarpetLib/src/operator_prototypes.hh
@@ -128,6 +128,22 @@ namespace CarpetLib {
                       ibbox3 const & restrict dstbbox,
                       ibbox3 const & restrict regbbox);
   
+  template <typename T>
+  void
+  interpolate_eno_3d_3tl (T const * restrict const src1,
+                          CCTK_REAL const t1,
+                          T const * restrict const src2,
+                          CCTK_REAL const t2,
+                          T const * restrict const src3,
+                          CCTK_REAL const t3,
+                          ivect3 const & restrict srcext,
+                          T * restrict const dst,
+                          CCTK_REAL const t,
+                          ivect3 const & restrict dstext,
+                          ibbox3 const & restrict srcbbox,
+                          ibbox3 const & restrict dstbbox,
+                          ibbox3 const & restrict regbbox);
+  
   
   
   template <typename T>