# A Navier-Stokes equations solver with tracers A fourth-order accurate solver for the solution to:

\displaystyle \frac{\partial \mathbf{u}}{\partial t} + \left(\mathbf{u} \cdot \mathbf{\nabla}\right)\mathbf{u} = -\mathbf{\nabla} p + \nu \nabla^2 \mathbf{u} + \mathbf{a},

with the constraint that,

\displaystyle \mathbf{\nabla} \cdot \mathbf{u} = 0.

Furthermore, a scalar s can be advected and diffused,

\displaystyle \frac{\partial s}{\partial t} = -\mathbf{u} \cdot \mathbf{\nabla}s + \kappa \nabla^2 s.

The velocity components are represented discretely as face averages. Their tendencies due to advection and diffusion are computed at vertices wheareas the projection operators acts on the face-averaged quantities.

#include "higher-order.h" // Higher-order functions and definitions
#include "poisson4b.h"    // 4th-order Projection scheme
#include "my_vertex.h"    // Vertex functions and definitions
#include "run.h"          // Time loop

The global variables are,

face vector u[];          // Face averaged values
face vector df[];         // Tendency for velocity
scalar p[], p2[];         // Cell-averaged scalars for projection
(const) vector a;         // Acceleration: Vertex point values
(const) scalar nu, kappa; // Viscosity and diffusivity: Vertex point values
mgstats mgp, mgp2;        // MG-solver statistics
extern scalar * tracers;  // Mandatory vertex-based tracers (maybe NULL)
scalar * dsl = NULL;      // Tendencies for these tracers

#define freegrad (layer_nr_y == 1 ? 2.*val(_s,0,0,0) - val(_s,0,-1,0) : 3*val(_s,0,0,0) - 2*val(_s,0,-1,0))
#if NOSLIP_TOP
u.n[top] = dirichlet_vert_top4(0.);
df.n[top] = dirichlet_vert_top4(.0);
u.t[top] = dirichlet_top(0);
df.t[top] = dirichlet_top(0);
#endif
#if NOSLIP_BOTTOM
u.t[bottom] = dirichlet_bottom4(0);
df.t[bottom] = dirichlet_bottom4(0);
#endif

## Runge-Kutta Time integration

We use a low-storage time integrator

#ifndef RKORDER
#define RKORDER (4)
#endif
#if (RKORDER == 3)
// Williamson, J. H.: Low-Storage Runge-Kutta schemes, J.
// Comput.Phys., 35, 48–56, 1980.
#define STAGES (3)
double An[STAGES] = {0., -5./9., -153./128.};
double Bn[STAGES] = {1./3., 15./16., 8./15.};
#else
// Carpenter, M.  H.  and Kennedy, C.  A.: Fourth-order
// 2N-storageRunge-Kutta schemes, Tech. Rep. TM-109112, NASA
// LangleyResearch Center, 1994
#define STAGES (5)
double An[STAGES] = {0.,
-567301805773. /1357537059087.,
-2404267990393./2016746695238.,
-3550918686646./2091501179385.,
-1275806237668./842570457699.};
double Bn[STAGES] = {1432997174477./9575080441755. ,
5161836677717./13612068292357.,
1720146321549./2090206949498. ,
3134564353537./4481467310338. ,
2277821191437./14882151754819.};
#endif

The time stepper is implemented below. It delineates between tracers and the velocity components.

void A_Time_Step (double dt,
void (* Lu) (face vector uf, face vector du,
scalar * ul, scalar * dul)) {
if (dsl == NULL && tracers != NULL)
dsl = list_clone (tracers);
scalar * dsltmp = list_clone (dsl);
face vector dftmp[];
for (int Stp = 0; Stp < STAGES; Stp++) {
Lu (u, dftmp, tracers, dsltmp);
foreach_face() {
df.x[]  = An[Stp]*df.x[] + dftmp.x[];
u.x[]  += Bn[Stp]*df.x[]*dt;
}
foreach() {
scalar s, ds, dst;
for (s, ds, dst in tracers, dsl, dsltmp) {
ds[]  = An[Stp]*ds[] + dst[];
s[]  += Bn[Stp]*ds[]*dt;
}
}
scalar * bound = list_concat ((scalar*){u}, tracers);
boundary (bound);
free (bound);
}
delete (dsltmp); free (dsltmp); dsltmp = NULL;
}

