{'fraction_non_alphanumeric': 0.022014051522248244, 'fraction_numerical': 0.00817720530835285, 'mean_word_length': 5.046135693215339, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hybrid Parallelization of Euler-Lagrange Simulations Based on MPI-3 Shared Memory

The use of Euler-Lagrange methods on unstructured grids extends their application area to more versatile setups. However, the lack of a regular topology limits the scalability of distributed parallel methods, especially for routines that perform a physical search in space. One of the most prominent slowdowns is the search for halo elements in physical space for the purpose of runtime communication avoidance. In this work, we present a new communication-free halo element search algorithm utilizing the MPI-3 shared memory model. This novel method eliminates the severe performance bottleneck of many-to-many communication during initialization compared to the distributed parallelization approach and extends the possible applications beyond those achievable with the previous approach. Building on these data structures, we then present methods for efficient particle emission, scalable deposition schemes for particle-field coupling, and latency hiding approaches. The scaling performance of the proposed algorithms is validated through plasma dynamics simulations of an open-source framework on a massively parallel system, demonstrating an efficiency of up to 80% on 131000 cores.


Introduction
The initialization time, i.e., the time from beginning of the code execution until the first computation step, plays a critical role in Euler-Lagrangian solvers in a high-performance computing context as it is closely linked with adequate load balancing. Ideally, each processor should receive equal load to achieve maximum overall simulation efficiency.
Accurate load estimation is already challenging in the case of a pure Euler solver as various cell sizes, time steps, local models and boundary conditions must be considered. Nonetheless, there are effective techniques to determine the local load and thus the grid distribution a priori, also known as static load balancing [1,2,3]. For these cases, a prolonged initialization period is acceptable if it results in improved runtime performance.
The presence of the Lagrangian phase adds substantial complexity as discrete particles introduce additional load, which is only weakly correlated with the local element sizes. Furthermore, particle concentrations may shift during the simulation, with high fluid and particle loads often occurring at the same mesh location, especially in fluid simulations. Load balancing approaches in Euler-Lagrangian solvers must adapt to these changes during runtime, which is referred to as dynamic load balancing. Over time, various load distribution strategies have evolved, which can generally be classified into two categories: 1) task parallelization and 2) domain partitioning. Task parallelization splits the work along the phase interface, distributing the fluid work and the particle work to different processors. The advantage is that due to the same nature of work within a phase, both groups of processors can internally subdivide the overall task in an optimal way. Implementations of this approach have been presented e.g., in Refs. [4,5,6].
The downside of task parallelization is the loss of any local connectivity, resulting in large communication effort, great memory requirements or both, rendering it inadequate for massively parallel computation where memory and interconnect bandwidth are scant [7]. Domain decomposition keeps the locality between the two phases intact but requires the load distribution to be performed on the combined work with examples of this approach published in Refs. [8,9,10,11]. Additionally, communication patterns become unpredictable as a processor may receive elements of discrete phase but not necessarily transmit them and vice versa [12]. Nevertheless, most massively parallel codes use domain decomposition, an approach we also follow.
As the focus of this work is on the discrete phase, an efficient solver for the continuous phase is presumed. Modern CFD solutions require high scalar performance and preferably minimal communication as storage and communication resources are unable to keep up with the steadily increase of available computing power [13,14]. High-order codes based on the Discontinuous Galerkin Spectral Element Method (DGSEM) have emerged as a well-suited approach as the fluid phase requires only the exchange of flux information at an element face, leading to a highly efficient numerical scheme while at the same time having dense, local operations [15]. Implementation as a solver for unstructured grids with possibly curved elements facilitates the creation of body-fitted domains even for complex geometries while retaining the high-order accuracy [16]. The presence of the Lagrangian phase necessitates dynamic load balancing, which should be regularly performed as particle loads can heavily shift during the simulation. Hereby, time spent on the load balancing step has to be kept minimal in order not to counteract benefits in application performance [17].
However, this task is non-trivial. While the DGSEM leads to a highly local scheme for the continuous phase, following the distributed memory approach means that each processor contains only local information on both the solution, associated quantities and the mesh without ready access to adjacent grid information in the case of distributed I/O.
The absence of this mesh information prevents full tracking of a particle in the event of it crossing a partition boundary if not remedied.
One method that allows the tracking of particles across different partitions is based on the idea of a shared layer of elements surrounding each partition, an idea akin to the ghost cell approach in e.g., finite volume methods. These halo regions contain the geometric information of neighboring cells within a given distance, here referred to as the halo distance, from the local domain and enable the completion of particle tracking on the initial processor, thus delaying the need for communication and requiring only the exchange of particles after having crossed the domain boundary [18]. While the process of halo element identification is straightforward for structured grids, it becomes significantly more complex for unstructured approaches as given within the present framework. Here, the search must be performed in physical space as the grid cells within a spatial region are not trivially mapped to locations in the mesh file [19].
Performing the search exclusively on processor-local mesh information, i.e., an inward search from the processor MPI border, leads to a severe performance bottleneck as the grid information for each cell residing on a single processor needs to be sent to a multitude of other processors within the halo distance, requiring many-to-many or in the worst case all-to-all point communication. This congests the interconnect infrastructure with the potential to stall code execution for minutes or even hours. As the information which needs to be exchanged grows both with grid size and number of processors involved in the simulation, the issue only becomes more urgent on modern massively parallel architectures where memory and communication constraints are all the more apparent. Moreover, as runtime load balancing is generally desirable to counter the shifting of particle loads during the simulation, the identification of the halo region must be performed multiple times during a given simulation, thus prompting the need for an efficient scheme which avoids detrimental effects on overall simulation performance. While this constitutes a new challenge for unstructured approaches that is absent in structured grids, this is easily outweighed by the advantages of such unstructured approaches for body fitted grids in domains of practical relevance.
Towards this goal, we present in this work a novel approach to unstructured Euler-Lagrangian simulations based on MPI-3 shared memory. Within this approach, we store information with compute node granularity and perform a multi-step communication-free parallel search on the compute node to identify the elements in the halo region, thereby retaining excellent scaling properties on today's massively parallel supercomputing architectures. Building on these data structures, we present methods for efficient particle emission, scalable deposition schemes for particlefield coupling, and latency hiding approaches. These efforts then give us the ability to conduct high-order simulations of rarefied gas flows at an industrial scale beyond generic test cases. To the best of our knowledge, this is the first unstructured framework enabling massively-parallel Euler-Lagrangian simulations at this problem size. The implementation considered in the present work is open-source and available on GitHub 1,2 .
The outline of this paper is as follows: The governing equations for non-equilibrium gas flows followed by the DGSEM scheme as well as the theory for particle motion and tracking are given in section 2. A high-level overview on the parallelization strategy is given in section 3. In section 4, we present the shared-memory approach for halo region determination and distribution, thereby shifting the communication load from the processor to the compute node level and alleviating the aforementioned scaling restrictions. The remainder of this section is designated to methods for 1 https://github.com/flexi-framework/flexi 2 https://github.com/piclas-framework/piclas emission, deposition on latency hiding based on the same shared-memory approach. The test case of an adiabatic box representing the optimum for the parallelization concept as well as near-application cases of a supersonic flow and a gyrotron resonator are discussed in section 5, and the scaling results are presented in section section 6. We conclude with a brief summary and give an outlook on further developments in section 7.

