sandbox/Antoonvh/inertial-particles.h

    Particle movements

    This header file implements a generic particle-path solving definition:

    \displaystyle \frac{\mathrm{d}\mathbf{x}_p}{\mathrm{d}t} = \mathbf{v}_p.

    Unlike flow-tracer particles, the intertial particle velocity is governed by the definition of its acceleration vector \mathbf{a}_p.

    \displaystyle \frac{\mathrm{d}\mathbf{v}_p}{\mathrm{d}t} = \mathbf{a}_p.

    For this purpose, a function can be specified following the prototype for p_acc(). Furthermore, particles must atleast have some additonal members to facilitate 2nd-order acurate advection and we aim to achieve compatibility with tracer-particles.h, which should then be included first.

    #include "run.h"
#ifndef ADD_PART_MEM 
#define ADD_PART_MEM coord u; long unsigned int tag; // u -> v_p
#endif

#include "particle.h"

#ifndef PA_INP
#define PA_INP //empty
#define INP (PA_inp){p()}
#endif

typedef struct PA_inp {
  particle p;
  PA_INP
} PA_inp;

coord p_acc (PA_inp inp);

Particles * inertial_particles = NULL;

extern scalar * _automatics_;

    Time integration

    A version of the second-order accurate Verlet method is used. (a.k.a. leaffrogging). Since p_acc may depend on p.u, we advance the particles in two solver timesteps.

    \displaystyle \mathbf{a}_{p,i} = \mathtt{p_{acc}}\left(\mathtt{p_i \ ADD\_ARGS}\right),

    \displaystyle \mathbf{v}_{p, i + 1} = \mathbf{v}_{p, i - 1} + \mathbf{a}_{p,i}\times \ 2 \mathtt{dt_i},

    \displaystyle \mathbf{x}_{p, i + 2} = \mathbf{x}_{p, i} + \mathbf{v}_{p,i + 1} \times \left( \mathtt{dt_i + dt_{i+1}}\right).

    The following events take care of the time integration.

    double dtprev;

event defaults (t = 0);

event init (t = 0);

event set_dtmax (i++);

event velocity (i++);

void ip_step1 (Particles P, int i) {
  double Dt = i > 0 ? 2*dt : dt;
  boundary (_automatics_);
  foreach_particle_in(P) {
    coord a = p_acc (INP);
    foreach_dimension()
      p().u.x += a.x*Dt;
  }
}

void ip_step2 (Particles P) {
  foreach_particle_in(P) {
    foreach_dimension()
      p().x += p().u.x*(dt + dtprev);
  }
}

event intertial_particles_step1 (i += 2, last) {
  dtprev = dt;
  foreach_P_in_list(inertial_particles) {
    particle_boundary (P);
    ip_step1 (P, i);
  }
}

event intertial_particles_step2 (i = 1; i += 2) {
  foreach_P_in_list(inertial_particles) 
    ip_step2 (P);
}

event free_intertial_particles (t = end) {
  if (inertial_particles)
    free (inertial_particles);
  inertial_particles = NULL;
}

    Utilities

    From tracer-particles.h, a bunch of functions are copied to provide some user interface for particle initializaiton.

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

void tag_ip_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 + offset;
}

void set_ip_particle_attributes (Particles p) {
  tag_ip_particles (p);
}

    Inertial-particle (ip) initialization

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

    Particles init_ip_cells(void) {
  int np = 0;
  foreach()
    np++;
  Particles p = new_inertial_particles (np);
  place_in_cells(p);
  particle_boundary (p); //?
  set_ip_particle_attributes (p);
  return p;
}

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

Particles init_ip_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_inertial_particles (inp.n);
    place_in_circle (p, inp);
  } else { 
    p = new_inertial_particles (0);
  }
  particle_boundary (p);
  set_ip_particle_attributes (p);
  return p;
}

    Todo

    • More Tests
    • Implement usages

    Tests

    Usage