# sandbox/Antoonvh/particles.h

n_part Partciles can be seeded in an array of coords called loc. The formulations on this page aim to advect it in a time loop by an external vector field u. There is an option for a third-order-accurate time-marhcing scheme, 2nd order is the default.

#include "run.h"
long unsigned int n_part;  //Number of particles (foreach MPI thread)
coord * loc;               //coordinates of particles
bool P_RK3 = false;        //Switch for RK-3 scheme

## Particle methods

First, some macros are defined to create a foreach_particle() iterator. It will give access to the positions (x,y,z), their writeable counterparts (p_x(), etc), velocity components (k1_x(), etc) and other stuff. It can even do simple reductions (well done qcc!).

foreach_dimension() {
@define ArrIP_x (j*dimension + cind.x)
@define p_x()  loc[j].x
@define pn_x() loc_n[j].x
@define k1_x() k1[j*dimension + cind.x]
@define k2_x() k2[j*dimension + cind.x]
}

@def PARTICLE_VARIABLES
bool LN = false;
if (loc_n != NULL)
LN = true;
double x = loc[j].x; NOT_UNUSED(x);
double xn = 0;
if (LN)
xn = loc_n[j].x;
NOT_UNUSED(xn);
double y = loc[j].y; NOT_UNUSED(y);
double yn = 0;
if (LN)
yn = loc_n[j].y;
NOT_UNUSED(yn);
double z = loc[j].z; NOT_UNUSED(z);
double zn = 0;
if (LN)
zn = loc_n[j].z;
NOT_UNUSED(zn);
@

@def foreach_particle() {
for (int j = 0; j < n_part; j++) {
PARTICLE_VARIABLES
@
@def end_foreach_particle()
}
}
@

In order to facilitate the higher-order advection, some storage arrays are declared.

long unsigned int n_part_a;                //allocated size of arrays
coord * loc_n;                             //Location storage
double *k1, *k2, dtf[2];                   //Velocities and timesteps
bool part_linear = true, start = true;     //Linear interpolation and RK3 first iteration
double tol = 1.e-2;                        //RK3: Tolerance on alpha != 2/3

typedef struct {
int x, y, z;
} coordi;
coordi cind = {0,1,2};

The above data arrays should be little concern to the user as the code tries to manage these scratch arrays in the background.

event init (t = 0) {
if (loc_n == NULL) {//We have not already been here
n_part_a = n_part + 1;
#if _MPI
n_part_a *= 2;
loc = realloc (loc, n_part_a*sizeof(coord));
#endif
loc_n = malloc (n_part_a*sizeof(coord));
k1 = malloc (n_part_a*dimension*sizeof(double));
k2 = malloc (n_part_a*dimension*sizeof(double));
}
}

event set_dtmax (i++);

void part_boundaries () {
coord mind = {X0, Y0, Z0};
foreach_particle() {
foreach_dimension() {
if (p_x() < mind.x)
p_x() += L0;
else if (p_x() > (mind.x + L0))
p_x() -= L0;
}
}
}

#if _MPI

Interpolate_array is convinient. We re-implement it, assuming all local particles.

void update_mpi (int step); //particle exchange function prototype
#endif
trace
void interpolate_array_local (scalar * list, coord * a, int n, double * v, bool linear) {
int j = 0;
for (int i = 0; i < n; i++) {
Point point = locate (a[i].x, a[i].y, a[i].z);
if (point.level < 0) {
fprintf (stderr, "#pid: %d found a non-local particle...\n", pid());
for (scalar s in list)
v[j++] = 0;
} else
for (scalar s in list) {
v[j++] =interpolate_linear (point,
(struct _interpolate){s, a[i].x, a[i].y, a[i].z});
}
}
}
#define interpolate_array interpolate_array_local

To facilitate the general case of time-dependend flow, particles are advected in two iterations between t_{n-2} and t_n. The formulation follows a variable two-stage Runge-Kutta scheme:

\displaystyle \begin{array}{c|cc} 0 & 0 & \\ \alpha & \alpha & \\ \hline & 1 - \frac{1}{2\alpha} & \frac{1}{2\alpha} \end{array}

Which can be extended to the RK-3 scheme of Sanderse and Veldman (2019) for \alpha \neq \frac{2}{3} \pm \mathtt{tol}:

