Particle advection test

    On this page we test particle advection for a non-trivial 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.

    We test the forward Euler scheme and the x/y-split semi-implicit method for time advancement.

    vector u[];
    #include "view.h"
    #include "poisson.h"
    #define BVIEW 1
    #include "particles.h"
    scalar omega[], psi[];
    int main(){
      L0 = 15;
      X0 = Y0 = -L0/2;
        periodic (left);
      init_grid (128);
      ssi = true; //Switch for the x/y-split semi-implicit method
    event init (t = 0) { //Init the dipole
      double k = 3.83170597;
      TOLERANCE = 10E-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 = pow(pow((x),2)+(pow((y),2)),0.5);
          double s = x/r;
          omega[] =  ((r<1)*((-2*j1(k*r)*s/(k*j0(k)))))*sq(k);
        boundary ({omega});
        poisson (psi, omega);
        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);


    We render two movies. Noting that drawing the vorticity field is the most expensive function on this page, followed by the save() command.

    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 (!ssi)
        save ("movie.mp4");
        save ("moviessi.mp4");

    We can compare the results from both methods,

    Forward Euler advection

    x/y-split semi-implicit advection

    It appears the tracers are not properly entrained by the dipole when using the forward method. The x/y split semi-implicit method appears to do better. We quantify this and log the number of tracers inside the dipole’s atmosphere:

    event log_tracers (i += 5) {
      int ni = 0;
      char fname[99];
      if (!ssi)
        sprintf (fname, "ni");
        sprintf (fname, "nissi");
      static FILE * fp1 = fopen (fname, "w");
      for (int j = 0; j < n_part; j++) {
        if (sq(loc[j].x) + sq(loc[j].y) <= 1.)
      fprintf(fp1, "%g\t%d\n", t, ni);

    This is the resulting plot:

       set xlabel 'time'
       set ylabel 'number of tracers'
       set key outside
       plot 'ni' u 1:2 w l lw 3 t 'Forward Euler', \
            'nissi' u 1:2 w l lw 3 t 'Split semi imp.'


    Quite convincing result for the split semi-implicit method. It is particularly of interest when you are interested in the statistics of particle paths over time. The x/y-split semi-implicit method works so well because, there exists a streamfunction (\psi) that serves as a Hamiltonian for the tracer paths. Combined with the fact the the split method is symplectic makes that, even with the presence of errors, particles remain ‘attached’ to their corresponding streamline over time.