sandbox/Antoonvh/parttest.c
Particle advection test
On this page we test particle advection for a non-trivial, steady flow field. We use the famous Lamb-Chaplygin vortex dipole soliton. It is well known that this structure is characterized by an enclosed stream line around the dipole’s atmosphere. You may inspect its structure here.
vector u[];
#include "view.h"
#include "poisson.h"
#define BVIEW 1
#include "particles.h"
scalar omega[];
int main() {
L0 = 15;
X0 = Y0 = -L0/2;
foreach_dimension()
periodic (left);
init_grid (256);
run();
P_RK3 = true;
run();
}
event init (t = 0) { //Init the dipole
scalar psi[];
double k = 3.83170597;
TOLERANCE = 1e-6;
The dipolar vortex is initialized on a grid with a maximum resolution that corresponds to 8192 \times 8192 equidistant cells such that the error in the velocity estimation is not more than 0.25% of the free stream velocity (U).
int maxlevel = 13;
while (depth() != maxlevel) {
adapt_wavelet ({u.x, u.y}, (double[]){0.0025, 0.0025}, maxlevel);
foreach() {
double r = sqrt(sq(x) + sq(y));
double s = x/r;
omega[] = ((r<1)*((-2*j1(k*r)*s/(k*j0(k)))))*sq(k);
}
boundary ({omega});
poisson (psi, omega);
boundary ({psi});
foreach() {
u.x[] = -((psi[0, 1] - psi[0, -1])/(2*Delta));
u.y[] = (psi[1, 0] - psi[-1, 0])/(2*Delta) + 1.;
}
boundary ({u.x, u.y});
}
DT = 0.01;
dt = dtnext (DT);
init_particles_2D_square_grid (40, 0, 0, 5.);// Add 40 x 40 = 1600 particles
}
event set_dtmax (i++)
dt = dtnext (DT);
Output
We render a movie.
event movie (t += 0.1; t <= L0*4.) {
view(width = 500, height = 500);
squares("omega", min = -10, max = 10, map = cool_warm);
scatter(loc = loc, s = 15);
if (!P_RK3)
save ("movie.mp4");
else
save ("movie_rk3.mp4");
}
Well done particles.h
!
event log_tracers (i += 5) {
int ni = 0;
char fname[99];
if (!P_RK3)
sprintf (fname, "particles");
else
sprintf (fname, "particles_rk3");
static FILE * fp1 = fopen (fname, "w");
for (int j = 0; j < n_part; j++) {
if (sq(loc[j].x) + sq(loc[j].y) <= 1.)
ni++;
}
fprintf(fp1, "%g\t%d\n", t, ni);
}
set xlabel 'time'
set ylabel 'number of tracers'
set key outside
plot 'particles' u 1:2 w l lw 3 t 'RK2',\
'particles_rk3' u 1:2 w l lw 3 t 'RK3'