From 0b07e99df3d194924a2e5904f5e13bd0c1dac0b8 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Mon, 19 Feb 2007 20:23:00 +0000 Subject: CarpetLib: Add ENO timer interpolator darcs-hash:20070219202316-dae7b-441bdad39e75963a86bedbcaea54bbf7ced0f8a7.gz --- Carpet/CarpetLib/src/data.cc | 128 ++++++++---- Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc | 272 +++++++++++++++++++++++++ Carpet/CarpetLib/src/make.code.defn | 1 + Carpet/CarpetLib/src/operator_prototypes.hh | 16 ++ 4 files changed, 383 insertions(+), 34 deletions(-) create mode 100644 Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc (limited to 'Carpet') 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 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 (srcs.AT(0)->storage()), - times.AT(0), - static_cast (srcs.AT(1)->storage()), - times.AT(1), - srcs.AT(0)->shape(), - static_cast (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 (srcs.AT(0)->storage()), + times.AT(0), + static_cast (srcs.AT(1)->storage()), + times.AT(1), + srcs.AT(0)->shape(), + static_cast (this->storage()), + time, + this->shape(), + srcs.AT(0)->extent(), + this->extent(), + box); + break; + + case 2: + assert (times.size() >= 3); + interpolate_3d_3tl (static_cast (srcs.AT(0)->storage()), + times.AT(0), + static_cast (srcs.AT(1)->storage()), + times.AT(1), + static_cast (srcs.AT(2)->storage()), + times.AT(2), + srcs.AT(0)->shape(), + static_cast (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 (srcs.AT(0)->storage()), - times.AT(0), - static_cast (srcs.AT(1)->storage()), - times.AT(1), - static_cast (srcs.AT(2)->storage()), - times.AT(2), - srcs.AT(0)->shape(), - static_cast (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 (srcs.AT(0)->storage()), + times.AT(0), + static_cast (srcs.AT(1)->storage()), + times.AT(1), + srcs.AT(0)->shape(), + static_cast (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 (srcs.AT(0)->storage()), + times.AT(0), + static_cast (srcs.AT(1)->storage()), + times.AT(1), + static_cast (srcs.AT(2)->storage()), + times.AT(2), + srcs.AT(0)->shape(), + static_cast (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 +#include +#include +#include + +#include +#include + +#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 + static inline + T + min3 (T const & x, T const & y, T const & z) + { + return min (x, min (y, z)); + } + + template + static inline + T + max3 (T const & x, T const & y, T const & z) + { + return max (x, max (y, z)); + } + + + + 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) + { + typedef typename typeprops::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 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 + 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 -- cgit v1.2.3