reset particlecontainer dm and ba

[MaxThevenet]

This PR proposes to define the PlasmaParticleContainer on the slice domain, instead of on the whole domain, which makes more sense for it.

Old implementation

• The slice multifab is defined on a custom BoxArray and DistributionMapping, defined in Fields::AllocData.
• The PlasmaParticleContainer is defined on the whole simulation domain, we just pay attention that all particles have the same z (taken to be the center of the last cell in z)

New implementation

• The slice Geometry, DistributionMapping and BoxArray are defined by function Hipace::DefineSliceGDB called in Hipace::MakeNewLevelFromScratch, and store as Hipace member variables.
• This BA and DM are passed to Fields::AllocData to initialized the slice MF
• This BA and DM and Geom are passed to Hipace::InitData, where we call the plasmaparticlecontainer Set functions for Geom, BA, DM before calling PlasmaParticleContainer::InitData.

This PR works OK for most tests, only the plasma z component checksum fails for the linear wake. But the z positions of plasma particles are meaningless, so this is not a problem (I would be curious to understand the change though).

However, when running in Debug mode, the code crashes at the AMREX_ASSERT(OK()); assertion at the end of PlasmaParticleContainer::InitParticles, and putting print statements couldn't help me understand this issue. @atmyers could you have a look, and let me know if you see something off in this PR?

• Small enough (< few 100s of lines), otherwise it should probably be split into smaller PRs
• Tested (describe the tests in the PR description)
• Runs on GPU (basic: the code compiles and run well with the new module)
• Contains an automated test (checksum and/or comparison with theory)
• Documented: all elements (classes and their members, functions, namespaces, etc.) are documented
• Constified (All that can be const is const)
• Code is clean (no unwanted comments, )
• Style and code conventions are respected at the bottom of https://github.com/Hi-PACE/hipace
• Proper label and GitHub project, if applicable

[atmyers]

I think that maybe there is an inconsistency between the z positions set during InitParticles and what amrex expects. If the BoxArray used to define the slice has a z-index of n, then the z position of the particle must be set so that this function:

https://github.com/AMReX-Codes/amrex/blob/development/Src/Particle/AMReX_ParticleUtil.H#L236

maps to n in the z direction. I'm but sure that's the case in the current PR. This coordinate is ignored in Hipace (as I understand it), but it is used in amrex functions like OK().

[MaxThevenet]

Thanks for your reply! So in the way the code is written, the slice domain should match exactly the last cell (in z) of the 3D domain. This means that the BoxArray should have only 1 cell (indexed 0) in the z direction, the prob_domain should be the left and right edge of the last cell, etc. Note that the particles are still given a z position consistent with that. In the example I am running, I would expect, for amrex::Math::floor((p.pos(2)-plo[2])*dxi[2]))

dx[2] = 0.6e-6
plo[2] = 5.94e-6 // 1or wherever this left edge is
p.pos[2] = 5.97e-6 // in the centre of this cell

so amrex::Math::floor((p.pos(2)-plo[2])*dxi[2])) should return 0, which is the only cell in the slice.

I checked that these things are right in ParticleContainerInit. But I think you are right, and another point is: when commenting out the line SetPlasmaGeometry the code PASSES THE OK.

[MaxThevenet]

These are the lines I am talking about.

[MaxThevenet]

@atmyers I think I see what is happening here. When I add some print statements to the lines where we reset the plasma container's geometry, BA and DM, replacing

    m_plasma_container.SetParticleBoxArray(lev, m_slice_ba);
    m_plasma_container.SetParticleDistributionMap(lev, m_slice_dm);
    m_plasma_container.SetParticleGeometry(lev, m_slice_geom);

by

    std::cout<<"m_slice_ba = \n"<<m_slice_ba<<'\n';
    std::cout<<"before set: m_plasma_container.ParticleBoxArray(lev) = \n"<<m_plasma_container.ParticleBoxArray(lev)<<'\n';
    m_plasma_container.SetParticleBoxArray(lev, m_slice_ba);
    std::cout<<"after  set: m_plasma_container.ParticleBoxArray(lev) = \n"<<m_plasma_container.ParticleBoxArray(lev)<<'\n';

    std::cout<<"m_slice_dm = \n"<<m_slice_dm<<'\n';
    std::cout<<"before set: m_plasma_container.ParticleDistributionMap(lev) = \n"<<m_plasma_container.ParticleDistributionMap(lev)<<'\n';
    m_plasma_container.SetParticleDistributionMap(lev, m_slice_dm);
    std::cout<<"before set: m_plasma_container.ParticleDistributionMap(lev) = \n"<<m_plasma_container.ParticleDistributionMap(lev)<<'\n';

    std::cout<<"m_slice_geom = \n"<<m_slice_geom<<"\n\n";
    std::cout<<"before set: m_plasma_container.Geom(lev) = \n"<<m_plasma_container.Geom(lev)<<"\n\n";
    m_plasma_container.SetParticleGeometry(lev, m_slice_geom);
    std::cout<<"after  set: m_plasma_container.Geom(lev) = \n"<<m_plasma_container.Geom(lev)<<"\n\n";

