sandbox/Antoonvh/tracer-particles.h

    Langrangian particles

    Tracer-particle coordinates \mathbf{p_i} are advected by a vector field \mathbf{u}:

    \displaystyle \frac{\mathrm{d}\mathbf{p_i}}{\mathrm{d}t} = \mathrm{u}|_{\mathbf{x} = \mathbf{p_i}}

    Definititions and types

    The generic particle class is used and tuned: For the time-advancement scheme, additional storage for velocities and locations is required. Furthermore, particles may be tagged. The tracer particles are stored in an array of groups of particles. Note that Particles reference a list of particles.

    extern vector u;
#ifndef ADD_PART_MEM
#define ADD_PART_MEM  coord u; coord u2; coord locn; long unsigned int tag; 
#endif
#include "particle.h"
Particles * tracer_particles = NULL;
bool P_RK2 = false;

    Time integration

    The used particle-advection scheme is the RK-3 variant 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..

    This scheme is used to advance the solution in two solver iterations between t_{n-2} and t_n. In case \alpha \approx \frac{2}{3}, or P_RK2 == true, we skip the last stage and switch to a generic second-order-accurate scheme.

    double tol = 3.e-2;
double dtf[2]; //timesteps

    The user is supposed to add particles in some relevant event, perhaps called init.

    event defaults (t = 0);

event init (t = 0);

event set_dtmax (i++);

event velocity (i++);

    By default, second-order-in-space “linear” interpolation is used to approximate \mathrm{u}|_{\mathbf{x} = \mathbf{p_i}}. This could be overloaded with higher-order methods.

    e.g for a 3rd order accurate method:

    ...
#include "higher-order.h" //This is in /sandbox/Antoonvh/
#define interpolate_linear interpolate_quadratic
#include "tracer-particles.h" 
...

    Velocities may found and stored in p().u or p().u2.

    void interpolate_vel (Particles p, int swtch) {
  if (swtch == 1) {
    foreach_particle_in(p) {
      if (locate (x, y, z).level < 0) {//freeze
	fprintf (stderr, "#pid: %d could not locate p.tag: %ld.\n",
		 pid(), p().tag); //This should not happen
	foreach_dimension()
	  p().u.x = 0;
      } else {
	foreach_dimension()
	  p().u.x =interpolate_linear (locate (x, y, z),
				       u.x, x, y, z);
      }
    }
  } else {
    foreach_particle_in(p) {
      if (locate (x, y, z).level < 0) {//freeze
	fprintf (stderr, "#pid: %d could not locate p.tag: %ld.\n",
		 pid(), p().tag); //This should not happen
	foreach_dimension()
	  p().u2.x = 0;
      } else {
	foreach_dimension()
	  p().u2.x =interpolate_linear (locate (x, y, z),
					u.x, x, y, z);
      }
    }
  }
}

    The three Runge-Kutta stages are defined below:

    void RK_step1 (Particles p, double dtf[2]) {
  interpolate_vel(p, 1);
  foreach_particle_in(p) { 
    foreach_dimension() {
      p().locn.x = p().x; //store the locations at t_n 
      p().x += dtf[0]*p().u.x;
    }
  }
}

void RK_step2 (Particles p, double dtf[2]) {
  interpolate_vel(p, 2);
  double a1 = -1, a2 = 2., h = dtf[1] + dtf[0];
  if (dtf[1] != dtf[0] || P_RK2) {
    double c = dtf[0]/h;
    if (fabs (c - 2./3.) > tol && !P_RK2)
      a2 = (c - 1.)/(c*(3.*c - 2.));
    else // Raltson's 2nd order method
      a2 = 1./(2*c);
    a1 = 1 - a2;
  }
  foreach_particle_in(p) {
    foreach_dimension() 
      p().x = p().locn.x + h*(a1*p().u.x + a2*p().u2.x);
  }
}

void RK_step3 (Particles p, 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);
    foreach_particle_in(p) {
      foreach_dimension() 
        p().x = p().locn.x + h*(b1*p().u.x + b2*p().u2.x +
				b3*interpolate_linear (locate (x, y, z),
				 u.x, x, y, z));
    }
  }
}

    The particles are advanced in time in various RK-stage events. Note that the particle coordinates (p().x) only contain the proper approximated location inbetween step3 and step1. These locations are otherwise stored in p().locn.

    event tracer_particles_step3 (i += 2, last);