Some default settings that should work for most scenarios.

event defaults (i = 0) {
#if TREE
for (scalar s in tracers) {
s.restriction = s.coarsen = restriction_vert;
s.refine = s.prolongation = refine_vert5;
}
u.x.refine = refine_face_solenoidal;
p.prolongation = refine_4th;
p2.prolongation = refine_4th;

foreach_dimension()
u.x.prolongation = refine_face_4_x;
#endif
CFL = 1.3;
compact_iters = 5;
}

event init (t = 0);

event call_timestep (t = 0) {
event ("timestep");
}

## Choosing the timestep size

Apart from the CFL condition, a stability criterion for the viscous term is included.

\displaystyle \mathrm{DI} < \frac{\mathrm{dt}\nu}{\Delta^2}

double DI = STAGES == 5 ? 0.2 : 0.1; //Maximum "Cell Diffusion" number
event timestep (i++, last) {
double dtm = HUGE;
foreach_face(reduction(min:dtm)) {
if (kappa.i)
if (DI*sq(Delta)/kappa[] < dtm)
dtm = DI*sq(Delta)/kappa[];
if (nu.i)
if (DI*sq(Delta)/nu[] < dtm)
dtm = DI*sq(Delta)/nu[];
if (fabs(u.x[]) > 0)
if (CFL*Delta/fabs(u.x[]) < dtm)
dtm = CFL*Delta/fabs(u.x[]);
}
dt = dtnext (min(DT, dtm));
}

### Diffusion

A 4th-order accurate second-derivative scheme is used for the viscous and diffusive terms.

#define D2SDX2 (-(s[-2] + s)/12. + 4.*(s + s[-1])/3. - 5.*s[]/2.)

## Computing the tendency fields

The tendency is computed from the field values for u and the tracers.

void adv_diff (face vector du, scalar * dsl) {
// Allocate some vertex vectors
vector v[], dv[];
v.n[top] = dirichlet_vert_top(0);
dv.n[top] = dirichlet_vert_top(a.y.i ? a.y[0,1] : 0);
v.n[bottom] = dirichlet_vert_bottom(0);
dv.n[bottom] = dirichlet_vert_bottom(0);
#if NOSLIP_TOP
v.t[top] = dirichlet_vert_top4(0);
dv.t[top] = dirichlet_vert_top4(0);
#endif
#if NOSLIP_BOTTOM
v.t[bottom] = dirichlet_vert_bottom4(0);
dv.t[bottom] = dirichlet_vert_bottom4(0);
#endif
scalar * trcrs  = list_concat ((scalar*){v}, tracers);
scalar * dtrcrs = list_concat ((scalar*){dv}, dsl);
for (scalar s in trcrs) {
}
s.prolongation = refine_vert5;
s.restriction  = restriction_vert;
foreach_dimension() {
}
}

The velocity tendency on vertices also requires boundary conditions

  foreach_dimension() {
dv.x.prolongation = refine_vert5;
dv.x.restriction = restriction_vert;
}

The face velocity field u is interpolated to the vertex-point values stored in v.

  foreach() {
foreach_dimension()
v.x[] = FACE_TO_VERTEX_4(u.x);
}
boundary ((scalar*){v});

We use a compact fourth-order upwind scheme to compute the gradients of all trcrs.

  compact_upwind (trcrs, grads, v);

The tendency is computed for each vertex;

  foreach() {
scalar ds; vector dsd;
for (ds, dsd in dtrcrs, grads) {
ds[] = 0;
foreach_dimension()
ds[] -= v.x[]*dsd.x[];
}
// Viscous term:
if (nu.i) {
foreach_dimension() {
scalar s = v.x, ds = dv.x;
foreach_dimension()
ds[] += nu[]*D2SDX2/sq(Delta);
}
}
// Diffusion term
if (kappa.i) {
scalar s, ds;
for (s, ds in tracers, dsl) {
foreach_dimension()
ds[] += kappa[]*D2SDX2/sq(Delta);
}
}
// Acceleration term:
if (a.x.i)
foreach_dimension()
dv.x[] += a.x[];
}