I get the following output:

m_slice_ba = 
(BoxArray maxbox(1)
       m_ref->m_hash_sig(0)
       ((0,0,0) (63,63,0) (0,0,0)) )

before set: m_plasma_container.ParticleBoxArray(lev) = 
(BoxArray maxbox(1)
       m_ref->m_hash_sig(0)
       ((0,0,0) (63,63,399) (0,0,0)) )

after  set: m_plasma_container.ParticleBoxArray(lev) = 
(BoxArray maxbox(1)
       m_ref->m_hash_sig(0)
       ((0,0,0) (63,63,0) (0,0,0)) )

m_slice_dm = 
(DistributionMapping
m_pmap[0] = 0
)

before set: m_plasma_container.ParticleDistributionMap(lev) = 
(DistributionMapping
m_pmap[0] = 0
)

before set: m_plasma_container.ParticleDistributionMap(lev) = 
(DistributionMapping
m_pmap[0] = 0
)

m_slice_geom = 
(0 (-0.0001,-0.0001,5.97e-05)(3.125e-06,3.125e-06,3e-07) 1)
(RealBox -0.0001 0.0001 -0.0001 0.0001 5.97e-05 6e-05 )((0,0,0) (63,63,0) (0,0,0))P(1,1,0)

before set: m_plasma_container.Geom(lev) = 
(0 (-0.0001,-0.0001,-6e-05)(3.125e-06,3.125e-06,3e-07) 1)
(RealBox -0.0001 0.0001 -0.0001 0.0001 -6e-05 6e-05 )((0,0,0) (63,63,399) (0,0,0))P(1,1,0)

after  set: m_plasma_container.Geom(lev) = 
(0 (-0.0001,-0.0001,-6e-05)(3.125e-06,3.125e-06,3e-07) 1)
(RealBox -0.0001 0.0001 -0.0001 0.0001 -6e-05 6e-05 )((0,0,0) (63,63,399) (0,0,0))P(1,1,0)

so it seems that the box array (and probably the DM, but that's hard to say) is properly reset, but the Geometry is not. @atmyers is this expected? Do you see something wrong in this code, or is this on the AMReX side? Note that, printing the particles longitudinal positions in numParticlesOutOfRange gave the right value for these tests (between 59.7e-6 and 60.0e-6).

[atmyers]

I think it's on the AMReX side - the problem seems to be that pc.Geom refers to the mesh geometry, but pc.ParticleGeom' points to the particle one. But internally in many places, including Redistribute, pc.Geomis used. The fixes are either to makepc.Geomalways refer to the particle one, or find and replace all the placespc.Geomis used withpc.ParticleGeom`. Probably the first is easier / more intuitive?

[MaxThevenet]

Both would work for me, but indeed the first one seems much easier to do. Should we try this?

[MaxThevenet]

Is this a change in AMReX, or should we override pc.Geom in Hipace++?

[atmyers]

Yes could you try this branch: AMReX-Codes/amrex#1470

[atmyers]

Is this a change in AMReX, or should we override pc.Geom in Hipace++?

This is a change in AMReX, no changes to Hipace++ should be needed.

[atmyers]

Hi @MaxThevenet,

Could you try AMReX-Codes/amrex#1532 and see if it fixes this issue?

Thanks,
Andrew

[MaxThevenet]

@atmyers this works like a charm! More precisely, after orthogonal changes in other PRs in Hipace++, this PR did not even run anymore, with an error in locate_particle. When running on your AMReX branch instead (just adding some cmake options for our CI, see last commit), all tests pass again \o/. I also tested 1 example in debug mode locally, and it runs fine, so I think we're on the right track here, thanks a lot!

[atmyers]

Cool, that has been merged into AMReX now.

[MaxThevenet]

Yep, thanks, I just updated this PR accordingly. I also removed the WIP tag, this is ready for review. Once this is merged, we can implement the longitudinal parallelisation for plasma particles.

[MaxThevenet]

Thanks for your review :)