src/test/gaussian.c

    Transcritical flow over a bump with multiple layers

    This test case is similar to transcritical.c but using simpler (and more physical) boundary conditions. Both hydrostatic and non-hydrostatic solutions are obtained. The non-hydrostatic solution is compared with a DNS solution using the Navier-Stokes/VOF solver.

    More details are given in Popinet (2019).

    #include "grid/multigrid1D.h"
#if STVT
# include "saint-venant.h"
#else // !STVT
# include "layered/hydro.h"
# if !HYDRO
#   include "layered/nh.h"
# endif
# include "layered/remap.h"
# include "layered/perfs.h"
#endif // !STVT

    The primary parameters are the flow rate, water level, viscosity and bump amplitude.

    #define QL 1.
#define HR 0.6
#define NU 1e-2
#define BA 0.4

int main()
{

    The domain is 30 meters long.

      L0 = 30.;
  G = 9.81;
  N = 512; // less damping with 1024

    The viscosity is set to NU and 20 layers are used to discretise vertically.

      nu = NU;
  nl = 20; // going to 30 changes very little
  run();
}

    Initialisation and boundary conditions

    We create a field hc to check convergence on h.

    scalar hc[];

    The inflow is a parabolic profile with a total flow rate Q. The function below computes the height zc of the middle of the layer and returns the corresponding velocity.

    #if STVT
double uleft (Point point, scalar s, double Q)
{
  double zc = zb[];
  for (int l = 0; l < s.l; l++)
    zc += layer[l]*h[];
  zc += layer[s.l]*h[]/2.;
  return 3./2.*Q/h[]*(1. - sq(zc/h[] - 1.));
}
#else
double uleft (Point point, scalar s, double Q)
{
  double H = 0.;
  int l = 0;
  double zc = zb[];
  for (scalar h in hl) {
    H += h[];
    if (l++ < s.l)
      zc += h[];
  }
  scalar h = hl[s.l];
  zc += h[]/2.;
  return 3./2.*Q/H*(1. - sq(zc/H - 1.));
}
#endif

    We initialise the topography and the initial thickness of each layer h.

    event init (i = 0) {
  foreach() {
    zb[] = BA*exp(- sq(x - 10.)/5.);
    hc[] = HR - zb[];
#if STVT
    h[] = hc[];
#else
    for (scalar h in hl)
      h[] = hc[]/nl;
#endif
  }

    The height of the free-surface \eta is imposed on the right boundary.

      eta[right] = dirichlet(HR);

    Boundary conditions are set for the inflow velocity and the outflow of each layer.

      for (vector u in ul) {
    u.n[left] = dirichlet (uleft (point, _s, QL*(t < 10. ? t/10. : 1.)));
    u.n[right] = neumann(0.);
  }
#if STVT
  h[right] = dirichlet(HR);
#else
  for (scalar h in hl) {
    h[right]  = dirichlet (HR/nl);
  }
#endif

    In the non-hydrostatic case, a boundary condition is required for the non-hydrostatic pressure \phi of each layer.

    #if !HYDRO && !STVT
  for (scalar phi in phil) {
    phi[right] = dirichlet(0.);
  }
#endif
}

    We can optionally add horizontal viscosity.

    #if 0
event viscous_term (i++)
{
  // add horizontal viscosity (small influence)
  vector u;
  scalar w;
  for (u,w in ul,wl) {
    scalar d2u[];
    foreach()
      d2u[] = (u.x[1] + u.x[-1] - 2.*u.x[])/sq(Delta);
    foreach()
      u.x[] += dt*nu*d2u[];
    foreach()
      d2u[] = (w[1] + w[-1] - 2.*w[])/sq(Delta);
    foreach()
      w[] += dt*nu*d2u[];
  }
  boundary ((scalar *)ul);
  boundary ((scalar *)wl);
}
#endif

    We check for convergence.

    event logfile (t += 0.1; i <= 100000) {
#if STVT
  double dh = change (h, hc);
#else
  scalar H[];
  foreach() {
    H[] = 0.;
    for (scalar h in hl)
      H[] += h[];
  }
  double dh = change (H, hc);
#endif
  if (i > 0 && dh < 1e-5)
    return 1;
}

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

    #if 0