\displaystyle \begin{array}{c|ccc} 0 & 0 & & \\ \alpha & \alpha & & \\ 1 &1+\frac{1- \alpha}{\alpha (3\alpha -2)} & -\frac{1- \alpha}{\alpha (3\alpha -2)} &\\ \hline & \frac{1}{2}-\frac{1}{6\alpha} & \frac{1}{6\alpha(1-\alpha)} & \frac{2-3\alpha}{6(1-\alpha)} \\ \end{array}

see:

B. Sanderse and AEP Veldman, Constraint-consistent Runge-Kutta methods for one-dimensional incompressible multiphase flow, J. Comp. Phys. (2019) link to JCP..

void RK_step1 (coord * loc, coord * loc_n, double * k1, double dtf[2]) {
interpolate_array ((scalar*){u}, loc, n_part, k1, part_linear);
foreach_particle()
foreach_dimension() {
pn_x() = p_x(); //store the locations at t_n
p_x() += dtf[0]*k1_x();
}
}

void RK_step2 (coord * loc, coord * loc_n, double * k1, double * k2, double dtf[2]) {
interpolate_array ((scalar*){u}, loc, n_part, k2, part_linear);
double a1 = -1, a2 = 2, h = dtf[1] + dtf[0];
if (dtf[1] != dtf[0] || !P_RK3) {
double c = dtf[0]/h;
if (fabs (c - 2./3.) > tol && P_RK3)
a2 = (c - 1.)/(c*(3.*c - 2.));
else //Raltson's 2nd order method
a2 = 1./(2*c);
a1 = 1 - a2;
}
foreach_particle()
foreach_dimension()
p_x() = pn_x() + h*(a1*k1_x() + a2*k2_x());
}

void RK_step3 (coord * loc, coord * loc_n, double * k1, double * k2, double dtf[2]) {
double h = dtf[1] + dtf[0];
double c = dtf[0]/h;
if (fabs(c - 2./3.) > tol) {// RK-3
double b1 = 0.5 - 1./(6.*c);
double b2 = 1./(6.*c*(1. - c));
double b3 = 1. - (b1 + b2);
double V[dimension*n_part];
interpolate_array ((scalar*){u}, loc, n_part, V, true);
foreach_particle()
foreach_dimension()
p_x() = pn_x() + h*(b1*k1_x() + b2*k2_x() + b3*V[ArrIP_x]);
}
}

Particle advection is performed in the advance_particles event. Varying between even and uneven iterations.

event advance_particles (i++, last) {
part_boundaries();
if (i%2 == 0) {
if (i > 0 && P_RK3) {
#if _MPI
update_mpi(3);
#endif
RK_step3 (loc, loc_n, k1, k2, dtf);
part_boundaries();
}
#if _MPI
update_mpi(1);
#endif
dtf[0] = dt;
RK_step1 (loc, loc_n, k1, dtf);
} else {
#if _MPI
update_mpi(2);
#endif
dtf[1] = dt;
RK_step2 (loc, loc_n, k1, k2, dtf);
}
}

event free_particles (t = end, last) {
free (loc);
free (loc_n);
free (k1);
free (k2);
loc_n = NULL;
}

## User functions for particle seeding

The following function initializes a particle at the centre of each grid cell.

void init_particles_in_cells(){
n_part = 0;
foreach()
n_part++;
loc = malloc (n_part*sizeof(coord));
int n = 0;
foreach() {
coord cc = {x, y, z};
foreach_dimension()
loc[n].x = cc.x;
n++;
}
}

Initialize particles in a 2D l\timesl grid centered at {xo, yo} with n\timesn particles:

void init_particles_2D_square_grid (int n, double xo, double yo, double l){
n_part = 0;
if (pid() == 0) {
n_part = sq(n);
loc = malloc (n_part*sizeof(coord));
int i = 0;
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
loc[i].x = xo - l/2. + (double)j*(l/((double)n - 1.));
loc[i].y = yo - l/2. + (double)k*(l/((double)n - 1.));
i++;
}
}
} else
loc = malloc (sizeof(coord));
}

The following function places n particles randomly in a circle with radius R at location {xo, yo}.

void init_particles_random_circle (int n, double xo, double yo, double R) {
n_part = 0;
if (pid() == 0) {
n_part = n;
loc = malloc (n_part*sizeof(coord));
int j = 0;
while (j < n_part) {
double a = noise();
double b = noise();
if (sq(a) + sq(b) < R) {
loc[j].x = a + xo;
loc[j].y = b + yo;
j++;
}
}
} else
loc = malloc (sizeof(coord));
}

## Visualization

For visualization with bview, the scatter() function can be used when BVIEW \neq 0.

use something like:

