add MultiFabRegister / ParticleContainerRegister classes

[BenWibking]

Describe the proposal
We often need to do the same set of operations over all MultiFabs and ParticleContainers, but C++ doesn't provide a capability to iterate over all of the member variables in a class, so we have to hard code them. However, can instead dynamically allocate MultiFabs and keep track of a list of all of the MultiFabs ourselves.

WarpX does this with its MultiFabRegister. This creates an object that behaves like a dictionary (e.g., in Python), where we can access (or create) the MultiFabs with a string, and also loop over all of them with an iterator.

We also need to do this with particle containers, since we may have many particle containers, and we don't want to hard code everywhere we need to redistribute particles, write them to disk, etc., for each particle container individually.

If we need to use the MultiFabs inside a amrex::ParallelFor kernel, we just get a reference to the MultiFab from the MultiFabManager while outside the ParallelFor, and from then on, use it as we normally would.

Describe alternatives you've considered
Keep as-is. If we need to loop over all MultiFabs or ParticleContainers, we do it manually.

Additional context
https://github.com/ECP-WarpX/WarpX/blob/fb9e7dd1a2946c1f293cd25f2061d78cf5d1d71f/Source/ablastr/fields/MultiFabRegister.H#L156

[BenWibking]

This will allow us to split state_old_cc_ and state_new_cc_ into separate hydro and radiation multifabs, which should have a small positive performance impact (since the radiation multifab does not need to have its ghost cells filled during the hydro update).

[BenWibking]

Also, this will allow us to split the state_old_fc_ and state_new_fc_ multifabs into divergence-free and non-divergence-free components. This is needed for AMR MHD, since we have to use a different AMR interpolater for the magnetic field and the face-centered velocity.

[BenWibking]

For reference, this is the WarpX PR for its equivalent MultiFabRegister class: ECP-WarpX/WarpX#5230

[markkrumholz]

Just following up on our call, my idea is that a particle type is defined by something like the following class:

class PhysicsParticleDescriptor {
   protected:
      *NeighborParticleContainer;  // pointer to the NeighborParticleContainer that contains particles of this type; needs to be a pointer so that it can work with containers with different template parameters
      int massIndex; // index that contains the mass for gravity; set to -1 for particle types that do not interact with gravity
      int lumIndex;  // index that contains the luminosity in radiation group 0; set to -1 for particle types that do not radiate
      bool interactsWithHydro;   // true if particle type interacts with hydro; the class must then defined the method hydroInteract
      virtual void hydroInteract(...) { } // method to implement interaction with hydro; defaults to a no-op, but can be overridden by derived classes
}

Then we have PhysicsParticleRegister, which contains a std::map containing (string, PhysicsParticleDescriptor) pairs, and which defines the depositGravity and deposityRadiation methods to add gravity and radiation source terms to the grid, and which is called on any member of the map that sets massIndex or lumIndex to a value other than -1.

[BenWibking]

This makes sense to me.

• I would prefer to avoid raw pointers if possible. We could use std::unique_ptr<NeighborParticleContainer> instead of *NeighborParticleContainer to do so. It is possible to assign a derived class type to a unique_ptr of the base class type, e.g.:

  std::unique_ptr<base> derived = std::make_unique<Derived>();
  

• There is a potential complication with putting this in a collection like std::map<std::string, quokka::PhysicsParticleDescriptor>. If this is a virtual base class, we would have to dynamic_cast from the derived class to the base class to put this in the same container.

  This is okay, but unnecessarily complicated. We could instead use std::function, which provides a safe function pointer: https://en.cppreference.com/w/cpp/utility/functional/function. std::function objects can be assigned named functions or lambda functions, e.g.:

  std::function<void(const Foo&, int)> f_add_display = &Foo::print_add;
  

  where Foo::print_add is a function that takes const Foo&, int as arguments and returns void.

---

So my preferred version of this object is:

