sandbox/geoffroy/testcase/morerivers.c

    A more real case for imposing flow rates

    In this example, we use the discharge() function in order to impose different flow rates on different rivers situated on the same boundary of a Saint-Venant simulation. Wa also use adaptative quadtrees.

    #include "saint-venant.h"
#include "discharge.h"

    The domain is 400 metres squared, centered on the origin. Time is in seconds. The minimum level of refinement is 2^4 and the maximum one is 2^8. We fix the maximum error on the water elevation to a low value (1e^{-3}). Therefore, the grid is refined to its maximum value where it has presence of fluid.

    #define MINLEVEL 4
#define MAXLEVEL 8
#define HMAXE 1e-3

int main(){
  L0 = 20.;
  X0 = - L0/2.;
  Y0 = - L0/2.;
  G = 9.81;
  N = 1 << MAXLEVEL; 
  run();
}

    Initial conditions

    We chose a complicated river bed with height “valleys” flowing in four lakes (see figure below). We set the field river[] with the values 1 and 2 in order to distinguish the rivers on the same boundary.

    event init (i = 0){

    We start with a dry riverbed, so that the problem does not have a natural timescale the Saint-Venant solver can use. We set a maximum timestep to set this timescale.

      DT = 1e-2;
  foreach() {

    We define a topography with two valleys at each side of the domain and one lake at each quarter of the map. See the map below.

        zb[] = 5e-4*(x*y*x*y - x*x - y*y + 2000*(cos(x/10.*2*3.14)+cos(y/10.*2*3.14)));
    reset
set isosample 20,20
set pm3d depthorder
f(x,y) = 5e-4*(x*y*x*y - x*x - y*y + 2000*(cos(x/10.*2*3.14)+cos(y/10.*2*3.14))) 
splot [-10:10][-10:10] f(x,y) with pm3d t'topography'
    Topography (script)

    Topography (script)

    The scalar river[] is defined by default in the discharge.h library. We use it to store the location of each rivers. Its value is 1 for the river at north-east and south-west of the map and 2 for the rivers at north-west and south-east.

    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
        river[] = x*y > 0 ? 1 : 2;
    reset
set isosample 20,20
set pm3d map
f(x,y) = (x*y < 0)+1
splot [-10:10][-10:10] f(x,y) t'river[]'
    river[] (script)

    river[] (script)

      }
    expected expression before ‘{’ token 
      boundary ({zb, river});
}

    Adaptation

    We use the adap_wavelet function to adapt the grid only on the wet cells. You can find a full explaination of this process in the tsunami example

    void adapt_H() {
#if QUADTREE
  scalar h1[];
  foreach(){
    h1[] = h[];
  }
  boundary ({h1});
  adapt_wavelet ({h1}, (double[]){HMAXE} ,MAXLEVEL,MINLEVEL);
#endif
}

    We refine at each time step.

    event do_adapt_H (i++) adapt_H();

    Boundary conditions

    We impose inflow on all boundaries. In addition, on each boundary, the tangential velocity is set to zero and a neumann condition is fixed on the normal velocity.

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

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

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

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

    We impose different flow rates : 1 m^3.s^{-1} for the rivers in the north-east quarter, 1.5 m^3.s^{-1} for the ones at north-west, 2 m^3.s^{-1} for the ones at south-west, and 0 at south-east. The process is more detailed in the multi-river test case.

    double eta1, eta2, eta3, eta4, eta5, eta6;

event inflow (i++) {
  // Top boundary
    implicit declaration of function ‘discharge’ [-Werror=implicit-function-declaration] 
      eta1 = discharge (0.5, top, 1); // N-E
    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
      eta2 = discharge (1, top, 2); // N-W
    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
      h[top] = max ( (river[] == 1 ? eta1 : eta2) - zb[] , 0) ;
  eta[top] = max ((river[] == 1 ? eta1 : eta2), zb[]);

  // Left boundary
  eta4 = discharge (1.5, left, 1); // S-W
    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
      eta3 = discharge (1, left, 2); // N-W
    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
      h[left] = max ((river[] == 1 ? eta4 : eta3) - zb[] , 0);
  eta[left] = max ((river[] == 1 ? eta4 : eta3), zb[]);

  // Bottom boundary
    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
      eta5 = discharge (1.5, bottom, 1); //S-W
    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
      h[bottom] = max ((river[] == 1 ? eta5 : zb[]) - zb[] , 0);
  eta[bottom] = max ((river[] == 1 ? eta5 : zb[]), zb[]);

  // Right boundary
    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
      eta6 = discharge (0.5, right, 1); //N-E
    ‘river’ undeclared (first use in this function)
    expected expression before ‘]’ token 
      h[right] = max ((river[] == 1 ? eta6 : zb[]) - zb[] , 0);
  eta[right] = max ((river[] == 1 ? eta6 : zb[]), zb[]);
}

    Outputs

    We compute the evolution of the water volumes in each quarter of the map.

    event volume (i += 10) {
  double volume1 = 0, volume2 = 0, volume3 = 0, volume4 = 0;
  foreach() {
    double dv = h[]*sq(Delta); 
    if (y > 0 && x > 0)
	volume1 += dv; // N-E
    else if(y > 0 && x < 0 )
      volume2 += dv; // N-W
    else if( y < 0 && x < 0 )
      volume3 += dv; // S-W
    else  volume4 += dv; // S-E-
  }
  fprintf (stderr, "%g %g %g %g %g\n",
	   t, volume1, volume2, volume3, volume4);
}

    We use gnuplot to produce an animation of the water surface.

    event animation (t <= 1; t += 0.01) {
  double dx = 2.*L0/(1 << MAXLEVEL), dy = dx;
  printf ("set title 't = %.3f'\n"
	  "sp [%g:%g][%g:%g][-5:5] '-'"
	  " u 1:2:($3+$4-.05) t 'free surface' w l lt 3,"
	  " '' u 1:2:4 t 'topography' w l lt 2\n",
	  t, X0, -X0, Y0, -Y0);
  for (double x = X0;  x <= X0 + L0; x += dx) {
    for (double y = Y0; y <= Y0 + L0; y += dy)
      printf ("%g %g %g %g\n",
	      x, y, interpolate (h, x, y),  interpolate (zb, x, y));
    putchar ('\n');
  }
  printf ("e\n"
	  "pause %.5lf \n\n", 0.);
}

    We record a movie for the field level which is storing the level of refinement.

    event movies (t+=0.01) {
  static FILE * fp1 = popen ("ppm2mpeg > level.mpg", "w");
  scalar l[];
  foreach()
    l[] = level;
  output_ppm (l, fp1, min = MINLEVEL, max = MAXLEVEL, n = 512);
}

    Results

    reset
