#include "b-flood/saint-venant-topo.h"
#include "b-flood/manning.h"
#include "b-flood/inflow.h"
#include "b-flood/hydro.h"
#include "terrain.h"
#include "input.h"

#define MAXLEVEL 10 // 1024 Volume
#define MINLEVEL 6
#define EMAXE 5e-4

int nf, nc, ncell;
timer t0;

int main()
{
  n = 0.0162;
  L0 = 43.4;
  X0 = 0.1;
  Y0 = 0;
  G = 9.81;
  dry = 1e-6;
  N = 1 << MAXLEVEL; // Which is equivalent to N = 2^(maxlevel)
  run();
}

// The case lasts 180 seconds
event stop (t = 180)
  fprintf (stderr, "# The end\n");

event init (i = 0)
{
  // We record the topography in a special way (k-tree) so that we can
  // use adaptive refinement on it.  
  terrain (zb, "./Topo/topo-house", NULL);
  
  // We have to set a time step because there is no water at the beginning.
  DT = 1;
  
  t0 = timer_start();
  double dx = L0/N;

  // We set the interval through which the water flow will come in.
  foreach()
    if (y <= 8.4 && y >= 5 && x <= 8*dx)
      river[] = 1;
  boundary ({river});
}

u.t[left] = dirichlet(0);
u.n[left] = neumann(0);

double l1, eta1, vola = 0;
event bound_left( i++ )
{
  l1 =  hydrograph ("hydro.asc", 1);
  if (l1 < 0)
    l1 = 0;
  eta1 = eta_b(l1,left,1);
  
  h[left] = dirichlet( max( river[] == 1 ? eta1 - zb[] : 0 , 0) );
  eta[left] = dirichlet( max( river[] == 1 ? eta1 : zb[]  , zb[]) );

  double vol = 0;
  foreach(reduction(+:vol)){
    if( h[] > dry )
      vol += h[] * sq(Delta);
  }
  double Q = 0;
  if( dt > 0 )
    Q =  (vol- vola)/dt;
  vola = vol;
 
  static FILE * fpp = fopen ("eta.dat", "w");
   if (i == 0)
     fprintf (fpp, "#t Q Qread EtaImp\n");
  fprintf (fpp, "%lf %lf %lf %lf \n", t, Q, l1, eta1);
  fflush (fpp);  
}

u.t[right] = neumann(0);
u.n[right] = neumann(0);
h[right] = neumann(0);
eta[right] = neumann(0);

event timestepreduce (t = 17) {
  DT = 0.01;
}

event adapt (t = 17; i++)
{
  scalar nh[];
  foreach()
    nh[] = h[] > dry ? h[] : 0;
  boundary ({nh});
  astats s = adapt_wavelet ({nh}, (double[]){EMAXE}, MAXLEVEL, MINLEVEL);
  nc += s.nc;
  nf += s.nf;
}

event logfile (i += 10)
{
  timing tr = timer_timing (t0, i, 0, NULL);
  stats s = statsf (h);
  scalar v[];
  foreach()
    v[] = norm(u);
  
  stats no = statsf (v);
  if( i == 0 )
    fprintf (stderr, "#1t 2treal 3i 4h.min 5h.max 6h.sum 7u.min 8u.max 9dt\n");  
  fprintf (stderr, "%g %g %d %g %g %g %g %g %g\n",
	   t,tr.real, i, s.min, s.max, s.sum,  no.min, no.max, dt);
}

