sandbox/M1EMN/BASIC/advecte.c

    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.tst;make advecte/plots    
 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)

    (script)

    Links

    • advecte1.c explains the notions of advection, testing the flux, coded with Basilisk

    • advecte1c.c the same than the previous one but with standard C

    • advecte.c advection with Basilisk, compares centered and decentered

    Bibliography

    • LeVeque chapitre 4, Finite Volume methods

    ready for new site 09/05/19