Theory
The approaches presented here are applicable to any unstructured Euler-Lagrange code and are already implemented in the two high-order open-source frameworks FLEXI 3 [20,21] and PICLas [18,22]. Both are actively developed at University of Stuttgart and share a common code basis for the DGSEM solver with FLEXI solving a continuous fluid phase prescribed by the compressible Navier-Stokes-Fourier equations while PICLas uses Maxwell's equations for the electromagnetic fields. FLEXI currently focuses on inertial particles in turbomachinery applications [23] with initial performance for the current methods presented in [24]. The present work focuses on solutions to non-equilibrium gas and plasma flows within the PICLas framework using the Particle-In-Cell (PIC) approach [25,26] as well as a particle based Bhatnagar-Gross-Krook (BGK) solver [27,28].

Non-equilibrium gas flows
Non-equilibrium gas and plasma flows are generally characterized by possibly charged particles that interact with an electromagnetic field where the statistical distribution is described by Boltzmann's equation [29,25] Here, f = f (x, v, t) represents the probability distribution function, i.e., the expected particle density at position x with velocity v. In the methods used here, the distribution function is approximated by particles. These can move freely and represent the phase space from a Lagrangian point of view through their location and velocity. The particles interact with each other through the right hand side described by a collision operator whereas the interaction for charged particles occurs through the Lorenz force F = F L . The Lorentz forces are calculated from electromagnetic fields, which are solved on a fixed grid in an Eulerian fashion. Typically, the Boltzmann collision integral is used as collision operator which effectively gives the change of the particle probability density function caused by binary particle-particle collisions.
However, since the Boltzmann collision integral is numerically difficult and time-consuming for various reasons [29], it is often approximated, e.g., by the Fokker-Planck (FP) solution algorithm [30,31] or the BGK approximation [27,28,32], see section 2.3.2. Beyond particle-particle collisions, charged particles experience the Lorentz force 3 https://www.flexi-project.org Here, q is the electric charge of a given particle whereas the electric field E and magnetic field B obey Maxwell's equations [33] ∂D ∂t with D being the electric displacement field, H the magnetic field strength whereas ρ and j represent the charge and current density, respectively. The field equations from Eq. (6) are thereby solved on a fixed grid, i.e. the Eulerian view. The coupling between the Eulerian view and the Lagrangian view arises on the one hand through the Lorentz force Eq. (2). Here, the fields are interpolated from the solution on the Euler grid to the Lagrangian particles in order to calculate the forces acting on these particles. On the other hand, the charge densities and current densities as source terms of the Maxwell equations on the Euler grid correspond to the zeroth and first moment of the distribution function, i.e., they are obtained by interpolating the particle data to the fixed Euler grid.

