diff --git a/Src/Particle/AMReX_ParticleContainer.H b/Src/Particle/AMReX_ParticleContainer.H index df71f2eeec1..966ca0a8441 100644 --- a/Src/Particle/AMReX_ParticleContainer.H +++ b/Src/Particle/AMReX_ParticleContainer.H @@ -1180,11 +1180,6 @@ public: Vector >& mf_to_be_filled, int lev_min, int ncomp, int finest_level, int ngrow=2) const; - void Interpolate (Vector >& mesh_data, - int lev_min, int lev_max); - - void InterpolateSingleLevel (MultiFab& mesh_data, int lev); - void AssignCellDensitySingleLevel (int rho_index, MultiFab& mf, int level, int ncomp=1, int particle_lvl_offset = 0) const; diff --git a/Src/Particle/AMReX_ParticleContainerI.H b/Src/Particle/AMReX_ParticleContainerI.H index 9dd1d39b953..fe83ce5ba1c 100644 --- a/Src/Particle/AMReX_ParticleContainerI.H +++ b/Src/Particle/AMReX_ParticleContainerI.H @@ -2481,56 +2481,6 @@ AssignCellDensitySingleLevel (int rho_index, } } -template class Allocator, class CellAssignor> -void -ParticleContainer_impl::Interpolate (Vector >& mesh_data, - int lev_min, int lev_max) -{ - BL_PROFILE("ParticleContainer::Interpolate()"); - for (int lev = lev_min; lev <= lev_max; ++lev) { - InterpolateSingleLevel(*mesh_data[lev], lev); - } -} - -template class Allocator, class CellAssignor> -void -ParticleContainer_impl:: -InterpolateSingleLevel (MultiFab& mesh_data, int lev) -{ - BL_PROFILE("ParticleContainer::InterpolateSingleLevel()"); - - if (mesh_data.nGrow() < 1) { - amrex::Error("Must have at least one ghost cell when in InterpolateSingleLevel"); - } - - const Geometry& gm = Geom(lev); - const auto plo = gm.ProbLoArray(); - const auto dxi = gm.InvCellSizeArray(); - - using ParIter = ParIter_impl; - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIter pti(*this, lev); pti.isValid(); ++pti) - { - auto& particles = pti.GetArrayOfStructs(); - auto ptd = pti.GetParticleTile().getParticleTileData(); - FArrayBox& fab = mesh_data[pti]; - const auto fabarr = fab.array(); - const Long np = particles.numParticles(); - - int nComp = fab.nComp(); - AMREX_FOR_1D( np, i, - { - auto p = make_particle{}(ptd,i); - amrex_interpolate_cic(p, nComp, fabarr, plo, dxi); - }); - } -} - template class Allocator, class CellAssignor> void diff --git a/Src/Particle/AMReX_Particle_mod_K.H b/Src/Particle/AMReX_Particle_mod_K.H index 2d1f91b4275..ac9254bd62c 100644 --- a/Src/Particle/AMReX_Particle_mod_K.H +++ b/Src/Particle/AMReX_Particle_mod_K.H @@ -233,80 +233,6 @@ void amrex_deposit_particle_dx_cic (P const& p, int nc, amrex::Array4 -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void amrex_interpolate_cic (P const& p, int nc, amrex::Array4 const& acc, - amrex::GpuArray const& plo, - amrex::GpuArray const& dxi) -{ -#if (AMREX_SPACEDIM == 1) - amrex::Real lx = (p.pos(0) - plo[0]) * dxi[0] + Real(0.5); - - int i = static_cast(amrex::Math::floor(lx)); - - amrex::Real xint = lx - static_cast(i); - - amrex::Real sx[] = {Real(1.0)-xint, xint}; - - for (int comp=0; comp < nc; ++comp) { - for (int ii = 0; ii <= 1; ++ii) { - amrex::Real acceleration = sx[ii]*acc(i+ii-1,0,0,comp); - amrex::ignore_unused(acceleration); - } - } -#elif (AMREX_SPACEDIM == 2) - amrex::Real lx = (p.pos(0) - plo[0]) * dxi[0] + Real(0.5); - amrex::Real ly = (p.pos(1) - plo[1]) * dxi[1] + Real(0.5); - - int i = static_cast(amrex::Math::floor(lx)); - int j = static_cast(amrex::Math::floor(ly)); - - amrex::Real xint = lx - static_cast(i); - amrex::Real yint = ly - static_cast(j); - - amrex::Real sx[] = {Real(1.0)-xint, xint}; - amrex::Real sy[] = {Real(1.0)-yint, yint}; - - for (int comp=0; comp < nc; ++comp) { - for (int jj = 0; jj <= 1; ++jj) { - for (int ii = 0; ii <= 1; ++ii) { - amrex::Real acceleration = sx[ii]*sy[jj]*acc(i+ii-1,j+jj-1,0,comp); - amrex::ignore_unused(acceleration); - } - } - } -#elif (AMREX_SPACEDIM == 3) - amrex::Real lx = (p.pos(0) - plo[0]) * dxi[0] + Real(0.5); - amrex::Real ly = (p.pos(1) - plo[1]) * dxi[1] + Real(0.5); - amrex::Real lz = (p.pos(2) - plo[2]) * dxi[2] + Real(0.5); - - int i = static_cast(amrex::Math::floor(lx)); - int j = static_cast(amrex::Math::floor(ly)); - int k = static_cast(amrex::Math::floor(lz)); - - amrex::Real xint = lx - static_cast(i); - amrex::Real yint = ly - static_cast(j); - amrex::Real zint = lz - static_cast(k); - - amrex::Real sx[] = {Real(1.0)-xint, xint}; - amrex::Real sy[] = {Real(1.0)-yint, yint}; - amrex::Real sz[] = {Real(1.0)-zint, zint}; - - for (int comp=0; comp < nc; ++comp) { - for (int kk = 0; kk <= 1; ++kk) { - for (int jj = 0; jj <= 1; ++jj) { - for (int ii = 0; ii <= 1; ++ii) { - amrex::Real acceleration = sx[ii]*sy[jj]*sz[kk]*acc(i+ii-1,j+jj-1,k+kk-1,comp); - amrex::ignore_unused(acceleration); - } - } - } - } -#else - amrex::Abort("Not implemented."); -#endif -} - } #endif diff --git a/Src/Particle/AMReX_TracerParticle_mod_K.H b/Src/Particle/AMReX_TracerParticle_mod_K.H index 433ad864e58..9df0864d33c 100644 --- a/Src/Particle/AMReX_TracerParticle_mod_K.H +++ b/Src/Particle/AMReX_TracerParticle_mod_K.H @@ -1,7 +1,7 @@ #ifndef AMREX_TRACERPARTICLE_MOD_K_H #define AMREX_TRACERPARTICLE_MOD_K_H -#include +#include #include #include #include @@ -9,200 +9,131 @@ #include #include #include - #include -namespace amrex{ +namespace amrex { +/** + \brief Linearly interpolates the mesh data to the particle position from cell-centered data. +*/ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cic_interpolate (const P& p, amrex::GpuArray const& plo, amrex::GpuArray const& dxi, - const amrex::Array4 & uccarr, - amrex::ParticleReal * val, int M) + const amrex::Array4& data_arr, + amrex::ParticleReal * val, int M = AMREX_SPACEDIM) { - AMREX_ASSERT(val != nullptr); - -#if (AMREX_SPACEDIM == 1) - - amrex::Real lx = (Real(p.pos(0)) - plo[0]) * dxi[0] - Real(0.5); //len - - int const i = static_cast(amrex::Math::floor(lx)); //cell - - amrex::Real xint = lx - static_cast(i); //frac - - amrex::Real sx[] = {Real(1.0) - xint, xint}; - - for (int d=0; d < M; ++d) - { - val[d] = ParticleReal(0.0); - for (int ii = 0; ii<=1; ++ii) - { - val[d] += static_cast(sx[ii]*uccarr(i+ii,0,0,d)); - } - } - - -#elif (AMREX_SPACEDIM == 2) - - amrex::Real lx = (Real(p.pos(0)) - plo[0]) * dxi[0] - Real(0.5); - amrex::Real ly = (Real(p.pos(1)) - plo[1]) * dxi[1] - Real(0.5); - - int const i = static_cast(amrex::Math::floor(lx)); - int const j = static_cast(amrex::Math::floor(ly)); - - amrex::Real xint = lx - static_cast(i); - amrex::Real yint = ly - static_cast(j); - - amrex::Real sx[] = {Real(1.0) - xint, xint}; - amrex::Real sy[] = {Real(1.0) - yint, yint}; - - for (int d=0; d < M; ++d) - { - val[d] = ParticleReal(0.0); - for (int jj = 0; jj <= 1; ++jj) - { - for (int ii = 0; ii <= 1; ++ii) - { - val[d] += static_cast(sx[ii]*sy[jj]*uccarr(i+ii,j+jj,0,d)); - } - } - } - -#elif (AMREX_SPACEDIM == 3) - - amrex::Real lx = (Real(p.pos(0)) - plo[0]) * dxi[0] - Real(0.5); - amrex::Real ly = (Real(p.pos(1)) - plo[1]) * dxi[1] - Real(0.5); - amrex::Real lz = (Real(p.pos(2)) - plo[2]) * dxi[2] - Real(0.5); - - int const i = static_cast(amrex::Math::floor(lx)); - int const j = static_cast(amrex::Math::floor(ly)); - int const k = static_cast(amrex::Math::floor(lz)); - - amrex::Real const xint = lx - static_cast(i); - amrex::Real const yint = ly - static_cast(j); - amrex::Real const zint = lz - static_cast(k); - - amrex::Real sx[] = {Real(1.0) - xint, xint}; - amrex::Real sy[] = {Real(1.0) - yint, yint}; - amrex::Real sz[] = {Real(1.0) - zint, zint}; - - for (int d=0; d < M; ++d) - { - val[d] = ParticleReal(0.0); - for (int kk = 0; kk<=1; ++kk) - { - for (int jj = 0; jj <= 1; ++jj) - { - for (int ii = 0; ii <= 1; ++ii) - { - val[d] += static_cast(sx[ii]*sy[jj]*sz[kk]*uccarr(i+ii,j+jj,k+kk,d)); - } - } - } - } -#endif + int start_comp = 0; + int ncomp_per_array = M; + int num_arrays = 1; + IntVect is_nodal = amrex::IntVect::TheZeroVector(); + linear_interpolate_to_particle (p, plo, dxi, &data_arr, val, &is_nodal, start_comp, ncomp_per_array, num_arrays); } +/** + \brief Linearly interpolates the mesh data to the particle position from node-centered data. +*/ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void cic_interpolate (const P& p, - amrex::GpuArray const& plo, - amrex::GpuArray const& dxi, - const amrex::Array4 & uccarr, - amrex::ParticleReal * val) +void cic_interpolate_nd (const P& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + const amrex::Array4& data_arr, + amrex::ParticleReal * val, int M = AMREX_SPACEDIM) { - cic_interpolate(p, plo, dxi, uccarr, val, AMREX_SPACEDIM); + int start_comp = 0; + int ncomp_per_array = M; + int num_arrays = 1; + IntVect is_nodal = amrex::IntVect::TheUnitVector(); + linear_interpolate_to_particle (p, plo, dxi, &data_arr, val, &is_nodal, start_comp, ncomp_per_array, num_arrays); } +/** + \brief Linearly interpolates the mesh data to the particle position from face-centered data. + The nth component of the data_arr array is nodal in the nth direction, and cell-centered in the others. +*/ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mac_interpolate (const P& p, amrex::GpuArray const& plo, amrex::GpuArray const& dxi, - amrex::GpuArray,AMREX_SPACEDIM> const& p_uccarr, + amrex::GpuArray,AMREX_SPACEDIM> const& data_arr, amrex::ParticleReal * val) { + int start_comp = 0; + int ncomp_per_array = 1; + int num_arrays = AMREX_SPACEDIM; + IntVect is_nodal[AMREX_SPACEDIM]; + for (int d=0; d < AMREX_SPACEDIM; ++d) { + is_nodal[d] = amrex::IntVect::TheZeroVector(); + is_nodal[d][d] = 1; + } + linear_interpolate_to_particle (p, plo, dxi, data_arr.data(), val, &is_nodal[0], start_comp, ncomp_per_array, num_arrays); +} -#if (AMREX_SPACEDIM == 1) - for (int d=0; d < AMREX_SPACEDIM; ++d) - { - amrex::Real lx = (Real(p.pos(0))-plo[0])*dxi[0] - static_cast(d != 0)*Real(0.5); - int const i = static_cast(amrex::Math::floor(lx)); - amrex::Real const xint = lx - static_cast(i); +/** + \brief Linearly interpolates the mesh data to the particle position from mesh data. + This general form can handle an arbitrary number of Array4s, each with different staggerings. +*/ +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void linear_interpolate_to_particle (const P& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + const Array4* data_arr, + amrex::ParticleReal * val, + const IntVect* is_nodal, + int start_comp, int ncomp, int num_arrays) +{ + AMREX_ASSERT(val != nullptr); - amrex::Real sx[] = {Real(1.0) - xint, xint}; + int ctr = 0; - val[d] = ParticleReal(0.0); - for (int ii = 0; ii <= 1; ++ii) - { - val[d] += static_cast((p_uccarr[d])(i+ii, 0, 0, 0)*sx[ii]); - } - } + for (int d = 0; d < num_arrays; d++) + { + AMREX_D_TERM(amrex::Real lx = (Real(p.pos(0))-plo[0])*dxi[0] - static_cast(!is_nodal[d][0])*Real(0.5);, + amrex::Real ly = (Real(p.pos(1))-plo[1])*dxi[1] - static_cast(!is_nodal[d][1])*Real(0.5);, + amrex::Real lz = (Real(p.pos(2))-plo[2])*dxi[2] - static_cast(!is_nodal[d][2])*Real(0.5)); + + // (i0,j0,k0) is the lower corner of the box needed for interpolation + // i0 = (i-1) if particle is lower than center of cell i + // i0 = (i ) if particle is higher than center of cell i + AMREX_D_TERM(int const i0 = static_cast(amrex::Math::floor(lx));, + int const j0 = static_cast(amrex::Math::floor(ly));, + int const k0 = static_cast(amrex::Math::floor(lz))); + + AMREX_D_TERM(amrex::Real const xint = lx - static_cast(i0);, + amrex::Real const yint = ly - static_cast(j0);, + amrex::Real const zint = lz - static_cast(k0)); + + amrex::Real sx[] = {amrex::Real(1.0) - xint, xint}; +#if (AMREX_SPACEDIM > 1) + amrex::Real sy[] = {amrex::Real(1.0) - yint, yint}; +#endif +#if (AMREX_SPACEDIM > 2) + amrex::Real sz[] = {amrex::Real(1.0) - zint, zint}; +#endif -#elif (AMREX_SPACEDIM == 2) - - for (int d=0; d < AMREX_SPACEDIM; ++d) - { - amrex::Real lx = (Real(p.pos(0))-plo[0])*dxi[0] - static_cast(d != 0)*Real(0.5); - amrex::Real ly = (Real(p.pos(1))-plo[1])*dxi[1] - static_cast(d != 1)*Real(0.5); - - int const i = static_cast(amrex::Math::floor(lx)); - int const j = static_cast(amrex::Math::floor(ly)); - - amrex::Real const xint = lx - static_cast(i); - amrex::Real const yint = ly - static_cast(j); - - amrex::Real sx[] = {Real(1.0) - xint, xint}; - amrex::Real sy[] = {Real(1.0) - yint, yint}; - - val[d] = ParticleReal(0.0); - for (int jj = 0; jj <= 1; ++jj) - { - for (int ii = 0; ii <= 1; ++ii) - { - val[d] += static_cast((p_uccarr[d])(i+ii, j+jj, 0, 0)*sx[ii]*sy[jj]); - } - } - } - - -#elif (AMREX_SPACEDIM == 3) - - for (int d=0; d < AMREX_SPACEDIM; ++d) - { - amrex::Real lx = (Real(p.pos(0))-plo[0])*dxi[0] - static_cast(d != 0)*Real(0.5); - amrex::Real ly = (Real(p.pos(1))-plo[1])*dxi[1] - static_cast(d != 1)*Real(0.5); - amrex::Real lz = (Real(p.pos(2))-plo[2])*dxi[2] - static_cast(d != 2)*Real(0.5); - - int const i = static_cast(amrex::Math::floor(lx)); - int const j = static_cast(amrex::Math::floor(ly)); - int const k = static_cast(amrex::Math::floor(lz)); - - amrex::Real const xint = lx - static_cast(i); - amrex::Real const yint = ly - static_cast(j); - amrex::Real const zint = lz - static_cast(k); - - amrex::Real sx[] = {Real(1.0) - xint, xint}; - amrex::Real sy[] = {Real(1.0) - yint, yint}; - amrex::Real sz[] = {Real(1.0) - zint, zint}; - - val[d] = ParticleReal(0.0); - for (int kk = 0; kk <=1; ++kk) - { - for (int jj = 0; jj <= 1; ++jj) - { - for (int ii = 0; ii <= 1; ++ii) - { - val[d] += static_cast((p_uccarr[d])(i+ii, j+jj, k+kk ,0)*sx[ii]*sy[jj]*sz[kk]); - } - } - } - } + for (int comp = start_comp; comp < ncomp; ++comp) { + val[ctr] = ParticleReal(0.0); +#if (AMREX_SPACEDIM > 2) + for (int kk = 0; kk <=1; ++kk) { #endif -} -} + +#if (AMREX_SPACEDIM > 1) + for (int jj = 0; jj <= 1; ++jj) { #endif + for (int ii = 0; ii <= 1; ++ii) { + val[ctr] += static_cast((data_arr[d])(IntVect(AMREX_D_DECL(i0+ii, j0+jj, k0+kk)), comp) * + AMREX_D_TERM(sx[ii],*sy[jj],*sz[kk])); + AMREX_D_TERM(},},}); + ctr++; + } // ncomp + } // d +} + +} // namespace amrex +#endif // include guard diff --git a/Tests/Particles/AssignDensity/main.cpp b/Tests/Particles/AssignDensity/main.cpp index fb4a0fb970b..d8caa32870e 100644 --- a/Tests/Particles/AssignDensity/main.cpp +++ b/Tests/Particles/AssignDensity/main.cpp @@ -68,10 +68,6 @@ void test_assign_density(TestParams& parms) myPC.InitRandom(num_particles, iseed, pdata, serialize); myPC.AssignCellDensitySingleLevel(0, partMF, 0, 1 + AMREX_SPACEDIM, 0); - // myPC.AssignDensitySingleLevel(0, partMF, 0, 4, 0); - - // myPC.InterpolateSingleLevel(acceleration, 0); - MultiFab::Copy(density, partMF, 0, 0, 1, 0); WriteSingleLevelPlotfile("plt00000", partMF, diff --git a/Tests/Particles/AssignMultiLevelDensity/main.cpp b/Tests/Particles/AssignMultiLevelDensity/main.cpp index 0ffa37dd695..21ae50725d2 100644 --- a/Tests/Particles/AssignMultiLevelDensity/main.cpp +++ b/Tests/Particles/AssignMultiLevelDensity/main.cpp @@ -98,14 +98,10 @@ void test_assign_density(TestParams& parms) double mass = 10.0; MyParticleContainer::ParticleInitData pdata = {{mass},{},{},{}}; - // myPC.InitRandom(num_particles, iseed, pdata, serialize, fine_box); myPC.InitRandom(num_particles, iseed, pdata, serialize); - //myPC.AssignDensity(0, true, partMF, 0, 1, 1); myPC.AssignDensity(0, partMF, 0, 1, nlevs-1); - myPC.Interpolate(acceleration, 0, nlevs-1); - for (int lev = 0; lev < nlevs; ++lev) { MultiFab::Copy(*density[lev], *partMF[lev], 0, 0, 1, 0); }