ExaParty always runs through the same sequence of setup steps prior to the actual calculation:

Create a regular mesh that is just as fine that it can accommodate the smallest predicted search radius of all particles.
Insert the particles (initial condition).
Run a set of initialisation steps, where particles can define some "warm up" calculation steps and are also sorted once more, should they have been inserted into the wrong cells.

The initial (regular) mesh

The first step is important, even when particles are allowed to change their search radius. It is important to tell the code base right from the start what the biggest and smallest radius might be, but the radius can change freely in-between those thresholds. Actually the radius might even become smaller, but the mesh width you pass in initially is used to constrain the AMR throughout the simulation. See initial_smoothing_lengths_and_search_radii for details on how to tweak the initial search radii.

The dynamic mesh refinement might decide to introduce an extremely fine mesh in some places and to erase the mesh somewhere else. However, we have to start somewhere, as we do not want to run the initial insertion of particles into a single cell, then create this initial grid (which is obviously guided by the particle distribution) and finally create the adaptive grid. Instead, the code creates some initial grid and load balances this one. After that, it inserts the particles (in parallel) and might still change both the load balancing and the mesh structure. But we have to start somewhere.

By default, the initial grid is a regular mesh controlled by the potential minimal particle search radius. As pointed out: this search radius can change over time and thus is solely an initial indicator. You can alter the initial mesh and hence start with an adaptive one right from the start. For this, you can take the step algorithm_step_create_grid which is an attribute of your exaparty.api.Project and add further action sets which refine your mesh according to your needs. However, the initial grid construction will try to establish a stationary grid, so ensure that you define a stationary grid, i.e. not a grid where you sometimes erase parts and sometimes add parts.

We have chosen the default initial grid to be fully refined so that there will be a sufficient number of cells to distribute among sub-trees to occupy many cores from the very beginning, i.e. without waiting for the adaptive refinement to kick in later. You can instead only create a coarse grid by passing the

        initial_grid_coarse=True

argument to

        Project.set_global_simulation_parameters()

in your python script. See src/exaparty/api/Project.py for more.

Altering the mesh size

To alter the initial mesh size, you have to tell each individual particle species that you are using what its maximum mesh size will be. Most codes set this quantity in the particle definition in the Python script, i.e. through the Python constructor via the

          min_h

and

          max_h

parameters in the particle constructor (see e.g. src/exaparty/api/particle/SPHParticle.py and exaparty.api.particle.SPHParticle()).

As Peano uses three-partitioning, the mesh size choices can seem to appear weird. Each Python script run produces a README-xxxx.md file which provides details on the arising mesh sizes. Alternatively, you can study the implementation in exaparty.api.Project.__compute_global_max_h().

Implementation

The ExaParty project object will export the arising global maximum mesh size into the Constants.h file as a simulation-global constant. It will also feed this quantity into the action set exaparty.api.Project.action_set_create_regular_grid:

action_set_create_regular_grid = peano4.toolbox.CreateRegularGrid(
    self._global_max_h
)

Here is what happens from hereon:

The global main file (cmp SWIFTMain.template.cpp) enters its phase 1. This means setting up the (initial) grid.
It then runs over the mesh in a while loop.
Per mesh traversal, the action set from action_set_create_regular_grid will be invoked and asked whether to decide the mesh further.
Once the global main finds out that the global maximum mesh size constraint is met, it sets its flag gridConstructed and consequently switches to phase 2 which is the actual particle insertion.

Load balancing

At this moment, we enable any load balancing only throughout the grid construction by default. That means: by default, we work with a static domain decomposition. As soon as we insert the particles, we disable any mesh load balancing. More sophisticated code variants switch on the load balancing every few iterations or enable it to rebalance it after each time step, but this is not enabled in the baseline version.

It makes sense to play around with different load balancing schemes and initial mesh construction strategies. Switching the load balancing scheme is done in the Python script vi the function exaparty.api.Project.set_load_balancing().

Playing around with the maximum particle radius

From the grid's perspective, there is no direct reason why you cannot pick the maximum particle radius relatively small, even though the particle radii eventually might become bigger. If you do this, it can happen that the original mesh is too fine and particles are inserted into a fine mesh where a search radius overlaps more than one neighbouring cell. In this case, the code will have to resort to a multiscale sorting approach, i.e. will have to move particles up within the grid hierarchy. So you work with an overly refined mesh and introduce some (initial) sorting overhead, but it might nevertheless be reasonable as it gives the load balancing more granularity.

