sandbox/M1EMN/Exemples/svdb.c
Rupture de barrage/ Dam Break
An example in simple C, this is not a Basilisk
example, this is the example of the lecture.
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,
definition of the fluxes, here a simple Rusanov \displaystyle F_{i-1}=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`
A note on the flux
| |
|->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)
so, in practice for the height, \displaystyle h_i^{n+1}= h_i^{n} + {\Delta t}\dfrac{F^n_{i+1/2}-F^n_{i-1/2}}{\Delta x},
for h[i]
the first component of F_{i-1} is the Fi
of the figure and is the subroutine FR1
, this gives
hn[i]=h[i]- dt*(fp[i+1]-fp[i])/dx;
for flow rate, we save the velocity and re compute the flow rate with depth:
q=h[i]*u[i]-dt*(fd[i+1]-fd[i])/dx ;
un[i]=q/hn[i];
Program/ code
you can extract this part from here
#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;
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.;}
}
// put friction here
// boundary conditions
// 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;
}
this is the end of the program. You copy and paste in your favorite editor (an ASCII editor). You save as svdb.c
in a specific directory somewhere in your hard drive.
You open two terminals (xterm
), one in wich you compile, the other in which you open`gnuplot
. then gnuplot
will genrate a figure in a fourth windon (the X11
window).
You never close the editor, you never close the two terminals. You play and circle between them.
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). 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 en première colonne x
, en seconde h
, en troisième Q
et en quatrième t
est créé. Pour le tracer, on lance gnuplot:
p 'solxhQt.OUT' w l
this gives a view i with h function of x at the different t superposed.
set ylabel "h"; set xlabel "x";
p 'solxhQt.OUT' w l
par défaut, 1a première colonne est l’abcisse, et la deuxième l’ordonnée, on rajoute Q
troisième variable en fonction de la première
reset; set xlabel "x";
p 'solxhQt.OUT' w l,'' u 1:3 w l
reset; set xlabel "x";
p 'solxhQt.OUT' w l,'' u 1:3 w l
Si on veut une représntation 3D (h,x) avec le temps t vers l’arrière
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
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)
for use of gnuplot
see http://basilisk.fr/sandbox/M1EMN/BASIC/gnuplot_examples.c
Go further
- 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
- test the simple tsunami model
h[i]=1+.2*exp(-x[i]*x[i]);
and see the right/left waves…
- You can add some viscosity in the code to test the example of viscous collapse
first change for a heap (two dams breaking!):
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+1;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
- You can put turbulent friction (this is the Dressler Problem See Chanson page 357 and see dam break Dressler )
turbulent friction
\displaystyle
\frac{d\mathbf{u}}{dt} = - C_f|\mathbf{u}|\frac{\mathbf{u}}{h}
with C_f=.05.
for(i=1;i<nx+1;i++) { if(h[i]>0.) un[i]= un[i]/(1+0.5*dt*fabs(u[i])/h[i]); }
- You can put granular friction…
ex http://basilisk.fr/sandbox/M1EMN/Exemples/savagestaron.c
etc.
- 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