src/axistream.h

    Axisymmetric stream function

    The axisymmetric or Stokes stream function \psi verifies \displaystyle u_r = -\frac{1}{r}\partial_z\psi \displaystyle 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 \displaystyle \omega = \partial_zu_r - \partial_ru_z \psi verifies the variable-coefficient Poisson equation \displaystyle \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 \displaystyle \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);
      boundary_flux ({g});
      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
      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 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);
      boundary ({psi});
      return mg_solve ({psi}, {omega}, residual_psi, relax_psi, NULL, 0, NULL, 1);
    }

    See also

    Usage

    Examples