sandbox/WMW/tophat_vary_time_off.c
// 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();