src/test/porous1.c

    Stokes flow through a complex porous medium, randomly refined

    This is similar to porous.c but tougher since the mesh refinement is now completely arbitrary (i.e. independent from the solution). This further tests robustness of the treatment of arbitrary embedded boundaries, with arbitrary levels of refinement.

    set size ratio -1
plot 'inter' w l t ''
    Randomly refined mesh and embedded boundary (script)

    Randomly refined mesh and embedded boundary (script)

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

int maxlevel = 5;

    The porous medium is the same as in porous.c.

    void porous (scalar cs, face vector fs)
{
  int ns = 800;
  double xc[ns], yc[ns], R[ns];
  srand (0);
  for (int i = 0; i < ns; i++)
    xc[i] = 0.5*noise(), yc[i] = 0.5*noise(), R[i] = 0.01 + 0.02*fabs(noise());

    Once we have defined the random centers and radii, we can compute the levelset function \phi representing the embedded boundary.

      vertex scalar phi[];
  foreach_vertex() {
    phi[] = HUGE;

    Since the medium is periodic, we need to take into account all the disk images using periodic symmetries.

        for (double xp = -L0; xp <= L0; xp += L0)
      for (double yp = -L0; yp <= L0; yp += L0)
	for (int i = 0; i < ns; i++)
	  phi[] = intersection (phi[], (sq(x + xp - xc[i]) +
					sq(y + yp - yc[i]) - sq(R[i])));
    phi[] = -phi[];
  }
  fractions (phi, cs, fs);
  fractions_cleanup (cs, fs);
}

    The domain is the periodic unit square centered on the origin.

    int main()
{

    Space and time are dimensionless. This is necessary to be able to use the ‘mu = fm’ trick.

      size (1. [0]);
  
  origin (-0.5, -0.5);
  periodic (right);
  periodic (top);

    We turn off the advection term. The choice of the maximum timestep and of the tolerance on the Poisson and viscous solves is not trivial. This was adjusted by trial and error to minimize (possibly) splitting errors and optimize convergence speed.

      stokes = true;
  DT = 2e-5 [0];
#if 1
  TOLERANCE = HUGE;
  NITERMIN = 2;
#else
  TOLERANCE = 1e-3;
  NITERMIN = 2;
#endif
  N = 1 << maxlevel;

  run();
}

scalar un[];

event init (t = 0)
{

    The mesh is refined randomly along the embedded boundary.

      porous (cs, fs);
  refine (cs[] && cs[] < 1. &&
	  level < 6 && (int)(100.*rand()/(double)RAND_MAX) == 0);
  porous (cs, fs);
  refine (cs[] && cs[] < 1. &&
	  level < 7 && (int)(200.*rand()/(double)RAND_MAX) == 0);
  porous (cs, fs);
  refine (cs[] && cs[] < 1. &&
	  level < 8 && (int)(400.*rand()/(double)RAND_MAX) == 0);

    We define the porous embedded geometry.

      porous (cs, fs);
  dump();

  FILE * fp = fopen ("inter", "w");
  output_cells (fp);
  output_facets (cs, fp, fs);
  foreach()
    fprintf (fp, "cs %g %g %g\n", x, y, cs[]);
  foreach_face()
    fprintf (fp, "fs %g %g %g\n", x, y, fs.x[]);
  fclose (fp);

    The gravity vector is aligned with the channel and viscosity is unity.

      const face vector g[] = {1.,0.};
  a = g;
  mu = fm;

    The boundary condition is zero velocity on the embedded boundary.

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

    We initialize the reference velocity.

      foreach()
    un[] = u.x[];
}

    We check for convergence of the solution.

    event logfile (i++; i <= 400)
{
  double avg = normf(u.x).avg, du = change (u.x, un)/(avg + SEPS);
  fprintf (stderr, "%d %d %d %d %d %d %d %d %.3g %.3g %.3g %.3g %.3g\n",
	   maxlevel, i,
	   mgp.i, mgp.nrelax, mgp.minlevel,
	   mgu.i, mgu.nrelax, mgu.minlevel,
	   du, mgp.resa*dt, mgu.resa, statsf(u.x).sum, normf(p).max);
}
    reset
set xlabel 'Iterations'
set logscale y
set ytics format '%.0e'
set key bottom left
plot '../porous1.ref' u 2:9 w l t '', '' u 2:10 w l t '', \
    '' u 2:11 w l t '', '' u 2:12 w l t '', '' u 2:13 w l t '', \
    'log' u 2:9 w p t 'du', '' u 2:10 w p t 'resp', \
    '' u 2:11 w p t 'resu', '' u 2:12 w p t 'u.x.sum', '' u 2:13 w p t 'p.max'
    Convergence history (script)

    Convergence history (script)

    See also