set key top left
set xlabel 'Time'
set ylabel 'Volume'
f(x)=a*x + b
fit [0.1:0.8] f(x) './log' u 1:2 via a,b
g(x)=c*x + d
fit [0.1:0.8] g(x) './log' u 1:3 via c,d
h(x)=k*x + l
fit [0.1:0.8] h(x) './log' u 1:4 via k,l
i(x)=m*x
fit [0.1:0.8] i(x) './log' u 1:5 via m
title1=sprintf("f(x) = %1.3f*t + %1.3f", a,b)
title2=sprintf("g(x) = %1.3f*t + %1.3f", c,d)
title3=sprintf("h(x) = %1.3f*t + %1.3f", k,l)
title4=sprintf("i(x) = %1.3f*t", m)
plot [0:0.8] './log' u 1:2 t 'North-east (q=1)', \
     f(x) t title1, \
     './log' u 1:3 t 'North-west (q=2)', \
     g(x) t title2, \
      './log' u 1:4 t 'South-west (q=3)', \
     h(x) t title3, \
     './log' u 1:5 t 'South-east (q=0)', \
     i(x) t title4
    Evolution of water volumes in both rivers (script)

    Evolution of water volumes in both rivers (script)

    reset
set view 40,15
set xlabel 'x'
set ylabel 'y'
set hidden3d
set term gif animate
set output 'movie.gif'
load './out'
    Animation of the free surface. (script)

    Animation of the free surface. (script)

    Animation of the level of refinement. Dark blue is 4 and dark red is 8.

    Animation of the level of refinement. Dark blue is 4 and dark red is 8.