struct PhysicsParticleDescriptor {
      std::unique_ptr<NeighborParticleContainer>;  // pointer to the NeighborParticleContainer that contains particles of this type
      int massIndex; // index that contains the mass for gravity; set to -1 for particle types that do not interact with gravity
      int lumIndex;  // index that contains the luminosity in radiation group 0; set to -1 for particle types that do not radiate
      bool interactsWithHydro;   // true if particle type interacts with hydro; the class must then defined the method hydroInteract
      std::function<void(...)> hydroInteract; // method to implement interaction with hydro
};

Then, the PhysicsParticleRegister can be a simple wrapper around a std::map<std::string, quokka::PhysicsParticleDescriptor>.

[BenWibking]

Note that function calls to neither virtual functions nor std::functions can be inlined by the compiler. So for performance reasons, they should not be called within an inner loop.

[markkrumholz]

Agreed: the idea is that these functions would themselves contain the on-GPU loops, rather than being called from within them. The general API should probably provide access to the full state vector. We may also need to think a bit about how to generalize this with respect to time stepping. I was thinking that these functions would be called in the Strang split stage, but for some things -- particularly tracer particles -- you might want to call them during the hydro update itself and have access to intermediate quantities like the outputs of the Riemann solvers. Not sure there is a really general way to do that.

[BenWibking]

Tracer particles could actually be handled in a split way. For second order accuracy, they only need to be advected with the Riemann solver velocity averaged over the hydro step. Currently, this is handled within the hydro update itself, because we don't average the Riemann solver velocity over multiple timesteps when the retry mechanism is triggered.

If we propertly average the face-centered velocities over time (including retry substeps), then we can advect the tracer particles in a separate function. Strictly speaking, we need to do that anyway in order to provide the correct boundary conditions for the tracer particle velocities on refined AMR levels (since the coarser levels could have had retries triggered).

Edit: I think the tracer velocity boundary conditions are not correct currently when we are doing AMR subcycling in time. This is a somewhat tricky issue to fix.

[BenWibking]

I think all of the particle interactions can be handled in a split way.

[markkrumholz]

Great, then no problem at will the design. We can just put a loop in the StrangSplitSourceTerm routine to loop over the elements of PhysicsParticleRegister, for each one check if interactsWithHydro is true (though maybe we should change the name to something more general, say doStrangSplitUpdate) and then calls the appropriate function.

[chongchonghe]

I finished a first version of this design: #820 . It's working great for CICParticles and RadParticles. I did this before even reading your comments on Dec 12, so I haven't thought about time stepping yet.

[chongchonghe]

This makes sense to me.

• I would prefer to avoid raw pointers if possible. We could use std::unique_ptr<NeighborParticleContainer> instead of *NeighborParticleContainer to do so. It is possible to assign a derived class type to a unique_ptr of the base class type, e.g.:

  std::unique_ptr<base> derived = std::make_unique<Derived>();
  

• There is a potential complication with putting this in a collection like std::map<std::string, quokka::PhysicsParticleDescriptor>. If this is a virtual base class, we would have to dynamic_cast from the derived class to the base class to put this in the same container.
  This is okay, but unnecessarily complicated. We could instead use std::function, which provides a safe function pointer: https://en.cppreference.com/w/cpp/utility/functional/function. std::function objects can be assigned named functions or lambda functions, e.g.:

  std::function<void(const Foo&, int)> f_add_display = &Foo::print_add;
  

  where Foo::print_add is a function that takes const Foo&, int as arguments and returns void.

So my preferred version of this object is:

struct PhysicsParticleDescriptor {
      std::unique_ptr<NeighborParticleContainer>;  // pointer to the NeighborParticleContainer that contains particles of this type
      int massIndex; // index that contains the mass for gravity; set to -1 for particle types that do not interact with gravity
      int lumIndex;  // index that contains the luminosity in radiation group 0; set to -1 for particle types that do not radiate
      bool interactsWithHydro;   // true if particle type interacts with hydro; the class must then defined the method hydroInteract
      std::function<void(...)> hydroInteract; // method to implement interaction with hydro
};

Then, the PhysicsParticleRegister can be a simple wrapper around a std::map<std::string, quokka::PhysicsParticleDescriptor>.

