/** # 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. ~~~gnuplot Randomly refined mesh and embedded boundary set size ratio -1 plot 'inter' w l t '' ~~~ */ #include "embed.h" #include "navier-stokes/centered.h" #include "view.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[]; } boundary ({phi}); fractions (phi, cs, fs); fractions_cleanup (cs, fs); } /** The domain is the periodic unit square centered on the origin. */ int main() { 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; #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); } /** ~~~gnuplot Convergence history 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' ~~~ ## See also * [Stokes flow through a complex porous medium](porous.c) */