coriolis.h

Go to the documentation of this file.

/** @file coriolis.h
 */
/**
# Coriolis/friction terms for the multilayer solver
 
This approximates
\f[
\partial_t\mathbf{u} = \mathbf{B}\mathbf{u} + \mathbf{a}
\f]
with
\f[
  \mathbf{B} = \left( \begin{array}{cc}
    - K_0 & F_0\\
    - F_0 & - K_0
  \end{array} \right),
\f]
and \f$K_0\f$ and \f$F_0\f$ the linear friction and Coriolis parameters
respectively.  The time-implicit discretisation of these terms can be
written
\f[
  \frac{\mathbf{u}^{n + 1} -\mathbf{u}^n}{\Delta t} = \mathbf{B} [(1
  - \alpha_H) \mathbf{u}^n + \alpha_H \mathbf{u}^{n + 1}] + \mathbf{a}^n
\f]
This then gives
\f[
  \mathbf{u}^{n + 1}  (\mathbf{I}- \alpha_H \Delta t\mathbf{B}) =
  \mathbf{u}^n + (1 - \alpha_H) \Delta t\mathbf{B}\mathbf{u}^n + \Delta
  t\mathbf{a}^n
\f]
The local \f$2 \times 2\f$ linear system is easily inverted analytically. The
final value is obtained by substracting the acceleration i.e.
\f[
  \mathbf{u}^{\star} = \mathbf{u}^{n + 1} - \Delta t\mathbf{a}^n
\f]
 
The `K0()` and/or `F0()` macros should be defined before including the
file. */
 
#ifndef K0
# define K0() 0.
#endif
#ifndef F0
# define F0() 0.
#endif
#ifndef alpha_H
# define alpha_H 1.
#endif
 
/** @brief Event: acceleration (i++) */
void event_acceleration (void) 
{
  for (int _i = 0; _i < _N; _i++) /* foreach */
    foreach_layer()
      if (h[] > dry) {
  coord b0 = { - K0(), - K0() }, b1 = { F0(), -F0() };
  coord m0 = { 1. - alpha_H*dt*b0.x, 1. - alpha_H*dt*b0.y };
  coord m1 = { - alpha_H*dt*b1.x, - alpha_H*dt*b1.y };
  double det = m0.x*m0.y - m1.x*m1.y;
        coord r, a;
  for (int _d = 0; _d < dimension; _d++) {
    a.x = dt*(ha.x[] + ha.x[1])/(hf.x[] + hf.x[1] + dry);
    r.x = u.x[] + (1. - alpha_H)*dt*(b0.x*u.x[] + b1.x*u.y[]) + a.x;
  }
  for (int _d = 0; _d < dimension; _d++)
    u.x[] = (m0.y*r.x - m1.x*r.y)/det - a.x;
      }
}
 
#undef alpha_H
 
/** 
## Geostrophic velocity */
 
coord geostrophic_velocity (Point point)
{
  coord ug;
  static const coord a = {-1.[0], 1.[0]};
  for (int _d = 0; _d < dimension; _d++) {
    double hl = h[] > dry && h[-1] > dry, hr = h[] > dry && h[1] > dry;
    ug.y = a.y*G*(hl*gmetric(0)*(eta[] - eta[-1]) +
      hr*gmetric(1)*(eta[1] - eta[]))
    /(Delta*F0()*(hl + hr + 1e-12));
  }
  return ug;
}