AMReX-Codes
MeshToParticle Limited to: Nearest node-centered | Linear cell-centered
Discussed in https://github.com/AMReX-Codes/amrex/discussions/4810
Originally posted by iljah November 20, 2025 This code
#define AMREX_SPACEDIM 3
#include "AMReX.H"
#include "AMReX_ParticleInterpolators.H"
#include "AMReX_Particles.H"
enum Particle_Real_Extras {Vx_i, Vy_i, Vz_i};
static constexpr int nr_extra_float_vars = Vz_i + 1;
enum Particle_Int_Extras {};
static constexpr int nr_extra_int_vars = 0;
using Particle = amrex::Particle<
nr_extra_float_vars, nr_extra_int_vars>;
using ParticleContainer = amrex::ParticleContainer<
nr_extra_float_vars, nr_extra_int_vars, 0, 0>;
int main(int argc, char* argv[]) {
using amrex::IntVect;
amrex::Initialize(argc, argv); {
amrex::Box domain{IntVect{0, 0, 0}, IntVect{2, 2, 1}};
amrex::BoxArray ba;
ba.define(domain);
ba.maxSize(2147483647);
amrex::Geometry geom; geom.define(
domain, amrex::RealBox({-1, -1, -1}, {2, 2, 1}),
amrex::CoordSys::cartesian, amrex::Array<int, 3>{1, 1, 1});
const auto grid_start = geom.ProbLoArray();
const auto dr = geom.CellSizeArray();
const auto inv_dr = geom.InvCellSizeArray();
amrex::DistributionMapping dm{ba};
amrex::MultiFab V{ba, dm, 3, 0};
for (amrex::MFIter iter(V); iter.isValid(); ++iter) {
const auto& box = iter.validbox();
const auto& V_ = V.array(iter);
amrex::ParallelFor(box, [=](int i, int j, int k){
const auto
x = grid_start[0] + dr[0]*(i + 0.5),
y = grid_start[1] + dr[1]*(j + 0.5),
z = grid_start[2] + dr[2]*(k + 0.5);
V_(i,j,k,2) = 0;
if (x < 0 and y < 0) {
V_(i,j,k,0) = -1;
V_(i,j,k,1) = +1;
}
if (x > 0 and y < 0) {
V_(i,j,k,0) = -1;
V_(i,j,k,1) = -1;
}
if (x < 0 and y > 0) {
V_(i,j,k,0) = +1;
V_(i,j,k,1) = +1;
}
if (x > 0 and y > 0) {
V_(i,j,k,0) = +1;
V_(i,j,k,1) = -1;
}
if (x < 1 and y < 1 and z < 0) {
std::cout
<< "Initialized velocity at ("
<<x<<","<<y<</*","<<z<<*/"): ("
<<V_(i,j,k,0)<<","<<V_(i,j,k,1)<<")"<<std::endl;
}
});
}
ParticleContainer particle_container(geom, dm, ba);
for (auto iter = particle_container.MakeMFIter(0); iter.isValid(); ++iter) {
auto& particles = particle_container.GetParticles(0)[
std::make_pair(iter.index(), iter.LocalTileIndex())];
for (int x = -1; x <= 1; x += 2)
for (int y = -1; y <= 1; y += 2) {
Particle particle;
particle.id() = x + 1;
particle.cpu() = y + 1;
particle.pos(0) = x / 2.0;
particle.pos(1) = y / 2.0;
particle.pos(2) = -0.5;
for (int i = 0; i < nr_extra_float_vars; i++) {
particle.rdata(i) = 0;
}
std::cout
<< "Created particle at ("
<< particle.pos(0) << ","
<< particle.pos(1) << /*","
<< particle.pos(2) << */")" << std::endl;
particles.push_back(particle);
}
}
particle_container.Redistribute();
amrex::MeshToParticle(particle_container, V, 0, [=](auto& particle1, const auto& data1){
amrex::ParticleInterpolator::Nearest interp(particle1, grid_start, inv_dr);
interp.MeshToParticle(particle1, data1, 0, Vx_i, 3,
[=](const auto& data2, int i, int j, int k, int comp){
return data2(i, j, k, comp);
},
[=](auto& particle2, int comp, auto value){
particle2.rdata(comp) += value;
});
});
for (ParticleContainer::ParIterType iter(particle_container, 0); iter.isValid(); ++iter) {
const auto& particles = iter.GetArrayOfStructs();
for (const auto& particle: particles) {
std::cout
<< "Interpolated velocity at ("
<< particle.pos(0) << ","
<< particle.pos(1) << /*","
<< particle.pos(2) << */"): ("
<< particle.rdata(Vx_i) << ","
<< particle.rdata(Vy_i) << ")" << std::endl;
}
}
} amrex::Finalize();
}
prints
Initializing AMReX (25.11)...
AMReX (25.11) initialized
Initialized velocity at (-0.5,-0.5): (-1,1)
Initialized velocity at (0.5,-0.5): (-1,-1)
Initialized velocity at (-0.5,0.5): (1,1)
Initialized velocity at (0.5,0.5): (1,-1)
Created particle at (-0.5,-0.5)
Created particle at (-0.5,0.5)
Created particle at (0.5,-0.5)
Created particle at (0.5,0.5)
Interpolated velocity at (-0.5,-0.5): (1,-1)
Interpolated velocity at (-0.5,0.5): (1,-1)
Interpolated velocity at (0.5,-0.5): (1,-1)
Interpolated velocity at (0.5,0.5): (1,-1)
AMReX (25.11) finalized
even though as far as I can tell it should print
Initializing AMReX (25.11)...
...
Interpolated velocity at (-0.5,-0.5): (-1,1)
Interpolated velocity at (-0.5,0.5): (1,1)
Interpolated velocity at (0.5,-0.5): (-1,-1)
Interpolated velocity at (0.5,0.5): (1,-1)
AMReX (25.11) finalized
How do I get MeshToParticle to interpolate as intended? Am I calculating x,y,z incorrectly around line 45?
https://github.com/AMReX-Codes/amrex/discussions/4810#discussioncomment-15026700
I think @AlexanderSinn comment is the reason we do not use this interpolator in WarpX, not all fields are nodal, but additionally we support B-splines on particle shapes.
@iljah in WarpX we use this nodal field gather: https://github.com/BLAST-WarpX/warpx/blob/25.11/Source/ablastr/particles/NodalFieldGather.H
And this staggered field gather: https://github.com/BLAST-WarpX/warpx/blob/25.11/Source/Particles/Gather/FieldGather.H
WarpX particles have B-spline shapes (e.g., Hockney). This could be written more effectively and with shared memory.
changed the title [-]Cannot get MeshToParticle to interpolate correctly[/-] [+]MeshToParticle Only Works for Nodal Fields
That is not accurate. It's
the Nearest only works if the grid is node-centered data, while the Linear interpolator only works for cell-centered data.
Fixed the title, thanks for the clarification!
Should we start with adding asserts, so people do not use it wrong as it is right now?
In #4822, I proposed
I think we have a few options. I would go with 2.
- Keep the current behavior. Add comments to the functions, including the hack above.
- Add a required index type to the constructors. This will break people's code.
- Add an optional argument for index type.
Note that we don't actually have index information that we can assert because the function takes Array4.
I prefer 2 because the functions are probably not widely used and it would be easy to fix.