spherical.h
📄 View in API Reference (Doxygen) · View on basilisk.fr
Test cases (1): lonlat
Examples (4): global-tides, held-suarez, tohoku, tsunami
Spherical coordinates
The default radius of the sphere is set to one.
double Radius = 1.;Rather than the classical mathematical convention, we use the geographic convention: *x* is the longitude within $[-180:180]$ degrees and *y* is the latitude within $[-90:90]$ degrees. $\Delta$ is the characteristic cell length expressed in the same units as *Radius*.
macro VARIABLES (Point point = point, int _ig = ig, int _jg = jg, int _kg = kg)
{
VARIABLES (point, _ig, _jg, _kg);
x = x < -180. ? x + 360. : x > 180. ? x - 360. : x;
Delta_x = Delta_y = Delta;
Delta *= Radius*pi/180.;
}On trees we need to define consistent refinement functions. The default restriction functions (averages) are already consistent (i.e. they preserve volume and length integrals).
Mathematically we have
with $\phi$ the latitude in radians. This is the definition we use for the length scale factors *fm*.
For the volume scale factor *cm*, more care needs to be taken to guarantee discrete volume conservation i.e.
This is achieved by computing the exact area of a cell i.e.
#if TREE
static void refine_cm_lonlat (Point point, scalar cm)
{
double phi = y*pi/180., dphi = Delta/(2.*Radius);
double sinphi = sin(phi);
fine(cm,0,0) = fine(cm,1,0) = (sinphi - sin(phi - dphi))/dphi;
fine(cm,0,1) = fine(cm,1,1) = (sin(phi + dphi) - sinphi)/dphi;
}
static void refine_face_x_lonlat (Point point, scalar fm)
{
if (!is_refined(neighbor(-1)))
fine(fm,0,0) = fine(fm,0,1) = 1.;
if (!is_refined(neighbor(1)) && neighbor(1).neighbors)
fine(fm,2,0) = fine(fm,2,1) = 1.;
fine(fm,1,0) = fine(fm,1,1) = 1.;
}
static void refine_face_y_lonlat (Point point, scalar fm)
{
double phi = y*pi/180., dphi = Delta/(2.*Radius);
if (!is_refined(neighbor(0,-1)))
fine(fm,0,0) = fine(fm,1,0) = cos(phi - dphi);
if (!is_refined(neighbor(0,1)) && neighbor(0,1).neighbors)
fine(fm,0,2) = fine(fm,1,2) = cos(phi + dphi);
fine(fm,0,1) = fine(fm,1,1) = cos(phi);
}
#endif
event metric (i = 0) {We initialise the scale factors, taking care to first allocate the fields if they are still constant.
if (is_constant(cm)) {
scalar * l = list_copy (all);
cm = new scalar;
free (all);
all = list_concat ({cm}, l);
free (l);
}
scalar cmv = cm;
foreach() {
double phi = y*pi/180., dphi = Delta/(2.*Radius);
cmv[] = (sin(phi + dphi) - sin(phi - dphi))/(2.*dphi);
}
if (is_constant(fm.x)) {
scalar * l = list_copy (all);
fm = new face vector;
free (all);
all = list_concat ((scalar *){fm}, l);
free (l);
}
face vector fmv = fm;
foreach_face(x)
fmv.x[] = 1.;
foreach_face(y)
fmv.y[] = cos(y*pi/180.);We set our refinement/prolongation functions on trees.
#if TREE
cm.refine = cm.prolongation = refine_cm_lonlat;
fm.x.prolongation = refine_face_x_lonlat;
fm.y.prolongation = refine_face_y_lonlat;
#endif
}