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 an 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 phi[];
      double eps = L0/(1 << 7)/1000.;
      foreach_vertex()
        phi[] = union (y - L0/4. + eps, - L0/4 + eps - y);
      boundary ({phi});
      fractions (phi, cs, fs);
    
      mu = fm;

    The boundary condition is zero velocity on the embedded boundary.

      u.n[embed] = dirichlet(0);
      u.t[embed] = dirichlet(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 %g\n", y, u.y[], p[], g.y[], cs[]);
      printf ("\n");
    }
    
    event profile (t = end)
    {
      p.nodump = false;
      dump();
      assert (normf(u.y).max < 1e-3);
    }