axistream.h

Go to the documentation of this file.

/** @file axistream.h
 */
/**
# Axisymmetric stream function
 
The axisymmetric or [Stokes stream
function](https://en.wikipedia.org/wiki/Stokes_stream_function) \f$\psi\f$
verifies
\f[
u_r = -\frac{1}{r}\partial_z\psi
\f]
\f[
u_z = +\frac{1}{r}\partial_r\psi
\f]
with \f$u_r\f$ and \f$u_z\f$ the radial and longitudinal velocity components
of an incompressible axisymmetric flow.
 
Given the definition of the vorticity \f$\omega\f$
\f[
\omega = \partial_zu_r - \partial_ru_z
\f]
\f$\psi\f$ verifies the variable-coefficient Poisson equation
\f[
\partial_z\left(\frac{1}{r}\partial_z\psi\right) +
\partial_r\left(\frac{1}{r}\partial_r\psi\right) = - \omega
\f]
This is not well-conditioned due to the divergence of the coefficients
at \f$r = 0\f$ and can be rewritten instead as
\f[
\frac{\partial^2\psi}{\partial z^2} + 
\frac{\partial^2\psi}{\partial r^2} - \frac{1}{r}\partial_r\psi 
= - \omega r
\f]
which is better behaved.
 
This equation can be inverted using the multigrid solver combined with
the relaxation and residual functions defined below. */
 
#include "poisson.h"
 
static void relax_psi (scalar * al, scalar * bl, int l, void * data)
{
  scalar a = al[0], b = bl[0];
  for (int _i = 0; _i < l; _i++)
    a[] = (- sq(Delta)*b[] - Delta*(a[0,1] - a[0,-1])/(2.*y) +
     a[1] + a[-1] + a[0,1] + a[0,-1])/4.;
}
 
static double residual_psi (scalar * al, scalar * bl, scalar * resl, void * data)
{
  scalar a = al[0], b = bl[0], res = resl[0];
  double maxres = 0.;
#if TREE
  /* conservative coarse/fine discretisation (2nd order) */
  vector g[];
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    g.x[] = face_gradient_x (a, 0);
  foreach (reduction(max:maxres)) {
    res[] = b[] + (a[0,1] - a[0,-1])/(2.*y*Delta);
    for (int _d = 0; _d < dimension; _d++)
      res[] -= (g.x[1] - g.x[])/Delta;
    if (fabs (res[]) > maxres)
      maxres = fabs (res[]);
  }
#else // !TREE
  /* "naive" discretisation (only 1st order on trees) */
  foreach (reduction(max:maxres)) {
    res[] = b[] + (a[0,1] - a[0,-1])/(2.*y*Delta);
    for (int _d = 0; _d < dimension; _d++)
      res[] += (face_gradient_x (a, 0) -
    face_gradient_x (a, 1))/Delta;  
    if (fabs (res[]) > maxres)
      maxres = fabs (res[]);
  }
#endif // !TREE
  return maxres;
}
 
/**
The function *axistream()* takes as input the \f$u\f$ velocity field (with
\f$x=z\f$ and \f$y=r\f$) and returns the corresponding axisymmetric stream
function. */
 
mgstats axistream (vector u, scalar psi)
{
  scalar omega[];
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    omega[] = y*(u.y[1] - u.y[-1] - u.x[0,1] + u.x[0,-1])/(2.*Delta);
    psi[] = 0.;
  }
  /* BC: psi[bottom] = dirichlet(0) */
  return mg_solve ({psi}, {omega}, residual_psi, relax_psi, NULL, 0, NULL, 1);
}
 
/**
## See also
 
* [Axisymmetric metric](axi.h)
*/