diffusion.h

Go to the documentation of this file.

/** @file diffusion.h
 */
/**
# Vertical diffusion
 
We consider the vertical diffusion of a tracer \f$s\f$ with a diffusion
coefficient \f$D\f$ for the multilayer solver.
 
For stability, we discretise the vertical diffusion equation implicitly as
\f[
\frac{(hs_l)^{n + 1} - (hs_l)^{\star}}{\Delta t} =
D \left( \frac{s_{l + 1} - s_l}{h_{l + 1 / 2}} -
\frac{s_l - s_{l - 1}}{h_{l - 1 / 2}} \right)^{n + 1}
\f]
which can be expressed as the linear system
\f[
\mathbf{Ms}^{n + 1} = \mathrm{rhs}
\f]
where \f$\mathbf{M}\f$ is a 
[tridiagonal matrix](https://en.wikipedia.org/wiki/Tridiagonal_matrix). 
The lower, principal and upper diagonals are *a*, *b* and *c* respectively. 
 
Boundary conditions on the top and bottom layers need to be added to close the
system. We chose to impose a Neumann condition on the free-surface i.e.
\f[
\partial_z s |_t = \dot{s}_t
\f]
and a Navier slip condition on the bottom i.e.
\f[
s|_b = s_b + \lambda_b \partial_z s|_b
\f] */
 
void vertical_diffusion (Point point, scalar h, scalar s, double dt, double D,
       double dst, double s_b, double lambda_b)
{
  double a[nl], b[nl], c[nl], rhs[nl];
 
  /**
  The *rhs* of the tridiagonal system is \f$h_l s_l\f$. */
      
  foreach_layer()
    rhs[point.l] = s[]*h[];
 
  /**
  The lower, principal and upper diagonals \f$a\f$, \f$b\f$ and \f$c\f$ are given by
  \f[
  a_{l > 0} = - \left( \frac{D \Delta t}{h_{l - 1 / 2}} \right)^{n + 1}
  \f]
  \f[
  c_{l < \mathrm{nl} - 1} = - \left( \frac{D \Delta t}{h_{l + 1 / 2}}
  \right)^{n + 1}
  \f]
  \f[
  b_{0 < l < \mathrm{nl} - 1} = h_l^{n + 1} - a_l - c_l
  \f]
  */
 
  for (int l = 1; l < nl - 1; l++) {
    a[l] = - 2.*D*dt/(h[0,0,l-1] + h[0,0,l]);
    c[l] = - 2.*D*dt/(h[0,0,l] + h[0,0,l+1]);
    b[l] = h[0,0,l] - a[l] - c[l];
  }
    
  /**
  For the top layer the boundary conditions give the (ghost)
  boundary value
  \f[
  s_{\mathrm{nl}} = s_{\mathrm{nl} - 1} + \dot{s}_t h_{\mathrm{nl} - 1},
  \f]
  which gives the diagonal coefficient and right-hand-side
  \f[
  b_{\mathrm{nl} - 1} = h_{\mathrm{nl} - 1}^{n + 1}
  - a_{\mathrm{nl} - 1}
  \f]
  \f[
  \mathrm{rhs}_{\mathrm{nl} - 1} = 
  (hs)_{\mathrm{nl} - 1}^{\star} + D \Delta t \dot{s}_t
  \f]
  */
 
  a[nl-1] = - 2.*D*dt/(h[0,0,nl-2] + h[0,0,nl-1]);
  b[nl-1] = h[0,0,nl-1] - a[nl-1];
  rhs[nl-1] += D*dt*dst;
 
  /**
  For the bottom layer a third-order discretisation of the Navier slip
  condition gives
  \f[
  \begin{aligned}
  b_0 & = h_0 + 2 \Delta t D \left( \frac{1}{h_0 + h_1} + \frac{h^2_1 + 3
  h_0 h_1 + 3 h^2_0}{\det} \right),\\
  c_0 & = - 2 \Delta t D \left( \frac{1}{h_0 + h_1} + \frac{h^2_0}{\det}
  \right),\\
  \text{rhs}_0 & = (hs_0)^{\star} + 2 \Delta t D s_b  \frac{h^2_1 + 3 h_0
  h_1 + 2 h^2_0}{\det},\\
  \det & = h_0 (h_0 + h_1)^2  + 2\lambda (3\,h_0 h_1 + 2\,h_0^2 + h_1^2),
  \end{aligned}
  \f]
  */
 
  double den = h[]*sq(h[] + h[0,0,1]) 
    + 2.*lambda_b*(3.*h[]*h[0,0,1] + 2.*sq(h[]) + sq(h[0,0,1]));
  b[0] = h[] + 2.*dt*D*(1./(h[] + h[0,0,1]) +
        (sq(h[0,0,1]) + 3.*h[]*h[0,0,1] + 3.*sq(h[]))/den);
  c[0] = - 2.*dt*D*(1./(h[] + h[0,0,1]) + sq(h[])/den);
  rhs[0] += 2.*dt*D*s_b*(sq(h[0,0,1]) + 3.*h[]*h[0,0,1] + 2.*sq(h[0]))/den;
 
  if (nl == 1) {
    b[0] += c[0];
    rhs[0] += (- c[0]*h[] - D*dt) * dst;
  }
    
  /**
  We can now solve the tridiagonal system using the [Thomas
  algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm). */
  
  for (int l = 1; l < nl; l++) {
    b[l] -= a[l]*c[l-1]/b[l-1];
    rhs[l] -= a[l]*rhs[l-1]/b[l-1];
  }
  a[nl-1] = rhs[nl-1]/b[nl-1];
  s[0,0,nl-1] = a[nl-1];
  for (int l = nl - 2; l >= 0; l--)
    s[0,0,l] = a[l] = (rhs[l] - c[l]*a[l+1])/b[l];
}
 
