From 2de0ecd32acad32a49cdd781330f6b2a49506f5d Mon Sep 17 00:00:00 2001 From: Roland Haas Date: Mon, 4 Jun 2012 17:50:53 -0700 Subject: CarpetLib: add fifth order cc restriction operator --- Carpet/CarpetLib/param.ccl | 1 + Carpet/CarpetLib/src/data.cc | 10 + Carpet/CarpetLib/src/make.code.defn | 1 + Carpet/CarpetLib/src/operator_prototypes_3d.hh | 14 ++ Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc | 252 +++++++++++++++++++++++++ 5 files changed, 278 insertions(+) create mode 100644 Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl index a259fff2d..8b3adcdd2 100644 --- a/Carpet/CarpetLib/param.ccl +++ b/Carpet/CarpetLib/param.ccl @@ -41,6 +41,7 @@ INT restriction_order_space "Order of restriction operator to use with use_highe { 1 :: "linear interpolation, this is Carpet's original implementation" 3 :: "third order accurate restriction for grid functions where prolongation is not (W)ENO" + 5 :: "fifth order accurate restriction for grid functions where prolongation is not (W)ENO" } 3 diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index ba1a438ce..6268a6e64 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -1165,6 +1165,16 @@ transfer_restrict (data const * const src, dstbox, srcregbox, dstregbox, NULL); break; + case 5: + // Don't use call_operator, because we parallelise ourselves + restrict_3d_cc_o5_rf2(static_cast (src->storage()), + src->shape(), + static_cast (this->storage()), + this->shape(), + srcbox, + dstbox, + srcregbox, dstregbox, NULL); + break; default: CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, "There is no restriction stencil with restriction_order_space==%d", restriction_order_space); diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index 89bdca337..d1d623685 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -35,6 +35,7 @@ SRCS = backtrace.cc \ interpolate_eno_3d_3tl.cc \ restrict_3d_cc_rf2.cc \ restrict_3d_cc_o3_rf2.cc \ + restrict_3d_cc_o5_rf2.cc \ restrict_3d_rf2.cc \ restrict_3d_vc_rf2.cc \ restrict_3d_dgfe_rf2.cc \ diff --git a/Carpet/CarpetLib/src/operator_prototypes_3d.hh b/Carpet/CarpetLib/src/operator_prototypes_3d.hh index dd5eb14d9..0d03f8cb9 100644 --- a/Carpet/CarpetLib/src/operator_prototypes_3d.hh +++ b/Carpet/CarpetLib/src/operator_prototypes_3d.hh @@ -396,6 +396,20 @@ namespace CarpetLib { + template + void + restrict_3d_cc_o5_rf2 (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 srcregbbox, + ibbox3 const & restrict dstregbbox, + void * extraargs); + + + template void restrict_3d_vc_rf2 (T const * restrict const src, diff --git a/Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc new file mode 100644 index 000000000..c9067c9d9 --- /dev/null +++ b/Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc @@ -0,0 +1,252 @@ +#include +#include + +#include +#include +#include + +#include "gdata.hh" +#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) +#define SRCOFF3(i,j,k) \ + offset3 (srcioff + (i), srcjoff + (j), srckoff + (k), \ + srciext, srcjext, srckext) +#define DSTOFF3(i,j,k) \ + offset3 (dstioff + (i), dstjoff + (j), dstkoff + (k), \ + dstiext, dstjext, dstkext) + + +// This operator offers fifth-order accurate restriction operators for cell +// centered grids when use_higher_order_restriction is set. This interpolation +// is done for samples at the cell centres (so this is not a ppm scheme or +// anything like that). It really only fits a polynomial of the form +// f(x,y,z) = \sum_{i,j,k=0}^5 a_{i,j,k} x^i y^j z^k +// to the fine cells and evaluates at x=y=z=0. So it is good for the metric, +// but bad for matter (since it will destroy the conservation). Because of +// this it should not be used for grid functions whose transport operator is +// not WENO or ENO which hopefully excludes all matter variables. +// TODO: use prolongate_3d_rf2 instead of writing new operators, its the same +// interpolation. + + // 1D restriction + template + static inline T restrict1 (T const * restrict const p, + size_t const d1) + { + typedef typename typeprops::real RT; + RT const den = 256; + RT const f2 = 3/den, f1 = 25/den, f0 = 150/den; + T const res = + + f2 * p[- 2] + - f1 * p[- 1] + + f0 * p[- 0] + + f0 * p[+ 1] + - f1 * p[+ 2] + + f2 * p[+ 3]; + return res; + } + + // 2D restriction + template + static inline T restrict2 (T const * restrict const p, + size_t const d1, size_t const d2) + { + typedef typename typeprops::real RT; + RT const den = 256; + RT const f2 = 3/den, f1 = 25/den, f0 = 150/den; + T const res = + + f2 * restrict1 (p - 2*d2, d1) + - f1 * restrict1 (p - 1*d2, d1) + + f0 * restrict1 (p - 0*d2, d1) + + f0 * restrict1 (p + 1*d2, d1) + - f1 * restrict1 (p + 2*d2, d1) + + f2 * restrict1 (p + 3*d2, d1); + return res; + } + + // 3D restriction + template + static inline T restrict3 (T const * restrict const p, + size_t const d1, size_t const d2, size_t const d3) + { + typedef typename typeprops::real RT; + RT const den = 256; + RT const f2 = 3/den, f1 = 25/den, f0 = 150/den; + T const res = + + f2 * restrict2 (p - 2*d3, d1, d2) + - f1 * restrict2 (p - 1*d3, d1, d2) + + f0 * restrict2 (p - 0*d3, d1, d2) + + f0 * restrict2 (p + 1*d3, d1, d2) + - f1 * restrict2 (p + 2*d3, d1, d2) + + f2 * restrict2 (p + 3*d3, d1, d2); + return res; + } + + template + void + restrict_3d_cc_o5_rf2 (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, + ibbox3 const & restrict regbbox, + void * extraargs) + { + assert (not extraargs); + + DECLARE_CCTK_PARAMETERS; + + 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 (reffact2 * srcbbox.stride() != dstbbox.stride())) { + CCTK_WARN (0, "Internal error: destination strides are not twice the source strides"); + } + + // 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.expanded_for(srcbbox).expand(1,1).is_contained_in(srcbbox) or + not regbbox.is_contained_in(dstbbox)) + { + cerr << "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 != gdata::allocated_memory_shape(srcbbox.shape() / srcbbox.stride()) or + dstext != gdata::allocated_memory_shape(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 (srcbbox.stride() % 2 == 0)); + assert (all ((regbbox.lower() - srcbbox.lower() - srcbbox.stride() / 2) % + srcbbox.stride() == 0)); + ivect3 const srcoff = + (regbbox.lower() - srcbbox.lower() - srcbbox.stride() / 2) / + srcbbox.stride(); + assert (all ((regbbox.lower() - dstbbox.lower()) % dstbbox.stride() == 0)); + ivect3 const dstoff = + (regbbox.lower() - dstbbox.lower()) / dstbbox.stride(); + + + + int const srciext = srcext[0]; + int const srcjext = srcext[1]; + int const srckext = srcext[2]; + + int const dstiext = dstext[0]; + int const dstjext = dstext[1]; + int const dstkext = dstext[2]; + + int const regiext = regext[0]; + int const regjext = regext[1]; + int const regkext = regext[2]; + + int const srcioff = srcoff[0]; + int const srcjoff = srcoff[1]; + int const srckoff = srcoff[2]; + + int const dstioff = dstoff[0]; + int const dstjoff = dstoff[1]; + int const dstkoff = dstoff[2]; + + // size_t const srcdi == SRCOFF3(1,0,0) - SRCOFF3(0,0,0); + size_t const srcdi = 1; + assert (srcdi == SRCOFF3(1,0,0) - SRCOFF3(0,0,0)); + size_t const srcdj = SRCOFF3(0,1,0) - SRCOFF3(0,0,0); + size_t const srcdk = SRCOFF3(0,0,1) - SRCOFF3(0,0,0); + + + + // Loop over coarse region +#pragma omp parallel for //collapse(3) + for (int k=0; k= 0 and + 2 * j - 2 + srcjoff >= 0 and + 2 * i - 2 + srcioff >= 0)) + { + cout << "restrict_3d_cc_o3_rf2.cc\n"; + cout << "regext " << regext << "\n"; + cout << "srcext " << srcext << "\n"; + cout << "srcbbox=" << srcbbox << "\n"; + cout << "dstbbox=" << dstbbox << "\n"; + cout << "regbbox=" << regbbox << "\n"; + cout << "i,j,k=" << i << " " << j << " " << k << "\n"; + assert(2 * k + 2 + srckoff < srckext); + assert(2 * j + 2 + srcjoff < srcjext); + assert(2 * i + 2 + srcioff < srciext); + assert(2 * k - 1 + srckoff >= 0); + assert(2 * j - 1 + srcjoff >= 0); + assert(2 * i - 1 + srcioff >= 0); + } +#endif + dst [DSTIND3(i, j, k)] = + restrict3 + (& src[SRCIND3(2*i, 2*j, 2*k)], srcdi, srcdj, srcdk); + + } + } + } + + } + + + +#define TYPECASE(N,T) \ + template \ + void \ + restrict_3d_cc_o5_rf2 (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, \ + ibbox3 const & restrict regbbox, \ + void * extraargs); +#include "typecase.hh" +#undef TYPECASE + + + +} // namespace CarpetLib -- cgit v1.2.3