sandbox/M1EMN/BASIC/kpz.c

Resolution of KPZ equation in 1D

A classic model for the evolution of the profile of a growing interface is the Kardar–Parisi–Zhang (KPZ) equation : \displaystyle \frac{\partial h}{\partial t} = -v_0 \frac{\partial h}{\partial x} + \frac{\lambda}{2} (\frac{\partial h}{\partial x})^2 + \nu \frac{\partial^2 h}{\partial x^2} +v_c without here the stochastic source. v_c is a constant feeding, \lambda an effect of curvature, \nu a diffusion effect, v_0 an extra slope effect. We propose here an explicit 1D resolution of KPZ equation compared with an analytical solution.

Code

mandatory declarations:

#include "grid/cartesian1D.h"
#include "run.h"

definition of the height of interface its O(Delta) derivative and it O(Delta^2) derivative, time step

scalar h[];
scalar hp[];
scalar hp2[]; 
double dt;

Main with definition of parameters

int main() {
  L0 = 10.;
  X0 = -L0/2;
  N = 200;
#define nu 1.    
  DT = (L0/N)*(L0/N)/nu/10 ;
#define lambda 1. 
#define EPS 0.1  
#define v0 1.
#define vc 2.
  run();
}

initial elevation: a “triangle”

event init (t = 0) {
  foreach()
    h[] =   (1-fabs(x))*(fabs(x)<1);
  boundary ({h});
  }

print data

event printdata (t += 0.1; t < 1) {
  foreach()
    fprintf (stdout, "%g %g %g\n", x, h[], t);
  fprintf (stdout, "\n\n");
}

integration

event integration (i++) {
  double dt = DT;
  scalar dh[];

finding the good next time step

  dt = dtnext (dt);

O(\Delta) derivative

  foreach()
    hp[] =  ( h[1,0] - h[0,0] )/Delta;
  boundary ({hp});

centered derivative O(\Delta^2)

  foreach()
    hp2[] =  ( h[1,0] - h[-1,0] )/2/Delta;
  boundary ({hp2});

explicit step \displaystyle \frac{ (h(x+\Delta x) - h(x))/\Delta -(h(x) - h(x-\Delta x))/\Delta }{\Delta } \simeq \frac{\partial^2 h}{\partial x^2}

  foreach()
    dh[] =  ( hp[] - hp[-1,0] )/Delta;
  boundary ({dh});

update viscous term

  foreach()
    h[] += dt*nu*dh[];
  boundary ({h});

update uniform growth

  foreach()
    h[] += dt*vc;
  boundary ({h});

update normal effect

  foreach()
    h[] += -dt*v0*hp2[];
  boundary ({h});

update curvature effect

  foreach()
    h[] += dt*(lambda/2.*hp2[]*hp2[]);
  boundary ({h});
  
}

Run

Then compile and run:

qcc  -g -O2 -DTRASH=1 -Wall  kpz.c -o kpz
./kpz > v.out

or with make

 make kpz.tst;make kpz/plots    
 make kpz.c.html ; open kpz.c.html 

Results

The analytical (for v_0=0) solution is \displaystyle h(x,t) = 2 \log \left(-\frac{1}{2} e^{-x/2} \left(e^{\frac{t}{4}+\frac{1}{2}} \text{erf}\left(\frac{t-x}{2 \sqrt{t}}\right)-e^{\frac{t}{4}+\frac{1}{2}} \text{erf}\left(\frac{t-x+1}{2 \sqrt{t}}\right)-e^{x/2} \text{erf}\left(\frac{x-1}{2 \sqrt{t}}\right)+e^{\frac{t}{4}+x+\frac{1}{2}} \text{erf}\left(\frac{t+x}{2 \sqrt{t}}\right)-e^{\frac{t}{4}+x+\frac{1}{2}} \text{erf}\left(\frac{t+x+1}{2 \sqrt{t}}\right)-e^{x/2} \text{erfc}\left(\frac{x+1}{2 \sqrt{t}}\right)-e^{x/2}\right)\right) in gnuplot type

 h(x,t)=2*log(-(-exp(x/2.) + exp(0.5 + t/4.)*erf((t - x)/(2.*sqrt(t))) -  exp(0.5 + t/4.)*erf((1 + t - x)/(2.*sqrt(t))) -  exp(x/2.)*erf((-1 + x)/(2.*sqrt(t))) +   exp(0.5 + t/4. + x)*erf((t + x)/(2.*sqrt(t))) -   exp(0.5 + t/4. + x)*erf((1 + t + x)/(2.*sqrt(t))) - exp(x/2.)*erfc((1 + x)/(2.*sqrt(t))))/(2.*exp(x/2.)))
  p'out' t'num'w lp,'' u 1:(h($1-$3,$3)+2*$3) t'exact' w l

 h(x,t)=2*log(-(-exp(x/2.) + exp(0.5 + t/4.)*erf((t - x)/(2.*sqrt(t))) -  exp(0.5 + t/4.)*erf((1 + t - x)/(2.*sqrt(t))) -  exp(x/2.)*erf((-1 + x)/(2.*sqrt(t))) +   exp(0.5 + t/4. + x)*erf((t + x)/(2.*sqrt(t))) -   exp(0.5 + t/4. + x)*erf((1 + t + x)/(2.*sqrt(t))) - exp(x/2.)*erfc((1 + x)/(2.*sqrt(t))))/(2.*exp(x/2.)))
 p'out' t'num'w lp,'' u 1:(h($1-$3,$3)+2*$3) t'exact' w l

(script)

Bibliography

M. Kardar, G. Parisi, Yi-C. Zhang, “Dynamic Scaling of Growing Interfaces” Phys. Rev. Lett. 56 (1986) 889.

M.T. Batchelor, R.V. Burne, B.I. Henryc; S.D. Watt “Deterministic KPZ model for stromatolite laminae” Physica A 282 (2000) 123–136

ready for new site 09/05/19