/** # Lagrangian particle advection 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 /** ## Advection scheme 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: $$\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}$: $$\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.](https://www.sciencedirect.com/science/article/pii/S0021999119300683). */ 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$\times$l grid centered at {xo, yo} with n$\times$n 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: ~~~literatec ... #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 /** ## Tests * [A quality test for the advection schemes](parttest.c) * [Quantative tests for the advection schemes](tp.c) * [Test for the MPI particle exchange](test_manyparticles.c) * [Particles, vof and tracer field in a vortex comparison](reversed.c) ## Usage * [All pages using particles.h](http://basilisk.fr/_search?patterns=particles.h) * [Tag a portion of a fluid](splash.c) * [Settling of volcanic ash](ash.c) * [Flow in a thermosyphon](tube.c) * [Planetary core convection](core.c) * [LES of isotropic turbulence](isotropicLES.c) * [Laminar mixing in 2D](laminarmixing.c) * [Axisymmetric mixing](coffee.c) ## Todo * Tag and trace particles with a constant unique number * Sort particles along grid iterator curve * Implement a proper Basilisk particles coordinates field */