/**
# Flow in a rotating bottom-driven cylindrical container

A cylindrical container is partially filled with liquid. The bottom
disk spins at constant speed. A steady state is eventually reached for
which the hydrostatic pressure due to the deformation of the
free-surface balances the centrifugal force. This is analogous to the
solid body rotation in "Newton's bucket" but more complex due to the
development of boundary layers and recirculations generated by
the discontinuous velocity boundary condition between the rotating
bottom disk and the stationary walls.

This example reproduces one of the cases studied by [Wen Yang,
2018](#yang2018) in his PhD.

We use the two-phase axisymmetric Navier--Stokes solver with swirl. */

#include "grid/multigrid.h"
#include "axi.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/swirl.h"
#include "two-phase.h"
#include "tension.h"
#include "axistream.h"
#include "view.h"

/**
## Boundary conditions

The left boundary is the "bottom" of the container. The azimuthal
velocity component *w* on the bottom disk is $\Omega r$ with the
angular velocity $\Omega = 1$ and $r = y$. */

double Omega = 1.;
u.t[left] = dirichlet(0);
w[left]   = dirichlet(Omega*y);

/**
The "top" boundary is a no-slip stationary wall. */

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

/**
The cylinder wall is also no-slip and stationary. */

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

/**
## Parameters

The independent parameters are the aspect ratio of the cylindrical
container $H/R$ (unity here), the aspect ratio of the fluid layer
$$
G = h/R
$$
the Froude number
$$
Fr = \frac{R\Omega^2}{g}
$$
with $g$ the acceleration of gravity, the Reynolds number of the
liquid
$$
Re = \frac{R^2\Omega}{\nu_l}
$$
the Weber number
$$
We = \frac{\rho_l\Omega^2R^3}{\sigma}
$$
the ratio of gas to liquid densities and dynamic viscosities
$$
\rho_r = \rho_g/\rho_l
$$
$$
\mu_r = \mu_g/\mu_l
$$
The values correspond to those used by Yang for figure 5.9 in his PhD
(see also p. 106). */

double R = 1 [1],
       G = 0.25 [0],
       Fr = 0.88 [0],
       Re = 3063 [0],
       We = 3153 [0],
       rhor = 1.205/1.2107e3 [0],
       mur = 18.2e-6/6.09e-2 [0];

/**
## Setup

The domain is the square box of unit size (i.e. $R = 1$), and the
angular velocity $\Omega$ is one (see the boundary condition for *w*
above). The gives the values for viscosities, densities and surface
tension below. Note that surface tension is commented out since it has
very little influence but slows down the calculation. */

int main()
{
  N = 128;
  mu1 = sq(R)*Omega/Re;
  rho2 = rho1*rhor;
  mu2 = mu1*mur;
  //  f.sigma = 1./We;
  DT = 0.1 [0,1];
  run();
}

/**
We add gravity (given by the Froude number). */

event acceleration (i++)
{
  face vector av = a;
  foreach_face(x)
    av.x[] -= R*sq(Omega)/Fr;
}

/**
The initial interface is flat and positioned at $z = RG$. */

event init (i = 0) {
  fraction (f, R*G - x);
}

#if 0
event snapshots (i += 100) {
  scalar psi[];
  psi[left] = dirichlet (0);
  psi[right] = dirichlet (0);
  psi[top] = dirichlet (0);
  axistream (u, psi);
  dump();
}
#endif

/**
## Results

We output the timestep, the total kinetic energy and the total
mass. */

event logfile (i += 10)
{
  double ke = 0.;
  foreach (reduction(+:ke))
    ke += dv()*rho(f[])*(sq(u.x[]) + sq(u.y[]) + sq(w[]));
  fprintf (stderr, "%g %g %g %g\n", t, dt, ke, statsf (f).sum);
}

/**
We make a movie of the streamlines and isolines of the azimuthal
velocity component. They can be directly compared with figure 5.9 of
[Yang, 2018](#yang2018), reproduced below. Results are close but not
identical.

![Streamlines (left) and azimuthal velocity (right). Interface in
 red.](yang/movie.mp4)

![Streamlines and interface reproduced from [Yang, 2018](#yang2018),
 Figure 5.9.](yang/fig-5-9.svg)
*/

event movie (t += 1; t <= 300)
{
  view (fov = 22.7232, quat = {0,0,-0.707,0.707},
	tx = -0.12824, ty = -0.446981,
	width = 830, height = 432);
  box();

  /**
  We compute the axisymmetric streamfunction $\psi$ from the velocity
  components. */
  
  scalar psi[];
  psi[left] = dirichlet (0);
  psi[right] = dirichlet (0);
  psi[top] = dirichlet (0);
  axistream (u, psi);

  /**
  The values for the isolines are the same as in [Yang,
  2018](#yang2018), caption of Figure 5.9, p. 107. */
  
  squares ("psi", linear = true, spread = -1);
  isoline ("psi", n = 21, min = 0, max = 0.006);
  isoline ("psi", -0.0012, lc = {0,1,0});
  draw_vof ("f", lc = {1,0,0}, lw = 2);
  translate (y = -1.15) {
    box();
    squares ("w", linear = true, spread = -1);
    isoline ("w", n = 21, min = 0, max = 1);
    draw_vof ("f", lc = {1,0,0}, lw = 2);
  }
  save ("movie.mp4");
}

#if 0 // TREE
event adapt (i++) {
  adapt_wavelet ({u,w}, (double[]){1e-3,1e-3,1e-3}, minlevel = 5, maxlevel = 8);
}
#endif

/**
## References

~~~bib
@phdthesis{yang2018,
  TITLE = {Études expérimentales et numériques d’écoulements
           tournants à surface libre déformable},
  AUTHOR = {Wen Yang},
  SCHOOL = {{Sorbonne Universit\'e}},
  YEAR = {2018},
  MONTH = Dec
}
~~~
*/