Possible to use 2x2x1 grid (with particles) in 3d?

[iljah]

Following 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 std::signbit;
	using amrex::IntVect;

	amrex::Initialize(argc,argv); {

	amrex::Box domain{IntVect{0, 0, 0}, IntVect{1, 1, 0}};
	amrex::BoxArray ba;
    ba.define(domain);
    ba.maxSize(2147483647);

	amrex::Geometry geom; geom.define(
		domain, amrex::RealBox({-1, -1, -0.5}, {1, 1, 0.5}),
		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};
std::cout<<__FILE__"("<<__LINE__<<")"<<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 / 10.0;
			particle.pos(1) = y / 10.0;
			particle.pos(2) = 0.2;
			for (int i = 0; i < nr_extra_float_vars; i++) {
				particle.rdata(i) = 0;
			}
			particles.push_back(particle);
		}
	}
	particle_container.Redistribute();
std::cout<<__FILE__"("<<__LINE__<<")"<<std::endl;

	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;
			});
	});
std::cout<<__FILE__"("<<__LINE__<<")"<<std::endl;
	} amrex::Finalize();
}

prints

Initializing AMReX (25.11)...
AMReX (25.11) initialized
main.cpp(63)
main.cpp(84)
amrex::Abort::0:: (1,1,1,0) is out of bound (0:1,0:1,0:0,0:2) !!!

which I can fix by using IntVect{1, 1, 1} instead of IntVect{1, 1, 0} but then I get 2 cells in z dimension. Is it possible to have 1 cell in any dimension with (3d?) AMReX?

[iljah]

If I change particle.pos(2) = 0.2; to particle.pos(2) = -0.2; it works, which is strange since particles should still be in same (and only) cell as realbox extends from -0.5 to +0.5 in z.

[AlexanderSinn]

You are using AMReX with Spacedim==3 for effectively a 2D problem. This is possible to do; however, some functions don't quite work correctly. In this MeshToParticle does interpolation in all 3 dimensions, including the one you don't have in the array. For performance reasons, bounds checking is usually omitted in particle interpolation, and it is assumed that the grids have enough ghost cells, which is not the case here.

You can use MeshToParticle with only two dimensions; however, you need to supply your own ParticleInterpolator that does not do any interpolation in the 3rd dimension.
See the bottom of https://github.com/AMReX-Codes/amrex/blob/development/Src/Particle/AMReX_ParticleInterpolators.H for the definitions of the existing two interpolators.

Basically you need to copy the implementation but replace AMREX_SPACEDIM with 2 (note: untested):

struct NearestInterp2D : public amrex::ParticleInterpolator::Base<NearestInterp2D, int>
{
    static constexpr int stencil_width = 1;
    int weights[3*stencil_width];
    static constexpr int dim = 2;

    static constexpr int nx = (dim >= 1) ? stencil_width - 1 : 0;
    static constexpr int ny = (dim >= 2) ? stencil_width - 1 : 0;
    static constexpr int nz = (dim >= 3) ? stencil_width - 1 : 0;

    template <typename P>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    NearestInterp2D (const P& p,
             amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
             amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi)
    {
        w = &weights[0];
        for (int i = 0; i < dim; ++i) {
            amrex::Real l = (p.pos(i) - plo[i]) * dxi[i] + 0.5;
            index[i] = static_cast<int>(amrex::Math::floor(l));
            w[i] = 1;
        }
        for (int i = dim; i < 3; ++i) {
            index[i] = 0;
            w[i] = 1;
        }
    }
};

And then replace amrex::ParticleInterpolator::Nearest with NearestInterp2D

[AlexanderSinn]

And possibly replace (p.pos(i) - plo[i]) * dxi[i] + 0.5; with (p.pos(i) - plo[i]) * dxi[i];

[iljah]

Thanks but looks like problem concerns all dimensions and probably only manifests if particle is located in positive half of last cell in particular dimension. If I change

particle.pos(0) = x / 10.0;
particle.pos(1) = y / 10.0;
particle.pos(2) = 0.2;

to

particle.pos(0) = x / 2.0;
particle.pos(1) = y / 2.0;
particle.pos(2) = -0.2; // note - not +

then I get amrex::Abort::0:: (1,2,0,0) is out of bound (0:1,0:1,0:0,0:2) !!!. If grid extents are [-1,+1] then it seems that this abort is triggered if particle(s) are in [0,+1] with 1 cell, [+0.5,+1] with 2 cells, etc.

[iljah]

Turns out that nearest neighbor interpolator assumes vertex-centered data so would fail whenever particle existed in last half of last cell in grid: #4810

Linear interpolator always needs two cells so how should I prepare my grid to make that work in one cell periodic grid? Should I use amrex::convert() somewhere or FillBoundary or some boundary condition?..

[iljah]

Solved in another question: changing amrex::MultiFab V{ba, dm, 3, 0}; to amrex::MultiFab V{ba, dm, 3, 1}; removes above abort.