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. Δ 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 = Δ;
  Δ *= Radius*π/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 cm=fm.y=cos(ϕ) fm.x=1 with ϕ 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. childrencmchild=4cmparent This is achieved by computing the exact area of a cell i.e. Δ2ϕdϕ/2ϕ+dϕ/2cosϕdϕ=Δ2(sin(ϕ+dϕ/2)sin(ϕdϕ/2))=Δ2cm

#if TREE
static void refine_cm_lonlat (Point point, scalar cm)
{
  double φ = y*π/180., dphi = Δ/(2.*Radius);
  double sinphi = sin(φ);
  fine(cm,0,0) = fine(cm,1,0) = (sinphi - sin(φ - dphi))/dphi;
  fine(cm,0,1) = fine(cm,1,1) = (sin(φ + 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 φ = y*π/180., dphi = Δ/(2.*Radius);
  if (!is_refined(neighbor(0,-1)))
    fine(fm,0,0) = fine(fm,1,0) = cos(φ - dphi);
  if (!is_refined(neighbor(0,1)) && neighbor(0,1).neighbors)
    fine(fm,0,2) = fine(fm,1,2) = cos(φ + dphi);
  fine(fm,0,1) = fine(fm,1,1) = cos(φ);
}
#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 φ = y*π/180., dphi = Δ/(2.*Radius);
    cmv[] = (sin(φ + dphi) - sin(φ - 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*π/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

  boundary ({cm, fm});
}

Usage

Examples

Tests