Gauge gauges[] = {  
  // gauge   x         y       zb
  {"S2",0.280560,7.585584, "#S2"},
  {"S4",4.255610,7.745400, "#S4"},
  {"S6S",6.753537,9.831888, "#S6S"},
  {"S6D",7.552486,7.696863, "#S6D"},
  {"S8D",10.743772,7.788068, "#S8D"},
  {"P1",2.113886,7.866009, "#P1"},
  {"P2",3.931183,8.893101, "#P2"},
  {"P3",4.527122,6.587255, "#P3"},
  {"P4",7.086251,8.768313, "#P4"},
  {"P5",10.163515,10.218714, "#P5"},
  {"P8",14.950365,11.804039, "#P8"},
  {"P9",15.684258,13.802726, "#P9"},
  {"P10",18.190768,10.429219, "#P10"},
  {"P13",19.814693,11.984270, "#P13"},
  {"P18",23.848138,14.744612, "#P18"},
  {"P19",25.469928,16.186758, "#P19"},
  {"P21",30.442363,18.357176, "#P21"},
  {"P23",36.466741,18.750427, "#P23"},
  {"P24",39.024308,22.827062, "#P24"},
  {"P25",40.443587,21.238121, "#P25"},
  {"P26",40.965033,26.182320, "#P26"},
  {NULL}
};
event gauges1 (t += 1) output_gauges (gauges, {h});

event movies ( t= 17; t += 0.1 )
{
  scalar m[], v[], fr[];
  
  foreach(){
    v[] = norm(u);
    m[] = (h[] > dry ) - 0.5;
    fr[] = h[] > dry ? v[]/sqrt(G*h[]) : 0;
  }    
  boundary ({m,v,fr});
  
  output_ppm (h, mask = m, min = 0, max = 0.4, n = N, linear = true, file = "height.mp4");
  output_ppm (h, mask = m, min = 0, max = 0.4, n = N, linear = true,
	      box = {{0,5},{12.5,12.5}}, file = "height-town1.mp4"); 
  output_ppm (h, mask = m, min = 0, max = 0.4, n = N, linear = true,
	      box = {{15,9},{25,19}}, file = "height-town2.mp4");
  
  

  output_ppm (v,  mask = m, min = 0, max = 3, n = N, linear = true, file = "vel.mp4");

  output_ppm (fr, mask = m, min = 0, max = 2, n = N, linear = true, file = "froude.mp4");
  
  scalar l = m;
  foreach()
    l[] = level;
  
  output_ppm (l, min = MINLEVEL, max = MAXLEVEL, n = N, file = "level.mp4");
  output_ppm (l, min = MINLEVEL, max = MAXLEVEL, n = N,
	      box = {{0,5},{12.5,12.5}}, file = "level-town1.mp4");
  output_ppm (l, min = MINLEVEL, max = MAXLEVEL, n = N,
	      box = {{15,9},{25,19}}, file = "level-town2.mp4");
}


event snap (t = 65)
{
  static FILE * rph = fopen ("raster_h-t65.asc", "w");
  output_grd (h, rph);

  scalar m[], v[];
    
  foreach(){
    v[] = norm(u);
    m[] = (h[] > dry ) - 0.5;
  }    
  boundary ({m,v});
  
  static FILE * rpu = fopen ("raster_v-t65.asc", "w");
  output_grd (v, rpu);
}

event ending (t = end)
{
  static FILE * rph = fopen ("raster_h.asc", "w");
  output_grd (h, rph);

  scalar m[], v[];
    
  foreach(){
    v[] = norm(u);
    m[] = (h[] > dry ) - 0.5;
  }    
  boundary ({m,v});
  
  static FILE * rpu = fopen ("raster_v.asc", "w");
  output_grd (v, rpu);
}

event figures (t = 17; t <= 67; t += 2)
{
  scalar m[], etam[];
  foreach() {
    etam[] = eta[]*(h[] > dry);
    m[] = etam[] - zb[];
  }
  boundary ({m, etam});
  char name[80];
  sprintf (name, "h-%g.png", t);
  output_ppm (h, min = 0, max = 0.4, mask = m, file = name, n = 1024,
	      linear = true,
	      opt = "-fill white -opaque black");
    
  sprintf (name, "level-%g.png", t);
  scalar l[];
  foreach()
    l[] = level;
  output_ppm (l, min = MINLEVEL, max = MAXLEVEL, file = name, n = 1024,
	      linear = false,);
}