Interpolation from finer grid to particle?

[iljah]

How should I interpolate from finer grid to particle in following program? Now nothing seems to happen. If I change mf_f, 1 to mf_c, 0 etc. in MeshToParticle() then it interpolates from coarser grid correctly.

#define AMREX_SPACEDIM 3
#include "AMReX.H"
#include "AMReX_MultiFab.H"
#include "AMReX_ParticleContainer.H"
#include "AMReX_ParticleInterpolators.H"
using namespace std;
using namespace amrex;
using ParticleContainer_ = ParticleContainer<3, 0, 0, 0>;
using Particle_ = ParticleContainer_::ParticleType;
int main(int argc, char* argv[]) {
	Initialize(argc, argv); {
	Box
		box_coarse{IntVect{0, 0, 0}, IntVect{2, 2, 2}},
		box_fine{IntVect{2, 2, 2}, IntVect{3, 3, 3}};
	Geometry
		geom_c{box_coarse, RealBox({-1.5,-1.5,-1.5}, {1.5,1.5,1.5}),
			CoordSys::cartesian, Array<int, 3>{0, 0, 0}},
		geom_f{box_fine, RealBox({-0.5,-0.5,-0.5}, {0.5,0.5,0.5}),
			CoordSys::cartesian, Array<int, 3>{0, 0, 0}};
	geom_f.refine(IntVect{2});
	BoxArray ba_c, ba_f;
	ba_c.define(box_coarse); ba_f.define(box_fine);
	DistributionMapping dm_c{ba_c}, dm_f{ba_f};
	MultiFab mf_c{ba_c, dm_c, 3, 0}, mf_f{ba_f, dm_f, 3, 0};

	for (amrex::MFIter iter(mf_c); iter.isValid(); ++iter) {
		const auto& box = iter.validbox();
		const auto& arr = mf_c.array(iter);
		amrex::ParallelFor(box, [&](int i, int j, int k){
			arr(i, j, k, 0) = i;
			arr(i, j, k, 1) = j;
			arr(i, j, k, 2) = k;
		});
	}
	for (amrex::MFIter iter(mf_f); iter.isValid(); ++iter) {
		const auto& box = iter.validbox();
		const auto& arr = mf_f.array(iter);
		amrex::ParallelFor(box, [&](int i, int j, int k){
			arr(i, j, k, 0) = i;
			arr(i, j, k, 1) = j;
			arr(i, j, k, 2) = k;
		});
	}

	ParticleContainer_ particle_container(
		{geom_c, geom_f}, {dm_c, dm_f},
		{ba_c, ba_f}, Vector<int>{2});
	for (
		auto iter = particle_container.MakeMFIter(0);
		iter.isValid(); ++iter
	) {
		auto& particles = particle_container.GetParticles(0)[
			std::make_pair(iter.index(), iter.LocalTileIndex())];
		Particle_ particle;
		particle.id() = 1;
		for (int i = 0; i < 3; i++) {
			particle.pos(i) = particle.rdata(i) = 0;
		}
		particles.push_back(particle);
	}
	particle_container.Redistribute();

	for (
		ParticleContainer_::ParIterType iter(particle_container, 0);
		iter.isValid(); ++iter
	) {
		const auto& particles = iter.GetArrayOfStructs();
		for (const auto& particle: particles) {
			cout << "before:";
			for (int i = 0; i < 3; i++) {
				cout<<particle.pos(i)<<","<<particle.rdata(i)<<",";
			} cout<<endl;
		}
	}

	amrex::MeshToParticle(particle_container, mf_f, 1, [&](auto& particle1, const auto& data1){
		amrex::ParticleInterpolator::Linear interp(
			particle1, geom_f.ProbLoArray(), geom_f.InvCellSizeArray());
		interp.MeshToParticle(particle1, data1, 0, 0, 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) {
			cout << "after:";
			for (int i = 0; i < 3; i++) {
				cout<<particle.pos(i)<<","<<particle.rdata(i)<<",";
			} cout<<endl;
		}
	}
	} Finalize(); return 0;
}

[WeiqunZhang]

The way how fine geometry is created is incorrect. The domain passed to the constructor should be the entire domain as if it's refined everywhere.

It is almost always the case that the fine Geometry is just amrex::refine(coarse_geometry, refinement_ratio). Here coarse_geometry is truly the coarse geometry. For example,

	Box
		box_coarse{IntVect{0, 0, 0}, IntVect{2, 2, 2}},
	Geometry
		geom_c{box_coarse, RealBox({-1.5,-1.5,-1.5}, {1.5,1.5,1.5}),
			CoordSys::cartesian, Array<int, 3>{0, 0, 0}},
	Geometry geom_f = amrex::refine(geom_c, IntVect{2});

See #4907.

[iljah]

Thanks but after refining as suggested the program prints neither before nor after:

-               geom_f{box_fine, RealBox({-0.5,-0.5,-0.5}, {0.5,0.5,0.5}),
-                       CoordSys::cartesian, Array<int, 3>{0, 0, 0}};
-       geom_f.refine(IntVect{2});
+               geom_f = refine(geom_c,IntVect(2));

No idea why that is.

[WeiqunZhang]

The next issue may be that particles are not on the fine level, because your fine boxarray comes from box_fine{IntVect{2, 2, 2}, IntVect{3, 3, 3}}. x=0.0 and y=0.0 would be exactly at the lower corner of the fine box. So maybe because of roundoff error, they are redistributed to coarse levels. You could try pos(i) = 0.0001.

Sorry for my curiosity. May I ask what's the purpose of these toy examples?

[WeiqunZhang]

Looks like the answer to your last question might be it does not print anything because the Particle iterator you use is on level 0, not level 1. AI is usually very good at catching these bugs. Here is what it suggests.

Summary of Recommended Changes

    Ghost Cells: Define mf_f with 1 ghost cell.
    Consistency: Define geom_f using amrex::refine(geom_c, 2).
    Output: Iterate over level 1 in your final print loop to see the particles that were moved there by Redistribute().

[WeiqunZhang]

This is the corrected code provided by an AI at LBNL.

#define AMREX_SPACEDIM 3
#include "AMReX.H"
#include "AMReX_MultiFab.H"
#include "AMReX_ParticleContainer.H"
#include "AMReX_ParticleInterpolators.H"

using namespace amrex;

using ParticleContainer_ = ParticleContainer<3, 0, 0, 0>;
using Particle_ = ParticleContainer_::ParticleType;

int main(int argc, char* argv[]) {
    Initialize(argc, argv); {
        // Coarse: 3 cells, indices 0, 1, 2. RealBox [-1.5, 1.5] -> dx=1.0
        Box box_coarse{IntVect{0, 0, 0}, IntVect{2, 2, 2}};
        // Fine: covers coarse index 1. With ref_ratio 2, indices are 2 and 3.
        Box box_fine{IntVect{2, 2, 2}, IntVect{3, 3, 3}};

        Geometry geom_c{box_coarse, RealBox({-1.5,-1.5,-1.5}, {1.5,1.5,1.5}),
                        CoordSys::cartesian, Array<int, 3>{0, 0, 0}};
        
        // Define geom_f directly on the refined index space. 
        // RealBox must match the physical footprint of indices 2-3 (which is [-0.5, 0.5])
        Geometry geom_f{box_fine, RealBox({-0.5,-0.5,-0.5}, {0.5,0.5,0.5}),
                        CoordSys::cartesian, Array<int, 3>{0, 0, 0}};

        BoxArray ba_c(box_coarse), ba_f(box_fine);
        DistributionMapping dm_c{ba_c}, dm_f{ba_f};

        // Use 1 ghost cell for interpolation safety
        MultiFab mf_c{ba_c, dm_c, 3, 1}, mf_f{ba_f, dm_f, 3, 1};

        // Initialize mf_f with values
        for (MFIter mfi(mf_f); mfi.isValid(); ++mfi) {
            const auto& box = mfi.validbox();
            const auto& arr = mf_f.array(mfi);
            ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                arr(i, j, k, 0) = i * 10.0; // Distinct values to verify
                arr(i, j, k, 1) = j * 10.0;
                arr(i, j, k, 2) = k * 10.0;
            });
        }
        mf_f.FillBoundary(geom_f.periodicity());

        ParticleContainer_ pc({geom_c, geom_f}, {dm_c, dm_f}, {ba_c, ba_f}, {2});
        
        // Add particle at level 0; Redistribute will move it to level 1 
        // because (0,0,0) is inside geom_f's RealBox
        if (ParallelDescriptor::MyProc() == 0) {
            auto& particles = pc.GetParticles(0)[std::make_pair(0, 0)];
            Particle_ p;
            p.id()  = Particle_::NextID();
            p.cpu() = ParallelDescriptor::MyProc();
            p.pos(0) = 0.0; p.pos(1) = 0.0; p.pos(2) = 0.0;
            for(int i=0; i<3; ++i) p.rdata(i) = 0.0;
            particles.push_back(p);
        }
        pc.Redistribute();

        // MeshToParticle on Level 1
        amrex::MeshToParticle(pc, mf_f, 1, [&](auto& p, const Array4<const Real>& data) {
            // Use local geometry info for the current level
            amrex::ParticleInterpolator::Linear interp(p, geom_f.ProbLoArray(), geom_f.InvCellSizeArray());
            
            interp.MeshToParticle(p, data, 0, 0, 3,
                [=](const Array4<const Real>& d, int i, int j, int k, int c) { return d(i, j, k, c); },
                [=](Particle_& p_ref, int c, Real val) { p_ref.rdata(c) += val; }
            );
        });

        // Print results from level 1
        for (ParticleContainer_::ParIterType pti(pc, 1); pti.isValid(); ++pti) {
            const auto& particles = pti.GetArrayOfStructs();
            for (int i = 0; i < pti.numParticles(); ++i) {
                const auto& p = particles[i];
                amrex::Print() << "After Interp: " << p.rdata(0) << ", " << p.rdata(1) << ", " << p.rdata(2) << std::endl;
            }
        }
    } Finalize(); return 0;
}

[WeiqunZhang]

For whatever reasons, it did not apply its own suggestion of Consistency: Define geom_f using amrex::refine(geom_c, 2)., probably because in this case, both work.