diffusion.h

Go to the documentation of this file.

/** @file diffusion.h
 */
/**
# Time-implicit discretisation of reaction--diffusion equations
 
We want to discretise implicitly the reaction--diffusion equation
\f[
\theta\partial_tf = \nabla\cdot(D\nabla f) + \beta f + r
\f] 
where \f$\beta f + r\f$ is a reactive term,  \f$D\f$ is the diffusion
coefficient and \f$\theta\f$ can be a density term.
 
Using a time-implicit backward Euler discretisation, this can be
written
\f[
\theta\frac{f^{n+1} - f^{n}}{dt} = \nabla\cdot(D\nabla f^{n+1}) + \beta
f^{n+1} + r
\f]
Rearranging the terms we get
\f[
\nabla\cdot(D\nabla f^{n+1}) + (\beta - \frac{\theta}{dt}) f^{n+1} =
- \frac{\theta}{dt}f^{n} - r
\f]
This is a Poisson--Helmholtz problem which can be solved with a
multigrid solver. */
 
#include "poisson.h"
 
/**
The parameters of the `diffusion()` function are a scalar field `f`,
scalar fields `r` and \f$\beta\f$ defining the reactive term, the timestep
`dt` and a vector field containing the diffusion coefficient
`D`. If `D` or \f$\theta\f$ are omitted they are set to one. If \f$\beta\f$ is
omitted it is set to zero. Both `D` and \f$\beta\f$ may be constant
fields.
 
Note that the `r`, \f$\beta\f$ and \f$\theta\f$ fields will be modified by the solver.
 
The function returns the statistics of the Poisson solver. */
 
trace
mgstats diffusion (scalar f, double dt,
       vector D = {{-1}},  // default 1
       scalar r = {-1},         // default 0
       scalar beta = {-1},      // default 0
       scalar theta = {-1})     // default 1
{
 
  /**
  If *dt* is zero we don't do anything. */
 
  if (dt == 0.) {
    mgstats s = {0};
    return s;
  }
 
  /**
  We define \f$f\f$ and \f$r\f$ for convenience. */
 
  scalar ar = automatic (r);
 
  /**
  We define a (possibly constant) field equal to \f$\theta/dt\f$. */
 
  const scalar idt[] = - 1./dt;
  const scalar theta_idt = theta.i >= 0 ? theta : idt;
  
  if (theta.i >= 0)
    for (int _i = 0; _i < _N; _i++) /* foreach */
      theta[] *= idt[];
 
  /**
  We use `r` to store the r.h.s. of the Poisson--Helmholtz solver. */
 
  if (r.i >= 0)
    for (int _i = 0; _i < _N; _i++) /* foreach */
      ar[] = theta_idt[]*f[] - ar[];
  else // r was not passed by the user
    for (int _i = 0; _i < _N; _i++) /* foreach */
      ar[] = theta_idt[]*f[];
 
  /**
  If \f$\beta\f$ is provided, we use it to store the diagonal term \f$\lambda\f$. */
 
  scalar lambda = theta_idt;
  if (beta.i >= 0) {
    for (int _i = 0; _i < _N; _i++) /* foreach */
      beta[] += theta_idt[];
    lambda = beta;
  }
 
  /**
  Finally we solve the system. */
 
  return poisson (f, ar, D, lambda);
}