Discontinuous Galerkin Spectral Element Method (DGSEM)
In order to enforce charge conservation, Maxwell's equations are cast into the purely hyperbolic Maxwell (PHM) form [34] which can then be solved using the high-order Discontinuous Galerkin Spectral Element Method (DGSEM) [35,36]. DG methods operate on a weak formulation of the conservation equations following the method of lines approach by projecting them onto a space of polynomial test functions in reference space. Collocation of interpolation and integration points yields a highly efficient scheme which is advanced by an explicit Runge-Kutta scheme in time. PICLas is designed as a solver for unstructured grids, thereby allowing the straightforward creation of body-fitted grids with curved boundaries even for complex geometries. The approximation of the solution by a high order polynomial in each element ensures that high-order accuracy is retained [16].

Particle Behavior
Particles contribute to the electrical and magnetic field (deposition, see section 4.5) but are simultaneously influenced by the field through the Lorentz force and optionally through particle-particle collisions.

Equation of Motion
Based on eq. (2), the change in position and velocity of each particle is given by the relativistic equation of motion, with γ being the Lorentz factor given as Here, x and v represent the position and velocity of a given particle in physical space, q and m its charge and mass, and c the speed of light.

Bhatnagar-Gross-Krook (BGK)
The BGK operator approximates the collision term in Eq. (1) to a simple relaxation form where the distribution function relaxes towards a target distribution function f t with a certain relaxation frequency ν: The original BGK model assumes that the target velocity distribution function is the Maxwellian velocity distribution with the particle density n, particle mass m, temperature T and the thermal particle velocity c = v − u from the particle velocity v and the average flow velocity u [27]. In order to obtain the correct Prandtl number in the flow, more complex target distribution functions must be used. A very frequently used target distribution function that is also applied here is the ellipsoidal statistical BGK target function, see details in Holway Jr [37]. In PICLas, the solution of the BGK equation is performed in a pure Lagrangian and stochastic manner with particles using the stochastic particle Bhatnagar-Gross-Krook (SP-BGK) method as described in [38,39]. It offers an efficient approach to simulate non-equilibrium flows in smaller Knudsen number regimes. All particles within a cell interact with each other by the relaxation process described in Eq. (10). This is similar to the collision process in the well-known Direct Simulation Monte Carlo (DSMC) method [29], which also happens only within a cell. Subsequently, the convection due to the free movement of the particles is modeled together with boundary conditions in the computational domain in order to perform a relaxation process again.

Localization and Tracking
Particles are tracked in physical space for all cases considered within this work. The first step within this approach is the solution of eqs. (7) and (8) to obtain the new particle position. In order to identify a particle crossing an element boundary and its recipient, all faces of the previous element are checked for intersections with the particle path, see fig. 1. This procedure is performed iteratively until no more intersections are found and thus the final element is determined. If an element face is representing a boundary face, the corresponding boundary conditions can intuitively be incorporated by adjusting the remaining particle path. For more details and an alternative tracking approach based on localization in the reference space, see [19].

Parallelization Strategy
High-performance clusters are almost exclusively constructed as distributed systems, connecting separate nodes through an interconnect. In general, any interconnect imposes bandwidth starvation compared to local memory while simultaneously incurring latency costs. Thus, efficient parallelization approaches need to employ two strategies: 1) Communication avoidance and 2) latency hiding. By reducing the amount of transferred data, congestion on the interconnect can be alleviated. Performing the communication in a non-blocking manner allows local work to continue, thereby obfuscating the additional latency of the interconnect.

Continuous Phase
Relying on the DGSEM allows PICLas to make extensive use of both strategies. By enforcing a basis with local support, the volume integral becomes a purely local operation and only the surface flux information has to be exchanged on the cell boundaries. To ensure fast initialization times, the unstructured fluid elements are pre-sorted along a space-filling curve (SFC) during mesh generation. SFC have shown their suitability for efficient calculation of new distributions during runtime for load balancing purposes in PIC simulations [40,41]. Furthermore, the SFC allows for highly parallel, non-overlapping disk storage access with an arbitrary number of processors [42]. Details on the implementation are given in [21].

