nsf4t.h

A Navier-Stokes equations solver with tracers

A Navier-Stokes equations solver with tracers

drawing

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;
    u.x.interpolant = interpolant_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[2])/12. + 4.*(s[1] + 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);
  vector * grads = NULL;                
  for (scalar s in trcrs) {
    vector dsd = new_vector ("gradient");
    grads = vectors_add (grads, dsd);
  }
  vector grad; scalar s;
  for (grad, s in grads, trcrs) {
    s.prolongation = refine_vert5;
    s.restriction  = restriction_vert;
    grad.t[top] = freegrad;
    grad.n[top] = freegrad;
    foreach_dimension() {
      grad.x.prolongation = refine_vert5;
      grad.x.restriction = restriction_vert;
    }
  }

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() {
    // Advection:
    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 
  delete ((scalar*)grads); free (grads);
  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) {
  adv_diff (du, 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);
}

event adapt (i++, last) ;

// 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 wavelet-based grid-adaptation function
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)
#include "adapt_field.h"
#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[1])/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[1]) - u.x[2])/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[1] - uc.y[-1]) + uc.y[-2] - uc.y[2]) -
                 (8*(uc.x[0,1] - uc.x[0,-1]) + uc.x[0,-2] - uc.x[0,2]))/(12.*Delta);
  }
  omega[top] = freegrad;
  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