layer.c

Singularities in the thin 2D film equation
- Numerical method
- Time integration
Outputs

Singularities in the thin 2D film equation

In a coming work of Dallaston et al., 2021 (under revision in JFM) it is studied the motion of a thin film of liquid under surface forces with origin in the solid substrate.

The nonlinear lubrication model used in the work is formed by the following equations:

\displaystyle h_t = \nabla \cdot \left[h^n \nabla p - h^{n-m-1} \nabla h \right] \displaystyle p = -\nabla^2 h

where h is the height of the film and p is the 2D pressure distribution. n and m are arbitrary exponents.

Numerical method

The discretization is of second order in space and of first order in time (although the code is prepared to increase it to second order). The resulting nonlinear set of equations is solved iteratively with a Newton-Raphson (N-R) scheme,

\displaystyle p = p^* + \delta p \quad h = h^* + \delta h

where (p^*, h^*), (\delta p, \delta h) and (p, h) are the trial, correction and converged solution, respectively.

#include "poisson.h"
#include "run.h"
#include "view.h"

An struct is created to facilitate the back and forth of data through the different functions.

struct Layer {
  scalar f;
  double dt;
  double m;
  double n;
  face vector D;
  face vector beta;
  face vector gamma;
  scalar lambda;
  double tolerance;
};

The linearization of the equations leads to the following set,

\displaystyle p + h_{xx} + h_{yy} = 0 \displaystyle (\alpha p_x)_x + (\alpha p_y)_y + (\beta h_x)_x + (\beta h_y)_y + (\gamma_x h)_x + (\gamma_y h)_y = h_t + (\gamma_x h^*)_x + (\gamma_y h^*)_y

with \displaystyle \alpha = (h^*)^n, \: \beta = -(h^*)^{n-m-1} \: \text{and} \: \gamma_{x,y} = n(h^*)^{n-1} p^*_{x,y} -(n-m-1)(h^*)^{n-m-2} h^*_{x,y}.

which is solved with the iterative Multigrid (MG) scheme integrated in Basilisk. The MG scheme requires the computation of the residual of the linear equation…

static double residual_layer (scalar * al, scalar * bl,
			      scalar * resl, void * data)
{
  scalar h = al[0], b = bl[0], res = resl[0];
  struct Layer * p = (struct Layer *) data;
  (const) face vector alpha = p->D;
  (const) face vector beta = p->beta;
  (const) face vector gamma = p->gamma;
  (const) scalar lambda = p->lambda;
  double maxres = 0.;
  
  face vector g[], hf[];

  foreach_face()
    hf.x[] = face_gradient_x (h, 0);
  boundary_flux ({hf});

  scalar pres[];
  foreach() {
    pres[] = 0.;
    foreach_dimension()
      pres[] -= (hf.x[1]- hf.x[0])/Delta;
  }
  boundary({pres});
  
  foreach_face() {
    g.x[] = alpha.x[]*face_gradient_x (pres, 0);
    hf.x[] = beta.x[]*face_gradient_x (h, 0) +
      gamma.x[]*(h[]+h[-1])/2.;
  }
  boundary_flux ({g, hf});
  foreach (reduction(max:maxres)) {
    res[] = b[] - lambda[]*h[];
    foreach_dimension()
      res[] -= (g.x[1] - g.x[])/Delta + (hf.x[1] - hf.x[])/Delta;

    if (fabs (res[]) > maxres)
      maxres = fabs (res[]);
  }
  boundary (resl);
  return maxres;
}

…together with the subsequent relaxation toward the solution of the linear system.

static void relax_layer (scalar * al, scalar * bl, int l, void * data)
{
  scalar h = al[0], b = bl[0];
  struct Layer * p = (struct Layer *) data;
  (const) face vector alpha = p->D;
  (const) face vector beta = p->beta;
  (const) face vector gamma = p->gamma;
  (const) scalar lambda = p->lambda;

  
#if JACOBI
  scalar c[];
#else
  scalar c = h;
#endif

  foreach_level_or_leaf (l) {
    double n =  sq(sq(Delta))*b[], d =  lambda[]*sq(sq(Delta));
    d -= 5*(alpha.x[] + alpha.x[1] + alpha.y[] + alpha.y[1]);
    d += 0.5*Delta*sq(Delta)*(gamma.x[1]-gamma.x[] + gamma.y[1]-gamma.y[]);
    foreach_dimension() {
      n -= (beta.x[1]*h[1] + beta.x[]*h[-1])*sq(Delta);
      d -= (beta.x[1] + beta.x[])*sq(Delta);
    }

    n -= alpha.x[]*(-h[-2] - h[-1,-1] + 5*h[-1] - h[-1,1] + h[0,-1] + h[0,1] + h[1]) +
      alpha.x[1]*(h[-1] + h[0,-1] + h[0,1] - h[1,-1] + 5*h[1] - h[1,1] - h[2]) +
      alpha.y[]*(-h[-1,-1] + h[-1] - h[0,-2] + 5*h[0,-1] + h[0,1] - h[1,-1] + h[1]) +
      alpha.y[1]*(h[-1] - h[-1,1] + h[0,-1] + 5*h[0,1] - h[0,2] + h[1] - h[1,1]);

    n -= 0.5*Delta*sq(Delta)*(h[1]*gamma.x[1] - h[-1]*gamma.x[] +
			      h[0,1]*gamma.y[1] - h[0,-1]*gamma.y[]); 
   
    c[] = n/d;
  }
  
#if JACOBI
  foreach_level_or_leaf (l)
    h[] = (h[] + 2.*c[])/3.;
#endif
  
#if TRASH
  scalar a1[];
  foreach_level_or_leaf (l)
    a1[] = h[];
  trash ({h});
  foreach_level_or_leaf (l)
    h[] = a1[];
#endif
}

