src/test/uf.c

Stability of the embedded face velocity interpolation

If a “naive” face interpolation is used (by the face_value() macro) to compute face velocities in the centered Navier–Stokes solver, an instability can appear, due to the amplifications of velocity perturbations by the third-order Dirichlet interpolation used to compute viscous fluxes.

This problem is avoided by using a embedded-fraction-weighted interpolation of the face velocities (see /src/embed.h).

#include "embed.h"
#include "navier-stokes/centered.h"
#include "view.h"

int main()
{
  origin (-0.5, -0.5);
  periodic (right);
  periodic (top);
  
  stokes = true;
  DT = 2e-5;
  TOLERANCE = HUGE;
  NITERMIN = 10;
  N = 32;

  run();
}

event init (t = 0)
{
  vertex scalar φ[];
  double eps = L0/(1 << 7)/1000.;
  foreach_vertex()
    φ[] = union (y - L0/4. + eps, - L0/4 + eps - y);
  fractions (φ, cs, fs);

  μ = fm;

The boundary condition is zero velocity on the embedded boundary.

  u.n[embed] = dirichlet_embed(0);
  u.t[embed] = dirichlet_embed(0);

  foreach()
    u.y[] = 1.;
  boundary ((scalar *){u});
}

event logfile (i++; i <= 100)
{
  fprintf (ferr, "%d %d %d %d %d %d %d %.3g %.3g %.3g %.3g\n",
	   i,
	   mgp.i, mgp.nrelax, mgp.minlevel,
	   mgu.i, mgu.nrelax, mgu.minlevel,
	   mgp.resa*dt, mgu.resa, normf(u.y).max, normf(p).max);

  foreach()
    if (x < -L0/2. + L0/N)
      printf ("%g %g %g %g\n", y, u.y[], p[], g.y[]);
  printf ("\n");
}

event profile (t = end)
{
  assert (normf(u.y).max < 1e-8);
  
  scalar p1[];
  foreach()
    p1[] = p[];
  dump();
}