// SCRIPT WHICH TAKES AN INITIAL SURFACE ELEVATION GENERATED BY THE DUFFY (1992) MODEL AT A CERTAIN TIME (CONST) AND SOLVES THE SAINT VENANT EQUATIONS. 
// OUTPUTS THE WAVE HEIGHT AS A FUNCTION OF TIME FOR R = 6. 
// The non dimensionalised duration of the eruption is T = 1, the water depth is 1, the source height is 0.4 and the source width is 0.5. 

#include "saint-venant.h"
double default_sea_level=0.;
#include "cgd_read2D.h"
#include "input.h"

/**


We then define a few useful macros and constants. */

#define MAXLEVEL 9
#define MINLEVEL 1
#define ETAE     1e-3 // error on free surface elevation (1 cm)
#define myconst CONST

int main()
{
  #if QUADTREE
  // 32^2 grid points to start with
  N = 1 << MINLEVEL;
#else // Cartesian
  // 1024^2 grid points
  N = 1 << MAXLEVEL;
#endif
  foreach_dimension() { // Free outflow boundary conditions
    u.n[right]= neumann(0);
    u.n[left]= neumann(0);
  }
   size (16.);
  origin (-8.,-8.);
  N = 1 << MAXLEVEL;
  run();
}

int adapt() {
#if QUADTREE
  scalar eta[];
  foreach()
    eta[] = h[] > dry ? h[] + zb[] : 0;
  boundary ({eta});

  astats s = adapt_wavelet ({eta}, (double[]){ETAE},
			    MAXLEVEL, MINLEVEL);
  fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
  return s.nf;
#else // Cartesian
  return 0;
#endif
}

event init (i = 0)
{

  foreach(){
    zb[]=-1.;
    h[] = max(0., - zb[]);}
  boundary ({h,zb});

/* Open the file. */
  char *fname = "/home/battershilll/basilisk/src/tests/output_tophat_t_CONST.cgd"; // This cgd file is generated from the Duffy_t_2D_spatial.m code
  FILE *fidin=fopen(fname,"r");
  deformation_cgd_read (x=0., y=0., fid = fidin, iterate = adapt );
  fprintf(stderr,"deformation\n");
  fclose ( fidin );

}


Gauge gauges[] = {
  {"WG_minlevel1_maxlevel9_error03_6_0", 6, 0}, //outputs at 6,0
  {NULL}
};

event output (t += 0.1; t <= 21) // NOTE: t = 0 is actually t = 3.1 (in terms of overall simulation). This needs to be changed in the output file. 
  output_gauges (gauges, {eta});

event movies (t += 0.1; t <= 21) {
  static FILE * fp = NULL;
  if (!fp) fp = popen ("ppm2mpeg > eta.mpg", "w");
  scalar m[], etam[];
  foreach() {
    etam[] = eta[]*(h[] > dry);
    m[] = etam[] - zb[];
  }
  boundary ({m, etam});
  output_ppm (etam, fp, n = 512, linear = true);

}

/**
## Adaptivity

And finally we apply our *adapt()* function at every timestep. */

event do_adapt (i++) adapt();