/**
# Transcritical flow over a bump

This is the Navier-Stokes/VOF version of [this test
case](/src/test/gaussian.c). A viscous fluid is injected to the left
of a Gaussian bump and exits the domain with a fixed water depth. This
creates a subcritical flow at inflow and outflow (i.e. the phase speed
of gravity waves is larger than the fluid velocity) and a
supercritical flow over the bump. In the absence of non-hydrostatic
terms this would create a viscous hydraulic jump (see
[/src/test/gaussian.c]()). When non-hydrostatic terms are included
(i.e. the full free-surface Navier-Stokes are solved), this leads to a
steady wavetrain of dispersive waves.

![Horizontal velocity field](gaussian-ns/u.x.png){ width=100% }

![Vertical velocity field](gaussian-ns/u.y.png){ width=100% }

See [Popinet (2020)](/Bibliography#popinet2020) for a more detailed
discussion. */

#include "embed.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "reduced.h"
#include "view.h"
#include "navier-stokes/perfs.h"

const double H0 = 0.6, QL = 1., BA = 0.4, NU = 1e-2, T0 = 10.;

double hl = H0;

u.n[left] = dirichlet((t < T0 ? t/T0 : 1.)*3./2.*QL/hl*(y < hl ? 1. - sq(y/hl - 1.) : 1.));
p[left] = neumann(0);
pf[left] = neumann(0);

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

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

int main()
{
  size (30 [1]);
  N = 32;
  rho1 = 0.001 [0];
  rho2 = 1. [0];
  mu2 = NU;
  mu1 = mu2*rho1/rho2/10.;
  G.y = - 9.81;
  Z.y = H0;
  run();
}

event init (i = 0)
{
  refine (y < 1.5 && level < 10);
  
  fraction (f, y - H0);

  double a = 5.;
  solid (cs, fs, y - BA*exp(- sq(x - 10.)/a) - 1e-3);
}

event update_hl (i++)
{
  hl = 0.;
  foreach_boundary (left, reduction(+:hl))
    hl += Delta*f[];
  hl = L0 - hl;
  printf ("%g %g\n", t, hl);
}

event snapshot (i += 10; t <= 70) {
  p.nodump = false;
  dump();
}

event maxdt (t <= 10.; t += 0.05);

event profiles (t += 5)
{
  fprintf (stderr, "# prof%g\n", t);
  output_facets (f, stderr);
}

event pictures (t = end)
{
  view (fov = 4.04484, tx = -0.498476, ty = -0.0923365, sy = 5,
	bg = {1,1,1},
	width = 1869, height = 390);
  draw_vof ("cs", filled = -1, fc = {1,1,1});
  draw_vof ("f", filled = 1, fc = {1,1,1});
  squares ("u.x", min = -0.5, max = 4, linear = true);
  isoline ("u.x", 0., lc = {1,1,1}, lw = 2);
  save ("u.x.png");

  draw_vof ("cs", filled = -1, fc = {1,1,1});
  draw_vof ("f", filled = 1, fc = {1,1,1});
  squares ("u.y", min = -0.8, max = 0.8, linear = true);
  save ("u.y.png");  
}