I'm using void *neighborParticleContainer_{}. I use radParticleDesc->neighborParticleContainer_ = RadParticles.get(); to add Particles into the containder, and use auto *container = static_cast<RadParticleContainer<problem_t> *>(descriptor->neighborParticleContainer_); to get the container. These are all recommended by VSCode Copilot (it's 10x smarter than me). It seems like in this way, you can put any XXXParticleContainer types into neighborParticleContainer_.

[chongchonghe]

Both of you talked about NeighborParticleContainer; do you mean amrex::NeighborParticleContainer, or just a concept of a container that contains amrex::AmrParticleContainer? amrex::NeighborParticleContainer doesn't exist in AMReX. Search NeighborParticleContainer in extern/amrex results in nothing. I'm using CICParticleContainer = amrex::AmrParticleContainer<4> like in the original CICParticles.hpp, and I'm using neighborParticleContainer_ in PhysicsParticleDescriptor to contain particles, but the 'neighbor' doesn't have any specific meaning.

[markkrumholz]

This is the class we mean: https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1NeighborParticleContainer.html.

[chongchonghe]

OK, I found it. Turns out it's because my vscode is ignore amrex in earch: "search.exclude": {"extern/amrex": true}.

[chongchonghe]

Can't use Neighbor Particles right now because we need AmrNeighborParticleContainer and it's currently missing. I'll wait for Andrew Myers's response in Slack.

[BenWibking]

This makes sense to me.

• I would prefer to avoid raw pointers if possible. We could use std::unique_ptr<NeighborParticleContainer> instead of *NeighborParticleContainer to do so. It is possible to assign a derived class type to a unique_ptr of the base class type, e.g.:

  std::unique_ptr<base> derived = std::make_unique<Derived>();
  

• There is a potential complication with putting this in a collection like std::map<std::string, quokka::PhysicsParticleDescriptor>. If this is a virtual base class, we would have to dynamic_cast from the derived class to the base class to put this in the same container.
  This is okay, but unnecessarily complicated. We could instead use std::function, which provides a safe function pointer: https://en.cppreference.com/w/cpp/utility/functional/function. std::function objects can be assigned named functions or lambda functions, e.g.:

  std::function<void(const Foo&, int)> f_add_display = &Foo::print_add;
  

  where Foo::print_add is a function that takes const Foo&, int as arguments and returns void.

So my preferred version of this object is:

struct PhysicsParticleDescriptor {
      std::unique_ptr<NeighborParticleContainer>;  // pointer to the NeighborParticleContainer that contains particles of this type
      int massIndex; // index that contains the mass for gravity; set to -1 for particle types that do not interact with gravity
      int lumIndex;  // index that contains the luminosity in radiation group 0; set to -1 for particle types that do not radiate
      bool interactsWithHydro;   // true if particle type interacts with hydro; the class must then defined the method hydroInteract
      std::function<void(...)> hydroInteract; // method to implement interaction with hydro
};

Then, the PhysicsParticleRegister can be a simple wrapper around a std::map<std::string, quokka::PhysicsParticleDescriptor>.

I'm using void *neighborParticleContainer_{}. I use radParticleDesc->neighborParticleContainer_ = RadParticles.get(); to add Particles into the containder, and use auto *container = static_cast<RadParticleContainer<problem_t> *>(descriptor->neighborParticleContainer_); to get the container. These are all recommended by VSCode Copilot (it's 10x smarter than me). It seems like in this way, you can put any XXXParticleContainer types into neighborParticleContainer_.

No, this is bad. Void pointers are unsafe and should never be used. You should use a std::unique_ptr to the base class. Casting to/from the derived to base class pointer type with dynamic_cast is safe.

[markkrumholz]

Void pointers are unsafe and should never be used.

Spoken like a modern-day C++ object-oriented coward! When I learned programming it was in C, and we didn't bother with "safety". We just accessed raw memory like Real Men! ;-)

[BenWibking]

The C++ Core Guidelines are worth reading. They are written by real C++ experts (many of whom wrote the C++ standard), not an LLM. Guideline C.146 says:

C.146: Use dynamic_cast where class hierarchy navigation is unavoidable
Reason dynamic_cast is checked at run time.

See https://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines#c146-use-dynamic_cast-where-class-hierarchy-navigation-is-unavoidable for more details.