sandbox/pressure_pool.c

Spatially Varying Pressure

Example of a constant spatially varying pressure over a square constant depth water body

A pressure field of Pa(x,y)-P_{ref}=Amp*sin(2*\pi*x/L)*cos(2*\pi*y/L) is applied to a square pool

If grid/cartesian1D.h is specified we are working only in 1D #include “grid/cartesian1D.h”

#include "saint-venant_ss.h"

#define MAXLEVEL 8
#define MINLEVEL 4
#define ETAE 1e-8

Inverse barometer conversion hPa to metres \displaystyle \eta = - \frac{\Delta p}{\rho_w g} Thus for pressure in hectoPascals the coefficient is \frac{100}{\rho_w g}.

double rho_inv = 1./1025.;
double amp = 200.;

double compare (scalar v, scalar vn)
{
  double maxC = 0.;
  foreach(reduction(maxC:max)) {
    double dv = fabs (v[] - vn[]);
    if (dv > maxC)
      maxC = dv;
  }
  return maxC;
}

Here we can set standard parameters in Basilisk

int main()
{
#if QUADTREE
  // 32^2 grid points to start with
  init_grid( 1 << MINLEVEL );
#else // Cartesian
  // 1024^2 grid points
  init_grid( 1 << MAXLEVEL );
#endif

Here we setup the domain geometry. For the moment Basilisk only supports square domains. This case uses metres east and north. We set the size of the box L0 and the coordinates of the lower-left corner (X0,Y0). In this case we are assuming a square ‘pool’ of length 1000 m.

  // the domain is
  size (1000.);
  origin(-L0/2.,-L0/2.);

G is the acceleration of gravity required by the Saint-Venant solver. This is the only dimensional parameter..

  G = 9.81;
  run();
}

Boundary conditions

We set the normal velocity component on the left, right bottom and top boundaries to a “radiation condition” with a reference sealevel of inverse barometer. Note that the sign is important and needs to reflect the orientation of the boundary.

u.x[left]   = - radiation(-Pa2m*(pa[]));
u.x[right]  = + radiation(-Pa2m*(pa[]));
u.y[bottom] = - radiation(-Pa2m*(pa[]));
u.y[top]    = + radiation(-Pa2m*(pa[]));

Adaptation

Here we define an auxilliary function which we will use several times in what follows. Again we have two #if...#else branches selecting whether the simulation is being run on an (adaptive) quadtree or a (static) Cartesian grid.

We want to adapt according to two criteria: an estimate of the error on the free surface position – to track the wave in time – and an estimate of the error on the maximum wave height hmax – to make sure that the final maximum wave height field is properly resolved.

We first define a temporary field (in the automatic variable ?) which we set to h+z_b but only for “wet” cells. If we used h+z_b everywhere (i.e. the default \eta provided by the Saint-Venant solver) we would also refine the dry topography, which is not useful.

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

h+z_b everywhere (i.e. the default \eta provided by the Saint-Venant solver) we would also refine the dry topography, which is not useful.

We can now use wavelet adaptation on the list of scalars {?,hmax} with thresholds {ETAE,HMAXE}. The compiler is not clever enough yet and needs to be told explicitly that this is a list of doubles, hence the (double[]) type casting .

The function then returns the number of cells refined.

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

event initiate(i=0)
{

The initial still water surface is at z=0 so that the water depth h is…

  double ramp = t < 120. ? t/120. : 1.;
  foreach() {
    zb[] = -2;
    h[] = max(0., - zb[]);
    pa[] = ramp*amp*sin(2.*pi*x/L0)*cos(2.*pi*y/L0);
  }
 boundary ({h,zb}); 
}

We want the simulation to stop when we are close to steady state. To do this we store the h field of the previous timestep in an auxilliary variable hn.

scalar hn[],etaref[];

event init_hn (i = 0) {
  foreach() {
    hn[] = h[];
    etaref[] = -Pa2m*pa[];
  }
}

Every 10 timesteps we check whether u has changed by more than 10^-8. If it has not, the event returns 1 which stops the simulation. We also output running statistics to the standard error.

event logfile (i+= 10; i <= 10000) {
  double dh = change (h, hn);
  if ( (i > 100 && dh < 1e-8 ) || i==10000) {
    double dhref = compare (eta,etaref);
    fprintf (stderr, "# Max difference is %.9g\n# Time to convergence %g",
	     dhref,t); 
    FILE * fp = fopen("diff", "w");
    fprintf(fp,"# 1:x 2:y 3:h 4:zb 5:eta 6:etaref 7:diff\n");
    for (double x = -500; x <= 500; x += 10) {
      for (double y = -500; y <= 500; y += 10)
	fprintf (fp, "%g %g %.9g %.9g %.9g %.9g %.9g\n", x, y, 
		 interpolate (h, x, y), 
		 interpolate (zb, x, y), 
		 interpolate (eta, x, y), 
		 interpolate (etaref, x, y), 
		 interpolate (etaref, x, y) - interpolate (eta, x, y));
      fprintf( fp, "\n");
    }
    fclose (fp);
    return 1; /* stop */
  }

  stats s = statsf (h);
  norm n = normf (u.x);
  double detaref = compare (eta,etaref);
  if (i == 0)
    fprintf (stderr, "t i h.min h.max h.sum u.x.rms u.x.max dh detaref dt\n");
  fprintf (stderr, "%g %d %g %g %g %g %g %g %g %g\n", t, i, s.min, s.max, 
	   s.sum, n.rms, n.max, dh, detaref, dt);
}

We also use a simple implicit scheme to implement quadratic bottom friction i.e. \displaystyle \frac{d\mathbf{u}}{dt} = - C_f|\mathbf{u}|\frac{\mathbf{u}}{h} with C_f=10^{-4}. We ramp up the pressure from zero

event source (i++) {
  double ramp = t < 120. ? t/120. : 1.;
  foreach() {
    pa[] = ramp*amp*sin(2.*pi*x/L0)*cos(2.*pi*y/L0);
    etaref[] = -Pa2m*pa[]; 
    eta[] = h[] > dry ? h[] + zb[] : 0.;
  }
  foreach() {
    double a_inv = h[] < dry ? 0. : h[]/(h[] + 1e-4*dt*norm(u));
    //double a = h[] < dry ? HUGE : 1. + 1e-2*dt*norm(u)/h[];
    foreach_dimension()
      u.x[] = u.x[]*a_inv;
  }
  boundary ({h,u,zb,pa,etaref});
}

Every 5 minutes, the h, z_b and hmax fields are interpolated bilinearly onto a n x n regular grid and written on standard output.

event snapshots (t += 300) {
  printf ("# file: t-%g\n", t);
  output_field ({eta}, stdout, n = 1 << MAXLEVEL, linear = true);
  scalar l[];
  foreach()
    l[]=level;
  printf ("# file: level-%g\n", t);
  output_field ({l});
}

event movies (t+=60) {
  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, min = -0.002, max = 0.002 , n = 512, linear = true);
  //  output_ppm (etam, fp, n = 512, linear = true);
#if QUADTREE
  static FILE * fp1 = NULL;
  if (!fp1) fp1 = popen ("ppm2mpeg > level.mpg", "w");
  scalar l = etam;
  foreach()
    l[] = level;
  output_ppm (l, fp1, min = MINLEVEL, max = MAXLEVEL, n = 512);
#endif
}

Adaptivity

We apply our adapt() function at every timestep.

event do_adapt (i++) adapt();

Results

After running and processing by gnuplot (using pressure_pool.plot) we get the following pictures and animations.

Evolution of the vorticity field with time.

Difference between calculated and reference surface level