diff --git a/CarpetX/src/interpolate.cxx b/CarpetX/src/interpolate.cxx index de484f68d..6df61fc1c 100644 --- a/CarpetX/src/interpolate.cxx +++ b/CarpetX/src/interpolate.cxx @@ -39,12 +39,34 @@ template struct interpolator { (centering & 0b001) != 0}; const GridDescBase &grid; +#ifdef CCTK_DEBUG + const int gi; +#endif const int vi; +#ifdef CCTK_DEBUG + const int patch; + const int level; +#endif const amrex::Array4 &vars; const vect &derivs; // Allow outer boundaries as interpolation sources const bool allow_boundaries; + interpolator(const GridDescBase &grid, CCTK_ATTRIBUTE_UNUSED const int gi, + const int vi, const int patch, const int level, + const amrex::Array4 &vars, const vect &derivs, + const bool allow_boundaries) + : grid(grid), +#ifdef CCTK_DEBUG + gi(gi), +#endif + vi(vi), +#ifdef CCTK_DEBUG + patch(patch), level(level), +#endif + vars(vars), derivs(derivs), allow_boundaries(allow_boundaries) { + } + static constexpr T eps() { using std::pow; return pow(std::numeric_limits::epsilon(), T(3) / 4); @@ -65,8 +87,9 @@ template struct interpolator { #ifdef CCTK_DEBUG using std::isfinite; if (!(isfinite(val))) { - std::cerr << "!isfinite i=" << i << " di=" << di << " val=" << val - << "\n"; + std::cerr << "!isfinite gi=" << gi + << " groupname=" << CCTK_FullGroupName(gi) << " vi=" << vi + << " i=" << i << " di=" << di << " val=" << val << "\n"; for (int c = -1; c <= +1; ++c) for (int b = -1; b <= +1; ++b) for (int a = -1; a <= +1; ++a) @@ -504,7 +527,10 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const int level = 0; const auto &restrict leveldata = patchdata.leveldata.at(level); const amrex::MFIter mfi(*leveldata.fab); - assert(mfi.isValid()); + // The mfi can be invalid if the number of processes does not evenly divide + // the number of blocks + if (!mfi.isValid()) + continue; ParticleTile *const particle_tile = &containers.at(patch).GetParticles( level)[make_pair(mfi.index(), mfi.LocalTileIndex())]; @@ -544,8 +570,8 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, ++pti) { const auto &particles = pti.GetArrayOfStructs(); const int component = MFPointer(pti).index(); - CCTK_VINFO("patch %d level %d component %d old_nparticles: %zu", - patch, lev, component, particles.size()); + // CCTK_VINFO("patch %d level %d component %d old_nparticles: %zu", + // patch, lev, component, particles.size()); for (const auto &particle : particles) { oldids.insert(particle.id()); // CCTK_VINFO(" id=%d proc=%d pos=[%g,%g,%g] locals=[%g,%g,%g] " @@ -578,8 +604,8 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, ++pti) { const int component = MFPointer(pti).index(); const auto &particles = pti.GetArrayOfStructs(); - CCTK_VINFO("patch %d level %d component %d new_nparticles: %zu", - patch, lev, component, particles.size()); + // CCTK_VINFO("patch %d level %d component %d new_nparticles: %zu", + // patch, lev, component, particles.size()); for (const auto &particle : particles) { newids.insert(particle.id()); // CCTK_VINFO(" id=%d proc=%d pos=[%g,%g,%g] locals=[%g,%g,%g] " @@ -600,10 +626,12 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, if (new_nparticles != old_nparticles) { for (const auto oldid : oldids) if (!newids.count(oldid)) - CCTK_VINFO("old id %d not present in new ids", oldid); + CCTK_VWARN(CCTK_WARN_ALERT, "old id %d not present in new ids", + oldid); for (const auto newid : newids) if (!oldids.count(newid)) - CCTK_VINFO("new id %d not present in old ids", newid); + CCTK_VWARN(CCTK_WARN_ALERT, "new id %d not present in old ids", + newid); CCTK_VERROR( "We lost interpolation points on patch %d. Before redistributing: " "%zu particles, after redistributing: %zu particles", @@ -637,6 +665,7 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const int patch = patchdata.patch; for (const auto &leveldata : patchdata.leveldata) { const int level = leveldata.level; + // TODO: use OpenMP for (amrex::ParIter<3, 2> pti(containers.at(patch), level); pti.isValid(); ++pti) { @@ -683,31 +712,36 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, switch (interpolation_order) { case 0: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 1: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 2: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 3: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 4: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } @@ -726,31 +760,36 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, switch (interpolation_order) { case 0: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 1: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 2: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 3: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 4: { const interpolator interp{ - grid, vi, vars, derivs, bool(allow_boundaries)}; + grid, gi, vi, patch, + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } @@ -850,6 +889,7 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, // Set result CCTK_REAL *const restrict *const restrict resultptrs = static_cast(resultptrs_); + assert(int(recvbuf.size()) == (nvars + 1) * npoints); for (int n = 0; n < npoints; ++n) { const int offset = (nvars + 1) * n; const int idx = int(recvbuf.at(offset));