sandbox/Antoonvh/force_geo.h
Pressure-gradient force in a rotating frame of reference
This file implements both a Coriolis force (F_{c} = \mathbf{f} \times \mathbf{v}) and and a large-scale pressure-gradient force \nabla_h p. Assuming that \mathbf{f} points in the vertial direction (zenith) with f_y= f_cor
>0. Then we can associate a so-called geostrophic velocity vector with the aforementionted focings.
\displaystyle \mathbf{U_g} = \frac{e_z}{\rho f}\times\nabla_h p
and \mathbf{U_g} = \{U_GEO
, V_GEO
, 0\}, which are two macros that should be #define
d before including this file, which should happen later than the inclusion of the centered slover.
note that the meteorological convention is rotated: \{u,v,w\} = `{u.z, u.x, u.y}.
#ifndef U_GEO
#define U_GEO 0.
#endif
#ifndef V_GEO
#define V_GEO 0.
#endif
double f_cor = 1e-4; //A non-zero constant
Implementation
in case of a two dimensional setup, we add a perpendicular velocity field for the above equations to make sense. The positive direction is e_z = e_x \times e_y (towards you in typical visualizations).
#if dimension == 2
#include "diffusion.h"
#include "tracer.h"
scalar uz[];
event tracer_diffusion (i++)
diffusion (uz, dt, mu);
#endif
Inspired by other pages, we possibly need to allocate the acceleration field in an early event.
event defaults (i = 0) {
if (is_constant(a.x)) {
a = new face vector;
foreach_face()
a.x[] = 0;
boundary ((scalar*){a});
#if dimension == 2
tracers = list_append (tracers, uz);
#endif
}
}
The “geostrophic velocity” gives a convinient basis to implement the pressure gradient force.
event acceleration (i++) {
face vector av = a;
boundary (all);
#if dimension == 3 // u.z should be available
foreach_face(x)
av.x[] += -f_cor*((u.z[] + u.z[-1])/2. - U_GEO);
foreach_face(z)
av.z[] += f_cor*((u.x[] + u.x[0,0,-1])/2. - V_GEO);
#else // we use uz
foreach_face(x)
av.x[] += -f_cor*((uz[] + uz[-1])/2. - U_GEO);
foreach()
uz[] += dt*f_cor*(u.x[] - V_GEO);
#endif
}