sandbox/M1EMN/BASIC/advecte1c.c

    Resolution of Advection equation in 1D

    it is exactly the same than advecte1.c but written with standard C

    Code

    mandatory declarations:

    #include <stdio.h>
    #include <stdlib.h> 
    #include <math.h>
    #include <string.h>

    definition of the field U, the flux F, time step

    double*x=NULL,*U=NULL,*F=NULL;
    double dt;  
    double cDelta;
    
    double L0,X0,Delta; 
    double t;
    int i,it=0;
    int N;

    Main with definition of parameters

    int main() {
      L0 = 12.;
      X0 = -L0/4;
      N = 128;    
      t=0;
      dt = (L0/N)/2;
      Delta = L0/N;

    dynamic allocation

      x= (double*)calloc(N,sizeof(double));
      U= (double*)calloc(N,sizeof(double));
      F= (double*)calloc(N,sizeof(double));

    Boundary condition

      U[0] = 0;
      U[N-1]=0; 
      
        for(i=0;i<N;i++)

    first cell between X0 and X0+Delta, centred in Delta/2 ith cell beween (i-1) Delta (left) and i Delta(right) centered in (i-1/2)Delta

    Delta=L0/N

        {  x[i]=X0+(i-1./2)*Delta;

    initial elevation: a “bump”, position

           U[i] = exp(-x[i]*x[i]);
      }

    begin the time loop

       while(t<=3){ 
         t = t + dt;
         it++;

    print data

    //event printdata (t += 1; t <=3)  
    if( (it == (int)(1/dt)) || (it == (int)(2/dt)) || (it == (int)(3/dt))  ){
     for(i=0;i<N;i++)
        fprintf (stdout, "%g %g %g  \n", x[i], U[i], t);
      fprintf (stdout, "\n\n");
    }

    flux

     for(i=1;i<N;i++)
       {
       //cDelta = Delta/dt;
       // cDelta = 0;
        cDelta = 1;
        F[i] = (U[i]+U[i-1])/2.  - cDelta *(U[i]-U[i-1])/2;
        }

    explicit step update \displaystyle U_i^{n+1}=U_i^{n} -{\Delta t} \dfrac{F(U_{i+1})-F(U_{i})}{\Delta x}

     for(i=0;i<N-1;i++)
        U[i] +=  - dt* ( F[i+1] - F[i] )/Delta;
     }
      
      free(U); 
      free(F);
      free(x);
    }

    Run

    Then compile and run:

     cc  -g -O2 -DTRASH=1 -Wall  advecte1c.c -o advecte1c ;./advecte1c > out

    or better

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

    Results

    One analytical solution is \displaystyle U(x,t) = exp(-(x-t)^2) in gnuplot type

     U(x,t)= exp(-(x-t)*(x-t))
     p'out' u ($1):($2)t'num'w l,'' u 1:(U($1,$3)) t'exact' w l
    which gives U(x,t) plotted here for t=0 1 2 3 and -3<x<9
     set xlabel "x"
     U(x,t)= exp(-(x-t)*(x-t))
     p'out' u ($1):($2)t'num'w p,'' u 1:(U($1,$3)) t'exact' w l
    (script)

    (script)

    Exercice

    Compare with the Basilisk advecte1.c Compare with advecte.c

    ready for new site 09/05/19