MeshToParticle Only Works for Nodal Fields

[WeiqunZhang]

Discussed in #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?

[WeiqunZhang]

#4810 (comment)

[ax3l]

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.

[WeiqunZhang]

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.

[WeiqunZhang]

See also #4822