/** # Resolution of flood wave in 1D This is a simple `C` code, not a `basilisk` one. # Model equations The flood wave or kinetic wave, see Whitham p 82, Fowler p 76, or Chanson, is a long wave approximation of shallow water equations were inertia and pressure gradient are neglected. There is only equilibrium between slope and friction. The waves go only downstream. We have to solve, with flux $\bar F=\bar h^{3/2}$: $$ \frac{\partial}{\partial \bar t} \bar h + \frac{\partial}{\partial \bar x} \bar F = 0 \text{ which is as well } \frac{\partial}{\partial \bar t} \bar h + \bar c \frac{\partial}{\partial \bar x} \bar h = 0 \text{ with } \bar c = \partial \bar F/\partial \bar h=\frac{3}{2}\sqrt{h}$$ ~~~gnuplot set samples 9 set label "U i-1" at 1.5,3.1 set label "U i" at 2.5,3.15 set label "U i+1" at 3.5,2.5 set xtics ("i-2" 0.5, "i-1" 1.5, "i" 2.5,"i+1" 3.5,"i+2" 4.5,"i+3" 5.5) set arrow from 2,1 to 2.5,1 set arrow from 3,1 to 3.5,1 set label "F i-1/2" at 2.1,1.25 set label "F i+1/2" at 3.1,1.25 set label "x i-1/2" at 1.5,0.25 set label "x i" at 2.4,0.25 set label "x i+1/2" at 3.,0.25 set label "x" at 0.5,2+sin(0) set label "x" at 1.5,2+sin(1) set label "x" at 2.5,2+sin(2) set label "x" at 3.5,2+sin(3) set label "x" at 4.5,2+sin(4) set label "x" at 5.5,2+sin(5) p[-1:7][0:4] 2+sin(x) w steps not,2+sin(x) w impulse not linec 1 ~~~ The numerical flux across face i+1/2 is denoted $F_{i+1/2}$ (or $F^n_{i+1/2}$ at time $n$), it is function (say $f$) of values before and after the face (i+1) which are $U_i$ and $U_{i+1}$ $$ F_{i+1/2}=f(U_i,U_{i+1}). $$ The position of the center of the cell is $x_{i}$. # Code mandatory declarations: */ #include #include #include #include /** definition of the field U, the flux F, time step */ double*x=NULL,*h=NULL,*F=NULL; double dt; double cDelta; double L0,X0,Delta; double t; int i,it=0; int N; /** Here we code the flux $h^{3/3}$ as a function of the values of the left ($h_g$) and right ($h_d$) values of the face, with an approximation for the velocity $\frac{3}{2}h^{1/2}$, $c= \frac{3}{2} \frac{h_g^{1/2} + h_d^{1/2}}2$ so that $$ F=\frac{h_g^{3/2} + h_d^{3/2}}2 - c \frac{h_d -h_g}2 $$ */ double FR1(double hg, double hd){ double c=1.5*(sqrt(hg)+sqrt(hd))/2; return (hg*sqrt(hg)+hd*sqrt(hd))*0.5-c*(hd-hg)*0.5;} /** Main with definition of parameters */ int main() { L0 = 12.; X0 = -L0/4; N = 256; t=0; dt = (L0/N)/2; Delta = L0/N; /** dynamic allocation */ x= (double*)calloc(N+2,sizeof(double)); h= (double*)calloc(N+2,sizeof(double)); F= (double*)calloc(N+2,sizeof(double)); /** loop for initial `h(x,0)`: initial elevation: a "bump" The celle `i=0` is a ghost cell/ The "first" (`i=1`) cell is 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` */ for(i=0;i out ~~~ ## Results in gnuplot type ~~~bash reset set xlabel "x" set ylabel "h" h(x)=.5+exp(-x*x) c(x)=3./2*sqrt(h(x)) set parametric set dummy x p[:][:]'out' w l , x,h(x) not w l lc -1,\ x+1*c(x),h(x) not w l lc -1,\ x+2*c(x),h(x) not w l lc -1,\ x+3*c(x),h(x) not w l lc -1,\ x+4*c(x),h(x) not w l lc -1,\ x ,0 not w l lc -1 ~~~ ~~~gnuplot reset set xlabel "x" set ylabel "h" h(x)=.5+exp(-x*x) c(x)=3./2*sqrt(h(x)) set parametric set dummy x p[:][:]'out' w l , x,h(x) not w l lc -1,\ x+1*c(x),h(x) not w l lc -1,\ x+2*c(x),h(x) not w l lc -1,\ x+3*c(x),h(x) not w l lc -1,\ x+4*c(x),h(x) not w l lc -1,\ x ,0 not w l lc -1 ~~~ which gives $h(x,t)$ numerical and exact plotted here for t=0 1 2 3 4 and $-3