Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

CarpetX: Add interpolator self-checks #326

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open

CarpetX: Add interpolator self-checks #326

Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
9db1121
CarpetX: Correct error check in interpolator
eschnett Nov 4, 2024
cf16d7a
CarpetX: Simplify error output in interpolator
eschnett Nov 4, 2024
2708b3f
CarpetX: Slight performance improvement in interpolator
eschnett Nov 4, 2024
757b12b
CarpetX: Improve safety/reduce warnings in interpolator
eschnett Nov 4, 2024
82db171
TestInterpolate: Test cell-centred interpolation
eschnett Nov 4, 2024
8395217
CarpetX: Run interpolator self-tests when built in debug mode
eschnett Dec 23, 2024
34e8e0b
Merge branch 'main' into eschnett/interpolator
eschnett Dec 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 59 additions & 19 deletions CarpetX/src/interpolate.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,12 +39,34 @@ template <typename T, int order, int centering> 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<const T> &vars;
const vect<int, dim> &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<const T> &vars, const vect<int, dim> &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<T>::epsilon(), T(3) / 4);
Expand All @@ -65,8 +87,9 @@ template <typename T, int order, int centering> 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)
Expand Down Expand Up @@ -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())];

Expand Down Expand Up @@ -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());
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd remove or move behind a #ifdef DEBUG this code.

for (const auto &particle : particles) {
oldids.insert(particle.id());
// CCTK_VINFO(" id=%d proc=%d pos=[%g,%g,%g] locals=[%g,%g,%g] "
Expand Down Expand Up @@ -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());
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd remove or move behind a #ifdef DEBUG this code.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(also similar in line 611ff below)

for (const auto &particle : particles) {
newids.insert(particle.id());
// CCTK_VINFO(" id=%d proc=%d pos=[%g,%g,%g] locals=[%g,%g,%g] "
Expand All @@ -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",
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -683,31 +712,36 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
switch (interpolation_order) {
case 0: {
const interpolator<CCTK_REAL, 0, 0b000> 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<CCTK_REAL, 1, 0b000> 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<CCTK_REAL, 2, 0b000> 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<CCTK_REAL, 3, 0b000> 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<CCTK_REAL, 4, 0b000> interp{
grid, vi, vars, derivs, bool(allow_boundaries)};
grid, gi, vi, patch,
level, vars, derivs, bool(allow_boundaries)};
interp.interpolate3d(particles, varresult);
break;
}
Expand All @@ -726,31 +760,36 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
switch (interpolation_order) {
case 0: {
const interpolator<CCTK_REAL, 0, 0b111> 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<CCTK_REAL, 1, 0b111> 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<CCTK_REAL, 2, 0b111> 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<CCTK_REAL, 3, 0b111> 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<CCTK_REAL, 4, 0b111> interp{
grid, vi, vars, derivs, bool(allow_boundaries)};
grid, gi, vi, patch,
level, vars, derivs, bool(allow_boundaries)};
interp.interpolate3d(particles, varresult);
break;
}
Expand Down Expand Up @@ -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<CCTK_REAL *const *>(resultptrs_);
assert(int(recvbuf.size()) == (nvars + 1) * npoints);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

assert is too late the code in CCTK_DEBUG in line 876 alreay touches all points.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This assert can trigger if there are points outside of the domain, can't it? 

In that case this assert does not adhere to the API set out in the ReferenceManual for CCTK InterpGridArrays which documents CCTK ERROR INTERP POINT OUTSIDE as a possible return value and does not indicate that the code would abort.

for (int n = 0; n < npoints; ++n) {
const int offset = (nvars + 1) * n;
const int idx = int(recvbuf.at(offset));
Expand Down
Loading