Halo Region
Halo regions bring these two strategies from continuous Eulerian to discrete Lagrangian phase. Since particle tracking is performed in unstructured physical space, geometric information along the considered particle path must be available at the time of tracking. By enriching the local DG domain with geometric information up to given physical distance from the domain boundaries, each processor can complete tracking to the final particle position. This halo distance is chosen as the maximum possible distance any particle can travel within a simulation time step, thereby ensuring a processor has all eligible elements accessible while simultaneously generating the minimal number of halo elements. Using this approach, communication is delayed until a particle is found to have left the processor domain after accounting for boundary conditions and only the minimum required information, i.e., the particle properties including the new particle position, must be communicated.
However, this shifts some work from the simulation time-stepping to the routines where geometric and neighboring information must be established or updated. This corresponds to the initialization and any load balancing step, which always also includes an update of this information. PICLas follows the commonly used restart-based load-balancing approach where the simulation is saved to disk and reloaded using an improved load distribution. As a result, this procedure relies on fast initialization times and is aided by the approaches outlined in the previous section. However, given the requirement to work on unstructured meshes, the required halo elements for the particulate phase can only be determined through a search in physical space even with the mesh elements already pre-sorted along the SFC.
As was previously outlined, performing the search and subsequent communication of mesh information using only processor-local information incurs severe performance penalties stemming from the differences in processor work from load distribution and the required many-to-many communication. The latter case in particular is exacerbated by modern many-core architectures, which consequently limits the scaling of the approach and necessitates the novel approach outlined in the subsequent section 4.

Load Balancing
Immediate benefit of the improved performance for the restart-based load balancing is the ability to increase the number of load evaluation and -if necessary -balancing steps. Following the classification by Watts and Taylor [43], the load evaluation can be based on the application, i.e., a priori using information on the algorithms involved; the system, i.e., at runtime using timing information, or a combination of both. While the application-based approach is simple to implement and successfully used for single-phase simulations with constant computation and communication time per degree of freedom [21], the determination of the correct weights becomes challenging for multi-phase flows. Hence, a hybrid approach is commonly considered more robust.
The approach utilized by PICLas relies on runtime measurements of the field solver and particle solver. Load evaluation steps are performed by comparing the runtime per rank at user-defined intervals. For this, total time spent in the field and particle solver are recorded prior to the load balancing step. Particles are assumed to remain in their element sufficiently long to assign their load to the element they are currently residing in. Thus, the load of a given element is estimated by the combination of the total time spent in field and particle solver, divided by the local number of elements and the element's contribution to the total tracking steps, respectively. The resulting time for a respective element is then given as t elem,tot = t field + δt particle · n particles,elem .
If an imbalance exceeding an acceptable threshold is detected, each rank gets assigned a new range of elements such that the load deviation becomes minimal, i.e., where i start and i end correspond to the respective element indices along the space-filling curve. More details on this implementation can be found in [44].

Implementation
This section describes the parallel implementation of the halo region search, emission and runtime deposition mechanisms. Following the MPI-3 shared memory paradigm, we store information with compute node granularity and perform subsequent routines in a communication-free way on the shared memory region. This chapter serves to illustrate the allocation of the shared memory window, the mesh distribution and aforementioned routines in section 4.1, section 4.2 and sections 4.3 to 4.5, respectively. Where applicable, we print the Fortran source code rather than pseudo-code to facilitate implementation in other scientific frameworks. Further information on how the stored mesh information is subsequently used for particle tracking is given in [19].

Shared Memory Allocation
Classical shared memory programming involves OpenMP. However, this approach is limited to single-node cases as OpenMP cannot handle distributed memory. Message Passing Interface (MPI)-3 introduces the concept of shared memory with the MPI Shared Memory (SHM) model. The resulting coding approach is also called "hybrid parallel programming" as it combines the shared memory approach of OpenMP with the distributed memory view of MPI. Memory regions allocated with MPI-3 SHM can be distributed arbitrarily between the processors while being accessible by all processors on the compute node. In our implementation for cache-coherent systems (MPI_WIN_UNIFIED), the shared memory window is allocated only on the compute node root to avoid offset calculations. The shared window is continuous in memory and thus can be directly read by each process. We do not employ RMA routines for store operations but ensure non-overlapping writes through data distribution along the SFC. MPI_WIN_SYNC calls ensure explicit synchronization and immediate availability of the written information on the compute node. The code for this approach is given in listing 1.

