src/test/layered.c

    Transcritical flow over a bump with multiple layers

    We want to reproduce the transcritical test case of Audusse et al, 2011, section 5.6.2.

    #include "grid/cartesian1D.h"
    #include "saint-venant.h"

    We need a field to store the variable bottom friction.

    scalar lambda[];
    
    int main() {
      X0 = 0.;
      L0 = 21.;
      G = 9.81;
      N = 256;

    The viscosity is set to \nu = 0.01 m^2/s and the bottom friction is variable.

      nu = 0.01;
      lambda_b = lambda;

    We vary the number of layers.

      nl = 2;  run();
      nl = 5;  run();
      nl = 15; run();
    }

    We impose the outlet water level.

    h[right]   = dirichlet(0.6);
    eta[right] = dirichlet(0.6);

    Initialisation

    We initialise the topography, the initial water depth h and we create a field hc to check convergence on h.

    scalar hc[];
    
    event init (i = 0) {
      foreach() {
        zb[] = max(0., 0.2*(1. - 1./sq(5.75/2.)*sq(x - 10.)));
        hc[] = h[]  = 0.6 - zb[];
      }

    We call the friction event (below) to initialize the bottom friction.

      event ("friction");

    Boundary conditions on velocity

    We impose a constant inflow of 1 m^2/s at the inlet and a Neumann condition at the outlet.

      for (vector u in ul) {
        u.n[left] = dirichlet(h[left] ? 1./h[left] : 0.);
        u.n[right] = neumann(0.);
      }
    }

    Bottom friction

    We use the Strickler relation: \displaystyle k(h,\mathbf{U}) = \frac{g}{S^2h^{1/3}}|\mathbf{U}| with S = 25 m^{1/3}/s the Strickler coefficient, h the water depth and \mathbf{U} the depth-averaged velocity. Note that we have to use a lower Strickler coefficient (i.e. larger friction) to get results comparable to those of Audusse et al, 2011.

    event friction (i++) {
      foreach() {
        double U = 0.;
        int l = 0;
        for (vector u in ul)
          U += u.x[]*layer[l++];
        double S = 25., k = G/(sq(S)*pow(h[],1./3.))*fabs(U);
        lambda[] = k > 0. ? nu/k : 0.;
      }
      boundary ({lambda});
    }

    We check for convergence.

    event logfile (t += 0.1; i <= 100000) {
      double dh = change (h, hc);
      if (i > 0 && dh < 1e-4)
        return 1;
    }

    Uncomment this part if you want on-the-fly animation.

    #if 0
    event output (i++) {
      static FILE * fp = popen ("gnuplot", "w");
      fprintf (fp,
               "set title 'nl=%d, t=%f'\n"
               "set xl 'x'\nset yl 'h'\n"
               "plot [0:21][] '-' u 1:2 w l t 'eta', '-' u 1:3 w l t 'zb'\n",
    	   nl, t); 
      foreach()
        fprintf (fp, "%g %g %g\n", x, eta[], zb[]);
      fprintf (fp, "e\n");
      fflush (fp);
    }
    #endif

    Outputs

    At the end of the simulation we save the profiles.

    event output (t = end) {
      char name[80];
      sprintf (name, "end-%d", nl);
      FILE * fp = nl == 15 ? stderr : fopen (name, "w");
      foreach() {
        fprintf (fp, "%g %g %g\n", x, eta[], zb[]);
        if (nl == 15) {
          double z = zb[];
          int l = 0;
          printf ("%g %g %g\n", x, z, u.x[]);
          for (vector u in ul) {
    	z += layer[l++]*h[];
    	printf ("%g %g %g\n", x, z, u.x[]);
          }
          printf ("\n");
        }
      }
    }

    Results

    set xr [0:21]
    set yr [0:1]
    set xlabel 'x'
    set ylabel 'z'
    plot 'end-2' u 1:3 w l t 'topography', \
         'end-2' u 1:2 w l t '2 layers', \
         'end-5' u 1:2 w l t '5 layers', \
         'log'   u 1:2 w l t '15 layers'
    Free surface and topography. This can be compared to figure 9 of Audusse et al, 2011. (script)

    Free surface and topography. This can be compared to figure 9 of Audusse et al, 2011. (script)

    set term PNG enhanced font ",10"
    set output 'vel.png'
    set pm3d
    set pm3d map interpolate 10,1
    unset key
    # jet colormap
    set palette defined ( 0 0 0 0.5647, 0.125 0 0.05882 1, 0.25 0 0.5647 1,	\
    0.375 0.05882 1 0.9333, 0.5 0.5647 1 0.4392, 0.625 1 0.9333 0, 0.75 1 0.4392 0,	\
    0.875 0.9333 0 0, 1 0.498 0 0 )
    
    splot 'out' u 1:2:3
    Horizontal velocity field (15 layers). (script)

    Horizontal velocity field (15 layers). (script)

    See also