diff --git a/CarpetX/src/interp.hxx b/CarpetX/src/interp.hxx index 6eca375a9..b710ddd3c 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,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 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); +}; + +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, diff --git a/CarpetX/src/interpolate.cxx b/CarpetX/src/interpolate.cxx index eee39a9a8..3a1d62124 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,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(cctkGH_); assert(in_global_mode(cctkGH)); @@ -458,40 +451,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 +506,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 +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 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 +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 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 +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(); @@ -587,10 +584,10 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, 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(); @@ -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 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(); @@ -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(); @@ -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(); @@ -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 idxs(npoints, false); - for (int n = 0; n < npoints; ++n) { + std::vector 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(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) @@ -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); @@ -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