Mesh Distribution
High-order meshes are created with the in-house preprocessor HOPR 4 [45]. The mesh elements are ordered along a space-filling curve and saved in binary HDF5 format for highly parallel access [21], together with likewise ordered side connectivity information and the grid coordinates. PICLas initially determines the number of elements per processor taking available load balancing information into account. Each processor only accesses mesh information for its region along the SFC and stores it in processor-local memory, thereby maximizing file system parallelism.
However, each compute node additionally allocates shared memory sufficient to hold the raw mesh information (as stored in the HDF5 file) and save its compute node information at the correct offset. Once every processor has finished the reading of the mesh, the compute node root processors perform an non-blocking IALLGATHERV operation on the interconnect, making use of available hardware offload capabilities. A graphical representation of this procedure is shown in fig. 2.
In addition to the unstructured computation grid, a Cartesian background mesh is created during runtime in order to reduce the eligible computational elements when performing particle localization procedures. The number of Cartesian elements in each direction is case-dependent and currently determined by the user. During runtime, particle intersection calculations are performed using either the physical element face corner coordinates in the case of a purely linear mesh or clipped Bézier surfaces for a curved grid. Additional particle information includes e.g., the face normal vectors, the distance of an element to the nearest boundary, the surrounding mesh node indices, and the tolerance of a curved element in reference space. Some of this information is only calculated given a specific tracking method and whether deposition is desired. Since part of this information varies depending on the mesh distribution and storing it would also inflate the size of the mesh file, these particle metrics are computed during the initialization phase.

Halo Element Search Algorithm
In order to minimize both computational effort and memory requirements, a two-step search algorithm to determine eligible halo elements is performed before calculating the particle mesh metrics. This approach not only alleviates computational effort during the initialization phase but also significantly reduces the memory footprint as the derived metrics are only stored for the local and actual halo elements. For the first step depicted in fig. 3a, a Cartesian bounding box around all mesh elements local to a compute node is calculated. This bounding box is then extended by the halo distance in each direction.Through projection of the bounding box onto the Cartesian background mesh (BGM), the corresponding limits for the required I,J,K indices are obtained.
Next, a similar Cartesian bounding box is created for each mesh element. This bounding box is again projected onto the background mesh, thus creating a mapping from each mesh element to the overlapping BGM cells. The BGM mapping of every element not located on the compute node is then compared against the BGM region previously extended by the halo distance. All elements whose BGM bounding boxes overlap with the extended bounding box are flagged as potential halo elements. Since the global mesh information is available in the MPI-3 shared memory array, this calculation is distributed among all compute node processors through slicing of the space-filling curve.
At this point, the potential halo elements need to be further reduced as the Cartesian bounding box arbitrarily extends beyond the compute node local mesh elements. However, comparing the distance of every potential halo elements again all available compute node local elements would measurably affect the initialization time. Thus, the number of elements to compare against is reduced by only considering the elements having at least one MPI boundary on the compute node circumference, therefore representing the MPI boundary of the compute node local mesh.
The second step shown in fig. 3b then calculates the radius of the convex hull of the elements on the MPI boundary and the radius of the potential halo elements. Each potential halo element is compared against each MPI-border element by subtracting the distance of the two barycenters from the sum of the two radii plus the halo distance. If the result is negative, the element is flagged as confirmed halo element. A positive value indicates that particles cannot reach the element within a given time increment and the element is subsequently discarded. This process is again distributed among all compute node processors through a uniform partitioning of the potential halo elements.  Once every compute node processor indicates that there are no more potential halo elements to check, a mapping containing first the compute node local elements followed by the compute node halo elements is built. This mapping allows for efficient looping when building the derived particle mesh metrics on the reduced mesh.

Emission
Since particles are tracked in physical space, this approach naturally has to extend to the particle emission as well.
In order to maintain good scalability, emission is performed in parallel with each processor calculating the initial Local element Halo element Periodic element Figure 4: Local, halo and periodic element within halo distance for the adiabatic box, see section 5.1.
particle positions within the complete emission region. Subsequently, the grid element corresponding to each position has to be identified and the particle sent to the respective processor. However, the emission region might extend beyond the compute node-local mesh even when including the halo region, meaning that the emitting processor cannot uniquely identify the target element. Worse yet, without knowledge of the elements outside the halo region, the emitting processor cannot distinguish between a valid position and one outside of the complete mesh. Yet, it is equally undesirable to retain all elements within the emission region.
The solution to the problem is again provided through the BGM. As mentioned in the previous paragraph, there exists a mapping from each BGM element to the overlapping grid elements located within the compute node and halo region. During initialization, each processor additionally provides the number of local elements per BGM cell. This number is summed up across all processors and stored with a mapping containing all processors which overlap with a given BGM cell. As the local elements are inherently distributed, this process is by design automatically scaling.
During the particle emission, after an emission position is computed, the processor calculates the associated BGM cell. Next, we compare the number of compute node grid elements mapped to this BGM cell with the total number of grid elements connected to this cell. If the numbers match, the processor flags the position as locally computable.
The other positions are gathered and sent to all processors associated with the BGM cell. Next, we perform the search algorithm of the locally computable positions, thereby acting as latency hiding. After identifying all local particle to element mappings, the search of communicated positions is performed on each processor. Since any position can only correspond to one single element, no further communication is required. All other processors silently discard the position.