...
#include "view.h"
#define BVIEW 1
#include "particles.h"
...
#if BVIEW
#include "scatter.h"
#endif

## Implementation of the MPI particle exchange:

The ugliest bit is saved for last. Local partciles outside the MPI domain are comunicated globally. In turn, each thread selects those inside their domain. This facilitates advection at large CFL numbers, grid adaptation and rebalancing.

#if _MPI
void update_mpi (int step) {
int outt, out = 0, in = 0, m = 0;
int psize = step < 2 ? 1 : step < 3 ? 3 : 4;
psize *= dimension;
//Count the number of outgoing particles per thread
foreach_particle()
if (locate (x, y, z).level < 0)
out++;
//get indices and outgoing data
int ind[out];
double senddata[out*psize];
foreach_particle() {
if (locate (x, y, z).level < 0) {
ind[m] = j;
int c = 0;
foreach_dimension()
senddata [m*psize + c++] = loc[j].x;
if (step > 1) {
foreach_dimension() {
senddata [m*psize + c++] = loc_n[j].x;
senddata [m*psize + c++] = k1[j*dimension + cind.x];
}
}
if (step > 2) {
foreach_dimension() {
senddata [m*psize + c++] = k2[j*dimension + cind.x];
}
}
m++;
}
}
//remove the senddata from arrays (shrink)
m = 0;
int j = 0;
fflush (stdout);
while (j < n_part - out) {
while (m < out ? j + m == ind[m] : 0)
m++;
while (m < out ? j < n_part - out && j + m != ind[m] : j < n_part - out) {
loc[j]   = loc[j + m];
if (step > 1) {
loc_n[j]   = loc_n[j + m];
foreach_dimension()
k1[j*dimension + cind.x]   =  k1[(j + m)*dimension + cind.x];
}
if (step > 2) {
foreach_dimension()
k2[j*dimension + cind.x]   =  k2[(j + m)*dimension + cind.x];
}
j++;
}
}
// Gather lost particles among threads:
// First, count all of them
int outa[npe()], outat[npe()];
outat[0] = 0;
MPI_Allgather (&out, 1, MPI_INT, &outa[0], 1, MPI_INT, MPI_COMM_WORLD);
//Compute displacements
for (int j = 1; j < npe(); j++)
outat[j] = outa[j - 1] + outat[j - 1];
//compute total
outt = outat[npe() - 1] + outa[npe() - 1];
// Allocate recieve buffer and gather
double recdata[outt*psize];
for (int j = 0; j < npe(); j++) {
outat[j] *= psize;
outa[j]  *= psize;
}
//send and recieve data
MPI_Allgatherv (senddata, outa[pid()], MPI_DOUBLE,
recdata, outa, outat, MPI_DOUBLE,
MPI_COMM_WORLD);
//count new particles
for (int j = 0; j < outt ; j++) {
coord a;
foreach_dimension()
a.x = recdata[j*psize + cind.x];
if (locate (a.x, a.y, a.z).level > 0)
in++;
}
int n_partn = n_part + in - out;
//Manage the memory if required...
if (n_partn > n_part_a || 2*(n_partn + 1) < n_part_a ) {
n_part_a = 2*(n_partn + 1);
loc = realloc   (loc  , n_part_a*sizeof(coord));
loc_n = realloc (loc_n, n_part_a*sizeof(coord));
k1 = realloc    (k1, n_part_a*sizeof(double)*dimension);
k2 = realloc    (k2, n_part_a*sizeof(double)*dimension);
}
//Collect new particles from recdata
if (in > 0) {
int indi[in];
m = 0;
for (int j = 0; j < outt; j++) {
coord a;
foreach_dimension()
a.x = recdata[j*psize + cind.x];
if (locate (a.x, a.y, a.z).level > 0) {
indi[m++] = j;
}
}
m = 0;
for (j = n_part - out; j < n_partn; j++) {
int c = 0;
foreach_dimension()
loc[j].x = recdata[indi[m]*psize + c++] ;
if (step > 1) {
foreach_dimension() {
loc_n[j].x = recdata [indi[m]*psize + c++];
k1[j*dimension + cind.x] = recdata [indi[m]*psize + c++];
}
}
if (step > 2) {
foreach_dimension() {
k2[j*dimension + cind.x] = recdata [indi[m]*psize + c++];
}
}
m++;
}
}
//Update n_part
n_part = n_partn;
}
#endif

## Todo

• Tag and trace particles with a constant unique number
• Sort particles along grid iterator curve
• Implement a proper Basilisk particles coordinates field