src/test/axi.c

    Convergence of axisymmetric viscous terms

    We wish to test the accuracy of the discretisation of viscous terms in axisymmetric coordinates. To do so, we “manufacture” a solution to the pure viscous diffusion problem \displaystyle \partial_t (yu) = \partial_x (y \tau_{x x}) + \partial_y (y \tau_{x y}) + s_x \displaystyle \partial_t (yv) = \partial_x (y \tau_{y x}) + \partial_y (y \tau_{y y}) + s_y - 2 \frac{v}{y} with \displaystyle \tau_{x x} = 2\partial_x u \displaystyle \tau_{x y} = \tau_{y x} = (\partial_y u + \partial_x v) \displaystyle \tau_{y y} = 2\partial_y v and where s_x and s_y are source terms to be determined. We look for solutions of the form \displaystyle u (x, y) = v (x, y) = y^a \cos x \sin y The viscous stress tensor is then (for \mu = 1) \displaystyle \tau_{x x} = - 2 y^a \sin x \sin y \displaystyle \tau_{x y} = \tau_{y x} = \cos x (y^a \cos y + ay^{a - 1} \sin y) - y^a \sin x \sin y \displaystyle \tau_{y y} = 2 \cos x (y^a \cos y + ay^{a - 1} \sin y) and the viscosity equation gives \displaystyle s_x = y^a [3 y \cos x \sin y - (2 a + 1) \cos x \cos y - a^2 y^{- 1} \cos x \sin y + y \sin x \cos y + (a + 1) \sin x \sin y] \displaystyle s_y = y^a [y \sin x \cos y + a \sin x \sin y + 3 y \cos x \sin y - 2 (2 a + 1) \cos x \cos y + 2 (1 - a^2) y^{- 1} \cos x \sin y] The game is then to solve the viscosity equation with the forcing terms defined above and recover the (stationary) solution for the velocity field.

    We use the axisymmetric metric, the viscous solver and the standard time loop.

    #include "grid/multigrid.h"
#include "axi.h"
#include "viscosity.h"
#include "run.h"

    We impose the correct boundary conditions for u (the default boundary conditions for v are already correct).

    vector u[];

u.t[top] = dirichlet(0);
u.t[bottom] = dirichlet(0);

int main()
{
  L0 = 2.*pi;
  periodic (right);
  for (N = 16; N <= 128; N *= 2)
    run();
}

    We solve only for the case a=1. The case a=0 is ill-posed because the forcing terms diverge near the axis of symmetry.

    double a = 1.;
mgstats mg;

event init (i = 0) {
  foreach()
    u.x[] = u.y[] = 0.;
}

event integration (i++)
{
  double dt = 1.;

    This is the time integration loop. We first add the forcing term…

      foreach() {
    u.x[] += dt*pow(y, a - 1.)*(3.*y*cos(x)*sin(y) - (2.*a + 1.)*cos(x)*cos(y)
				- sq(a)*cos(x)*sin(y)/y	
				+ y*sin(x)*cos(y) + (a + 1.)*sin(x)*sin(y));
    u.y[] += dt*pow(y, a - 1.)*(y*sin(x)*cos(y) + a*sin(x)*sin(y)
				+ 3.*y*cos(x)*sin(y)
				- 2.*(2.*a + 1.)*cos(x)*cos(y) +
				2.*(1. - sq(a))*cos(x)*sin(y)/y);
  }

    …and then solve for viscosity.

      mg = viscosity (u, fm, cm, dt);
}

event error (i = 20)
{

    Twenty timesteps are enough to converge toward the stationary solution.

      scalar e[];
  foreach() {
    printf ("%g %g %g %g %g\n", x, y, u.x[], u.y[], pow(y,a)*cos(x)*sin(y));
    e[] = u.x[] - pow(y,a)*cos(x)*sin(y);
  }
  norm n = normf (e);
  fprintf (stderr, "%d %g %g %g %d %d\n",
	   N, n.avg, n.rms, n.max, mg.i, mg.nrelax);
}

    As expected we get second-order convergence.

    ftitle(a,b) = sprintf("%.0f/x^{%4.2f}", exp(a), -b)
f(x)=a+b*x
fit f(x) 'log' u (log($1)):(log($4)) via a,b
f2(x)=a2+b2*x
fit f2(x) 'log' u (log($1)):(log($2)) via a2,b2
set xlabel 'Resolution'
set ylabel 'Error'
set logscale
set xrange [8:256]
set cbrange [1:2]
set xtics 8,2,256
set grid ytics
set yrange [1e-5:]
plot 'log' u 1:4 t 'max', exp(f(log(x))) t ftitle(a,b), \
     'log' u 1:2 t 'norm1', exp(f2(log(x))) t ftitle(a2,b2)
    Convergence (script)

    Convergence (script)