The intermediate tendency for the velocity components needs to be re-interpolated to face-averaged values.

  boundary ((scalar*){dv});
foreach_face()
du.x[] = VERTEX_TO_FACE_4(dv.x);
// cleanup
free (dtrcrs); free (trcrs);
}

## Time integration

Chorin’s operator-splitting method is employed. Furthermore, we keep track of the worst multigrid stratistics for all stages in mgp

void Navier_Stokes (face vector u, face vector du, scalar * sl, scalar * dsl) {
boundary_flux ({du});
mgstats mgt = project (du, p, dt = dt);
mgp.i      = max(mgp.i     , mgt.i);
mgp.nrelax = max(mgp.nrelax, mgt.nrelax);
mgp.resa   = max(mgp.resa  , mgt.resa);
mgp.resb   = max(mgp.resb  , mgt.resb);
mgp.sum    = max(mgp.sum   , mgt.sum);
}

In order to prevent the accumulation of the divergence’ residuals, the solution is projected after each itegration step.

event advance (i++, last) {
mgp = (mgstats){0}; // reset
A_Time_Step (dt, Navier_Stokes);
mgp2 = project (u, p2);
}

// Clean up tracer tendency
event rm_dfl (t = end) {
delete (dsl); free (dsl); dsl = NULL;
}

## Utilities

Utilities include,

• a function that computes a 2nd-order-accurate estimate of the vorticity (in the z-direction at cell centres.

• A log event prototype

#include "utils.h"

### A Wavelet-based grid-adaptation helper function

It can help to reduce the likelyhood of many small/narrow high-resolution islands.

#if (TREE)
#endif

The log event;

event logger (i++) {
fprintf (stderr, "%d %g %d %d %d %d %ld %d\n", i, t, mgp.i,
mgp.nrelax, mgp2.i, mgp2.nrelax, grid->tn, grid->maxdepth);
}

Funtions to convert between face and centered fields.

void vector_to_face (vector uc) {
foreach_face()
u.x[] = (-uc.x[-2] + 7*(uc.x[-1] + uc.x[]) - uc.x)/12.;
boundary ((scalar *){u});
}

void face_to_vector (vector uc) {
foreach_dimension()
uc.x.prolongation = refine_4th;
foreach() {
foreach_dimension()
uc.x[] = (-u.x[-1] + 13.*(u.x[] + u.x) - u.x)/24.;
}
#if NOSLIP_TOP
uc.t[top] = dirichlet_top4 (0);
uc.n[top] = dirichlet_top4 (0);
#endif
boundary ((scalar*){uc});
}

void vorticityf (face vector u, scalar omega) {
vector uc[];
face_to_vector (uc);
foreach() {
omega[] = ((8*(uc.y - uc.y[-1]) + uc.y[-2] - uc.y) -
(8*(uc.x[0,1] - uc.x[0,-1]) + uc.x[0,-2] - uc.x[0,2]))/(12.*Delta);
}
omega.prolongation = refine_4th;
boundary ((scalar*){omega});
}

#if dimension == 3
void vorticityf3 (face vector u, vector omega) {
vector uc[];
face_to_vector (uc);
foreach() {
foreach_dimension()
omega.x[] = ((8*(uc.z[0,1] - uc.z[0,-1]) + uc.z[0,-2] - uc.z[0,2]) -
(8*(uc.y[0,0,1] - uc.y[0,0,-1]) + uc.y[0,0,-2] - uc.y[0,0,2]))/(12*Delta);
}
foreach_dimension()
omega.x.prolongation = refine_5th;
boundary ((scalar*){omega});
}
#endif

## To do

• Proper box boundaries (stratified flows)
• Three dimensional simulations
• Critical evaluation