Deposition with Shape Function
The charged particles are responsible for the source terms of the field equations and therefore are coupled with the underlying grid on which the field equations are solved. The source terms themselves are determined from the respective moments as described by the distribution function. Here, this is achieved by mapping the particle position and velocity to the grid via shape functions that smoothly distribute the charge and current densities of the particles on the grid, which is referred to as deposition. The cut-off radius of the deposition is determined by the physical problem and can range across multiple elements of the grid as shown in fig. 5, hence, deposition may occur in processorlocal elements as well as elements that are of different processors or even different nodes. Thus, processors require the communication of either the deposited properties or the particle properties, which are in turn deposited by the receiving processor. Since the presented parallelization concept allows elements to be uniquely globally identified, the particles are deposited by the host process of the particle. Source terms that are possibly deposited in other processors are stored in a separate array for communication. Subsequently, a message is created from this array for all processors in whose elements were deposited and communicated regardless of which node they are located on. In order to communicate only with processors that can potentially exchange source terms, a list of all reachable processors in the halo region is initially created using the shape function radius for each processor. In order to avoid multiple communications between all processors in this initial process, this is done in a two-step communication procedure.
In a first step, the list with the processors to be communicated is sent to the MPI roots of the compute node for each processor. These node leaders on each compute node gather the information for distribution to other processors from all corresponding processors and store them in a shared array so that the information about the necessary exchange processors is available to all processors.

Latency Hiding
The main goal of latency hiding is to allow communication and computation to overlap completely, which means that during the time of communication, parts of the algorithm are already being carried out, so that there is no waiting time during communication. There are two basic problems with latency hiding for PICLas. PICLas is a modular toolbox, so e.g., PIC can be used as a module without or with collision term, within PIC there is the distinction whether electromagnetic or electrostatic simulations should be carried out or which interpolation method between particles and grid should be used. On the collision term side, there is also a wide choice of methods such as DSMC, BGK, FP and others. Each possible combination of modules with different types of time integration has different requirements as to which data must be available and when. Obviously, it is therefore not possible to use a latency hiding method that represents the optimum for all possible module combinations. Therefore, in the following we will only concentrate on the methods used in this work and already described in the theory section. The other problem is that many particle methods such as DSMC, BGK or FP have very sequential structures, so that only a few parts can already be calculated if the information from particles from other cores is not yet available. Most of the time, calculations that are carried out in parallel with communication require a large amount of additional memory, because quantities that would otherwise only be needed locally in the cell have to be stored for all cells so that they are available after communication.

Maxwell-PIC with Shape Functions
Maxwell-PIC within this context refers to the electrodynamic PIC method in which the complete set of Maxwell's equations are solved as described in section 2.1. The time integration scheme used in this example is a 5-stage 4 thorder low-storage Runge-Kutta (RK) method [46]. For simulations of this type, there are three MPI communications per time step which should require latency masking, i.e., the communication of the flux data for the Maxwell-DG solver, the current densities and charge densities deposited on the grid as source terms for the DG solver by the shape functions and the particle data that leave the processor after the movement. In general, the compute time of the discrete phase is dominating compared to the continuous phase, hence the focus of latency hiding is to avoid stalling of the particle routines.
A flow chart of one RK stage is given in fig. 6. Performing the deposition as outlined in section 4.5 allows to hide the costly exchange of volume data behind the particle operators. Interpolation, calculation of the Lorenz forces as well as particle tracking are purely local operations. At their end, particles are already assigned to their final processor, so the communication of particle data can start. As the continuous phase requires two communication steps, we start by extrapolating the field data to the element faces and communicate the surfaces data. This corresponds to (N−)D information and thus requires considerably less interconnect time than the previous field data. As was already outlined, the DG volume integral is a purely local operation and is thus performed on the first half of the local elements to hide the communication latency. After receiving the surface data, the numerical flux on the MPI sides is calculated and immediately sent back. Local operations on the inner sides as well as the remaining volume integrals can commence.
Once these routines return, the numerical flux on the MPI sides should be received and the surface integral on the MPI sides is calculated. The RK stage concludes by incorporating the particle data into the local arrays.

SP-BGK
The SP-BGK method requires the allocation of the particles to the cells in order to calculate the moments of the distribution function per cell. However, this is a purely discrete phase method, thus no field solver part or deposition Receive ρ,j