event tracer_particles_step3 (i = 2; i += 2) {
  if (!P_RK2) {
    foreach_P_in_list (tracer_particles) {
      particle_boundary (P);
      RK_step3 (P, dtf);
    }
  }
}

event tracer_particles_step1 (i += 2, last) {
  dtf[0] = dt;
  foreach_P_in_list (tracer_particles) {
    particle_boundary (P);
    RK_step1 (P, dtf);
  }
}

event tracer_particles_step2 (i += 2, last);
event tracer_particles_step2 (i = 1; i += 2) {
  dtf[1] = dt;
  foreach_P_in_list (tracer_particles) {
    particle_boundary (P);
    RK_step2 (P, dtf);
    
  }
}

event free_tracers (t = end) {
  free (tracer_particles);
  tracer_particles = NULL;
}

    Initializing tracer particles

    n New tracer particles can be declared by calling new_tracer_particles (n). This is a building-block of tracer-particle initialization.

    Particles new_tracer_particles (long unsigned int n) {
  Particles p = new_particles (n);
  int l = 0, t = 0;
  if (tracer_particles != NULL) {
    while (pn[l] != terminate_int) {
      if (l == tracer_particles[t]) 
	t++;
      l++;
    }
  }
  tracer_particles = realloc (tracer_particles, (t + 1)*sizeof(Particles));
  tracer_particles[t] = p;
  return p;
}

void tag_particles (Particles p) {
  long unsigned int offset = 0;
#if _MPI
  long unsigned int pntr[npe()], tag_start[npe()];
  MPI_Gather (&pn[p], 1, MPI_UNSIGNED_LONG,
	      pntr, 1, MPI_UNSIGNED_LONG,
	      0, MPI_COMM_WORLD);
  if (pid() == 0) {
    tag_start[0] = 0;
    for (int tr = 1; tr < npe(); tr++) 
      tag_start[tr] = tag_start[tr - 1] + pntr[tr - 1];
  }
  MPI_Scatter (tag_start, 1, MPI_UNSIGNED_LONG,
	       &offset, 1, MPI_UNSIGNED_LONG,
	       0, MPI_COMM_WORLD);
#endif
  foreach_particle_in(p)
    p().tag = _j_particle + offset;
}

    We need a ugly hack to access ADD_PART_MEMBERS, as they only become available after the #define macros are processed.

    @define part p 
void set_a_particle_attributes (particle * p) {
  double x = p->x, y = p->y, z = p->z;
  foreach_dimension() {
    part->locn.x = p->x;
    part->u.x = interpolate (u.x, x, y, z);
    part->u2.x = 0.;
  }
}

void set_particle_attributes (Particles p) {
  interpolate_vel (p, 1);
  foreach_particle_in(p) {
    foreach_dimension() {
      p().locn.x = p().x;
      //p().u2.x = p().u.x; //Do or dont?
    }
  }
  tag_particles (p);
}

    Tracer-particle (tp) initialization

    Tracer-particles may be initialized by calling the following functions.

    Particles init_tp_cells(void) {
  int np = 0;
  foreach(serial)
    np++;
  Particles p = new_tracer_particles (np);
  place_in_cells(p);
  particle_boundary (p); //?
  set_particle_attributes (p);
  return p;
}

Particles init_tp_square (struct Init_P inp) {
  if (!inp.n)
    inp.n = 10;
  if (!inp.l)
    inp.l = L0;
  Particles p;
  if (pid() == 0) {
    p = new_tracer_particles (sq(inp.n));
    place_in_square (p, inp);
  } else { 
    p = new_tracer_particles (0);
  }
  particle_boundary (p);
  set_particle_attributes (p);
  return p;
}

Particles init_tp_circle (struct Init_P inp) {
  Particles p;
   if (!inp.n)
     inp.n = 100;
  if (!inp.l)
    inp.l = L0/2.;
  if (pid() == 0) {
    p = new_tracer_particles (inp.n);
    place_in_circle (p, inp);
  } else { 
    p = new_tracer_particles (0);
  }
  particle_boundary (p);
  set_particle_attributes (p);
  return p;
}

    To do

    • Implement functions for useful tracer-particle diagnostics
    • Check if tracer-particles can exist alongside other Particles
    • Box boundary conditions for particles (other than periodic)
    • Tag particles (for MPI)
    • Bind tracer particles with vof interfaces

    Tests

    Usage

    See also