viscosity-embed.h

Go to the documentation of this file.

/** @file viscosity-embed.h
 */
/**
# Viscous solver with embedded boundaries
 
We consider the Stokes equations
\f[
\rho \mathbf{u}_t = \rho \mathbf{g} +
\nabla \cdot [\mu ( \nabla \mathbf{u} + \nabla^T \mathbf{u})]
\f]
with \f$\mathbf{g}\f$ the acceleration and \f$T\f$ the transpose.
 
In the case of incompressible flow and uniform viscosity, the viscous term
can be further simplified. Since
\f[
\nabla \cdot (\mu  \nabla \mathbf{u}) =
(\nabla \mu) \cdot \nabla \mathbf{u} + \mu \nabla (\nabla \cdot  \mathbf{u}) = 0
\f]
the viscous term reduces to
\f[
\nabla \cdot [\mu ( \nabla \mathbf{u} + \nabla^T \mathbf{u})] =
\nabla \cdot (\mu \nabla^T \mathbf{u}) = \nabla^2 (\mu \mathbf{u})
\f]
and the equations for each velocity component are decoupled. For each
component \f$v_i\f$, the following discrete implicit equation is solved
using the multigrid solver
\f[
\frac{\Delta t} \nabla (\mu \nabla v_i^{n+1}) - (\rho + \lambda_i) v_i^{n+1}
+ \underbrace{(\rho v_i^{n} + g_i\, Delta t)}_{r_i} = 0
\f]
\f$\lambda_i\f$ is a possible extra term due to the metric. */
 
#include "poisson.h"
 
struct Viscosity {
  vector mu;
  scalar rho;
  double dt;
  double (* embed_flux) (Point, scalar, vector, double *);
};
 
#if AXI
# define lambda ((coord){0, dt*(mu.x[] + mu.x[1]      \
        + mu.y[] + mu.y[0,1])*sq(cs[])    \
  /(fm.x[] + fm.x[1] + fm.y[] + fm.y[0,1] + SEPS)/(cm[] + SEPS)})
 
#else // not AXI
# define lambda ((coord){0.,0.,0.})
#endif
 
// Note how the relaxation function uses "naive" gradients i.e. not
// the face_gradient_* macros.
 
static void relax_diffusion (scalar * a, scalar * b, int l, void * data)
{
  struct Viscosity * p = (struct Viscosity *) data;
  const vector mu = p->mu;
  const scalar rho = p->rho;
  double dt = p->dt;
  vector u = vector(a[0]), r = vector(b[0]);
 
  double (* embed_flux) (Point, scalar, vector, double *) = p->embed_flux;
  for (int _i = 0; _i < l, nowarning; _i++) {
    double avgmu = 0.;
    for (int _d = 0; _d < dimension; _d++)
      avgmu += mu.x[] + mu.x[1];
    avgmu = dt*avgmu + SEPS;
    for (int _d = 0; _d < dimension; _d++) {
      double c = 0.;
      double d = embed_flux ? embed_flux (point, u.x, mu, &c) : 0.;
      scalar s = u.x;
      double a = 0.;
      for (int _d = 0; _d < dimension; _d++)
  a += mu.x[1]*s[1] + mu.x[]*s[-1];
      u.x[] = (dt*a + (r.x[] - dt*c)*sq(Delta))/
  (sq(Delta)*(rho[] + lambda.x + dt*d) + avgmu);
    }
  }
  
#if TRASH
  vector u1[];
  for (int _i = 0; _i < l; _i++)
    for (int _d = 0; _d < dimension; _d++)
      u1.x[] = u.x[];
  /* trash: u */;
  for (int _i = 0; _i < l; _i++)
    for (int _d = 0; _d < dimension; _d++)
      u.x[] = u1.x[];
#endif
}
 
static double residual_diffusion (scalar * a, scalar * b, scalar * resl, 
          void * data)
{
  struct Viscosity * p = (struct Viscosity *) data;
  const vector mu = p->mu;
  const scalar rho = p->rho;
  double dt = p->dt;
  double (* embed_flux) (Point, scalar, vector, double *) = p->embed_flux;
  vector u = vector(a[0]), r = vector(b[0]), res = vector(resl[0]);
  double maxres = 0.;
#if TREE
  /* conservative coarse/fine discretisation (2nd order) */
  for (int _d = 0; _d < dimension; _d++) {
    scalar s = u.x;
    vector g[];
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      g.x[] = mu.x[]*face_gradient_x (s, 0);
    foreach (reduction(max:maxres), nowarning) {
      double a = 0.;
      for (int _d = 0; _d < dimension; _d++)
  a += g.x[] - g.x[1];
      res.x[] = r.x[] - (rho[] + lambda.x)*u.x[] - dt*a/Delta;
      if (embed_flux) {
  double c, d = embed_flux (point, u.x, mu, &c);
  res.x[] -= dt*(c + d*u.x[]);
      }
      if (fabs (res.x[]) > maxres)
  maxres = fabs (res.x[]);
    }
  }
#else
  /* "naive" discretisation (only 1st order on trees) */
  foreach (reduction(max:maxres), nowarning)
    for (int _d = 0; _d < dimension; _d++) {
      scalar s = u.x;
      double a = 0.;
      for (int _d = 0; _d < dimension; _d++)
  a += mu.x[0]*face_gradient_x (s, 0) - mu.x[1]*face_gradient_x (s, 1);
      res.x[] = r.x[] - (rho[] + lambda.x)*u.x[] - dt*a/Delta;
      if (embed_flux) {
  double c, d = embed_flux (point, u.x, mu, &c);
  res.x[] -= dt*(c + d*u.x[]);
      }
      if (fabs (res.x[]) > maxres)
  maxres = fabs (res.x[]);
    }
#endif
  return maxres;
}
 
#undef lambda
 
double TOLERANCE_MU = 0.; // default to TOLERANCE
 
trace
mgstats viscosity (vector u, vector mu, scalar rho, double dt,
       int nrelax = 4, scalar * res = NULL)
{
  vector r[];
  for (int _i = 0; _i < _N; _i++) /* foreach */
    for (int _d = 0; _d < dimension; _d++)
      r.x[] = rho[]*u.x[];
 
  restriction ({mu, rho});
  struct Viscosity p = { mu, rho, dt };
  p.embed_flux = u.x.boundary[embed] != antisymmetry ? embed_flux : NULL;
  return mg_solve ((scalar *){u}, (scalar *){r},
       residual_diffusion, relax_diffusion, &p, nrelax, res,
       minlevel = 1, // fixme: because of root level
                                  // BGHOSTS = 2 bug on trees
       tolerance = TOLERANCE_MU ? TOLERANCE_MU : TOLERANCE);
}