```
# The strange choices for radii R1,R2 and eccentricity ECC come from
# the 'bipolar' variant
Define R1 (1./sinh(1.5))
Define R2 (1./sinh(1.))
Define X1 (1./tanh(1.5))
Define X2 (1./tanh(1.))
Define ECC (X2 - X1)

1 0 GfsSimulation GfsBox GfsGEdge {} {
PhysicalParams { L = 2.5 }
# Upper bound on time, it should converge much before this
Time { end = 100 }
Refine LEVEL
Solid (- ellipse (0.,ECC,R2,R2))
Solid (ellipse (0.,0.,R1,R1))
ApproxProjectionParams { tolerance = 1e-6 }
AdvectionParams { scheme = none }
SourceViscosity 1
# Fixed outer cylinder and rotating inner cylinder (tangential velocity unity)
SurfaceBc U Dirichlet (x*x + y*y > 1.5*R1*R1 ? 0. : - ay/R1)
SurfaceBc V Dirichlet (x*x + y*y > 1.5*R1*R1 ? 0. :   ax/R1)
# Stop when stationnary
EventStop { step = 1e-2 } U 1e-4 DU

Global {
#include "wannier.c"
double psi (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return p;
}
double ux (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return u;
}
double uy (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return v;
}
}

OutputScalarNorm { istep = 1 } du { v = DU }
OutputErrorNorm { start = end } { awk '{ print LEVEL,\$5,\$7,\$9 }' } { v = Velocity } {
s = {
double p, u, v;
psiuv (x, y - ECC, R1, R2, ECC, 1., 0., &u, &v, &p);
return sqrt (u*u + v*v);
}
v = EU
}
OutputSimulation { start = end } end-LEVEL.gfs
}
GfsBox {}
```