h and p are solved in each timestep with the interface below.

trace
mgstats layer (struct Layer p)
{
  if (p.dt == 0.) {
    mgstats s = {0};
    return s;
  }

  scalar f = p.f;
  scalar b[], fold[], fa[];
  const scalar lambda[] = - 1./p.dt;
  double m = p.m, n = p.n;
  face vector D[], beta[], gamma[];

The first trial solution is that of the previous time step, h^*=h(t- \delta t).

  foreach() 
    fold[] = f[];
  //  boundary({fa});

  double maxres;
  mgstats s;
  int niter = 0, nitermax = 3;

  double defaultol = TOLERANCE;
  if (p.tolerance)
    TOLERANCE = p.tolerance;

In each iteration step of the N-R scheme, the auxiliary functions \alpha(h^*), \beta(h^*) and \gamma_{x,y}(h^*) are first computed.

  do {
    foreach_face()
      beta.x[] = face_gradient_x (f, 0);
    boundary_flux ({beta});
    
    scalar pres[];
    foreach() {
      pres[] = 0.;
      foreach_dimension()
	pres[] -= (beta.x[1]- beta.x[0])/Delta;
    }
    boundary({pres});
    
    face vector g[];
    foreach_face () {
      double ff = (f[]+f[-1])/2.;
      D.x[] = pow(ff, n);
      beta.x[] = -pow(ff, n-m-1.);
      gamma.x[] = n*pow(ff, n-1.)*(pres[]-pres[-1])/Delta
	-(n-m-1.)*pow(ff, n-m-2.)*(f[]-f[-1])/Delta;
      g.x[] = gamma.x[]*ff;
    }
    boundary_flux({g});
  
    foreach() {
      b[] = -fold[]/p.dt;
      fa[] = f[];
      foreach_dimension()
	b[] += (g.x[1] - g.x[])/Delta;
    }
    
    restriction ({beta, D, lambda, gamma});
    p.lambda = lambda;
    p.D = D;
    p.beta = beta;
    p.gamma = gamma;

The MG linear solver is launched next,

    s = mg_solve ({f}, {b}, residual_layer, relax_layer, &p);

The convergence of the N-R is checked below. If the convergence is not reached and the number iteration is below a threshold value, a new iteration proceed.

    maxres = 0; 
    foreach (reduction(max:maxres)) {
      double res = fabs(f[]-fa[]);
      if (res > maxres)
	maxres = res;
    }
    boundary({fa});
    niter++;
    //    printf("maxres = %g, niter = %d\n", maxres, niter);
  } while(maxres > 1e-8 && niter < nitermax);

We restore the default tolerance.

  if (p.tolerance)
    TOLERANCE = defaultol;

  return s;
}

Time integration

The variables we wish follow are: the height of the film h, the inverse of the heigth hi and the level of refinement l. The problem is fully symmetric and the the domain is a [-1.5:1.5]x[-1.5:1.5] square.

#define T2OR 0 // active second order in time integration
#define pert1 0.05
#define pert2 0.03
#define HMIN A*(1-pert1)*(1-pert2)

scalar h[], hi[], l[];
double A = 0.05, n = 3., m = 1., DTMAX = 1e-3;
double dtold , hmin, hmin1, xmin = 0.5, ymin = 0.5;
int levelm;
int LEVEL= 12, lrejct = 0, LEVEL_MIN = 7;
double dtfactor = 0.995;

#if T2OR
scalar h2[], h1[], err[];
double errorstep = 0., DTMIN = 1e-6;
#endif

initial timestep can be given by screen.

