diff options
author | Roland Haas <roland.haas@physics.gatech.edu> | 2012-06-02 17:36:04 -0700 |
---|---|---|
committer | Barry Wardell <barry.wardell@gmail.com> | 2012-09-11 18:23:22 +0100 |
commit | 7ebd620a45f689e1e96831dc1963669d32fe44a3 (patch) | |
tree | 1d4f533d0c280d2ed7f14d461c40e8db827f89fb | |
parent | f8631186ea6e25bbe375d94dc8b78bb33f318267 (diff) |
CarpetLib: add preliminary support for higher order restriction
-rw-r--r-- | Carpet/CarpetLib/param.ccl | 6 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/data.cc | 13 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/dh.cc | 11 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/make.code.defn | 1 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/operator_prototypes_3d.hh | 14 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc | 246 |
6 files changed, 291 insertions, 0 deletions
diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl index 2ad563960..b0a825171 100644 --- a/Carpet/CarpetLib/param.ccl +++ b/Carpet/CarpetLib/param.ccl @@ -33,6 +33,12 @@ BOOLEAN use_dgfe "Use DGFE operators instead of Lagrange operators" +BOOLEAN use_cc_o3 "Use third order cell centered restriction operators instead of first order" +{ +} "no" + + + BOOLEAN output_bboxes "Output bounding box information to the screen" STEERABLE=always { } "no" diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 4a0c53e45..9a26eec8a 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -1141,6 +1141,19 @@ transfer_restrict (data const * const src, srcregbox, dstregbox, NULL); break; } + if (use_cc_o3 and + transport_operator != op_WENO and + transport_operator != op_ENO) { // HACK + // Don't use call_operator, because we parallelise ourselves + restrict_3d_cc_o3_rf2(static_cast <T const *> (src->storage()), + src->shape(), + static_cast <T *> (this->storage()), + this->shape(), + srcbox, + dstbox, + srcregbox, dstregbox, NULL); + break; + } // Don't use call_operator, because we parallelise ourselves restrict_3d_cc_rf2(static_cast <T const *> (src->storage()), src->shape(), diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 704ce8ee4..572fe2e75 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -1034,13 +1034,20 @@ regrid (bool const do_init) full_dboxes const& obox = full_olevel.AT(oc); ibset needrecv = allrestricted & obox.owned; + if(use_cc_o3) { + // NOTE: change in behaviour (affects only outer boundaries I think)!!!! + // NOTE: b/c of this we need a low-level sync after the restrict + needrecv = allrestricted & obox.interior; + } // Cannot restrict into buffer zones assert ((allrestricted & obox.buffers).empty()); for (int c = 0; c < h.components(rl); ++ c) { full_dboxes const& box = full_level.AT(c); + // HORRIBLE HACK ibbox const contracted_exterior = + use_cc_o3 ? box.interior.expand(ivect(int(h.refcent==cell_centered))).contracted_for(odomext) : box.exterior.contracted_for(odomext); ibset const ovlp = needrecv & contracted_exterior; @@ -1051,6 +1058,10 @@ regrid (bool const do_init) ibbox const send = recv.expanded_for(box.exterior); ASSERT_c (send <= box.exterior, "Refinement restriction: Send region must be contained in exterior"); + if(use_cc_o3) { + ASSERT_c (send <= box.interior.expand(ivect(int(h.refcent==cell_centered))), + "Refinement restriction: Send region must be contained in interior"); + } sendrecv_pseudoregion_t const preg (send, c, recv, oc); fast_olevel.fast_ref_rest_sendrecv.push_back(preg); diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index 67c93abe6..89bdca337 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -34,6 +34,7 @@ SRCS = backtrace.cc \ interpolate_3d_5tl.cc \ interpolate_eno_3d_3tl.cc \ restrict_3d_cc_rf2.cc \ + restrict_3d_cc_o3_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 dce5b5c51..b9cb7d343 100644 --- a/Carpet/CarpetLib/src/operator_prototypes_3d.hh +++ b/Carpet/CarpetLib/src/operator_prototypes_3d.hh @@ -372,6 +372,20 @@ namespace CarpetLib { + template <typename T> + void + restrict_3d_cc_o3_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 <typename T, int centi, int centj, int centk> void restrict_3d_vc_rf2 (T const * restrict const src, diff --git a/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc new file mode 100644 index 000000000..c8ecc040f --- /dev/null +++ b/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc @@ -0,0 +1,246 @@ +#include <cctk.h> +#include <cctk_Parameters.h> + +#include <algorithm> +#include <cassert> +#include <cmath> + +#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) + + + + template <typename T> + void + restrict_3d_cc_o3_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<T>::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]; + + + + RT const den = 4096; + + RT const f1 = 1/den; + RT const f2 = 9/den; + RT const f3 = 81/den; + RT const f4 = 729/den; + + + + // Loop over coarse region +#pragma omp parallel for //collapse(3) + for (int k=0; k<regkext; ++k) { + for (int j=0; j<regjext; ++j) { + for (int i=0; i<regiext; ++i) { + +#if(1 || defined(CARPET_DEBUG)) + if(not (2 * k + 2 + srckoff < srckext and + 2 * j + 2 + srcjoff < srcjext and + 2 * i + 2 + srcioff < srciext and + 2 * k - 1 + srckoff >= 0 and + 2 * j - 1 + srcjoff >= 0 and + 2 * i - 1 + 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)] = + - f1 * src [SRCIND3(2*i-1, 2*j-1, 2*k-1)] + + f2 * src [SRCIND3(2*i , 2*j-1, 2*k-1)] + + f2 * src [SRCIND3(2*i+1, 2*j-1, 2*k-1)] + - f1 * src [SRCIND3(2*i+2, 2*j-1, 2*k-1)] + + f2 * src [SRCIND3(2*i-1, 2*j , 2*k-1)] + - f3 * src [SRCIND3(2*i , 2*j , 2*k-1)] + - f3 * src [SRCIND3(2*i+1, 2*j , 2*k-1)] + + f2 * src [SRCIND3(2*i+2, 2*j , 2*k-1)] + + f2 * src [SRCIND3(2*i-1, 2*j+1, 2*k-1)] + - f3 * src [SRCIND3(2*i , 2*j+1, 2*k-1)] + - f3 * src [SRCIND3(2*i+1, 2*j+1, 2*k-1)] + + f2 * src [SRCIND3(2*i+2, 2*j+1, 2*k-1)] + - f1 * src [SRCIND3(2*i-1, 2*j+2, 2*k-1)] + + f2 * src [SRCIND3(2*i , 2*j+2, 2*k-1)] + + f2 * src [SRCIND3(2*i+1, 2*j+2, 2*k-1)] + - f1 * src [SRCIND3(2*i+2, 2*j+2, 2*k-1)] + + f2 * src [SRCIND3(2*i-1, 2*j-1, 2*k )] + - f3 * src [SRCIND3(2*i , 2*j-1, 2*k )] + - f3 * src [SRCIND3(2*i+1, 2*j-1, 2*k )] + + f2 * src [SRCIND3(2*i+2, 2*j-1, 2*k )] + - f3 * src [SRCIND3(2*i-1, 2*j , 2*k )] + + f4 * src [SRCIND3(2*i , 2*j , 2*k )] + + f4 * src [SRCIND3(2*i+1, 2*j , 2*k )] + - f3 * src [SRCIND3(2*i+2, 2*j , 2*k )] + - f3 * src [SRCIND3(2*i-1, 2*j+1, 2*k )] + + f4 * src [SRCIND3(2*i , 2*j+1, 2*k )] + + f4 * src [SRCIND3(2*i+1, 2*j+1, 2*k )] + - f3 * src [SRCIND3(2*i+2, 2*j+1, 2*k )] + + f2 * src [SRCIND3(2*i-1, 2*j+2, 2*k )] + - f3 * src [SRCIND3(2*i , 2*j+2, 2*k )] + - f3 * src [SRCIND3(2*i+1, 2*j+2, 2*k )] + + f2 * src [SRCIND3(2*i+2, 2*j+2, 2*k )] + + f2 * src [SRCIND3(2*i-1, 2*j-1, 2*k+1)] + - f3 * src [SRCIND3(2*i , 2*j-1, 2*k+1)] + - f3 * src [SRCIND3(2*i+1, 2*j-1, 2*k+1)] + + f2 * src [SRCIND3(2*i+2, 2*j-1, 2*k+1)] + - f3 * src [SRCIND3(2*i-1, 2*j , 2*k+1)] + + f4 * src [SRCIND3(2*i , 2*j , 2*k+1)] + + f4 * src [SRCIND3(2*i+1, 2*j , 2*k+1)] + - f3 * src [SRCIND3(2*i+2, 2*j , 2*k+1)] + - f3 * src [SRCIND3(2*i-1, 2*j+1, 2*k+1)] + + f4 * src [SRCIND3(2*i , 2*j+1, 2*k+1)] + + f4 * src [SRCIND3(2*i+1, 2*j+1, 2*k+1)] + - f3 * src [SRCIND3(2*i+2, 2*j+1, 2*k+1)] + + f2 * src [SRCIND3(2*i-1, 2*j+2, 2*k+1)] + - f3 * src [SRCIND3(2*i , 2*j+2, 2*k+1)] + - f3 * src [SRCIND3(2*i+1, 2*j+2, 2*k+1)] + + f2 * src [SRCIND3(2*i+2, 2*j+2, 2*k+1)] + - f1 * src [SRCIND3(2*i-1, 2*j-1, 2*k+2)] + + f2 * src [SRCIND3(2*i , 2*j-1, 2*k+2)] + + f2 * src [SRCIND3(2*i+1, 2*j-1, 2*k+2)] + - f1 * src [SRCIND3(2*i+2, 2*j-1, 2*k+2)] + + f2 * src [SRCIND3(2*i-1, 2*j , 2*k+2)] + - f3 * src [SRCIND3(2*i , 2*j , 2*k+2)] + - f3 * src [SRCIND3(2*i+1, 2*j , 2*k+2)] + + f2 * src [SRCIND3(2*i+2, 2*j , 2*k+2)] + + f2 * src [SRCIND3(2*i-1, 2*j+1, 2*k+2)] + - f3 * src [SRCIND3(2*i , 2*j+1, 2*k+2)] + - f3 * src [SRCIND3(2*i+1, 2*j+1, 2*k+2)] + + f2 * src [SRCIND3(2*i+2, 2*j+1, 2*k+2)] + - f1 * src [SRCIND3(2*i-1, 2*j+2, 2*k+2)] + + f2 * src [SRCIND3(2*i , 2*j+2, 2*k+2)] + + f2 * src [SRCIND3(2*i+1, 2*j+2, 2*k+2)] + - f1 * src [SRCIND3(2*i+2, 2*j+2, 2*k+2)]; + + } + } + } + + } + + + +#define TYPECASE(N,T) \ + template \ + void \ + restrict_3d_cc_o3_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 |