Implementation

Parallelisation

Since this library is built on a MPI-parallelised DNS solver, Eulerian domain is decomposed in \(y\) direction. Regarding particles, update procedures such as momentum exchange are conducted on the same Eulerian domain and thus the parallelisation is straightforward. The velocity and position of the gravity centers are stored in a Lagrangian way and shared among all processes. Since they are \(\mathcal{O} \left( N_p \right)\) operations, the computational cost is negligibly small compared to the fluid solver \(\mathcal{O} \left( N_x \times N_y \right)\): proportional to the number of grid points.

Collisions (judgements and computations of repellent forces), on the other hand, require relatively heavy computational loads since we need \(\mathcal{O} \left( N_p^2 \right)\) operations. Although it is still much smaller than \(N_x \times N_y\), it can be fairly heavy to be completed by a main process when the volume fraction or the domain size is large. Thus this library parallelises the procedure simply as the same way as the decomposing the Eulerian domain. Namely, there are following collision pairs to be considered:

	\(0\)	\(1\)	\(2\)	\(3\)	\(\cdots\)	\(N-2\)	\(N-1\)
\(0\)		\(0\)	\(1\)	\(2\)	\(\cdots\)	\(N-3\)	\(N-2\)
\(1\)			\(N-1\)	\(N\)	\(\cdots\)	\(2N-5\)	\(2N-4\)
\(2\)				\(2N-3\)	\(\cdots\)	\(3N-8\)	\(3N-7\)
\(\cdots\)					\(\cdots\)
\(N-1\)					\(\cdots\)		\(\frac{N \left( N-1 \right)}{2}-1\)

Here the left-most column \(n_0\) and the top-most row \(n_1\) show indices of the particles. The other numbers show the indices of pairs, which are separated to each processor and considered independently.

In practice, we need to translate the index of collision pair n to the corresponding particle indices \(n_0\) and \(n_1\), which are handled by get_particle_indices.