diff options
Diffstat (limited to 'Carpet/CarpetLib/src')
-rw-r--r-- | Carpet/CarpetLib/src/accumulate_3d.cc | 146 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/data.cc | 43 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/gdata.cc | 4 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/ggf.cc | 39 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/ggf.hh | 12 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/make.code.defn | 1 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/operator_prototypes_3d.hh | 10 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/operators.hh | 1 |
8 files changed, 229 insertions, 27 deletions
diff --git a/Carpet/CarpetLib/src/accumulate_3d.cc b/Carpet/CarpetLib/src/accumulate_3d.cc new file mode 100644 index 000000000..f57cc06c0 --- /dev/null +++ b/Carpet/CarpetLib/src/accumulate_3d.cc @@ -0,0 +1,146 @@ +#include <algorithm> +#include <cassert> +#include <cmath> +#include <cstdlib> +#include <iostream> + +#include <cctk.h> +#include <cctk_Parameters.h> + +#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 + accumulate_3d (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 regbbox) + { + if (any (srcbbox.stride() != regbbox.stride() or + dstbbox.stride() != regbbox.stride())) + { +#pragma omp critical + { + cout << "accumulate_3d.cc:" << endl + << "srcbbox=" << srcbbox << endl + << "dstbbox=" << dstbbox << endl + << "regbbox=" << regbbox << endl; + 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()) { +#pragma omp critical + CCTK_WARN (0, "Internal error: region extent is empty"); + } + + if (not regbbox.is_contained_in(srcbbox) or + not regbbox.is_contained_in(dstbbox)) + { +#pragma omp critical + { + cout << "accumulate_3d.cc:" << endl + << "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 != srcbbox.shape() / srcbbox.stride() or + dstext != dstbbox.shape() / dstbbox.stride())) + { +#pragma omp critical + 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(); + + + + ptrdiff_t const srciext = srcext[0]; + ptrdiff_t const srcjext = srcext[1]; + ptrdiff_t const srckext = srcext[2]; + + ptrdiff_t const dstiext = dstext[0]; + ptrdiff_t const dstjext = dstext[1]; + ptrdiff_t const dstkext = dstext[2]; + + ptrdiff_t const regiext = regext[0]; + ptrdiff_t const regjext = regext[1]; + ptrdiff_t const regkext = regext[2]; + + ptrdiff_t const srcioff = srcoff[0]; + ptrdiff_t const srcjoff = srcoff[1]; + ptrdiff_t const srckoff = srcoff[2]; + + ptrdiff_t const dstioff = dstoff[0]; + ptrdiff_t const dstjoff = dstoff[1]; + ptrdiff_t const dstkoff = dstoff[2]; + + + + // Loop over region + for (int k=0; k<regkext; ++k) { + for (int j=0; j<regjext; ++j) { + for (int i=0; i<regiext; ++i) { + + dst [DSTIND3(i, j, k)] += src [SRCIND3(i, j, k)]; + + } + } + } + + } + + + +#define TYPECASE(N,T) \ + template \ + void \ + accumulate_3d (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 regbbox); +#include "typecase.hh" +#undef TYPECASE + + + +} // namespace CarpetLib diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 49229008a..2c38fcd9a 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -429,27 +429,38 @@ copy_from_innerloop (gdata const * const gsrc, assert (0); } + if (transport_operator != op_accumulate) { #if CARPET_DIM == 3 - call_operator<T> (& copy_3d, - static_cast <T const *> (src->storage()), - src->shape(), - static_cast <T *> (this->storage()), - this->shape(), - srcbox, - dstbox, - box); + call_operator<T> (& copy_3d, + static_cast <T const *> (src->storage()), + src->shape(), + static_cast <T *> (this->storage()), + this->shape(), + srcbox, + dstbox, + box); #elif CARPET_DIM == 4 - call_operator<T> (& copy_4d, - static_cast <T const *> (src->storage()), - src->shape(), - static_cast <T *> (this->storage()), - this->shape(), - srcbox, - dstbox, - box); + call_operator<T> (& copy_4d, + static_cast <T const *> (src->storage()), + src->shape(), + static_cast <T *> (this->storage()), + this->shape(), + srcbox, + dstbox, + box); #else # error "Value for CARPET_DIM not supported" #endif + } else { + call_operator<T> (& accumulate_3d, + static_cast <T const *> (src->storage()), + src->shape(), + static_cast <T *> (this->storage()), + this->shape(), + srcbox, + dstbox, + box); + } } diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index c9f9f8ad2..3d3c06524 100644 --- a/Carpet/CarpetLib/src/gdata.cc +++ b/Carpet/CarpetLib/src/gdata.cc @@ -305,7 +305,7 @@ find_source_timelevel (vector <CCTK_REAL> const & times, CCTK_REAL const max_time = * max_element (times.begin(), times.end()); // TODO: Use a real delta-time from somewhere instead of 1.0 CCTK_REAL const some_time = abs (min_time) + abs (max_time) + 1.0; - if (op != op_copy) { + if (op != op_copy and op != op_accumulate) { if (time < min_time - eps * some_time or time > max_time + eps * some_time) { @@ -333,7 +333,7 @@ find_source_timelevel (vector <CCTK_REAL> const & times, } } if (timelevel == -1) { - if (op == op_copy) { + if (op == op_copy or op == op_accumulate) { timelevel = 0; } } diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index ce6fe3ee3..314a04843 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -346,6 +346,26 @@ sync_all (comm_state & state, +// Accumulate the boundaries of all components +void +ggf:: +accumulate_all (comm_state & state, + int const tl, int const rl, int const ml) +{ + assert (transport_operator == op_accumulate); + // Accumulate + static Timer timer ("accumulate_all"); + timer.start (); + transfer_from_all (state, + tl,rl,ml, + & dh::fast_dboxes::fast_sync_sendrecv, + tl,rl,ml, + NULL, true); + timer.stop (0); +} + + + // Prolongate the boundaries of all components void ggf:: @@ -361,7 +381,7 @@ ref_bnd_prolongate_all (comm_state & state, vector<int> tl2s; static Timer timer ("ref_bnd_prolongate_all"); timer.start (); - if (transport_operator != op_copy) { + if (transport_operator != op_copy and transport_operator != op_accumulate) { // Interpolation in time if (not (timelevels(ml,rl) >= prolongation_order_time+1)) { char * const fullname = CCTK_FullName (varindex); @@ -518,7 +538,8 @@ transfer_from_all (comm_state & state, srpvect const dh::fast_dboxes::* sendrecvs, vector<int> const & tl2s, int const rl2, int const ml2, CCTK_REAL const & time, - mdata * const srcstorage_) + mdata * const srcstorage_, + bool const flip_send_recv) { assert (rl1>=0 and rl1<h.reflevels()); assert (ml1>=0 and ml1<h.mglevels()); @@ -558,9 +579,11 @@ transfer_from_all (comm_state & state, // Interpolation orders assert (transport_operator != op_none); int const pos = - transport_operator == op_copy ? 0 : d.prolongation_orders_space.AT(rl2); + transport_operator == op_copy or transport_operator == op_accumulate + ? 0 : d.prolongation_orders_space.AT(rl2); int const pot = - transport_operator == op_copy ? 0 : prolongation_order_time; + transport_operator == op_copy or transport_operator == op_accumulate + ? 0 : prolongation_order_time; vector<const gdata*> gsrcs(tl2s.size()); @@ -568,8 +591,12 @@ transfer_from_all (comm_state & state, for (srpvect::const_iterator ipsendrecv = psendrecvs.begin(); ipsendrecv!=psendrecvs.end(); ++ ipsendrecv) { - pseudoregion_t const & psend = (* ipsendrecv).send; - pseudoregion_t const & precv = (* ipsendrecv).recv; + typedef sendrecv_pseudoregion_t srp; + typedef pseudoregion_t sendrecv_pseudoregion_t::* srp_pr; + srp_pr const send_field = flip_send_recv ? &srp::recv : &srp::send; + srp_pr const recv_field = flip_send_recv ? &srp::send : &srp::recv; + pseudoregion_t const & psend = (* ipsendrecv).*send_field; + pseudoregion_t const & precv = (* ipsendrecv).*recv_field; ibbox const & send = psend.extent; ibbox const & recv = precv.extent; assert (all (send.stride() == h.baseextent(ml2,rl2).stride())); diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index d9412f822..56263fa80 100644 --- a/Carpet/CarpetLib/src/ggf.hh +++ b/Carpet/CarpetLib/src/ggf.hh @@ -124,6 +124,9 @@ public: // Synchronise the boundaries of a component void sync_all (comm_state& state, int tl, int rl, int ml); + // Accumulate from boundaries of a component + void accumulate_all (comm_state& state, int tl, int rl, int ml); + // Prolongate the boundaries of a component void ref_bnd_prolongate_all (comm_state& state, int tl, int rl, int ml, CCTK_REAL time); @@ -168,14 +171,16 @@ protected: srpvect const dh::fast_dboxes::* sendrecvs, vector<int> const & tl2s, int rl2, int ml2, CCTK_REAL const & time, - mdata * srcstorage = 0); + mdata * srcstorage = 0, + bool flip_send_recv = false); void transfer_from_all (comm_state & state, int tl1, int rl1, int ml1, srpvect const dh::fast_dboxes::* sendrecvs, int tl2, int rl2, int ml2, - mdata * srcstorage = 0) + mdata * srcstorage = 0, + bool flip_send_recv = false) { vector <int> tl2s(1); tl2s.AT(0) = tl2; @@ -185,7 +190,8 @@ protected: sendrecvs, tl2s, rl2, ml2, time, - srcstorage); + srcstorage, + flip_send_recv); } diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index 96d621172..f63ce2b62 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -25,6 +25,7 @@ SRCS = balance.cc \ th.cc \ timestat.cc \ vect.cc \ + accumulate_3d.cc \ copy_3d.cc \ copy_4d.cc \ interpolate_3d_2tl.cc \ diff --git a/Carpet/CarpetLib/src/operator_prototypes_3d.hh b/Carpet/CarpetLib/src/operator_prototypes_3d.hh index 81a384d76..c4793f606 100644 --- a/Carpet/CarpetLib/src/operator_prototypes_3d.hh +++ b/Carpet/CarpetLib/src/operator_prototypes_3d.hh @@ -58,6 +58,16 @@ namespace CarpetLib { ibbox3 const & restrict dstbbox, ibbox3 const & restrict regbbox); + template <typename T> + void + accumulate_3d (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 regbbox); + template <typename T> diff --git a/Carpet/CarpetLib/src/operators.hh b/Carpet/CarpetLib/src/operators.hh index 07e4e5616..df3828d44 100644 --- a/Carpet/CarpetLib/src/operators.hh +++ b/Carpet/CarpetLib/src/operators.hh @@ -12,6 +12,7 @@ enum operator_type op_restrict, // restrict only, do not prolongate op_copy, // use simple copying for prolongation // (needs only one time level) + op_accumulate, // accumulate (sum) into the destination op_Lagrange, // Lagrange interpolation (standard) op_ENO, // use ENO stencils (for hydro) op_WENO, // use WENO stencils (for hydro) |