src/test/kh-ns.c

    Lock exchange (Kelvin–Helmoltz shear instability)

    This is the centered Navier–Stokes solver version of the multilayer lock-exchange test where a more detailed description can be found.

    Animation of the density field

    set term svg enhanced size 1000,170 font ",10"
    set size ratio -1
    unset key
    unset ytics
    plot 'log' u 1:2 w l lc rgbcolor "black"
    Ten equally-spaced contours of density at t= 21. (script)

    Ten equally-spaced contours of density at t= 21. (script)

    #include "grid/multigrid.h"
    #include "navier-stokes/centered.h"
    #include "tracer.h"
    #include "diffusion.h"
    #include "navier-stokes/perfs.h"

    We add a tracer and define the acceleration of gravity.

    scalar T[], * tracers = {T};
    double G = 9.81;

    The acceleration field is allocated to store the Boussinesq buoyancy term.

    face vector av[];
    
    int main()
    {

    We define a rectangular domain using (eight) parallel domains.

      dimensions (ny = 1);
      L0 = npe();
      X0 = -L0/2.;
      N = 64*npe();

    Acceleration field, maximum timestep and viscosity.

      a = av;
      DT = 0.03;
      const face vector muc[] = {2e-4,2e-4};
      mu = muc;

    This is required due to a problem with the automatic increase in the number of relaxations in poisson.h.

      // fixme
      TOLERANCE = HUGE;
      NITERMIN = 4;
      
      run();

    The contour lines are saved in ‘log’ to serve as reference solution.

      system ("gnuplot -e 'set table' kh-ns.plot | sed '/^#.*/d' > log");
    }

    We add diffusion for the tracer, with a Peclet number unity.

    event tracer_diffusion(i++) {
      diffusion (T, dt, mu);
    }

    This is the Bousinesq buoyancy vertical acceleration, with a relative density ratio \Delta\rho/\rho = 10^{-2}.

    event acceleration (i++)
    {
      double dr = 0.01;
      foreach_face(y)
        av.y[] = G*dr*(T[] + T[0,-1])/2.;
    }

    The initial conditions with a hyperbolic tangent initial profile, and limited tracer gradient.

    event init (i = 0)
    {
      T.gradient = minmod;
      foreach()
        T[] = - 0.5*tanh((x + 0.1*cos(pi*y/2.))/0.04);
      boundary ({T});
    }

    Density field at t= 21.

    event density (t = 21)
    {
      output_field ({T}, box = {{-L0/2.,0},{L0/2,1}});
    }

    We generate an animation of the density field.

    event movie (t += 1; t <= 30)
    {
      output_ppm (T, file = "T.mp4", spread = -1, linear = true, n = 128*npe(),
    	      box = {{-L0/2.,0},{L0/2,1}});
    }