Arbitrary Coordinates

For most general relativity applications, the provided coordinates should prove sufficient. In case they do not, Athena++ provides support for arbitrary coordinate systems. Here we detail how to make use of this advanced feature.

Enrolling new coordinates

Arbitrary coordinates require writing a new problem generator file. The only additional requirement is that the user must define and enroll a function to calculate the metric and its derivatives.

For example, suppose we were to recreate the built-in spherical Kerr-Schild coordinates. In the problem generator, we would declare the function

void SphericalKerrSchild(Real x1, Real x2, Real x3, ParameterInput *pin,
    AthenaArray<Real> &g, AthenaArray<Real> &g_inv, AthenaArray<Real> &dg_dx1,
    AthenaArray<Real> &dg_dx2, AthenaArray<Real> &dg_dx3);

This function would need to be enrolled, much like other user-defined functions:

void Mesh::InitUserMeshData(ParameterInput *pin)
{
  ...
  EnrollUserMetric(SphericalKerrSchild);
  ...
}

The function will take the spatial coordinates of a point, as well as a pointer to the parameter list, as inputs. It should populate the 5 given 1D arrays with the appropriate values. These arrays are the covariant metric components, the contravariant components, and the derivatives of the covariant components with respect to the three spatial coordinates in turn. Note that even arbitrary coordinates are limited to be stationary in Athena++, removing time dependence from consideration.

The indices should use the enum values I00, ..., I33, of which there should be NMETRIC = 10. These are meant for symmetric 4×4 matrices like the metric. The first number is less than or equal to the second number.

For spherical Kerr-Schild, the coordinates can be defined as follows:

void SphericalKerrSchild(Real x1, Real x2, Real x3, ParameterInput *pin,
    AthenaArray<Real> &g, AthenaArray<Real> &g_inv, AthenaArray<Real> &dg_dx1,
    AthenaArray<Real> &dg_dx2, AthenaArray<Real> &dg_dx3)
{
  // Extract inputs
  Real m = pin->GetReal("coord", "m");
  Real a = pin->GetReal("coord", "a");
  Real r = x1;
  Real theta = x2;
  Real phi = x3;

  // Calculate intermediate quantities
  Real s = std::sin(theta);
  Real c = std::cos(theta);
  Real delta = SQR(r) - 2.0*m*r + SQR(a);
  Real sigma = SQR(r) + SQR(a) * SQR(c);
  Real xi = SQR(r) - SQR(a) * SQR(c);

  // Set covariant components
  g(I00) = -(1.0 - 2.0*m*r/sigma);
  g(I01) = 2.0*m*r/sigma;
  g(I02) = 0.0;
  g(I03) = -2.0*m*a*r*SQR(s)/sigma;
  g(I11) = 1.0 + 2.0*m*r/sigma;
  g(I12) = 0.0;
  g(I13) = -(1.0 + 2.0*m*r/sigma) * a*SQR(s);
  g(I22) = sigma;
  g(I23) = 0.0;
  g(I33) = (SQR(r) + SQR(a) + 2.0*m*SQR(a)*r*SQR(s)/sigma) * SQR(s);

  // Set contravariant components
  g_inv(I00) = -(1.0 + 2.0*m*r/sigma);
  g_inv(I01) = 2.0*m*r/sigma;
  g_inv(I02) = 0.0;
  g_inv(I03) = 0.0;
  g_inv(I11) = delta/sigma;
  g_inv(I12) = 0.0;
  g_inv(I13) = a/sigma;
  g_inv(I22) = 1.0/sigma;
  g_inv(I23) = 0.0;
  g_inv(I33) = 1.0/(sigma*SQR(s));

  // Set r-derivatives of covariant components
  dg_dx1(I00) = -2.0*m*xi/SQR(sigma);
  dg_dx1(I01) = -2.0*m*xi/SQR(sigma);
  dg_dx1(I02) = 0.0;
  dg_dx1(I03) = 2.0*m*a*xi*SQR(s)/SQR(sigma);
  dg_dx1(I11) = -2.0*m*xi/SQR(sigma);
  dg_dx1(I12) = 0.0;
  dg_dx1(I13) = 2.0*m*a*xi*SQR(s)/SQR(sigma);
  dg_dx1(I22) = 2.0*r;
  dg_dx1(I23) = 0.0;
  dg_dx1(I33) = (2.0*r - 2.0*m*SQR(a)*xi*SQR(s)/SQR(sigma)) * SQR(s);

  // Set theta-derivatives of covariant components
  dg_dx2(I00) = 4.0*m*SQR(a)*r*s*c/SQR(sigma);
  dg_dx2(I01) = 4.0*m*SQR(a)*r*s*c/SQR(sigma);
  dg_dx2(I02) = 0.0;
  dg_dx2(I03) = -4.0*m*a*r*s*c/sigma * (1.0 + SQR(a)*SQR(s)/sigma);
  dg_dx2(I11) = 4.0*m*SQR(a)*r*s*c/SQR(sigma);
  dg_dx2(I12) = 0.0;
  dg_dx2(I13) = -4.0*m*a*r*s*c/sigma * (1.0 + SQR(a)*SQR(s)/sigma) - 2.0*a*s*c;
  dg_dx2(I22) = -2.0*SQR(a)*s*c;
  dg_dx2(I23) = 0.0;
  dg_dx2(I33) = 4.0*m*SQR(a)*r*s*SQR(s)*c/sigma * (2.0 + SQR(a)*SQR(s)/sigma)
      + 2.0 * (SQR(r) + SQR(a)) * s*c;

  // Set phi-derivatives of covariant components
  for (int n = 0; n < NMETRIC; ++n) {
    dg_dx3(n) = 0.0;
  }
  return;
}

Performance considerations

The built-in coordinate systems are internally optimized by consideration of the formulas for quantities such as the metric coefficients, the discretized areas and volumes, and the connection coefficients. Symmetries of the metric generally allow these quantities to be broken into parts that depend on 1 or at most 2 coordinates. Thus these pieces can be stored in small arrays and assembled when needed.

For arbitrary coordinates this is no longer feasible. Instead, all the needed quantities are computed and stored in 3D arrays. This large memory usage can adversely affect performance, especially on machines with limited cache size.

Since these computations are only done during initialization, the speed of the user-defined function is generally not a performance issue.