# 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.

```
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"
```

```
#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.

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);
}
```

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}});
}
```