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 λ[];

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

The viscosity is set to ν=0.01m2/s and the bottom friction is variable.

  ν = 0.01;
  lambda_b = λ;

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);
η[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: k(h,U)=gS2h1/3U with S=25m1/3/s the Strickler coefficient, h the water depth and 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);
    λ[] = k > 0. ? ν/k : 0.;
  }
  boundary ({λ});
}

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, η[], 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, η[], 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

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

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

Horizontal velocity field (15 layers).

Horizontal velocity field (15 layers).

See also