From 079f9bb97f1c25ff4b1788d2993afc02f3e3263e Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Sat, 29 Nov 2025 02:57:13 -0600 Subject: [PATCH 1/3] CarpetX: Split CarpetX_Interpolate into setup and interpolation routines --- CarpetX/src/interp.hxx | 39 ++++++++++ CarpetX/src/interpolate.cxx | 144 +++++++++++++++++++++--------------- 2 files changed, 122 insertions(+), 61 deletions(-) diff --git a/CarpetX/src/interp.hxx b/CarpetX/src/interp.hxx index 6eca375a9..4784bc94f 100644 --- a/CarpetX/src/interp.hxx +++ b/CarpetX/src/interp.hxx @@ -3,6 +3,13 @@ #include +#include +#include + +#include +#include +#include + // Scheduled functions extern "C" int CarpetX_InterpGridArrays( cGH const *const cGH, int const N_dims, int const local_interp_handle, @@ -23,6 +30,38 @@ 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}; + const CCTK_INT nvars{0}; + const CCTK_INT *restrict const varinds{nullptr}; + const CCTK_INT *restrict const operations{nullptr}; + const CCTK_INT allow_boundaries{0}; + + std::vector symmetry_reflected_z; + std::vector > 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, + const CCTK_INT nvars, + const CCTK_INT *restrict const varinds, + const CCTK_INT *restrict const operations, + const CCTK_INT allow_boundaries); +}; + +void InterpolateFromSetup(const InterpolationSetup &setup, + 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, diff --git a/CarpetX/src/interpolate.cxx b/CarpetX/src/interpolate.cxx index eee39a9a8..91feb3510 100644 --- a/CarpetX/src/interpolate.cxx +++ b/CarpetX/src/interpolate.cxx @@ -20,10 +20,7 @@ #include #include #include -#include #include -#include -#include #include #include @@ -433,16 +430,15 @@ 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, const CCTK_INT nvars, + const CCTK_INT *restrict const varinds, + const CCTK_INT *restrict const operations, const CCTK_INT allow_boundaries) + : npoints(npoints), nvars(nvars), varinds(varinds), operations(operations), + allow_boundaries(allow_boundaries) { DECLARE_CCTK_PARAMETERS; const cGH *restrict const cctkGH = static_cast(cctkGH_); assert(in_global_mode(cctkGH)); @@ -458,40 +454,49 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, std::vector localsx(npoints); std::vector localsy(npoints); std::vector 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 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]; + } } } @@ -504,7 +509,8 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, std::vector posx(npoints); std::vector posy(npoints); std::vector 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); @@ -520,11 +526,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 pinned_particle_tiles(ghext->num_patches()); for (int patch = 0; patch < ghext->num_patches(); ++patch) { PinnedParticleTile &pinned_particle_tile = pinned_particle_tiles.at(patch); @@ -552,14 +553,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 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(patchdata.amrcore.get()); const int level = 0; const auto &restrict leveldata = patchdata.leveldata.at(level); const amrex::MFIter mfi(*leveldata.fab); @@ -569,7 +569,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(); @@ -605,7 +605,7 @@ 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; @@ -643,6 +643,11 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, } #endif } +} + +void InterpolateFromSetup(const InterpolationSetup &setup, + const CCTK_POINTER resultptrs_) { + DECLARE_CCTK_PARAMETERS; // Define result variables const int nprocs = amrex::ParallelDescriptor::NProcs(); @@ -653,12 +658,12 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, struct givi_t { int gi, vi; }; - std::vector givis(nvars); - for (int v = 0; v < nvars; ++v) { - int gi = CCTK_GroupIndexFromVarI(varinds[v]); + std::vector givis(setup.nvars); + for (int v = 0; v < setup.nvars; ++v) { + int gi = CCTK_GroupIndexFromVarI(setup.varinds[v]); assert(gi >= 0); assert(gi < CCTK_NumGroups()); - int vi = varinds[v] - CCTK_FirstVarIndexI(gi); + int vi = setup.varinds[v] - CCTK_FirstVarIndexI(gi); assert(vi >= 0); assert(vi < CCTK_NumVarsInGroupI(gi)); givis.at(v) = {gi, vi}; @@ -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(); @@ -679,11 +684,11 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const int np = pti.numParticles(); const auto &particles = pti.GetArrayOfStructs(); - std::vector > varresults(nvars); + std::vector > varresults(setup.nvars); // TODO: Don't re-calculate interpolation coefficients for each // variable - for (int v = 0; v < nvars; ++v) { + for (int v = 0; v < setup.nvars; ++v) { const int gi = givis.at(v).gi; const int vi = givis.at(v).vi; const auto &restrict groupdata = *leveldata.groupdata.at(gi); @@ -694,7 +699,7 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const amrex::Array4 &vars = groupdata.mfab.at(tl)->array(pti); vect derivs; - int op = operations[v]; + int op = setup.operations[v]; while (op > 0) { const int dir = op % 10 - 1; if (dir >= 0) { @@ -714,35 +719,35 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, case 0: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 1: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 2: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 3: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 4: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } @@ -762,35 +767,35 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, case 0: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 1: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 2: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 3: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 4: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(allow_boundaries)}; + level, vars, derivs, bool(setup.allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } @@ -816,7 +821,7 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const int id = particles[n].idata(1); auto &result = results.at(proc); result.push_back(id); - for (int v = 0; v < nvars; ++v) + for (int v = 0; v < setup.nvars; ++v) result.push_back(varresults.at(v).at(n)); } } @@ -879,12 +884,12 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, // Set result CCTK_REAL *const restrict *const restrict resultptrs = static_cast(resultptrs_); - if (int(recvbuf.size()) != (nvars + 1) * npoints) + if (int(recvbuf.size()) != (setup.nvars + 1) * setup.npoints) CCTK_ERROR("Internal error"); - for (int n = 0; n < npoints; ++n) { - const int offset = (nvars + 1) * n; + for (int n = 0; n < setup.npoints; ++n) { + const int offset = (setup.nvars + 1) * n; const int idx = int(recvbuf.at(offset)); - for (int v = 0; v < nvars; ++v) + for (int v = 0; v < setup.nvars; ++v) resultptrs[v][idx] = recvbuf.at(offset + 1 + v); } @@ -894,21 +899,38 @@ 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); - assert(varinds[0] == CCTK_VarIndex("Weyl::Psi4re")); - assert(varinds[1] == CCTK_VarIndex("Weyl::Psi4im")); + assert(setup.nvars == 2); + assert(setup.varinds[0] == CCTK_VarIndex("Weyl::Psi4re")); + assert(setup.varinds[1] == CCTK_VarIndex("Weyl::Psi4im")); // l^a = et^a + er^a // 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, + nvars, varinds, operations, allow_boundaries); + InterpolateFromSetup(setup, resultptrs_); +} + } // namespace CarpetX From 5fd11a0e28b7a5ecdea53f8f81f2bf81f874b799 Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Sat, 29 Nov 2025 03:10:10 -0600 Subject: [PATCH 2/3] CarpetX: Fix debug mode compilation issues in interpolator changes --- CarpetX/src/interpolate.cxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/CarpetX/src/interpolate.cxx b/CarpetX/src/interpolate.cxx index 91feb3510..5619a4ea6 100644 --- a/CarpetX/src/interpolate.cxx +++ b/CarpetX/src/interpolate.cxx @@ -587,10 +587,10 @@ CarpetX::InterpolationSetup::InterpolationSetup( std::size_t old_nparticles = 0; std::set 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(); @@ -611,10 +611,10 @@ CarpetX::InterpolationSetup::InterpolationSetup( std::size_t new_nparticles = 0; std::set 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(); @@ -870,14 +870,14 @@ void InterpolateFromSetup(const InterpolationSetup &setup, comm); #ifdef CCTK_DEBUG // Check consistency of received ids - std::vector idxs(npoints, false); - for (int n = 0; n < npoints; ++n) { - const int offset = (nvars + 1) * n; + std::vector idxs(setup.npoints, false); + for (int n = 0; n < setup.npoints; ++n) { + const int offset = (setup.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 From 45a8aeb2dd46574c1266b0c7dbcce06cbd93f5c6 Mon Sep 17 00:00:00 2001 From: Lucas Timotheo Sanches Date: Sat, 29 Nov 2025 03:42:49 -0600 Subject: [PATCH 3/3] CarpetX: Don't save Cactus data that should not be saved in interpolation setup --- CarpetX/src/interp.hxx | 16 +++------ CarpetX/src/interpolate.cxx | 69 +++++++++++++++++++------------------ 2 files changed, 40 insertions(+), 45 deletions(-) diff --git a/CarpetX/src/interp.hxx b/CarpetX/src/interp.hxx index 4784bc94f..b710ddd3c 100644 --- a/CarpetX/src/interp.hxx +++ b/CarpetX/src/interp.hxx @@ -40,11 +40,6 @@ struct InterpolationSetup { using ParticleTile = Container::ParticleTileType; const CCTK_INT npoints{0}; - const CCTK_INT nvars{0}; - const CCTK_INT *restrict const varinds{nullptr}; - const CCTK_INT *restrict const operations{nullptr}; - const CCTK_INT allow_boundaries{0}; - std::vector symmetry_reflected_z; std::vector > containers{}; @@ -52,14 +47,13 @@ struct InterpolationSetup { 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_REAL *restrict const globalsz); }; -void InterpolateFromSetup(const InterpolationSetup &setup, +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 diff --git a/CarpetX/src/interpolate.cxx b/CarpetX/src/interpolate.cxx index 5619a4ea6..3a1d62124 100644 --- a/CarpetX/src/interpolate.cxx +++ b/CarpetX/src/interpolate.cxx @@ -434,11 +434,8 @@ 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, const CCTK_INT nvars, - const CCTK_INT *restrict const varinds, - const CCTK_INT *restrict const operations, const CCTK_INT allow_boundaries) - : npoints(npoints), nvars(nvars), varinds(varinds), operations(operations), - allow_boundaries(allow_boundaries) { + const CCTK_REAL *restrict const globalsz) + : npoints(npoints) { DECLARE_CCTK_PARAMETERS; const cGH *restrict const cctkGH = static_cast(cctkGH_); assert(in_global_mode(cctkGH)); @@ -645,7 +642,10 @@ CarpetX::InterpolationSetup::InterpolationSetup( } } -void InterpolateFromSetup(const InterpolationSetup &setup, +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; @@ -658,12 +658,12 @@ void InterpolateFromSetup(const InterpolationSetup &setup, struct givi_t { int gi, vi; }; - std::vector givis(setup.nvars); - for (int v = 0; v < setup.nvars; ++v) { - int gi = CCTK_GroupIndexFromVarI(setup.varinds[v]); + std::vector givis(nvars); + for (int v = 0; v < nvars; ++v) { + int gi = CCTK_GroupIndexFromVarI(varinds[v]); assert(gi >= 0); assert(gi < CCTK_NumGroups()); - int vi = setup.varinds[v] - CCTK_FirstVarIndexI(gi); + int vi = varinds[v] - CCTK_FirstVarIndexI(gi); assert(vi >= 0); assert(vi < CCTK_NumVarsInGroupI(gi)); givis.at(v) = {gi, vi}; @@ -684,11 +684,11 @@ void InterpolateFromSetup(const InterpolationSetup &setup, const int np = pti.numParticles(); const auto &particles = pti.GetArrayOfStructs(); - std::vector > varresults(setup.nvars); + std::vector > varresults(nvars); // TODO: Don't re-calculate interpolation coefficients for each // variable - for (int v = 0; v < setup.nvars; ++v) { + for (int v = 0; v < nvars; ++v) { const int gi = givis.at(v).gi; const int vi = givis.at(v).vi; const auto &restrict groupdata = *leveldata.groupdata.at(gi); @@ -699,7 +699,7 @@ void InterpolateFromSetup(const InterpolationSetup &setup, const amrex::Array4 &vars = groupdata.mfab.at(tl)->array(pti); vect derivs; - int op = setup.operations[v]; + int op = operations[v]; while (op > 0) { const int dir = op % 10 - 1; if (dir >= 0) { @@ -719,35 +719,35 @@ void InterpolateFromSetup(const InterpolationSetup &setup, case 0: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 1: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 2: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 3: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 4: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } @@ -767,35 +767,35 @@ void InterpolateFromSetup(const InterpolationSetup &setup, case 0: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 1: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 2: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 3: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } case 4: { const interpolator interp{ grid, gi, vi, patch, - level, vars, derivs, bool(setup.allow_boundaries)}; + level, vars, derivs, bool(allow_boundaries)}; interp.interpolate3d(particles, varresult); break; } @@ -821,7 +821,7 @@ void InterpolateFromSetup(const InterpolationSetup &setup, const int id = particles[n].idata(1); auto &result = results.at(proc); result.push_back(id); - for (int v = 0; v < setup.nvars; ++v) + for (int v = 0; v < nvars; ++v) result.push_back(varresults.at(v).at(n)); } } @@ -872,7 +872,7 @@ void InterpolateFromSetup(const InterpolationSetup &setup, // Check consistency of received ids std::vector idxs(setup.npoints, false); for (int n = 0; n < setup.npoints; ++n) { - const int offset = (setup.nvars + 1) * n; + const int offset = (nvars + 1) * n; const int idx = int(recvbuf.at(offset)); assert(!idxs.at(idx)); idxs.at(idx) = true; @@ -884,12 +884,12 @@ void InterpolateFromSetup(const InterpolationSetup &setup, // Set result CCTK_REAL *const restrict *const restrict resultptrs = static_cast(resultptrs_); - if (int(recvbuf.size()) != (setup.nvars + 1) * setup.npoints) + if (int(recvbuf.size()) != (nvars + 1) * setup.npoints) CCTK_ERROR("Internal error"); for (int n = 0; n < setup.npoints; ++n) { - const int offset = (setup.nvars + 1) * n; + const int offset = (nvars + 1) * n; const int idx = int(recvbuf.at(offset)); - for (int v = 0; v < setup.nvars; ++v) + for (int v = 0; v < nvars; ++v) resultptrs[v][idx] = recvbuf.at(offset + 1 + v); } @@ -902,9 +902,9 @@ void InterpolateFromSetup(const InterpolationSetup &setup, if (reflection_z) { // The code below is only valid for Psi4 - assert(setup.nvars == 2); - assert(setup.varinds[0] == CCTK_VarIndex("Weyl::Psi4re")); - assert(setup.varinds[1] == CCTK_VarIndex("Weyl::Psi4im")); + assert(nvars == 2); + assert(varinds[0] == CCTK_VarIndex("Weyl::Psi4re")); + assert(varinds[1] == CCTK_VarIndex("Weyl::Psi4im")); // l^a = et^a + er^a // n^a = et^a - er^a // m^a = etheta^a + i ephi^a @@ -928,9 +928,10 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_INT *restrict const operations, const CCTK_INT allow_boundaries, const CCTK_POINTER resultptrs_) { - const InterpolationSetup setup(cctkGH_, npoints, globalsx, globalsy, globalsz, - nvars, varinds, operations, allow_boundaries); - InterpolateFromSetup(setup, resultptrs_); + const InterpolationSetup setup(cctkGH_, npoints, globalsx, globalsy, + globalsz); + InterpolateFromSetup(setup, nvars, varinds, operations, allow_boundaries, + resultptrs_); } } // namespace CarpetX