/**
# Viscous friction between layers
 
By default the viscosity is zero and we impose free-slip on the
free-surface and no-slip on the bottom boundary
i.e. \f$\dot{\mathbf{u}}_t = 0\f$, \f$\mathbf{\lambda}_b = 0\f$, \f$\mathbf{u}_b
= 0\f$. */
 
double nu = 0.;
const vector lambda_b[] = {0,0,0}, dut[] = {0,0,0}, u_b[] = {0,0,0};
 
/**
In the [layered solver](hydro.h), vertical viscosity is applied to the
velocity field just after advection, but before the pressure
gradient/acceleration term is applied. To take the pressure gradient
into account, we first apply the acceleration of the previous
timestep, apply vertical viscosity and then substract the previous
acceleration. */
 
/** @brief Event: viscous_term (i++,last) */
void event_viscous_term (void) 
{
  if (nu > 0.) {
    for (int _i = 0; _i < _N; _i++) /* foreach */ {
      foreach_layer()
  for (int _d = 0; _d < dimension; _d++)
    u.x[] += dt*(ha.x[] + ha.x[1])/(hf.x[] + hf.x[1] + dry);
      for (int _d = 0; _d < dimension; _d++)
  vertical_diffusion (point, h, u.x, dt, nu,
          dut.x[], u_b.x[], lambda_b.x[]);
      foreach_layer()
  for (int _d = 0; _d < dimension; _d++)
    u.x[] -= dt*(ha.x[] + ha.x[1])/(hf.x[] + hf.x[1] + dry);
    }
  }
}
 
/**
## Horizontal diffusion
 
This approximates
\f[
h \partial_t s = D \nabla \cdot (h \nabla s)
\f]
with \f$D\f$ the diffusion coefficient. Note that metric terms linked to
the slope of the layers are not taken into account. Note also that the
time discretisation is explicit so that the timestep must be limited
(manually) by \f$\min(\Delta^2/D)\f$. */
 
void horizontal_diffusion (scalar * list, double D, double dt)
{
  if (D > 0.) {
    scalar * d2sl = list_clone (list);
    foreach_layer() {
      for (int _i = 0; _i < _N; _i++) /* foreach */ {
  scalar s, d2s;
  for (s,d2s in list,d2sl) {
    double a = 0.;
    for (int _d = 0; _d < dimension; _d++)
      a += (hf.x[]*fm.x[]/(cm[-1] + cm[])*(s[-1] - s[]) +
      hf.x[1]*fm.x[1]/(cm[1] + cm[])*(s[1] - s[]));
    d2s[] = 2.*a/(cm[]*sq(Delta));
        }
      }
      for (int _i = 0; _i < _N; _i++) /* foreach */
  if (h[] > dry) {
    scalar s, d2s;
    for (s,d2s in list,d2sl)
      s[] += dt*D*d2s[]/h[];
  }
    }
    delete (d2sl);
    free (d2sl);
  }
}
 
/**
### Time-implicit integration
 
When the timestep constraint due to explicit time integration is too
strong, one can use the time-implicit version below. The diffusion equation
\f[
h \partial_t s = D \nabla \cdot (h \nabla s)
\f]
is discretised implicitly as
\f[
h \frac{s^{n + 1} - s^n}{\Delta t} = \nabla \cdot (Dh \nabla s^{n + 1})
\f]
which can be reordered as
\f[
\nabla \cdot (\Delta tDh \nabla s^{n + 1}) - hs^{n + 1} = - hs^n
\f]
i.e. a Poisson--Helmoltz equation of the form
\f[
\nabla \cdot (\alpha \nabla s^{n + 1}) + \lambda s^{n + 1} = b
\f]
which is solved using the [multigrid Poisson--Helmholz solver](/src/poisson.h#user-interface).
 
The function returns convergence statistics for the multigrid solver. */
 
#include "poisson.h"
 
trace
mgstats implicit_horizontal_diffusion (scalar * list, double D, double dt)
{
  mgstats mg = {0};
  if (D > 0.) {
    scalar b[], lambda[];
    vector alpha[];
    foreach_layer() {
      for (int _i = 0; _i < _N; _i++) /* foreach_face */
  alpha.x[] = 2.*dt*D*hf.x[]*fm.x[]/(cm[-1] + cm[]);
      for (int _i = 0; _i < _N; _i++) /* foreach */
  lambda[] = - cm[]*h[];
      for (int _s = 0; _s < 1; _s++) /* scalar in list */ {
  for (int _i = 0; _i < _N; _i++) /* foreach */
    b[] = lambda[]*s[];
  mg = poisson (s, b, alpha, lambda);
      }
    }
  }
  return mg;
}
 
/**
## References
 
~~~bib
@hal{popinet2020, hal-02365730}
 
@hal{devita2019, hal-02295398}
~~~
*/