diff options
Diffstat (limited to 'Carpet/CarpetInterp2/src/fasterp.cc')
-rw-r--r-- | Carpet/CarpetInterp2/src/fasterp.cc | 60 |
1 files changed, 30 insertions, 30 deletions
diff --git a/Carpet/CarpetInterp2/src/fasterp.cc b/Carpet/CarpetInterp2/src/fasterp.cc index be0acfc74..54848d6dc 100644 --- a/Carpet/CarpetInterp2/src/fasterp.cc +++ b/Carpet/CarpetInterp2/src/fasterp.cc @@ -84,7 +84,7 @@ namespace CarpetInterp2 { } dist::mpi_struct_descr_t const descr[] = { ENTRY(int, mrc), -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK ENTRY(int, pn), ENTRY(int, ipos), ENTRY(int, ind), @@ -109,7 +109,7 @@ namespace CarpetInterp2 { { os << "fasterp_iloc_t{" << "mrc=" << mrc << "," -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK << "pn=" << pn << "," << "ipos=" << ipos << "," << "ind=" << ind << "," @@ -198,7 +198,7 @@ namespace CarpetInterp2 { assert (order <= max_order); CCTK_REAL const eps = 1.0e-12; -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK mrc = iloc.mrc; pn = iloc.pn; ipos = iloc.ipos; @@ -222,7 +222,7 @@ namespace CarpetInterp2 { exact[d] = true; } -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK ind = iloc.ind; #endif ind3d = iloc.ind3d; @@ -249,8 +249,7 @@ namespace CarpetInterp2 { // round is not available with PGI compilers // CCTK_REAL const rx = round(x); CCTK_REAL const rx = floor(x+0.5); -#warning "TODO: make eps a relative error my taking max(abs(xmin),abs(xmax)) into account" - if (abs(x - rx) < eps) { + if (abs(x - rx) < eps * (1.0 + abs(x))) { // The interpolation point coincides with a grid point; no // interpolation is necessary (this is a special case) iorigin[d] += int(rx); @@ -278,11 +277,12 @@ namespace CarpetInterp2 { } // Set 3D location of stencil anchor -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK ind = iloc.ind + iorigin; if (not (all (ind>=0 and ind+either(exact,0,order)<lsh))) { stringstream buf; - buf << "*this=" << *this << " order=" << order << " lsh=" << lsh; + buf << "*this=" << *this << " iloc=" << iloc << " " + << "lsh=" << lsh << " order=" << order; CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, "Array access out of bounds: %s. Most likely, the interpolation point is too close to an outer or to a symmetry/inter-patch boundary. More information follows...", buf.str().c_str()); @@ -290,7 +290,7 @@ namespace CarpetInterp2 { } #endif ind3d = iloc.ind3d + index(lsh, iorigin); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK assert (index(lsh,ind) == ind3d); #endif return 0; @@ -314,7 +314,7 @@ namespace CarpetInterp2 { CCTK_REAL * restrict const vals) const { -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK assert (all (lsh == saved_lsh)); #endif assert (O0 == 0 or not exact[0]); @@ -325,7 +325,7 @@ namespace CarpetInterp2 { size_t const dj = di * lsh[0]; size_t const dk = dj * lsh[1]; -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK assert (all (ind>=0 and ind+ivect(O0,O1,O2)<lsh)); assert (ind3d == index(lsh,ind)); assert (int(di*ind[0] + dj*ind[1] + dk*ind[2]) == index(lsh,ind)); @@ -352,7 +352,7 @@ namespace CarpetInterp2 { } } } -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK # if 0 for (size_t v=0; v<varptrs.size(); ++v) { vals[v] = ind[0] + 1000 * (ind[1] + 1000 * ind[2]); @@ -467,14 +467,14 @@ namespace CarpetInterp2 { } os << "],"; os << "exact=" << exact << ","; -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK os << "pn=" << pn << ","; os << "mrc=" << mrc << ","; os << "ipos=" << ipos << ","; os << "ind=" << ind << ","; #endif os << "ind3d=" << ind3d; -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK os << "," << "saved_lsh=" << saved_lsh; #endif os << "}"; @@ -710,7 +710,7 @@ namespace CarpetInterp2 { // Store result fasterp_iloc_t & iloc = ilocs.AT(n); iloc.mrc = mrc_t(m, rl, c); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK iloc.pn.p = dist::rank(); iloc.pn.n = n; iloc.ipos = ipos * ext.stride(); @@ -855,7 +855,7 @@ namespace CarpetInterp2 { fasterp_iloc_t::mpi_datatype(), comm_world); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK // Ensure that the ilocs were sent to the right processor for (int p=0; p<nprocs; ++p) { #pragma omp parallel for @@ -962,7 +962,7 @@ namespace CarpetInterp2 { int const ierr = sloc.calc_stencil (iloc, send_comp.lsh, order); if (ierr) { CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, - "Could not determind valid interpolation stencil for point %d on map %d, refinement level %d, component %d", + "Could not determine valid interpolation stencil for point %d on map %d, refinement level %d, component %d", n, iloc.mrc.m, iloc.mrc.rl, iloc.mrc.c); } send_comp.locs.push_back (sloc); @@ -986,7 +986,7 @@ namespace CarpetInterp2 { } #endif -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK for (int comp=0; comp<ncomps; ++comp) { send_comp_t const & send_comp = send_proc.comps.at(comp); assert (int(send_comp.locs.size()) == send_comp.npoints); @@ -1000,7 +1000,7 @@ namespace CarpetInterp2 { } // for pp -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK { if (verbose) CCTK_INFO ("Compare send and receive counts"); @@ -1125,7 +1125,7 @@ namespace CarpetInterp2 { vector<CCTK_REAL> recv_points (recv_descr.npoints * nvars); fill_with_poison (recv_points); vector<MPI_Request> recv_reqs (recv_descr.procs.size()); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK vector<pn_t> recv_pn (recv_descr.npoints); vector<MPI_Request> recv_reqs_pn (recv_descr.procs.size()); #endif @@ -1136,7 +1136,7 @@ namespace CarpetInterp2 { recv_proc.npoints * nvars, dist::mpi_datatype<CCTK_REAL>(), recv_proc.p, mpi_tag, comm_world, & recv_reqs.AT(pp)); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK MPI_Irecv (& recv_pn.AT(recv_proc.offset), recv_proc.npoints * 2, dist::mpi_datatype<int>(), recv_proc.p, mpi_tag, @@ -1150,7 +1150,7 @@ namespace CarpetInterp2 { vector<CCTK_REAL> send_points (send_descr.npoints * nvars); fill_with_poison (send_points); vector<MPI_Request> send_reqs (send_descr.procs.size()); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK vector<pn_t> send_pn (send_descr.npoints); vector<MPI_Request> send_reqs_pn (send_descr.procs.size()); #endif @@ -1159,7 +1159,7 @@ namespace CarpetInterp2 { vector<CCTK_REAL> computed_points (send_proc.npoints * nvars); fill_with_poison (computed_points); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK vector<pn_t> computed_pn (send_descr.npoints); #endif for (size_t comp=0; comp<send_proc.comps.size(); ++comp) { @@ -1192,7 +1192,7 @@ namespace CarpetInterp2 { size_t const ind = (send_comp.offset + n) * nvars; send_comp.locs.AT(n).interpolate (send_comp.lsh, order, varptrs, &computed_points.AT(ind)); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK computed_pn.AT(send_comp.offset + n) = send_comp.locs.AT(n).pn; #endif } @@ -1207,14 +1207,14 @@ namespace CarpetInterp2 { send_points.AT((send_proc.offset + n) * nvars + v) = computed_points.AT(nn * nvars + v); } -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK send_pn.AT(send_proc.offset + n) = computed_pn.AT(nn); #endif } // TODO: Use a scatter index list instead of a gather index // list, and combine the scattering with the calculation -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK #pragma omp parallel for for (int n=0; n<send_proc.npoints; ++n) { assert (send_pn.AT(send_proc.offset + n).p == send_proc.p); @@ -1229,7 +1229,7 @@ namespace CarpetInterp2 { send_proc.npoints * nvars, dist::mpi_datatype<CCTK_REAL>(), send_proc.p, mpi_tag, comm_world, & send_reqs.AT(pp)); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK MPI_Isend (& send_pn.AT(send_proc.offset), send_proc.npoints * 2, dist::mpi_datatype<int>(), send_proc.p, mpi_tag, @@ -1240,7 +1240,7 @@ namespace CarpetInterp2 { // Wait for Irecvs to complete if (verbose) CCTK_INFO ("Waiting for MPI_Irevcs to complete"); MPI_Waitall (recv_reqs.size(), & recv_reqs.front(), MPI_STATUSES_IGNORE); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK MPI_Waitall (recv_reqs.size(), & recv_reqs_pn.front(), MPI_STATUSES_IGNORE); #endif @@ -1255,7 +1255,7 @@ namespace CarpetInterp2 { recv_points.AT(nn * nvars + v) = poison; #endif } -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK assert (recv_pn.AT(nn).p == dist::rank()); assert (recv_pn.AT(nn).n == n); #endif @@ -1264,7 +1264,7 @@ namespace CarpetInterp2 { // Wait for Isends to complete if (verbose) CCTK_INFO ("Waiting for MPI_Isends to complete"); MPI_Waitall (send_reqs.size(), & send_reqs.front(), MPI_STATUSES_IGNORE); -#ifdef CARPETINTER2_CHECK +#ifdef CARPETINTERP2_CHECK MPI_Waitall (send_reqs.size(), & send_reqs_pn.front(), MPI_STATUSES_IGNORE); #endif |