src/axistream.h
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_zu_r - \partial_ru_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(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) */
face vector g[];
foreach_face()
g.x[] = face_gradient_x (a, 0);
foreach (reduction(max:maxres)) {
res[] = b[] + (a[0,1] - a[0,-1])/(2.*y*Delta);
foreach_dimension()
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);
foreach_dimension()
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 u velocity field (with x=z and y=r) and returns the corresponding axisymmetric stream function.
mgstats axistream (vector u, scalar psi)
{
scalar omega[];
foreach() {
omega[] = y*(u.y[1] - u.y[-1] - u.x[0,1] + u.x[0,-1])/(2.*Delta);
psi[] = 0.;
}
psi[bottom] = dirichlet(0);
return mg_solve ({psi}, {omega}, residual_psi, relax_psi, NULL, 0, NULL, 1);
}