Skip to content

CarpetX: Interpolation setup #369

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

Open
lucass-carneiro wants to merge 3 commits into
base: main
Choose a base branch
from
Open

CarpetX: Interpolation setup #369

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
079f9bb
CarpetX: Split CarpetX_Interpolate into setup and interpolation routines
lucass-carneiro Nov 29, 2025
5fd11a0
CarpetX: Fix debug mode compilation issues in interpolator changes
lucass-carneiro Nov 29, 2025
45a8aeb
CarpetX: Don't save Cactus data that should not be saved in interpola…
lucass-carneiro Nov 29, 2025
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
33 changes: 33 additions & 0 deletions CarpetX/src/interp.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,13 @@

#include <cctk.h>

#include <AMReX_AmrParticles.H>
#include <AMReX_Particles.H>

#include <vector>
#include <array>
#include <memory>

// Scheduled functions
extern "C" int CarpetX_InterpGridArrays(
cGH const *const cGH, int const N_dims, int const local_interp_handle,
Expand All @@ -23,6 +30,32 @@ extern "C" CCTK_INT CarpetX_DriverInterpolate(
CCTK_POINTER const output_arrays[]);

namespace CarpetX {

// Create a cache of data required for interpolation
struct InterpolationSetup {
using Container = amrex::AmrParticleContainer<3, 2>;
using Particle = Container::ParticleType;
using PinnedParticleTile = amrex::ParticleContainer_impl<
Particle, 0, 0, amrex::PinnedArenaAllocator>::ParticleTileType;
using ParticleTile = Container::ParticleTileType;

const CCTK_INT npoints{0};
std::vector<char> symmetry_reflected_z;
std::vector<std::shared_ptr<Container> > containers{};

InterpolationSetup(const CCTK_POINTER_TO_CONST cctkGH_,
const CCTK_INT npoints,
const CCTK_REAL *restrict const globalsx,
const CCTK_REAL *restrict const globalsy,
const CCTK_REAL *restrict const globalsz);
};

void InterpolateFromSetup(const InterpolationSetup &setup, const CCTK_INT nvars,
const CCTK_INT *restrict const varinds,
const CCTK_INT *restrict const operations,
const CCTK_INT allow_boundaries,
const CCTK_POINTER resultptrs_);

// a dummy routine for now
// TODO: implement this for actual local interpolation
int InterpLocalUniform(int N_dims, int param_table_handle,
Expand Down
113 changes: 68 additions & 45 deletions CarpetX/src/interpolate.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,7 @@
#include <array>
#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
#include <map>
#include <set>
#include <utility>
#include <vector>

Expand Down Expand Up @@ -433,16 +430,12 @@ extern "C" CCTK_INT CarpetX_DriverInterpolate(
return 0;
}

extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
const CCTK_INT npoints,
const CCTK_REAL *restrict const globalsx,
const CCTK_REAL *restrict const globalsy,
const CCTK_REAL *restrict const globalsz,
const CCTK_INT nvars,
const CCTK_INT *restrict const varinds,
const CCTK_INT *restrict const operations,
const CCTK_INT allow_boundaries,
const CCTK_POINTER resultptrs_) {
CarpetX::InterpolationSetup::InterpolationSetup(
const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT npoints,
const CCTK_REAL *restrict const globalsx,
const CCTK_REAL *restrict const globalsy,
const CCTK_REAL *restrict const globalsz)
: npoints(npoints) {
DECLARE_CCTK_PARAMETERS;
const cGH *restrict const cctkGH = static_cast<const cGH *>(cctkGH_);
assert(in_global_mode(cctkGH));
Expand All @@ -458,40 +451,49 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
std::vector<CCTK_REAL> localsx(npoints);
std::vector<CCTK_REAL> localsy(npoints);
std::vector<CCTK_REAL> localsz(npoints);

if (have_MultiPatch_GlobalToLocal2) {
MultiPatch_GlobalToLocal2(npoints, globalsx, globalsy, globalsz,
patches.data(), localsx.data(), localsy.data(),
localsz.data());
} else {
for (int n = 0; n < npoints; ++n)
for (int n = 0; n < npoints; ++n) {
patches.at(n) = 0;
for (int n = 0; n < npoints; ++n)
}

for (int n = 0; n < npoints; ++n) {
localsx.at(n) = globalsx[n];
for (int n = 0; n < npoints; ++n)
}

for (int n = 0; n < npoints; ++n) {
localsy.at(n) = globalsy[n];
for (int n = 0; n < npoints; ++n)
}

for (int n = 0; n < npoints; ++n) {
localsz.at(n) = globalsz[n];
}
}

// Apply symmetries to coordinates
std::vector<bool> symmetry_reflected_z;
assert(!reflection_x);
assert(!reflection_y);
assert(!reflection_upper_x);
assert(!reflection_upper_y);
assert(!reflection_upper_z);

if (reflection_z) {
symmetry_reflected_z.resize(npoints);
assert(ghext->num_patches() == 1);
constexpr int patch = 0;
const amrex::Geometry &geom = ghext->patchdata.at(patch).amrcore->Geom(0);
const CCTK_REAL *restrict const xmin = geom.ProbLo();
#pragma omp simd
#pragma omp parallel for simd
for (int n = 0; n < npoints; ++n) {
const bool refl = localsz[n] < xmin[2];
const auto refl = localsz[n] < xmin[2];
symmetry_reflected_z[n] = refl;
if (refl)
if (refl) {
localsz[n] = 2 * xmin[2] - localsz[n];
}
}
}

Expand All @@ -504,7 +506,8 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
std::vector<CCTK_REAL> posx(npoints);
std::vector<CCTK_REAL> posy(npoints);
std::vector<CCTK_REAL> posz(npoints);
#pragma omp simd

#pragma omp parallel for simd
for (int n = 0; n < npoints; ++n) {
const int patch = patches.at(n);
const amrex::Geometry &geom = ghext->patchdata.at(patch).amrcore->Geom(0);
Expand All @@ -520,11 +523,6 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
}

// Create particle containers
using Container = amrex::AmrParticleContainer<3, 2>;
using Particle = Container::ParticleType;

using PinnedParticleTile = typename amrex::ParticleContainer_impl<
Particle, 0, 0, amrex::PinnedArenaAllocator>::ParticleTileType;
std::vector<PinnedParticleTile> pinned_particle_tiles(ghext->num_patches());
for (int patch = 0; patch < ghext->num_patches(); ++patch) {
PinnedParticleTile &pinned_particle_tile = pinned_particle_tiles.at(patch);
Expand Down Expand Up @@ -552,14 +550,13 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
pinned_particle_tiles.at(patch).push_back(p);
}

using ParticleTile = Container::ParticleTileType;
std::vector<Container> containers(ghext->num_patches());
containers.resize(ghext->num_patches());
for (int patch = 0; patch < ghext->num_patches(); ++patch) {
const PinnedParticleTile &pinned_particle_tile =
pinned_particle_tiles.at(patch);

const auto &restrict patchdata = ghext->patchdata.at(patch);
containers.at(patch) = Container(patchdata.amrcore.get());
containers.at(patch) = std::make_shared<Container>(patchdata.amrcore.get());
const int level = 0;
const auto &restrict leveldata = patchdata.leveldata.at(level);
const amrex::MFIter mfi(*leveldata.fab);
Expand All @@ -569,7 +566,7 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
continue;
}

ParticleTile &particle_tile = containers.at(patch).GetParticles(
ParticleTile &particle_tile = containers.at(patch)->GetParticles(
level)[make_pair(mfi.index(), mfi.LocalTileIndex())];

const auto old_np = particle_tile.numParticles();
Expand All @@ -587,10 +584,10 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
std::size_t old_nparticles = 0;
std::set<int> oldids;
{
const auto &levels = container.GetParticles();
const auto &levels = container->GetParticles();
for (const auto &level : levels) {
const int lev = int(&level - levels.data());
for (amrex::ParConstIter<3, 2> pti(container, lev); pti.isValid();
for (amrex::ParConstIter<3, 2> pti(*container, lev); pti.isValid();
++pti) {
const auto &particles = pti.GetArrayOfStructs();
const int component = MFPointer(pti).index();
Expand All @@ -605,16 +602,16 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
}
#endif

container.Redistribute();
container->Redistribute();

#ifdef CCTK_DEBUG
std::size_t new_nparticles = 0;
std::set<int> newids;
{
const auto &levels = container.GetParticles();
const auto &levels = container->GetParticles();
for (const auto &level : levels) {
const int lev = int(&level - levels.data());
for (amrex::ParConstIter<3, 2> pti(container, lev); pti.isValid();
for (amrex::ParConstIter<3, 2> pti(*container, lev); pti.isValid();
++pti) {
const int component = MFPointer(pti).index();
const auto &particles = pti.GetArrayOfStructs();
Expand Down Expand Up @@ -643,6 +640,14 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
}
#endif
}
}

void InterpolateFromSetup(const InterpolationSetup &setup, const CCTK_INT nvars,
const CCTK_INT *restrict const varinds,
const CCTK_INT *restrict const operations,
const CCTK_INT allow_boundaries,
const CCTK_POINTER resultptrs_) {
DECLARE_CCTK_PARAMETERS;

// Define result variables
const int nprocs = amrex::ParallelDescriptor::NProcs();
Expand Down Expand Up @@ -670,8 +675,8 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
const int level = leveldata.level;

// TODO: use OpenMP
for (amrex::ParIter<3, 2> pti(containers.at(patch), level); pti.isValid();
++pti) {
for (amrex::ParIter<3, 2> pti(*setup.containers.at(patch), level);
pti.isValid(); ++pti) {
const MFPointer mfp(pti);
const GridDesc grid(leveldata, mfp);
// const int component = mfp.index();
Expand Down Expand Up @@ -865,23 +870,23 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
comm);
#ifdef CCTK_DEBUG
// Check consistency of received ids
std::vector<bool> idxs(npoints, false);
for (int n = 0; n < npoints; ++n) {
std::vector<bool> idxs(setup.npoints, false);
for (int n = 0; n < setup.npoints; ++n) {
const int offset = (nvars + 1) * n;
const int idx = int(recvbuf.at(offset));
assert(!idxs.at(idx));
idxs.at(idx) = true;
}
for (int n = 0; n < npoints; ++n)
for (int n = 0; n < setup.npoints; ++n)
assert(idxs.at(n));
#endif

// Set result
CCTK_REAL *const restrict *const restrict resultptrs =
static_cast<CCTK_REAL *const *>(resultptrs_);
if (int(recvbuf.size()) != (nvars + 1) * npoints)
if (int(recvbuf.size()) != (nvars + 1) * setup.npoints)
CCTK_ERROR("Internal error");
for (int n = 0; n < npoints; ++n) {
for (int n = 0; n < setup.npoints; ++n) {
const int offset = (nvars + 1) * n;
const int idx = int(recvbuf.at(offset));
for (int v = 0; v < nvars; ++v)
Expand All @@ -894,6 +899,7 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
assert(!reflection_upper_x);
assert(!reflection_upper_y);
assert(!reflection_upper_z);

if (reflection_z) {
// The code below is only valid for Psi4
assert(nvars == 2);
Expand All @@ -903,12 +909,29 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
// n^a = et^a - er^a
// m^a = etheta^a + i ephi^a
// Psi4 = C_abcd m-bar^b n^b m-bar^c n^d
for (int n = 0; n < npoints; ++n) {
if (symmetry_reflected_z[n]) {
for (int n = 0; n < setup.npoints; ++n) {
if (setup.symmetry_reflected_z[n]) {
resultptrs[0][n] = -resultptrs[0][n];
resultptrs[1][n] = +resultptrs[1][n];
}
}
}
}

extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
const CCTK_INT npoints,
const CCTK_REAL *restrict const globalsx,
const CCTK_REAL *restrict const globalsy,
const CCTK_REAL *restrict const globalsz,
const CCTK_INT nvars,
const CCTK_INT *restrict const varinds,
const CCTK_INT *restrict const operations,
const CCTK_INT allow_boundaries,
const CCTK_POINTER resultptrs_) {
const InterpolationSetup setup(cctkGH_, npoints, globalsx, globalsy,
globalsz);
InterpolateFromSetup(setup, nvars, varinds, operations, allow_boundaries,
resultptrs_);
}

} // namespace CarpetX