# Resolution of simple transport equation, centered and decentered

Explicit resolution of \displaystyle \frac{\partial h}{\partial t}+ v_0 \frac{\partial h}{\partial x}=0

with various schemes for comparison: centered derivative O(\Delta^2):

\displaystyle \frac{h(t+\Delta t,x) - h(t ,x)}{\Delta t}+ v_0 \frac{h(t,x+\Delta x) - h(t ,x-\Delta x)}{2 \Delta x}=0 it is not a good idea, it is unstable, down stream decentered derivative O(\Delta), is stable, but damped:

\displaystyle \frac{h(t+\Delta t,x) - h(t ,x)}{\Delta t}+ v_0 \frac{h(t,x) - h(t ,x-\Delta x)}{\Delta x}=0

## Code

mandatory declarations:

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

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

scalar h0[];
scalar hc[];
scalar hd[];
scalar hp[];
scalar hdp[];
scalar hp2[];
double dt;

Main with definition of parameters

int main() {
L0 = 8.;
X0 = -L0/2;
N = 128/2;
DT = (L0/N)/10 ;
#define v0 .5
run();
}

initial elevation: two “triangles”

event init (t = 0) {
foreach(){
h0[] =   -(1-fabs(x-2))*(fabs(x-2)<1) + (1-fabs(x+2))*(fabs(x+2)<1);
hc[]=h0[];
hd[]=h0[];
}
boundary ({hc,hd});
}

print data

event printdata (t += .5; t <= 2) {
foreach()
fprintf (stdout, "%g %g %g %g \n", x, hc[],hd[], t);
fprintf (stdout, "\n");
}

integration

event integration (i++) {
double dt = DT;

finding the good next time step

  dt = dtnext (dt);

O(\Delta) derivative

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

down stream decentered derivative O(\Delta), stable, but damped

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

centered derivative O(\Delta^2), it is not a good idea, it is unstable

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

update of the centered and decentered fields
\displaystyle h(t+\Delta t,x)= h(t ,x)- \Delta t v_0 \frac{h(t,x+\Delta x) - h(t ,x-\Delta x)}{2 \Delta x} and \displaystyle h(t+\Delta t,x)= h(t ,x)- \Delta t v_0 \frac{h(t,x) - h(t ,x-\Delta x)}{\Delta x}

  foreach(){
hc[] += -dt*v0*hp2[];
hd[] += -dt*v0*hdp[];
}
boundary ({hc,hd});
}

## Run

Then compile and run:

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

or better

 ln -s ../../Makefile Makefile
make advecte.c.html ; open advecte.c.html
 source ../c2html.sh advecte

## Results

The analytical (for v_0=0) solution is \displaystyle h(x,t) = h_0(x-v_0t) in gnuplot type

  v0=.5
h0(x,t)=-(1-abs(x-v0*t-2))*(abs(x-v0*t-2)<1)  + (1-abs(x-v0*t+2))*(abs(x-v0*t+2)<1)
p'out' t'cent.'w lp,'' u 1:3 t'dec.','' u 1:(h0($1,$4)) t'exact' w l

which gives h(x,t) plotted here for t=0 .5 … 2.0 and -8<x<8

 v0=.5
h0(x,t)=-(1-abs(x-v0*t-2))*(abs(x-v0*t-2)<1)  + (1-abs(x-v0*t+2))*(abs(x-v0*t+2)<1)
p'out' t'cent.'w lp,'' u 1:3 t'dec.','' u 1:(h0($1,$4)) t'exact' w l (script)