1 0 GfsSimulation GfsBox GfsGEdge {} {
  Time { iend = 100 dtmax = 1e-2 }
  Refine 6
  Solid (ellipse (0,0,0.25,0.25)) { flip = 1 scale = 1.9999 }
  Solid (ellipse (0,0,0.25,0.25))
  ApproxProjectionParams { tolerance = 1e-6 }
  AdvectionParams { scheme = none }
  SourceViscosity {} {
    double mu, ty, N, mumax = 1000.;
    double m;

    switch (MODEL) {
    case 0: /* Newtonian */
      mu = 1.; ty = 0.; N = 1.; break;
    case 1: /* Power-law (shear-thinning) */
      mu = 0.08; ty = 0.; N = 0.5; break;
    case 2: /* Herschel-Bulkley */
      mu = 0.0672; ty = 0.12; N = 0.5; break;
    case 3: /* Bingham */
      mu = 1.; ty = 10.; N = 1.; break;
    }
    if (D2 > 0.)
      m = ty/(2.*D2) + mu*exp ((N - 1.)*log (D2));
    else {
      if (ty > 0. || N < 1.) m = mumax;
      else m = N == 1. ? mu : 0.;
    }
    return MIN (m, mumax);
  } {
    # Crank-Nicholson does not converge for these cases, we need backward Euler
    # (beta = 0.5 -> Crank-Nicholson, beta = 1 -> backward Euler)
    beta = 1
  }
  SurfaceBc U Dirichlet (x*x + y*y > 0.140625 ? 0. : - ay)
  SurfaceBc V Dirichlet (x*x + y*y > 0.140625 ? 0. :   ax)
  EventStop { istep = 1 } U 1e-4 DU

  OutputScalarNorm { istep = 1 } du-MODEL { v = DU }
  OutputLocation { start = end } { awk '{if ($1 != "#") print $2,$8;}' > prof-MODEL } profile
}
GfsBox {}