Memory issues with larger numbers of particles and replicated meshes

[Waqar-ukaea]

Right now AEGIS relies upon a replicating the mesh for each MPI rank. For large enough meshes and enough MPI ranks this causes me to run out of memory on my local system (~32GB). I don't recall running into memory issues when working on HPC systems such as CSD3 though.

It might be worth investigating methods of reducing this memory overhead. Options to consider:

• Use valgrind-massif to investigate memory usage at runtime.
• Use valgrind-memcheck to check for memory leaks.
• Implementing a hybrid MPI + OpenMP parallel algorithm so that a single MPI rank is allocated per node and OpenMP threads are used to access the same mesh available to each MPI rank.
• MOAB has some parallel mesh related operations in its documentation - so another route could be to explore what options I have here. Could be a route in for a distributed mesh?

[Waqar-ukaea]

In hindsight, this might be a combination of two different problems - the replicated meshes incur a memory hit at the start but my naive implementation of a "root array" of particles on rank 0 is causing the memory usage to blow up with large enough meshes. I am relying on a std::vector<Particle> which uses std::vector::emplace_back() to expand the array of particles on rank 0. However this problem is exacerbated when working with large meshes, since AEGIS launches at least 1 particle per mesh element, a larger mesh does result in more particles being launched. This is resulting in performance issues when the array needs to have its size reallocated since emplace_back() calls to double the size of the array - so as the array gets larger the memory overhead increases massively.

Edit - Changed the title to reflect that this is more than just a replicated mesh or array of particles problem.

The line in question can be found in ParticleSimulation::setup_sources:556

[Waqar-ukaea]

I've had some more thoughts about the mesh replication side of things. Right now AEGIS reads in an entire MOAB .h5m mesh alongside the surfaces of interest provided in the aegis_settings.json. I will need to have a deeper dive into the code to see what exactly us happening but I think the full mesh is replicated on each MPI rank which is pretty wasteful memory wise. We don't need a fully meshed geometry on each rank, instead the only required things are:

• The full mesh of the target surfaces
• A coarse mesh of the shadowing geometry (including the target surfaces)

This does beg the question - would it actually be more beneficial to do what SMARDDA does and have two separate meshes? A coarse faceting of the entire geometry (target and shadow) and a fine faceted mesh of the target. The added convenience of using a single mesh for shadowing and target may not actually be worth it in this case.

[makeclean]

We should have a chat about this, my feeling is a different MPI parallelisation strategy using hybrid OpenMP & MPI, i.e. each MPI task reproduces the whole problem but then OpenMP threads handle the on-node parallelism, akin to how a traditional neutronic code would handle it.

I think its beneficial to have one geometry all at high resolution, since I think with the above change, we should be able to handle it. Might need to think about how to handle the condensing of results, just moving to OpenMP you will be reducing your memory usage by a factor of Ncpuspernode, for most problems I can recall, the memory usage of Double Down problems is typically less than 12 Gb for the whole mesh & BVH.

[Waqar-ukaea]

We should have a chat about this, my feeling is a different MPI parallelisation strategy using hybrid OpenMP & MPI, i.e. each MPI task reproduces the whole problem but then OpenMP threads handle the on-node parallelism, akin to how a traditional neutronic code would handle it.

I agree, this feels like probably the best approach to take.

I think its beneficial to have one geometry all at high resolution, since I think with the above change, we should be able to handle it.

I think for AEGIS currently only doing wall->plasma tracking I can make an argument for multiple geometries but you're right, changing parallelism strategy will dramatically reduce the memory usage, so there isn't really a reason against a single high resolution mesh.

• Gives us a coherent single source of truth for the geometry.
• It makes it easier for a user to setup their geometry, not having to produce multiple geometries and multiple meshes. And in fact ideally they could use the same tet mesh which is used for other analyses.
• If we want to implement a forwards (plasma->wall) tracking we wont need this concept of a "Target" and "Shadow" geometry and instead need to keep track of which elements are hit across the entire mesh so a single high resolution mesh would indeed be better.