Commit fd8a69e3 authored by Ann Almgren's avatar Ann Almgren
Browse files

1) Add to run-time inputs

2) Fix formatting of ParticlesOnGpus
parent 7b2de8f8
Loading
Loading
Loading
Loading
+19 −4
Changes for docs/source/Inputs.rst: 19 added lines, 4 removed lines.
Original line number Diff line number Diff line
@@ -6,7 +6,7 @@ Run-time Inputs
The following inputs must be preceded by "amr."

+-----------------+-----------------------------------------------------------------------+-------------+-----------+
| File            | Description                                                           |   Type      | Default   |
|                 | Description                                                           |   Type      | Default   |
+=================+=======================================================================+=============+===========+
| max_step        | Maximum number of time steps to take                                  |    Int      | None      |
+-----------------+-----------------------------------------------------------------------+-------------+-----------+
@@ -37,7 +37,7 @@ The following inputs must be preceded by "amr."
The following inputs must be preceded by "geometry."

+-----------------+-----------------------------------------------------------------------+-------------+-----------+
| File            | Description                                                           |   Type      | Default   |
|                 | Description                                                           |   Type      | Default   |
+=================+=======================================================================+=============+===========+
| coord_sys       | 0 for Cartesian                                                       |   Int       |   0       |
+-----------------+-----------------------------------------------------------------------+-------------+-----------+
@@ -48,14 +48,15 @@ The following inputs must be preceded by "geometry."
| prob_hi         | High corner of physical domain (physical not index space)             |   Reals     | None      |
+-----------------+-----------------------------------------------------------------------+-------------+-----------+


The following inputs must be preceded by "mfix."

+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| File               | Description                                                           |   Type      | Default   |           
|                    | Description                                                           |   Type      | Default   |           
+====================+=======================================================================+=============+===========+
| fixed_dt           | Should we use a fixed timestep?                                       |    Int      |   0       |
+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| dt_min             | Abort if dt gets smaller than this value                              |    Real     |  1.e-6    |
+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| dt_max             | Maximum value of dt if calculating with cfl                           |    Real     |  1.e14    |
+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| cfl                | CFL constraint (dt < cfl * dx / u) if fixed_dt not 1                  |    Real     |   0.5     |
@@ -64,5 +65,19 @@ The following inputs must be preceded by "mfix."
+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| explicit_diffusion | Should we use explicit or implicit diffusion?                         |   Int       |    1      | 
+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| eb_ho_dirichlet    | Should we use ray tracing to compute dphi/dn for Dirichlet EB faces?  |   Int       |    0      | 
+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| particle_init_type | How do we initialize the particles?   "Auto" vs AsciiFile             |   String    | AsciiFile |
+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| do_initial_proj    | Should we do the initial projection?                                  |    Bool     |  True     |
+--------------------+-----------------------------------------------------------------------+-------------+-----------+
| initial_iterations | How many pressure iterations before starting the first timestep       |  Int        | 3         |
+--------------------+-----------------------------------------------------------------------+-------------+-----------+

The following inputs must be preceded by "particles."

+--------------------+---------------------------------------------------------------------------+-------------+-----------+
|                    | Description                                                               |   Type      | Default   |           
+====================+===========================================================================+=============+===========+
| removeOutOfRange   |   Should we remove particles at initialization that are touching the wall |    Int      |   1       |
+--------------------+---------------------------------------------------------------------------+-------------+-----------+
+22 −27
Changes for docs/source/ParticlesOnGpus.rst: 22 added lines, 27 removed lines.
Original line number Diff line number Diff line
@@ -3,38 +3,35 @@ Particles on GPUs

The particle components of MFIX-Exa are a natural candidate for offloading to the GPU. 
The particle kernels are compute-intensive and can in principle be processed asynchronously with parts of the fluid advance.
Additionally, particle methods have been successfully offloaded using AMReX in the ECP code WarpX. 
The core components of the particle method in MFIX-Exa are:

\begin{enumerate}
    \item Neighbor List Construction
    \item Particle-Particle Collisions
    \item Particle-Wall Collisions
\end{enumerate}
- Neighbor List Construction
- Particle-Particle Collisions
- Particle-Wall Collisions

