sandbox/M1EMN/Exemples/svdb.c

    Rupture de barrage/ Dam Break

    An example in simple C (not a Basilisk example). We solve here the dambreak problem with a simple scheme with Rusanov flux. \displaystyle \dfrac{U_i^{n+1}-U_i^{n}}{\Delta t}+\dfrac{F^n_{i+1/2}-F^n_{i-1/2}}{\Delta x}=0,

    #include <stdio.h>
    #include <stdlib.h> 
    #include <math.h> 
    #include <string.h>
    
      double*x=NULL,*h=NULL,*u=NULL,*Q=NULL;
      double t,dt,tmax,dx,Z; 
      int nx;

    definition of the fluxes, here a simple Rusanov \displaystyle F_i=f(U_{i-1},U_i)= \dfrac{F(U_{i-1})+F(U_i)}{2}-c\dfrac{U_{i}-U_{i-1}}{2}, which is coded as \displaystyle f(U_G,U_D)= \dfrac{F(U_G)+F(U_D)}{2}-c\dfrac{U_D-U_G}{2} with c=\sup\limits_{U=U_G,U_D}({\sup\limits_{j\in\{1,2\}}} |\lambda_j(U)|), where \lambda_1(U) and \lambda_2(U) are the eigen values.

    A note on the number of points:

    D refeers to Domain, we compute the values for the nx points in the domain, but we need the index 0 and nx+1. Those two extra localisations are the ghost cells, G refers to Ghost.

     G   D   D   D   D    D    G
     |   |   |   |   |    |   |
     0 | 1 | 2 | 3 |...| nx| nx+1
     --|---|---|---|---|---|
        x1  x2             L
     <-------------------->
     <-->
     
     nx cells, +2 ghost cells, (nx+2) data
     
     first cell, index 1,  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`
     
     |           |
     |->Fi       |->Fi+1
     |  thetai   |
     xi+Delta/2
     xi   xi+Delta/2

    Fi is the flux across interface between hi-1 and hi (sometimes Fi-1/2)

    Fi+1 is the flux across interface between hi and hi+1 (sometimes Fi+1/2)

    double FR1(double ug,double ud,double hg,double hd)
     { double c;
       c=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
       return (hg*ug+hd*ud)*0.5-c*(hd-hg)*0.5;
     }
    double FR2(double ug,double ud,double hg,double hd)
     { double c;
       c=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
       return (ug*ug*hg + hg*hg/2. + ud*ud*hd + hd*hd/2.)*0.5 - c*(hd*ud-hg*ug)*0.5;
     }
      
    /*     -------------------------------------------------    */        
    int main (int argc, const char *argv[]) {
        int  i,it=0;
        FILE *g;
        double*fp=NULL,*fd=NULL,*un=NULL,*hn=NULL;
        double  q;
        
    // parametres --------------------------------------------------------------  
    dt=0.005;
    tmax=10;
    dx=0.125;
    nx=320;
    t=0;    
      
      fprintf(stderr,"  ------------------\\                          \n"); 
      fprintf(stderr,"                   <--\\                        \n");   
      fprintf(stderr,"                        \\--->                  \n");   
      fprintf(stderr,"  ________________________\\____________________\n");   
    /* ----------------------------------------------------------------------*/      
          x= (double*)calloc(nx+2,sizeof(double));
          h= (double*)calloc(nx+2,sizeof(double));
          Q= (double*)calloc(nx+2,sizeof(double));
          u= (double*)calloc(nx+2,sizeof(double));
          fp=(double*)calloc(nx+2,sizeof(double));
          fd=(double*)calloc(nx+2,sizeof(double));
          un=(double*)calloc(nx+2,sizeof(double));
          hn=(double*)calloc(nx+2,sizeof(double));
          
     fprintf(stderr,"tmax= %lf \n",tmax);
          
    // initialisation cond init ----------------------------      
      for(i=0;i<=nx+1;i++)
        {  x[i]=(i-nx/2)*dx;
       //dambreak
         Z=0;
         h[i]=1*( x[i]<0);
         u[i]=0;
         Q[i]=u[i]*h[i];
         }
           
    // initialisation du fichier de sortie ----------------------      
          g = fopen("solxhQt.OUT", "w");
          fclose(g); 
                      
        while(t<tmax){   // boucle en temps
          t=t+dt;
          it=it+1;      
     
       for(i=1;i<=nx+1;i++)
        {    
            fp[i]=FR1(u[i-1],u[i],h[i-1],h[i]);
            fd[i]=FR2(u[i-1],u[i],h[i-1],h[i]);  
        }
       
       for(i=1;i<nx+1;i++)
        {
          hn[i]=h[i]- dt*(fp[i+1]-fp[i])/dx;   //conservation de la masse
      if(h[i]>0.){                             //conservation qunatité de mouvement
          q=h[i]*u[i]-dt*(fd[i+1]-fd[i])/dx ;  
          un[i]=q/hn[i];}
          else{
          un[i]=0.;}
        }
        
        // flux nul en entree sortie 
        hn[0]=hn[1];
        un[0]=un[1];  
        hn[nx+1]=hn[nx];
        un[nx+1]=un[nx];
       
       //swap      
      for(i=0;i<=nx;i++)
        { 
          h[i]=hn[i];
          u[i]=un[i];
          Q[i]=hn[i]*un[i];
        }
    
      if(it%100==0){
        /* Saving the fields */ 
      g = fopen("solxhQt.OUT", "a");
      for (i=0; i<=nx;i++)
      {      fprintf(g,"%lf %lf %lf %lf \n",x[i],h[i],Q[i],t);}
       fprintf(g,"\n");
       fclose(g);
      }
    #ifdef gnuX
    if(it%10==0){
      double  y1=-.1,y2=1.25;
      printf("t=%lf\n",t);
      printf("set xlabel \"x\" ; set title \"t=%lf \"\n",t);
         y1=-.1,y2=1.25;  
      printf("set label \"+\" at 0,.296296296296296 \n"); 
      printf("h(x)=(((x-0)<-t)+((x-0)>-t)*(2./3*(1-(x-0)/(2*t)))**2)*(((x-0)<2*t)) \n");
    
       
      printf("p[%lf:%lf][%lf:%lf] h(x) t'h exact','-' u 1:2 t'Q'  w l,'-'  t'h' w lp,'-'  t'Z' w l linec -1\n ",
             -nx*dx/2,nx*dx/2,y1,y2); 
        for (i=0; i<=nx;i++)
        {
        printf("%lf %lf \n",x[i],Q[i]);}
        printf("e \n");
        for (i=0; i<=nx;i++)
        {
        printf("%lf %lf \n",x[i],h[i]+Z);}
        printf("e \n");
        for (i=0; i<=nx;i++)
        {
          printf("%lf %lf \n",x[i],Z);}
        printf("e \n");}
    #endif   
          }
          
         free(x);
         free(fp);
         free(fd); 
         free(un);
         free(hn);  
        return 0;
    }

    Run avec cc

    Programme en C simple fichier (en cours on a regardé un fichier plus compliqué). Cas de rupture de barrage flux Rusanov (essayer HLC) compilation on le compile et on crée l’exécutable db

    cc -O3 -ffast-math -std=c99 -lm svdb.c -o db

    pour lancer le programme:

     ./db 

    Alternatively, we can compile with gnuX option which allows to pipe in gnuplot

    cc -O3 -ffast-math -std=c99 -lm svdb.c -o db -DgnuX

    to run the code

     ./db | gnuplot -persist

    Run avec make

    Avec le makefile

     make svdb.tst
     make svdb/plots
     make svdb.c.html

    Results

    un fichier appelé solxhQt.OUT avec x h Q et t est créé, pour le tracer, on lance gnuplot:

     set ylabel "t"; set xlabel "x"; set hidden3d;
     sp[-20:20][0:10][0:2] 'solxhQt.OUT' u 1:4:2 not w l linec 3

    this gives a view in 3D with h function of x at the different t.

      set ylabel "t"; set xlabel "x";  set hidden3d;
      splot  'solxhQt.OUT' u 1:4:2 not w l
    (script)

    (script)

    An other plot, comparison exact and computed at two times

      set ylabel "h"; set xlabel "x";  
      set label "+" at 0,.296296296296296 
      h(x,t)=(((x-0)<-t)+((x-0)>-t)*(2./3*(1-(x-0)/(2*t)))**2)*(((x-0)<2*t))
      p  'solxhQt.OUT' u 1:(($4==5)||($4==10)? $2 :NaN)w p, h(x,5),h(x,10)
    (script)

    (script)

    Go further

    1. First you can change the initial boundary condition to mimick a dam break in a river
     h[i]=0.1+.9*( x[i]<0);

    observe the formation of a shock

    1. You can add some viscosity in the code to test the example of viscous collapse

    first change for a heap:

    h[i]=1*( fabs(x[i])<1);

    then add viscosity \frac{\partial u}{\partial t} = - 3*u /h with implicit discretization u(t+\Delta t,x)- u(t,x) = - 3 * u(t+\Delta t,x) \Delta t /h(x,t)^2 hence \displaystyle u(t+\Delta t,x)= \frac{u(t,x)}{1+ 3 \Delta t /h(t,x)^2}

    for(i=1;i<nx;i++) { if(h[i]>0.) un[i]= un[i]/(1+3*dt/h[i]/h[i]); } 

    check the selsimilar variable:

    p[-5:5]'solxhQt.OUT' u ($1/$4**.2):($4>10?$2*$4**.2:0) w l
    1. You can put turbulent friction (this is the Dressler Problem See Chanson page 357 and see dam break Dressler )

    2. You may change the boundary conditions,

    here BC are output BC, cells u[0]and u[nx+1] are ghost cells

    no flux, we copy the values

     hn[0]=hn[1];
     un[0]=un[1];
     hn[nx+1]=hn[nx];
     un[nx+1]=un[nx];

    It is changed to have no penetration at the walls: a zero velocity at the final face

    dirichlet 0: the final face value is the linear extrapolation of the nx and the nx+1 \displaystyle u= \frac{u_{nx}+u_{nx+1}}{2}

     hn[0]=hn[1];
     un[0]=-un[1];
     hn[nx+1]=hn[nx];
     un[nx+1]=-un[nx];

    Links

    the same example of dam break with Basilisk

    Biblio

    • A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows Emmanuel Audusse, Franccois Bouchut, Marie-Odile Bristeau, Rupert Klein and Benoıt Perthame
    • Olivier Delestre Thèse
    • cours PYL
    • cours PYL sur numérique
    • Hubert Chanson The hydraulics of open channel flow: an introduction ; basic principles" Elsevier Butterworth-Heinemann, Amsterdam, The Netherlands. 2004

    // version init mai 2012 version mars 2013 // PYL