# Incompressible Navier–Stokes solver (MAC formulation)

We wish to approximate numerically the incompressible Navier–Stokes equations \displaystyle \partial_t\mathbf{u}+\nabla\cdot(\mathbf{u}\otimes\mathbf{u}) = -\nabla p + \nabla\cdot(\nu\nabla\mathbf{u}) \displaystyle \nabla\cdot\mathbf{u} = 0

We will use the generic time loop, a CFL-limited timestep and we will need to solve a Poisson problem.

#include "run.h"
#include "timestep.h"
#include "poisson.h"

The Markers-And-Cells (MAC) formulation was first described in the pioneering paper of Harlow and Welch, 1965. It relies on a face discretisation of the velocity components u.x and u.y, relative to the (centered) pressure p. This guarantees the consistency of the discrete gradient, divergence and Laplacian operators and leads to a stable (mode-free) integration.

scalar p[];
face vector u[];

The only parameter is the viscosity coefficient \nu.

The statistics for the (multigrid) solution of the Poisson problem are stored in mgp.

double nu = 0.;
mgstats mgp;

We need to explicitly set the zero normal velocity condition.

u.n[left]   = 0;
u.n[right]  = 0;
u.n[top]    = 0;
u.n[bottom] = 0;

We apply boundary conditions after user initialisation.

event init (i = 0)
{
boundary ({p,u});
}

## Time integration

In a first step, we compute \mathbf{u}_* \displaystyle \frac{\mathbf{u}_* - \mathbf{u}_n}{dt} = \nabla\cdot\mathbf{S} with \mathbf{S} the symmetric tensor \displaystyle \mathbf{S} = - \mathbf{u}\otimes\mathbf{u} + \nu\nabla\mathbf{u} = \left(\begin{array}{cc} - u_x^2 + 2\nu\partial_xu_x & - u_xu_y + \nu(\partial_yu_x + \partial_xu_y)\\ \ldots & - u_y^2 + 2\nu\partial_yu_y \end{array}\right)

The timestep for this iteration is controlled by the CFL condition (and the timing of upcoming events).

double dtmax;

event set_dtmax (i++,last) dtmax = DT;

event stability (i++,last) {
dt = dtnext (timestep (u, dtmax));
}

{

We allocate a local symmetric tensor field. To be able to compute the divergence of the tensor at the face locations, we need to compute the diagonal components at the center of cells and the off-diagonal component at the vertices.

  symmetric tensor S[]; // fixme: this does not work on trees

We average the velocity components at the center to compute the diagonal components.

  foreach()
foreach_dimension()
S.x.x[] = - sq(u.x[] + u.x[1,0])/4. + 2.*nu*(u.x[1,0] - u.x[])/Delta;

We average horizontally and vertically to compute the off-diagonal component at the vertices.

  foreach_vertex()
S.x.y[] =
- (u.x[] + u.x[0,-1])*(u.y[] + u.y[-1,0])/4. +
nu*(u.x[] - u.x[0,-1] + u.y[] - u.y[-1,0])/Delta;

We only need to apply boundary conditions to the diagonal components.

  boundary ({S.x.x, S.y.y});

Finally we compute \displaystyle \mathbf{u}_* = \mathbf{u}_n + dt\nabla\cdot\mathbf{S}

  foreach_face()
u.x[] += dt*(S.x.x[] - S.x.x[-1,0] + S.x.y[0,1] - S.x.y[])/Delta;
}

### Projection

In a second step we compute \displaystyle \mathbf{u}_{n+1} = \mathbf{u}_* - \Delta t\nabla p with the condition \displaystyle \nabla\cdot\mathbf{u}_{n+1} = 0 This gives the Poisson equation for the pressure \displaystyle \nabla\cdot(\nabla p) = \frac{\nabla\cdot\mathbf{u}_*}{\Delta t}

event projection (i++,last)
{
boundary ((scalar *){u});
const face vector unity[] = {1.,1.};
mgp = project (u, p, unity, dt, mgp.nrelax);
}