/**
# 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, slope and pressure gradient are neglected.
 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 <stdio.h>
#include <stdlib.h> 
#include <math.h>
#include <string.h>
/** 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 = 128*2*2;
  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<N+2;i++)
    {  x[i]=X0+(i-1./2)*Delta;  
       h[i] = 0.5+exp(-x[i]*x[i]);
  }
/**
 begin the time loop
*/  
   while(t<=4){
     t = t + dt;
     it++;	
/** 
print data at t=1 t=2 t=3 
*/
       
if((it == 1)||(it == (int)(1/dt))
   || (it == (int)(2/dt))
   || (it == (int)(3/dt))
   || (it == (int)(4/dt)) ){
 for(i=0;i<=N;i++)
    fprintf (stdout, "%g %g %g  \n", x[i], h[i], t);
  fprintf (stdout, "\n\n");
}
/** 
flux 
*/    
 for(i=1;i<=N+1;i++)
   {
    F[i] = FR1(h[i-1],h[i]);
    }
/** 
explicit step
update 
$$
 h_i^{n+1}=h_i^{n} -{\Delta t} \dfrac{F_{i+1/2}-F_{i-1/2}}{\Delta x}
$$*/ 
 for(i=1;i<=N;i++)
    h[i] +=  - dt* ( F[i+1] - F[i] )/Delta;
       
/**
     Boundary condition
*/
   h[0] = 0.5;
   h[N+1]=h[N];
 }
    
  
  free(h);
  free(F);
  free(x);
}
/**
## Run
Then compile and run:

~~~bash
 cc  floodwaveC.c -lm ; ./a.out> 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<x<9$
 
 
## Note

the output is the standard one (terminal), we can save in a file and plot the generated file with `gnuplot`, see 
[http://basilisk.fr/sandbox/M1EMN/BASIC/gnuplot_examples.c]()


    FILE * fp = fopen ("out.txt", "w")
    for(i=0;i<=N;i++)
    fprintf (fp, "%g %g %g  \n", x[i], h[i], t);
    fprintf (fp, "\n\n");
    fclose(fp);
   

 
 
 
## Links
 
 * [http://basilisk.fr/sandbox/M1EMN/BASIC/advecte1C.c]()
 
 * [http://basilisk.fr/sandbox/M1EMN/BASIC/advecte1.c]()
 
 * [http://basilisk.fr/sandbox/M1EMN/Exemples/floodwave.c]()
 
 * [http://basilisk.fr/sandbox/M1EMN/BASIC/gnuplot_examples.c]()
 
## Bibliography
 
 * [Lagrée P-Y](http://www.lmm.jussieu.fr/~lagree/COURS/MFEnv/MFEnv.pdf)
 "Equations de Saint Venant et application, Ecoulements en milieux naturels" Cours MSF12, M1 UPMC
 
  * [Lagrée P-Y](http://www.lmm.jussieu.fr/~lagree/COURS/MFEnv/code_C_saintvenant.pdf)
  "Résolution numérique des équations de Saint-Venant,
mise en oeuvre en volumes finis par un solveur de Riemann bien balancé" Cours MSF12, M1 UPMC
 
 * G. B. Whitham "Linear and Nonlinear Waves" Wiley-Interscience, p 82
 
 * Fowler ["Mathematics of the environment"](https://people.maths.ox.ac.uk/fowler/courses/mathenvo/envonotes.pdf)
 
 
ready for new site  01/22
*/