/**
# Boundary layer on a rotating disk

[Von Kármán, 1921](#karman1921) showed that the steady flow of an
incompressible liquid of kinematic viscosity $\nu$ induced by an
infinite plane disk rotating at angular velocity $\Omega$ can be
described by a similarity solution. In effect, using
$\zeta=z\sqrt{\Omega/\nu}$ and setting the axial velocity $U$,
radial velocity $V$ and azimuthal velocity $W$ as
$$
U=\sqrt{\nu \Omega} F(\zeta) \quad V=\Omega r H(\zeta) 
\quad W = \Omega r G(\zeta)
$$
the Navier-Stokes equations reduce to a couple of ODEs:
$$
F'''-F\, F'' +F'^2/2 +2G^2  = 0 \quad \mathrm{and} \quad G''-F\,G'+G\,F' = 0
$$
with boundary conditions
$$
F(0)=F'(0)=0 \: G(0)=1. \quad \mathrm{and} \quad F'(\infty)=G(\infty)=0.
$$
where the prime denotes differentiation with respect to $\zeta$.

To reproduce this solution numerically, we use the axisymmetric
Navier--Stokes solver with azimuthal velocity (swirl). */

#include "grid/multigrid.h"
#include "axi.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/swirl.h"

/**
The left boundary is the rotating disk with $\Omega = 1$ and a no-slip
condition for the tangential velocity i.e. */
  
u.t[left] = dirichlet(0);
w[left]   = dirichlet(y);

/**
We use an open (outflow) boundary condition for the right boundary. */

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

/**
The top boundary condition is more tricky but the following seems to
work. */

u.n[top] = neumann(0);
p[top] = neumann(0);

/**
We use a constant viscosity but it needs to be weighted by the
(axisymmetric) metric. */

face vector muv[];

event properties (i++) {
  foreach_face()
    muv.x[] = 0.2*fm.x[];
}

/**
The computational domain is $12\times 12$ and we limit the
timestep. */

int main()
{
  size (12);
  N = 128;
  mu = muv;
  DT = 2e-2;
  run();
}

/**
We wait until the boundary layer is fully developed and
quasi-stationary. We only consider values close to the origin to
minimize the influence of boundaries (von Kármán's solution is valid
in an infinite domain). */

event end (t = 20)
{
  foreach()
    if (x*x + y*y < 8)
      fprintf (stderr, "%g %g %.4g %.4g\n", x, y, u.x[], w[]);
}

/**
~~~gnuplot Axial $F(\zeta)$ and azimuthal $G(\zeta)$ dimensionless velocity components
set xlabel '{/Symbol z}'
set key center right
nu = 0.2
plot [0:6]'analytical' u 1:2 w l t '-F({/Symbol z})',	     \
          'log' u ($1/sqrt(nu)):(-$3/sqrt(nu)) t 'Basilisk', \
          'analytical' u 1:3 w l t 'G({/Symbol z})',	     \
	  'log' u ($1/sqrt(nu)):($4/$2) t 'Basilisk'
~~~

## References

~~~bib
@article{karman1921,
  title={{\"U}ber laminare und turbulente Reibung},
  author={Karman, Th V},
  journal={ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift 
           f{\"u}r Angewandte Mathematik und Mechanik},
  volume={1},
  number={4},
  pages={233--252},
  year={1921},
  publisher={Wiley Online Library}
}
~~~

## See also

* [Same case with 
   Gerris](http://gerris.dalembert.upmc.fr/gerris/tests/tests/swirl.html)
*/