# src/test/taylor-green.c

# 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
event tracer_advection (i++)
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 32^{2} to 256^{2}.

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