# Axisymmetric stream function

The axisymmetric or Stokes stream function $\psi$ verifies ${u}_{r}=-\frac{1}{r}{\partial }_{z}\psi$ ${u}_{z}=+\frac{1}{r}{\partial }_{r}\psi$ with ${u}_{r}$ and ${u}_{z}$ the radial and longitudinal velocity components of an incompressible axisymmetric flow.

Given the definition of the vorticity $\omega$ $\omega ={\partial }_{z}{u}_{r}-{\partial }_{r}{u}_{z}$ $\psi$ verifies the variable-coefficient Poisson equation ${\partial }_{z}\left(\frac{1}{r}{\partial }_{z}\psi \right)+{\partial }_{r}\left(\frac{1}{r}{\partial }_{r}\psi \right)=-\omega$ This is not well-conditioned due to the divergence of the coefficients at $r=0$ and can be rewritten instead as $\frac{{\partial }^{2}\psi }{\partial {z}^{2}}+\frac{{\partial }^{2}\psi }{\partial {r}^{2}}-\frac{1}{r}{\partial }_{r}\psi =-\omega r$ 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];
foreach_level_or_leaf (l)
a[] = (- sq(Δ)*b[] - Δ*(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) */
face vector g[];
foreach_face()
boundary_flux ({g});
foreach (reduction(max:maxres)) {
res[] = b[] + (a[0,1] - a[0,-1])/(2.*y*Δ);
foreach_dimension()
res[] -= (g.x[1] - g.x[])/Δ;
#else // !TREE
/* "naive" discretisation (only 1st order on trees) */
foreach (reduction(max:maxres)) {
res[] = b[] + (a[0,1] - a[0,-1])/(2.*y*Δ);
foreach_dimension()
res[] += (face_gradient_x (a, 0) -
#endif // !TREE
if (fabs (res[]) > maxres)
maxres = fabs (res[]);
}
boundary (resl);
return maxres;
}``````

The function axistream() takes as input the $u$ velocity field (with $x=z$ and $y=r$) and returns the corresponding axisymmetric stream function.

``````mgstats axistream (vector u, scalar ψ)
{
scalar ω[];
foreach() {
ω[] = y*(u.y[1] - u.y[-1] - u.x[0,1] + u.x[0,-1])/(2.*Δ);
ψ[] = 0.;
}
ψ[bottom] = dirichlet(0);
boundary ({ψ});
return mg_solve ({ψ}, {ω}, residual_psi, relax_psi, NULL, 0, NULL, 1);
}``````