event gnuplot (i += 20) {
  static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
  if (i == 0)
    fprintf (fp, "set term x11\n");
  fprintf (fp,
	   "set title 'nl = %d, t = %.2f'\n"
	   "p [%g:%g][0:]'-' u 1:3:2 w filledcu lc 3 t '',"
	   " '' u 1:(-1):3 t '' w filledcu lc -1", nl, t,
	   X0, X0 + L0);
  int i = 4;
  for (scalar h in hl)
    fprintf (fp, ", '' u 1:%d w l lw 2 t ''", i++);
  fprintf (fp, "\n");
  foreach_leaf() {
    double H = 0.;
    for (scalar h in hl)
      H += h[];
    fprintf (fp, "%g %g %g", x, zb[] + H, zb[]);
    double z = zb[];
    for (scalar h in hl) {
      fprintf (fp, " %g", z);
      z += h[];
    }
    fprintf (fp, "\n");
  }
  fprintf (fp, "e\n\n");
  //  fprintf (fp, "pause 0.05\n");
  fflush (fp);
}
#endif

    Outputs

    At the end of the simulation we save the profiles.

    event profiles (t += 5)
{
  foreach_leaf() {
#if STVT
    double H = h[];
#else
    double H = 0.;
    for (scalar h in hl)
      H += h[];
#endif
    fprintf (stderr, "%g %g %g\n", x, zb[] + H, zb[]);
  }
  fprintf (stderr, "\n\n");
}

    We also save the velocity and non-hydrostatic pressure fields.

    event output (t = end)
{
#if HYDRO  
  scalar * wl = NULL;
  for (scalar h in hl) {
    scalar w = new scalar;
    wl = list_append (wl, w);
  }
  vertical_velocity (wl);
  foreach() {
    double wm = 0.;
    scalar h, w;
    for (h,w in hl,wl) {
      double w1 = w[];
      w[] = (w1 + wm)/2.;
      wm = w1;
    }
  }
  boundary (wl);
#endif // HYDRO
  foreach_leaf() {
    double z = zb[];
    vector u = ul[0];
#if HYDRO
    scalar w = wl[0], h = hl[0];
    printf ("%g %g %g %g\n", x, z, u.x[], w[]);
    for (u,h,w in ul,hl,wl) {
      z += h[];
      printf ("%g %g %g %g\n", x, z, u.x[], w[]);
    }
#elif STVT
    scalar w = wl[0];
    printf ("%g %g %g %g\n", x, z, u.x[], w[]);
    int l = 0;
    for (u,w in ul,wl) {
      z += layer[l++]*h[];
      printf ("%g %g %g %g\n", x, z, u.x[], w[]);
    }
#else // !STVT
    scalar w = wl[0], h = hl[0], phi = phil[0];
    printf ("%g %g %g %g %g\n", x, z, u.x[], w[], phi[]);
    for (u,h,w,phi in ul,hl,wl,phil) {
      z += h[];
      printf ("%g %g %g %g %g\n", x, z, u.x[], w[], phi[]);
    }
#endif // !STVT   
    printf ("\n");
  }
#if HYDRO
  delete (wl), free (wl);
#endif
  printf ("# end = %g\n", t);
}

    Results

    We compare the hydrostatic results obtained with the layered solver and those obtained with the Saint-Venant multilayer code.

    set yr [0.5:1.2]
set xlabel 'x'
set ylabel 'z'
plot '../gaussian-hydro/log' u 1:2 w l t 'Multilayer', \
     '../gaussian-stvt/log' u 1:2 w l t 'De Vita et al.'
    Evolution of the free surface (hydrostatic). (script)

    Evolution of the free surface (hydrostatic). (script)

    The results which follow are obtained with the non-hydrostatic solver.

    plot 'log' u 1:2 w l t ''
    Evolution of the free surface (non-hydrostatic) (script)

    Evolution of the free surface (non-hydrostatic) (script)

    We compare the results obtained with the layered solver to those obtained with the Navier-Stokes VOF solver (and without metric).

    plot 'log' w l t 'Multilayer', \
     '../../examples/gaussian-ns/log' w l t 'Navier-Stokes VOF'
    Evolution of the free surface for both solvers (script)

    Evolution of the free surface for both solvers (script)

    plot 'log' index 13 w l t 'Multilayer', \
     '../../examples/gaussian-ns/log' index 'prof70' w l t 'Navier-Stokes VOF', \
     'gaussian.nometric' w l t 'no metric'
    Final free surface profiles (script)

    Final free surface profiles (script)

    set term PNG enhanced font ",15" size 2048,512
set yr [0:1.2]
set output 'vel.png'
set pm3d
set pm3d map interpolate 10,4
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 )
set contour base
set cntrparam levels discrete 0
set cntrlabel onecolor
set cntrparam bspline
splot 'out' u 1:2:3 lt 3 lc rgb "#ffffff" lw 2
    Horizontal velocity field (script)

    Horizontal velocity field (script)

    set term PNG enhanced font ",15" size 2048,512
set output 'w.png'
# set cbrange [-0.1:0.1]
unset contour
splot 'out' u 1:2:4
    Vertical velocity field (script)

    Vertical velocity field (script)

    set term PNG enhanced font ",15" size 2048,512
set output 'phi.png'
# set cbrange [-0.2:0.35]
splot 'out' u 1:2:5
    Non-hydrostatic pressure (script)

    Non-hydrostatic pressure (script)

    See also