/** # 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 : $$\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 $$\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: ~~~bash qcc -g -O2 -DTRASH=1 -Wall kpz.c -o kpz ./kpz > v.out ~~~ or with make ~~~bash make kpz.tst;make kpz/plots make kpz.c.html ; open kpz.c.html ~~~ ## Results The analytical (for $v_0=0$) solution is $$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 ~~~bash 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 ~~~ which gives $h(x,t)$ plotted here for t=0 .1 .2 .3 ... .9 1.0 and \$-5