Initial conditions

The initial conditions are added by one mesh sweep modelled by the project's attribute algorithm_step_initial_conditions. It represents the grid sweep, i.e. mesh sweep, which sets up the initial condition. That is, if you want to insert particles, you have to add action sets to this algorithm step calling add_action_set(). The passed action sets have to throw all the particles you want to have into your mesh. There are three mainstream ways to implement initial conditions for the simulations:

Insert particle by coordinates (src/exaparty/api/input/InsertParticlesByCoordinates.py).
Insert particle along a regular mesh (src/exaparty/api/input/InsertParticlesAlongCartesianGrid.py).
Insert particles from an HDF5 file (src/exaparty/api/input/InsertParticlesFromHDF5File.py).

Nothing stops you from additing further action sets to algorithm_step_initial_conditions. Some scenarios, for example, hard code the initial condition of further particles. The planet orbit benchmark is such an example, where no input files are required.

The particle inserters as discussed above typically only set a subset of the particle attributes. They are generic tools and thus cannot know which attributes your particle actually hosts and requires to be set. To discuss how to initialise the attributes, we have to distinguish two types of attributes:

There are attributes which are domain-specific and have to be set by the initial conditions to meaningful values.
There are attributes which have to be calculated prior to the kick off of the simulation.

To set the intial conditions, it is worth looking into the documentation of your particle inserter of choice. Most of them allow you to add some further initialisation snippets or tailor the reader, such that you can set the particle attributes of interest. Alternatively, you might want to create a subclass of a generic file reader that is specific to your setup and, for example, grabs additional particle attributes from the input file.

If you don't want to mess with the particle reader, nothing stops you from adding an additional action set to the initial condition sweep. This set can run over all particles within a vertex - once they are inserted - and set meaningful attribute values.

If you have to calculate particle attributes after the particles are thrown into the domain, these calculations are best modelled as steps within the particle lifecycle.

Insert particles from an HDF5 file

A basic HDF5 reader that is compatible with IC files from the original Swift code is implemented in InsertParticlesFromHDF5File.py. In order to make use of the HDF5 Initial conditions reader, Peano must have been configured with HDF5 enabled and linked. Instructions are provided here. Also make sure to install the Python hdf5 libraries, which may be required for python scripts which generate the initial conditions.

General remarks on particle insertion

There are a few key policies to follow whenever you insert particles:

Particles have to be created on the heap, i.e. via an explicit new
Every particle has to be initialised. You do this by calling toolbox::particles::init(). As you need need this helper routine, you have to include toolbox/particles/particles.h.
Particles have to be inserted into cells, and you should not try to do this on your own. Instead use the helper function from particles.h:

toolbox::particles::insertParticleIntoCell(
  marker,
  particle,
  fineGridVerticesParticleSet
);

Obviously, the names for the parameters for toolbox::particles::insertParticleIntoCell() depend on the location from where you call the routine. As you insert particles into cells (even though the code will then assign them to vertices adjacent to the cells), you plug into touchCellFirstTime() or touchCellLastTime(). I recommend the former. These two events always are given a parameter

   const peano4::datamanagement::CellMarker& marker

and you can use peano4::datamanagement::CellMarker::isContained() to check if a particle should be fitted into a cell.

Initial Smoothing Lengths and Search Radii

The particles' initial smoothing lengths and search radii need special attention because those quantities don't only play a role with the physics and hydrodynamics equations to be solved, but are also instrumental in the core algorithms and grid construction.

The particles' smoothing lengths, interaction radii, and search radii are related quantities. See the discussion in Terminology for details.

The initial maximal search radii of particles can be set indirectly through the maximal and minimal smoothing lengths, which are set by the h_hydro_max and h_hydro_min parameters of the particle species (see SPHParticle.set_parameters() in src/exaparty/api/particle/SPHLeapfrogFixedSearchRadius.py). Note that those parameters are min/max particle smoothing lengths, not search radii. To convert between smoothing lengths and search radii, look into exaparty.api.sphtools.sph_kernel_get_smoothing_length_from_search_radius() and exaparty.api.sphtools.sph_kernel_get_interaction_radius().

Regarding the actual smoothing lengths of particles, we require the user to provide an initial guess. The better the guess is, the fewer iterations we need to do to determine it in the initial setup. A guess of zero will lead to problems. Usually, the guess \( a \cdot dx \), where \(a>0\) is a factor of something a few and \(dx\) is the mean interparticle distance (e.g. boxsize / particles per dimension) is good enough.