diff options
Diffstat (limited to 'Carpet')
-rw-r--r-- | Carpet/CarpetLib/src/make.code.defn | 1 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc | 314 |
2 files changed, 315 insertions, 0 deletions
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index 81a0352ad..0e7153bef 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -30,6 +30,7 @@ SRCS = bbox.cc \ interpolate_eno_3d_3tl.cc \ restrict_3d_cc_rf2.cc \ restrict_3d_rf2.cc \ + restrict_3d_vc_rf2.cc \ restrict_4d_rf2.cc \ prolongate_3d_cc_rf2.cc \ prolongate_3d_o1_rf2.cc \ diff --git a/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc new file mode 100644 index 000000000..b7b92e4e3 --- /dev/null +++ b/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc @@ -0,0 +1,314 @@ +#include <algorithm> +#include <cassert> +#include <cmath> + +#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) + + + + // 0D "restriction" + template <typename T> + struct restrict0 { + static inline T call (T const * restrict const p) + { + return * p; + } + }; + + // 1D restriction + template <typename T, int centi> + struct restrict1 { + static inline T call (T const * restrict const p, + size_t const d1); + }; + template <typename T> + struct restrict1<T,0> { + static inline T call (T const * restrict const p, + size_t const d1) + { + T const res = restrict0<T>::call (p); + return res; + } + }; + template <typename T> + struct restrict1<T,1> { + static inline T call (T const * restrict const p, + size_t const d1) + { + typedef typename typeprops<T>::real RT; + RT const one = 1; + RT const half = one/2; + T const res = + half * restrict0<T>::call (p ) + + half * restrict0<T>::call (p + d1); + return res; + } + }; + + // 2D restriction + template <typename T, int centi, int centj> + struct restrict2 { + static inline T call (T const * restrict const p, + size_t const d1, size_t const d2); + }; + template <typename T, int centi> + struct restrict2<T,centi,0> { + static inline T call (T const * restrict const p, + size_t const d1, size_t const d2) + { + T const res = restrict1<T,centi>::call (p, d1); + return res; + } + }; + template <typename T, int centi> + struct restrict2<T,centi,1> { + static inline T call (T const * restrict const p, + size_t const d1, size_t const d2) + { + typedef typename typeprops<T>::real RT; + RT const one = 1; + RT const half = one/2; + T const res = + half * restrict1<T,centi>::call (p , d1) + + half * restrict1<T,centi>::call (p + d2, d1); + return res; + } + }; + + // 3D restriction + template <typename T, int centi, int centj, int centk> + struct restrict3 { + static inline T call (T const * restrict const p, + size_t const d1, size_t const d2, size_t const d3); + }; + template <typename T, int centi, int centj> + struct restrict3<T,centi,centj,0> { + static inline T call (T const * restrict const p, + size_t const d1, size_t const d2, size_t const d3) + { + T const res = restrict2<T,centi,centj>::call (p, d1, d2); + return res; + } + }; + template <typename T, int centi, int centj> + struct restrict3<T,centi,centj,1> { + static inline T call (T const * restrict const p, + size_t const d1, size_t const d2, size_t const d3) + { + typedef typename typeprops<T>::real RT; + RT const one = 1; + RT const half = one/2; + T const res = + half * restrict2<T,centi,centj>::call (p , d1, d2) + + half * restrict2<T,centi,centj>::call (p + d3, d1, d2); + return res; + } + }; + + + + template <typename T, int centi, int centj, int centk> + void + restrict_3d_vc_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 regbbox) + { + DECLARE_CCTK_PARAMETERS; + + typedef typename typeprops<T>::real RT; + ivect const icent (centi, centj, centk); + assert (all (icent == 0 or icent == 1)); + bvect const cent (icent); + + + + 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).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 (either (cent, srcbbox.stride() % 2 == 0, true))); + assert (all ((regbbox.lower() - srcbbox.lower() - + either (cent, srcbbox.stride() / 2, 0)) % + srcbbox.stride() == 0)); + ivect3 const srcoff = + (regbbox.lower() - srcbbox.lower() - + either (cent, srcbbox.stride() / 2, 0)) / + 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 = 1; + size_t const srcdj = srcdi * srciext; + size_t const srcdk = srcdj * srcjext; + + + + // Loop over coarse 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)] = + restrict3<T,centi,centj,centk>::call + (& src[SRCIND3(2*i, 2*j, 2*k)], srcdi, srcdj, srcdk); + + } + } + } + + } + + + +#define TYPECASE(N,T) \ + template \ + void \ + restrict_3d_vc_rf2<T,0,0,0> (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 \ + void \ + restrict_3d_vc_rf2<T,0,0,1> (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 \ + void \ + restrict_3d_vc_rf2<T,0,1,0> (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 \ + void \ + restrict_3d_vc_rf2<T,0,1,1> (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 \ + void \ + restrict_3d_vc_rf2<T,1,0,0> (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 \ + void \ + restrict_3d_vc_rf2<T,1,0,1> (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 \ + void \ + restrict_3d_vc_rf2<T,1,1,0> (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 \ + void \ + restrict_3d_vc_rf2<T,1,1,1> (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 |