# Taylor–Green vortices

Taylor–Green vortices are one of the few exact non-trivial solutions of the incompressible Euler equations. In this test case, we use this solution as initial condition and check whether the numerical scheme can respect the balance between non-linear advection terms and pressure gradients. Numerical diffusion will in particular introduce dissipation. This dissipation can be quantified and is a useful measure of the accuracy of the numerical scheme.

We solve the incompressible Euler equations on a Cartesian (multi)grid either with the centered Navier-Stokes solver or with the “all Mach” solver (in incompressible mode).

``````#include "grid/multigrid.h"
#if ALL_MACH
# include "all-mach.h"
# include "bcg.h"
# define u q

advection ((scalar *){q}, uf, dt, (scalar *){g});
#else
# include "navier-stokes/centered.h"
#endif

int main() {``````

The domain is unity, centered on the origin and periodic in all directions.

``````  origin (-0.5,-0.5);
foreach_dimension()
periodic (right);``````

We check convergence with spatial resolution from 322 to 2562.

``````  for (N = 32; N <= 256; N *= 2)
run();
}

event init (i = 0) {``````

This is the initial Taylor–Green solution for velocity and pressure.

``````  foreach() {
u.x[] = - cos(2.*π*x)*sin(2.*π*y);
u.y[] =   sin(2.*π*x)*cos(2.*π*y);
p[]   = - (cos(4.*π*x) + cos(4.*π*y))/4.;
}``````

We also need to define the initial centered pressure gradient (this improves the accuracy of the initial conditions).

``````  boundary ({p});
foreach()
foreach_dimension()
g.x[] = - (p[1] - p[-1])/(2.*Δ);
boundary ((scalar *){g});
}

event logfile (i++) {``````

We log the evolution of the maximum divergence and of the total kinetic energy.

``````  scalar div[], ke[];
foreach() {
div[] = (u.x[1,0] - u.x[-1,0] + u.y[0,1] - u.y[0,-1])/(2.*Δ);
ke[] = sq(u.x[]) + sq(u.y[]);
}
printf ("%d %d %g %g %g\n", N, i, t, normf(div).max, statsf(ke).sum);
}

event error (t = 2) {``````

At $t=2$ we compute the error on the norm of the velocity.

``````  scalar e[];
foreach() {
double u0 = - cos(2.*π*x)*sin(2.*π*y);
double v0 =   sin(2.*π*x)*cos(2.*π*y);
e[] = norm(u) - sqrt(sq(u0) + sq(v0));
}
norm n = normf (e);
fprintf (stderr, "%d %g %g %g\n", N, n.avg, n.rms, n.max);
}``````

For this particular case, the Bell–Collela–Glaz advection scheme converges at third-order.

The divergence of the centered velocity field is well-behaved.

By fitting the decrease of the kinetic energy, we get an estimate of the numerical viscosity.