int main(int argc, char * argv[]) {
  if (argc > 1)
    DTMAX = atof (argv[1]);

  periodic(right);
  periodic(top);
  origin(0.,0.);
  init_grid (1 << LEVEL_MIN);
  L0 = 3;
  run();
}

Initially, a small sinusoidal pertubation is seeded.

event init (i = 0) {
#if T2OR
  dtold = DTMIN;
  h2.nodump = true;
  h1.nodump = true;
  err.nodump = true;
#endif
  hmin = HMIN;
  hmin1 = HMIN + 1e-9;
  foreach() 
    h[] = A*(1.-pert1*cos(2*M_PI*(x/L0-0.5)))*
    (1.-pert2*cos(2*M_PI*(y/L0-0.5)));
  boundary ({h});
}

Several timestep control procedures have been implemented. For example, a time step control based in the error in time integration is computed below. Has only sense if second order in time is used.

#if T2OR
double ckstp (double dtmin, double dtmax, double dt, double error)
{
         
double dtgrow= 1.189207115;     
    
if ( error<1e-4 )
  {
    lrejct = lrejct + 1;
    if( (lrejct > 10) & (error < 2.5e-5))
      dt = dt * dtgrow;
  }
 else  // optional
   {
     dt=0.5*dt;
     lrejct = 0;
   }
 return ( min( dtmax, max(dtmin, dt))) ;
}
#endif

Or the time step is proportional to the minimum heigth.

double mastep ( double hmin) {
  return 1.e-4*hmin;
}

Or the time step depends on the maximum level of refinement. This criteria although implemented has not been really used.

double steplevel ( double levelmax) {
return 1.e-4/pow(2.,levelmax);
}

In this case, the time step come after a linear interpolation. (Method proposed by M.A. Herrada)

double lininterp (double h1, double h2, double dt,
		  double alpha, double dtmax) {
  printf("%g %g %g %g %g \n",
	 h1,h2,dt,alpha, (alpha*h2-h1)/(h2-h1)*dt -dt);
  double step = min(((alpha*h2-h1)/(h2-h1)*dt -dt), dtmax);
    return (step <= 0. ? dtmax : step); 
}


event integration (i++) {
#if T2OR  
    dt = dtnext (ckstp (DTMIN, DTMAX, dtold, errorstep));
    dtold = dt;
    foreach() {
      h2[] = h[];
      h1[] = h[];
    }
  boundary({h2, h1});
  layer (h1, dt, m, n);  
  layer (h2, dt/2., m, n);  
  layer (h2, dt/2., m, n);

  double errormax = 0;
  foreach() {
    h[] = 2.*h2[] -h1[];
    err[] = fabs(h2[]-h1[]);
    if (err[] > errormax)
      errormax = err[];
  }
  boundary({h2, h1});
  stats s = statsf(err);
  errorstep = sqrt(s.sum);
#else
  dtold = dt;
  dt = dtnext (lininterp (hmin1, hmin, dtold, dtfactor, DTMAX));
  layer (h, dt, m, n);
#endif

Once h[] is computed, we calculate its inverse and the spatial distribution of the level of refinement. Beside, position and value of the minimum hmin is found.

  double minres = 1000; 
  levelm = -1;
  hmin1 = hmin;
  foreach (reduction(min:minres)) {
    hi[] = 1/h[];
    l[] = level;
    double res = fabs(h[]);
      if(levelm < level)
	levelm = level;
    if (res < minres) {
      xmin = x;
      ymin = y;
      hmin = res;
      minres = res;
    }
  }
}

In each timestep the grid is adapted accordingly to the error estimation of the inverse of h.

#if TREE
event adapt (i++) {
  adapt_wavelet ({hi}, (double[]){0.2}, LEVEL, LEVEL_MIN);
}
#endif

Outputs

initial initial

The plots of figure 12 of the paper are generated.

event safe_figure (t = {0, 3.8944}) {
  view (fov = 20., tx = -0.5, ty = -0.5,
	width = 600, height = 600, samples = 1);
 squares("hi", linear =1, spread = -1);
 char name[80];
 sprintf(name,"hi_n%g_m%g_ho%g_t%g.png", n, m, A, t);
 save (name);
}

We could save also the time evolution of the minimum h.

#if 0
event output (i += 5) {
  static FILE * fp = fopen("hmin", "w");
  if (i==0)
    fprintf (fp, "#m = %g n = %g ho = %g LEVEL = %d LEVEL_MIN = %d dt = %g alpha = %g  \n",
	     m, n, A, LEVEL, LEVEL_MIN, DTMAX, dtfactor);
  fprintf (fp, "%.15g %.15g %.15g %d\n", t, hmin, dt, levelm);
  fflush(fp);
}
#endif