How to make grid periodic (for particles [interpolation])?

[iljah]

If I understood correctly, following code should create a periodic 1x1x1 grid spanning [0,1] in all dimensions and interpolate grid data to one particle:

#define AMREX_SPACEDIM 3
#include "AMReX.H"
#include "AMReX_ParticleInterpolators.H"
#include "AMReX_Particles.H"
using Particle = amrex::Particle<1, 0>;
using ParticleContainer = amrex::ParticleContainer<1, 0, 0, 0>;

int main(int argc, char* argv[]) {
	using amrex::IntVect;
	amrex::Initialize(argc, argv); {
	amrex::Box domain{IntVect{0, 0, 0}, IntVect{0, 0, 0}};
	amrex::BoxArray ba;
    ba.define(domain);
    ba.maxSize(2147483647);
	amrex::Geometry geom; geom.define(
		domain, amrex::RealBox({0, 0, 0}, {1, 1, 1}),
		amrex::CoordSys::cartesian, amrex::Array<int, 3>{1, 1, 1});
	amrex::DistributionMapping dm{ba};
	amrex::MultiFab grid{ba, dm, 1, 0};

	ParticleContainer particle_container(geom, dm, ba);
	for (auto iter = particle_container.MakeMFIter(0); iter.isValid(); ++iter) {
		Particle particle;
		particle.id() = 0;
		particle.pos(0) = 0.1;
		particle.pos(1) = 0.1;
		particle.pos(2) = 0.1;
		auto& particles = particle_container.GetParticles(0)[
			std::make_pair(iter.index(), iter.LocalTileIndex())];
		particles.push_back(particle);
	}
	particle_container.Redistribute();

	const auto grid_start = geom.ProbLoArray();
	const auto inv_dr = geom.InvCellSizeArray();
	amrex::MeshToParticle(particle_container, grid, 0, [&](auto& particle1, const auto& data1){
		amrex::ParticleInterpolator::Nearest interp(particle1, grid_start, inv_dr);
		interp.MeshToParticle(particle1, data1, 0, 0, 1,
			[](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;});
	});
	} amrex::Finalize();
	return 0;
}

but if I change e.g. particle.pos(0) = 0.1; to particle.pos(0) = 0.6; I get amrex::Abort::0:: (1,0,0,0) is out of bound (0:0,0:0,0:0,0:0) !!!. How can I make above grid truly periodic (for particle <-> mesh interpolation if that's a separate thing)?

[WeiqunZhang]

The issue is not periodic boundary. The issue is #4810 (reply in thread)
As @AlexanderSinn pointed out,

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.

You could work around the issue by doing

auto grid_start = geom.ProbLoArray();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
    grid_start[idim] += geom.CellSize(idim) * Real(0.5);
}

@ax3l I think you wrote the original code, what's your take? I think we have a few options. I would go with 2.

1. Keep the current behavior. Add comments to the functions, including the hack above.
2. Add a required index type to the constructors. This will break people's code.
3. Add an optional argument for index type.

[iljah]

If I use Linear instead of Nearest it fails no matter where particle is located so doesn't seem like assumption of node or other centered data is the fundamental problem here.

I'd like above program to work regardless of any particular interpolation function etc. while having only real cell, is that possible? Should I introduce ghost/boundary/other cells? Something like this from heat equation test seems relevant phi_old.FillBoundary(geom.periodicity());, would it work with particles or should I use something else?

[AlexanderSinn]

The standard AMReX way would be to add ghost cells, call FillBoundary with the periodicity, and let particles interpolate from the ghost cells. Redistribute would be called on the particles after every push to periodically shift the particles positions so that they are inside the domain. The advantage of this approach is that this also works if the domain is distributed across multiple MPI ranks, as both FillBoundary and Redistribute will do the necessary MPI communications. However, if you only have a single process, it can work without ghost cells by adding a modulo operation to the index calculation for particles. For this, MeshToParticle needs to be replaced with a custom function. This version does nearest grid point interpolation and does not have the 0.5 offset.

    const auto grid_start = geom.ProbLoArray();
    const auto inv_dr = geom.InvCellSizeArray();
    const amrex::IntVect domain_lo = domain.smallEnd();
    const amrex::IntVect domain_len = domain.length();
    amrex::MeshToParticle(particle_container, grid, 0, [&](auto& particle1, const auto& data1){

        const amrex::Real xmid = (particle1.pos(0) - grid_start[0]) * inv_dr[0]; // + 0.5;
        const amrex::Real ymid = (particle1.pos(1) - grid_start[1]) * inv_dr[1]; // + 0.5;
        const amrex::Real zmid = (particle1.pos(2) - grid_start[2]) * inv_dr[2]; // + 0.5;

        const int idx_start_x = static_cast<int>(amrex::Math::floor(xmid));
        const int idx_start_y = static_cast<int>(amrex::Math::floor(ymid));
        const int idx_start_z = static_cast<int>(amrex::Math::floor(zmid));

        // periodically sift the index
        int i = (idx_start_x - domain_lo[0]) % domain_len[0];
        int j = (idx_start_y - domain_lo[1]) % domain_len[1];
        int k = (idx_start_z - domain_lo[2]) % domain_len[2];
        if (i < 0) {
            i += domain_len[0];
        }
        if (j < 0) {
            j += domain_len[1];
        }
        if (k < 0) {
            k += domain_len[2];
        }
        i += domain_lo[0];
        j += domain_lo[1];
        k += domain_lo[2];

        particle1.rdata(0) += data1(i, j, k, 0);

    });

[iljah]

Thanks, fix for Nearest and Linear was surprisingly simple: just had to change amrex::MultiFab grid{ba, dm, 1, 0}; to amrex::MultiFab grid{ba, dm, 1, 1};. Real code will use FillBoundary.