Send Particle Data
Receive Particle Data to hide the communication time of the particle communication as in the PIC method. Therefore, latency hiding is implemented here in two different approaches. In the first approach, the elements of each processor are divided into exchange elements in which particles from other processors can potentially move within a time step and purely local elements in which this is not possible. For this purpose, the halo region is extended from the MPI boundaries into each processor's own computational domain and the elements are flagged accordingly as shown in fig. 7. With this information, the BGK operator can already be applied to the purely local elements during particle communication, as this operator is cell local. After all particles have been received, the BGK operator is then applied to the exchange elements. This type of latency hiding is only effective if each processor has enough elements so that a minimum amount of purely local elements exists. However, this is often not the case for very high numbers of processors.
Therefore, the assignment of the particles to the elements is additionally divided. All particles that do not leave the processor's slice of the computational domain are already assigned to the elements during the particle communication.
Within the elements, an adaptive octree is used to create subcells in order to better capture any gradients [47,18]. As far as possible, parts of this assignment are also carried out during the communication of the particles. Subsequently, after the particle data is received, the element assignment is done for the received particles. This second step of latency hiding has the advantage that it can always be performed regardless of the number of purely local elements.

MPI boundary
Purely local elements Potential exchange elements

Test Cases
Three test cases were selected to investigate the efficiency of the proposed methods. The first test case consists of a generic setup representing a weak scaling case of an adiabatic periodic Cartesian box to assess specific performance metrics in a regular and adaptable setup. The second and third simulation showcase the strong scaling of practical application setups with fixed sizes on unstructured grids. This section presents the test cases themselves with the scaling results given in section 6.
All simulations are performed on the HPE

Adiabatic Box
The first test is an adiabatic box which represents the optimum for the parallelization. The domain itself is a fullyperiodic 3D box into which methane with a particle density of n = 1 · 10 23 1/m 3 is homogeneously inserted. The start temperature is T ∞ = 2000 K, resulting in a required time step of ∆t = 3 · 10 −10 s [28]. The simulation is carried out using the stochastic particle Bhatnagar-Gross-Krook (SP-BGK) method as described in section 2.3.2, with the periodic structure of the domain ensuring that the number of particles remains homogeneous in the entire computational domain and thus the computational load for all processors is equally distributed over the entire computation time.
Nevertheless, of course, the particles move at each time step and the collision operator is also executed. As

Supersonic Flow around 70°Sphere-Cone
The second test case represents the near-application setup of a supersonic flow. The geometry, dimensions and the test case itself are adapted from the paper by Hollis et al. [49]. The test vehicle is a 70°sphere-cone blunt-body in an high-enthalpy carbon-dioxide flow. The inflow at an angle of attack has a temperature of T ∞ = 126 K, a velocity of u ∞ = 2030 m/s and a density of ρ ∞ = 5.9 · 10 −3 kg/m 3 resulting in a Mach number of M ∞ = 11.4. The flow is again simulated with the SP-BGK method.
The 3D grid with a total of 1 638 637 hexahedral cells is shown in fig. 8. The particle number in the simulation is 1.25 · 10 8 , where the time step to resolve the stiff BGK relaxation term (see Pfeiffer [28]) is chosen to be ∆t = 3 · 10 −10 s. The average number of communication partners for each core ranges between 15 and 21, increasing with the number of cores and a total number of 1000 time steps is carried out during the test.

140 GHz Gyrotron Resonator
The third test case is a gyrotron resonator operating at 140 GHz. The details of the setup are found in [50], an adaption of the original setup found in [51]. The geometry resembles a tapered hollow cylinder, the diameter of which increases along the symmetry axis as depicted in fig. 10.  The resulting charge density within the domain that is created by the electron hollow beam is shown in fig. 10.
Only the elements that show a charge density contain simulation particles. Therefore, the majority of the simulation domain is empty, leading to a strong workload imbalance between empty elements and those containing simulation particles. This imbalance is addressed by the timer-based dynamic load balancing as described in section 3.3, which assigns weights to each element and partitions the complete domain into segments such that each segment has the same computational load. The average number of communication partners for each core ranges between 19 and 80, increasing with the number of cores and a total number of 522 time steps is carried out during the test.

Results
In the following, this paper will not interpret the results physically, but will examine the problematics of such flows for parallelization. First, the scaling of the initialization phase, i.e., the construction of the halo region and communicators, is examined. This is followed by an evaluation of the calculation phase, considering aspects such as field/particle operator and load imbalance. The total number of elements and particles as well as the initialization and execution times (without I/O and initialization) are summarized for all test cases in table 1.

Initialization
Within the context of this paper, the initialization comprises the complete code startup including the initial insertion of the particles or the corresponding restart routines in the case of a continued simulation. The initialization    6 50 · 10 6 100 · 10 6 200 · 10 6 400 · 10 6 800 · 10 6 1600 · 10 6 3200 · 10 6 6400 · 10 6 12 800 · 10 6   One exception is the run on a single node for both setups. For this case the initialization time of the gyrotron is shorter than with two nodes. The advantageous effect is that in the case of a single node, all grid cells are automatically compute node local elements and thus no halo region needs to be constructed, thus saving initialization time. Another point is visible for the supersonic test case. As the node count exceeds 16, the peculiarity of the test system becomes apparent. Since 16 nodes are directly connected to a switch, each doubling of the node numbers introduces an additional hop, leading to an increase in the initialization time as communication calls now have reduced bandwidth available. This yields an increased variation in initialization and calculation times and therefore, the minimum and maximum as well as the average values are depicted.

