Cannot get MeshToParticle to interpolate correctly

[iljah]

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?

[iljah]

x,y,z should be correct, added (i,j,k) to above output:

Initialized velocity at (-0.5,-0.5)(0,0,0): (-1,1)
Initialized velocity at (0.5,-0.5)(1,0,0): (-1,-1)
Initialized velocity at (-0.5,0.5)(0,1,0): (1,1)
Initialized velocity at (0.5,0.5)(1,1,0): (1,-1)

[iljah]

Same result with amrex::Array<int, 3>{0, 0, 0} instead of amrex::Array<int, 3>{1, 1, 1}.

[AlexanderSinn]

Looking at the index calculation for the Nearest interpolator, to me it looks like it is for node-centered field arrays, making it half a cell off cell-centered data.

amrex/Src/Particle/AMReX_ParticleInterpolators.H

Lines 197 to 199 in 5102120

	amrex::Real l = (p.pos(i) - plo[i]) * dxi[i] + 0.5;
	index[i] = static_cast<int>(amrex::Math::floor(l));
	w[i] = 1;

I believe for cell-centered data, the + 0.5 should be removed.

[iljah]

Thanks, do you mean that nearest interpolator only works with node-centered data even if I'd convert something to cell-centered or if something already is cell-centered? My x,y,z calculation is same as e.g. https://github.com/AMReX-Codes/amrex-tutorials/blob/0da557baac1736dbf99f4a273ebdba3be25dd6dd/ExampleCodes/Basic/HeatEquation_EX0_C/Source/main.cpp#L131 and as far as I can tell I'm also creating cell-centered data.

[iljah]

Hmm linear has same index calculation, I remember seeing many permutations of functions interpolating between cells, faces, edges and vertices in amrex, are there separate functions for those also for particles?

[AlexanderSinn]

If my understanding of ParticleInterpolator is correct, then the Nearest only works if the grid is node-centered data (but does not explicitly check/convert this), while the Linear interpolator only works for cell-centered data. It also has the offset of +0.5 but this is actually needed since it uses two cells per direction. Other codes based on AMReX typically implement custom interpolators that also support higher orders. 
https://github.com/Hi-PACE/hipace/blob/development/src/particles/particles_utils/ShapeFactors.H
https://github.com/BLAST-WarpX/warpx/blob/development/Source/Particles/ShapeFactors.H
Note how every second order has an offset of 0.5. The input for these functions already includes an offset of 0.5.

For testing, it would be good if you could have, for example, a 20^2 grid (from -10 to +10) and then use ParticleToMesh to interpolate a constant of 1 for a particle at a known location and plot the result. 

Yes, I think those interpolation functions are for mesh-to-mesh interpolations. 

[iljah]

Thanks, plot indeed shows interpolation from vertices in a grid with 3 cells spanning [-1,+2] in each dimension. I guess this also explains why interpolation doesn't work within [1.5, 2] since that data would be fetched from 4th cell which doesn't exist.