Merging of two vortices (centered Euler solver)

This test is similar to stream.c but uses the centered Navier–Stokes solver (without viscosity). It also shows how to convert a vorticity field into a velocity field.

#include "navier-stokes/centered.h"

The domain is centered on (0,0) and the maximum level of refinement is 8 i.e. the initial grid has 28=256 grid points per dimension.

#define MAXLEVEL 8

// This is necessary for convergence when lowering the tolerance
uf.n[left]   = 0.;
uf.n[right]  = 0.;
uf.n[top]    = 0.;
uf.n[bottom] = 0.;

int main()
  origin (-0.5, -0.5);
  init_grid (1 << MAXLEVEL);
  //  TOLERANCE = 1e-12;

For the centered Navier–Stokes solver, the primary variables are the velocity and pressure field. We need to convert the initial vorticity field into the velocity field. To do so we first declare the (temporary) streamfunction ψ and vorticity ω fields. We also set appropriate boundary conditions for the streamfunction.

event init (t = 0)
  scalar ψ[], ω[];

  ψ[left]   = dirichlet(0);
  ψ[right]  = dirichlet(0);
  ψ[top]    = dirichlet(0);
  ψ[bottom] = dirichlet(0);

We then initialise both fields, using the same initial condition for the vorticity as in stream.c, and apply boundary conditions. Note that it is necessary to initialise the streamfunction, as the Poisson solver requires an initial guess.

  double dd = 0.1;
  foreach() {
    ω[] = (exp(-(sq(x - dd) + sq(y))/(dd/10.)) +
	       exp(-(sq(x + dd) + sq(y))/(dd/10.)));
    ψ[] = 0.;
  boundary ({ψ,ω});

We then solve the Poisson equation 2ψ=ω and compute the centered velocity components by differentation of the streamfunction i.e. ux=yψ uy=xψ

  poisson (ψ, ω);
  coord f = {-1.,1.};
      u.x[] = f.x*(ψ[0,1] - ψ[0,-1])/(2.*Δ);
  boundary ((scalar *){u});

We output some statistics on the vorticity field and Poisson solver at the start and end of the simulation.

event logfile (t <= 30; t += 1) {
  scalar ω[];
  vorticity (u, ω);
  stats s = statsf (ω);
  fprintf (ferr, "%g %d %g %g %g %d\n", t, i, dt, s.sum, s.max, mgp.i);

We make animations of the vorticity and level of refinement.

event movie (t += 0.2; t <= 30) {
  static FILE * fp = popen ("ppm2mpeg > vort.mpg", "w");
  scalar ω[];
  vorticity (u, ω);
  output_ppm (ω, fp, linear = true);

  static FILE * fp1 = popen ("ppm2mpeg > level.mpg", "w");
    ω[] = level;
  output_ppm (ω, fp1, spread = 2);

  static FILE * fp2 = popen ("ppm2mpeg > pid.mpg", "w");
    ω[] = pid();
  output_ppm (ω, fp2, min = 0, max = npe() - 1);

We output the vorticity and level fields at regular intervals in a format compatible with gnuplot.

event output (t += 5) {
  static int nf = 0;
  scalar ω[];
  vorticity (u, ω);

  char name[80];
  sprintf (name, "omega-%d", nf);
  FILE * fp = fopen (name, "w");
  output_field ({ω}, fp, linear = true);
  fclose (fp);
  scalar l = ω;
    l[] = level;
  sprintf (name, "level-%d", nf);
  fp = fopen (name, "w");
  output_field ({l}, fp);
  fclose (fp);

If we are using a tree grid, it is adapted using wavelet error control on both components of the velocity field. Note that the error thresholds need to be specified twice (once for each component of vector u).

#if TREE
event adapt (i++) {
  adapt_wavelet ((scalar *){u}, (double[]){5e-5,5e-5}, MAXLEVEL);


After running and processing by gnuplot we get the following pictures and animations.