Simulation Performance
Weak Scaling. The parallel efficiency η N for the weak scaling was determined by where t 128 and t N are the computational time using 128 and N cores, respectively. Since the problem size increases linearly with the number of cores used, the parallel efficiency should ideally remain around one.
Already with low core counts, weak scaling of the adiabatic box in fig. 12 shows a decreasing efficiency with an increasing number of cores. As the setup is inherently ideally load balanced, the increase in computing time is assumed to be related to an increasing time for particle communication via the nodes. However, the parallel efficiency with 32 nodes (4096 cores) is still around 0.9, which is satisfactory. Reaching 64 nodes, the influence of the reduced bandwidth becomes visible due to the higher particle density in the domain. For even higher node numbers, the efficiency increases as the weak scaling results benefit from communication locality, thereby compensating for the interconnect penalty.  Strong Scaling. The strong scaling was calculated using the basis of one node corresponding to 128 cores by with the respective parallel efficiency η N determined by where t 128 and t N are the computational time using 128 and N cores, respectively. The strong scaling results of the supersonic flow and gyrotron test case are depicted in figs. 13a and 14a with the parallel efficiency shown in figs. 13b and 14b, respectively.
Due to the load balancing as described in section 3.3, whereby the grid cells are weighted by the number of particles they currently contain when they are divided among the processors and the latency hiding as described in section 4.6.2, the parallel efficiency of the supersonic flow test case remains around one up to 16 nodes, which corresponds to 2048 processors. For larger processor numbers, the parallel efficiency then drops due to the large    discrepancies in computational load between different regions of the computational domain. Some elements contain so many particles that they can no longer be better distributed by the load balancing and thus also the latency hiding can no longer work properly. In this case, it would be necessary to increase the number of mesh elements in the regions of high particle density, thereby allowing the computational load to be better distributed. Nevertheless, even with 64 nodes, which corresponds to 8192 processors, a parallel efficiency of 0.77 is still achieved for this case with a very nonuniform load distribution. This setup thus retains the better scaling properties than the adiabatic box, presumably because the lower particle number poses less requirements on the interconnect bandwidth. Nonetheless, it is still sensitive to an even further increase in latency as is seen when 64 nodes are exceeded.
Compared to the strong scaling of the supersonic flow case, the scaling behavior is distinctively worse in the gyrotron case. As visible in fig. 14b, the parallel efficiency drops significantly faster as the number of processors increases, so that with 64 nodes (8192 cores), a parallel efficiency of only 0.37 is achieved. The main reason for this is most likely the more complex load balancing for this case. On the one hand, the presence of the Maxwell solver introduces differing loads per cell between field solver and the particle solver. On the other hand, the number of particles per element differs by orders of magnitude, as shown in fig. 10b. Numerous cells contain no particles at all while only very few carry the particle beam shaped as a hollow cylinder. As a consequence, not only is an ideal load distribution increasingly difficult to achieve as the number of processors increases, but efficient latency hiding also becomes progressively unattainable. At any given time, there exists a considerable fraction of processors that partake in the particle communication without ever receiving particles in their corresponding cells. These processors thus stall the particle communication of the remaining cores, a phenomenon which cannot be hidden. Ultimately, in this test case, more and more cores have to wait for communication, resulting in a scaling which is not optimal.

Conclusions
As computers become increasingly parallel, code areas with previously negligible performance impact such as initialization become increasingly more relevant due to their influence on load balancing. In the work presented here, our aim was to contribute to this challenge by presenting a massively parallel, communication-free approach to build the halo region required for Euler-Lagrange simulations. The use of the MPI-3 shared memory model enabled us to utilize the ever increasing core count per socket without introducing additional interconnect load. Based on this programming model, we developed new methods for emission, deposition and latency hiding which were implemented in the open-source plasma dynamics framework PICLas. This framework was applied to a generic test setup as well as two practical application cases. In all setups, we were able to show respectable initialization times as long as no interconnect congestion occurred in other parts of the startup phase. Furthermore, we were able to retain good efficiency for both the weak and the strong scaling case of the BGK setup while also highlighting challenges inherent to the Euler-Lagrange setup including strong variations in particle density in the PIC setup.
In the future, we will extend our research to communication-minimizing decomposition approaches for the Euler-Lagrange codes. These pose additional challenges since the continuous and the disperse phase entail separate communication regions with machine-dependent costs, resulting in a multi-point optimization problem.