src/spherical.h
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.
map()
{
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 \displaystyle cm = fm.y = \cos(\phi) \displaystyle fm.x = 1 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. \displaystyle \sum_{children}cm_{child} = 4cm_{parent} This is achieved by computing the exact area of a cell i.e. \displaystyle \Delta^2 \int^{\phi + d\phi/2}_{\phi - d \phi/2} \cos\phi d\phi = \Delta^2 (\sin(\phi + d\phi/2) - \sin(\phi - d\phi/2)) = \Delta^2 cm
#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
}