Of these operations, the neighbor list construction requires the most care. A neighbor list is a pre-computed list of all the neighbors a given particle can interact with over the next $n$ timesteps. Neighbor lists are usually constructed by binning the particles by an interaction distance, and then performing the $N^2$ distance check only on the particles in neighboring bins. In detail, the CPU version of the neighbor list algorithm is as follows:
Of these operations, the neighbor list construction requires the most care. 
A neighbor list is a pre-computed list of all the neighbors a given particle can interact with over the next *n* timesteps. 
Neighbor lists are usually constructed by binning the particles by an interaction distance, 
and then performing the N\ :sup:`2` distance check only on the particles in neighboring bins. In detail, the CPU version of the neighbor list algorithm is as follows:

\begin{enumerate}
    \item For each tile on each level, loop over the particles, identifying the bin it belongs to.
    \item Add the particle to a linked-list for the cell that ``owns`` it.
    \item For each cell, loop over all the particles, and then loop over all potential collisions partners in the neighboring cells.
    \item If a collision partner is close enough, add it to that particle's neighbor list.
\end{enumerate}
- For each tile on each level, loop over the particles, identifying the bin it belongs to.
- Add the particle to a linked-list for the cell that `owns` it.
- For each cell, loop over all the particles, and then loop over all potential collisions partners in the neighboring cells.
- If a collision partner is close enough, add it to that particle's neighbor list.

To port this algorithm to the GPU, we use the parallel algorithms library Thrust, distributed as part of the CUDA Toolkit. Thrust provides parallel sorting, searching, and prefix summing algorithms that are particularly useful in porting particle algorithms. To construct the neighbor list on the GPU, we follow the basic approach used by Canaba, a product of the Particle Co-Design Center:

\begin{enumerate}
    \item Sort the particles on each grid by bin, using a parallel counting sort. We use Thrust's `exclusive\_scan` function to implement the prefix sum phase of the sort, and hand-coded kernels for the rest. This step does not actually involving rearranging the particle data - rather, we compute a permutation that would put the particles in order without actually reordering them.
    \item Once the particles are sorted by bin, we can loop over the particles in neighboring bins. We make two passes over the particles. First, we launch a kernel to count the number of collision partners for each particle.
    \item Then, we sum these numbers and allocate space for our neighbor list.
    \item Finally, we make a another pass over the particles, putting them into to list at the appropriate place.
\end{enumerate}
- Sort the particles on each grid by bin, using a parallel counting sort. We use Thrust's `exclusive\_scan` function to implement the prefix sum phase of the sort, and hand-coded kernels for the rest. This step does not actually involving rearranging the particle data - rather, we compute a permutation that would put the particles in order without actually reordering them.
- Once the particles are sorted by bin, we can loop over the particles in neighboring bins. We make two passes over the particles. First, we launch a kernel to count the number of collision partners for each particle.
- Then, we sum these numbers and allocate space for our neighbor list.
- Finally, we make a another pass over the particles, putting them into to list at the appropriate place.

Note that we build a \emph{full} neighbor list, meaning that if particle $i$ appears in particle $j$'s list, then particle $j$ also appears in particle $i$'s list. This simplifies the force-computation step when using these lists, since the forces and torques for a given particle can be updated without atomics.

The final on-grid neighbor list data structure consists of two arrays. First, we have the neighbor list itself, stored as a big, 1D array of particle indices. Then, we have an `offsets` array that stores, for each particle, where in the neighbor list array to look. The details of this data structure have been hidden inside an iterator, so that user code can look like:

\begin{minted}{cpp}
.. code-block:: c

       // now we loop over the neighbor list and compute the forces                                                                                                                                                
        AMREX_FOR_1D ( np, i,
        {
@@ -51,7 +48,6 @@ The final on-grid neighbor list data structure consists of two arrays. First, we

                ...
            }
\end{minted}

Note that, because of our use of managed memory to store the particle data and the neighbor list, the above code will work when compiled for either CPU or GPU.

@@ -60,12 +56,11 @@ The above algorithm deals with constructing a neighbor list for the particles on
Once the neighbor list has been constructed, collisions with both particles and walls can easily be processed. 

We have created a GPU branch of MFIX that is capable of running with GPU support. As of this writing, the following operations in MFIX have been offloaded:
\begin{enumerate}
    \item Neighbor particles / neighbor list construction
    \item Particle-particle collisions
    \item Particle-wall collisions
    \item PIC Deposition (used in putting the drag force and solids volume fraction on the grid)
\end{enumerate}

- Neighbor particles / neighbor list construction
- Particle-particle collisions
- Particle-wall collisions
- PIC Deposition (used in putting the drag force and solids volume fraction on the grid)