src/layered/coriolis.h
Coriolis/friction terms for the multilayer solver
This approximates \displaystyle \partial_t\mathbf{u} = \mathbf{B}\mathbf{u} + \mathbf{a} with \displaystyle \mathbf{B} = \left( \begin{array}{cc} - K_0 & F_0\\ - F_0 & - K_0 \end{array} \right), and K_0 and F_0 the linear friction and Coriolis parameters respectively. The time-implicit discretisation of these terms can be written \displaystyle \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 This then gives \displaystyle \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 The local 2 \times 2 linear system is easily inverted analytically. The final value is obtained by substracting the acceleration i.e. \displaystyle \mathbf{u}^{\star} = \mathbf{u}^{n + 1} - \Delta t\mathbf{a}^n
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
event acceleration (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;
foreach_dimension() {
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;
}
foreach_dimension()
u.x[] = (m0.y*r.x - m1.x*r.y)/det - a.x;
}
}
#undef alpha_H
coord geostrophic_velocity (Point point)
{
coord ug;
static const coord a = {-1., 1.};
foreach_dimension()
ug.y = a.y*G*(hf.x[]*gmetric(0)*(eta[] - eta[-1]) +
hf.x[1]*gmetric(1)*(eta[1] - eta[]))
/(Delta*F0()*(hf.x[] + hf.x[1] + dry));
return ug;
}