Coordinate Systems and Meshes

Select a Coordinate System

Currently the following coordinate systems are available:

• Cartesian (default; $x$, $y$, $z$)
• Cylindrical coordinates ($R$, $\phi$, $z$)
• Spherical-polar coordinates ($r$, $\theta$, $\phi$)
• Various general-relativistic coordinates:
  • Minkowski
  • Schwarzschild
  • Kerr-Schild
  • User-defined metric

To select a coordinate system, configure the code with the --coord [coord] option.

> python ./configure.py --prob sphtorus --coord spherical_polar

The names of coordinate systems are corresponding to the file names in src/coordinates, and are listed in the help message (use the -h option).

2D Polar Coordinates

Note, if a complete circular mesh in 2D is desired, the appropriate choice for the coordinate system is --coord=cylindrical, not spherical_polar. The x2 polar angle coordinate is strictly limited to the $[0,\ \pi]$ range in both 2D and 3D spherical-polar coordinates. The x2 azimuthal angle coordinate can take on any value in $[0,\ 2\pi]$ for the cylindrical coordinate system in 2D and 3D.

In Athena++, 2D spherical-like coordinates are used to simulate a very thin "orange slice" of the domain.

Mesh Spacing

By default, the cell size is simply set evenly: $\Delta x = L_x / N_x$. It is sometimes convenient to use geometric spacing, i.e. the cell size multiplies (or divides) by a constant factor $f_\mathrm{rat}$ as one moves to an adjacent cell. This is especially useful in spherical-polar coordinates so that it can maintain the cell aspect ratio. To enable this, set x1rat, x2rat, or x3rat in the <mesh> block in the input file.

    <mesh>
    ...
    x1rat = 1.01

While Athena++ uses formulae consistent with nonuniform mesh spacing, too large a mesh ratio should not be used; a non-fatal warning is raised if the value is not between 0.9 and 1.1 everywhere. The reconstruction algorithm reduces to first-order if the grid spacing changes too much from cell-to-cell, and moreover, the nonuniform truncation error can produce anomalous reflections of waves.

Note, the geometric spacing is equivalent to logarithmic spacing (uniform $\Delta x'$ in logspace given by the coordinate transformation $x' = \log(x)$) only for a particular choice of x1rat given x1min and x1max: $f_\mathrm{rat} = (x_\mathrm{max} / x_\mathrm{min})^{1/N_x}$.

Note that this ratio should be adjusted when resolution is manually changed. For example, when the number of cells in the direction is doubled, the ratio should be adjusted to square root of the original one to maintain consistency. This adjustment is done automatically whenever mesh refinement is invoked.

User-Defined Mesh Generator

More generally, any kind of mesh spacing can be used through the user-defined mesh generator feature. Similar to the user-defined boundary conditions, define the MeshGenerator function first:

    Real MyMeshSpacingX2(Real x, RegionSize rs);

And then enroll it in Mesh::InitUserMeshData:

    EnrollUserMeshGenerator(X2DIR, MyMeshSpacingX2);

Finally, the input variable mesh/x2rat=-1.0 must be set in The Input File or on the command line in order to run a simulation with a user-defined MeshGenerator function in that direction.

The MeshGenerator function in x1, for example, must be a monotonic, 1:1 mapping from $[0,\ 1]$ to [rs.x1min, rs.x1max], which returns a coordinate location between rs.x1min and rs.x1max given the logical location of the cell. For example, the first active cell will have its left and right $x$-faces logically located at $0$ and $1/N_x$, which in coordinates will be rs.x1min and the MeshGenerator's returned value given $1/N_x$.

For example, let us consider a mesh generator which produces high resolution near the midplane in the $\theta$ (x2) direction in spherical-polar coordinates.

    Real CompressedX2(Real x, RegionSize rs)
    {
      Real t=2.0*x-1.0;
      Real w=0.25*(t*(t*t+1.0))+0.5;
      return w*rs.x2max+(1.0-w) * rs.x2min;
    }

This function gives mesh spacing proportional to $(3 \psi^2 + 1)$, where $\psi$ is the angle from the midplane, and the mesh spacing near the pole is about 4 times larger than that near the midplane. The code will give a warning if the resulting mesh contains a steep change in the mesh size.

Creating a new Coordinates class

For almost every problem, one of the existing Newtonian/special-relativistic coordinates (Cartesian, cylindrical, spherical-polar) or general-relativistic coordinates (Minkowski, Schwarzschild, Kerr-Schild) will suffice. One can also easily specify an arbitrary metric in GR as long as a Cartesian topology (e.g., no polar axis singularity) is acceptable (see Arbitrary Coordinates). However, in rare cases, one might wish to extend the code to another coordinate system.

The following outlines the steps one can take to augment Athena++ with a new coordinates class, named MyCoord here for illustration:

• Add the definition of the class, derived from the Coordinates class, to src/coordinates/coordinates.hpp:
  • class MyCoord : public Coordinates {friend class HydroSourceTerms;...};
  • This can be as simple as what Cartesian does.
  • Any of the virtual functions declared in the Coordinates base class should be repeated here, without virtual, if and only if they are to be redefined from their Cartesian defaults.
  • Any newly required scratch arrays for storing internal data should be declared in the protected section of the Coordinates definition.
• Create a new file named something appropriate: src/coordinates/my_coord.cpp:
  • At a minimum, this must define a constructor that does most of what the Cartesian constructor in src/coordinates/cartesian.cpp does.
    • It must initialize the x1v (cell-centered position), dx1f (face-to-face change in coordinate across cell), and dx1v (center-to-center change in coordinate between adjacent cells) arrays given x1f (face coordinates). These three arrays are already allocated. The same must be done for x2 and x3.
    • It must similarly initialize the h* and dh* coefficients if there is a chance these coordinates will be used with diffusion operators or orbital advection.
    • It must similarly initialize the x*s* arrays to be compatible with MHD-AMR.
  • This file must also define the functions overridden in the class definition. See some of the non-Cartesian files for nontrivial examples. To be safe, any lengths, areas, or volumes should be exact. There are only some cases (e.g., cell widths in GR, which are only used for calculating safe timesteps) in which approximations will not harm the accuracy of the code.
  • Any existing or new scratch arrays can be allocated and initialized in the constructor, to be used for any computation. None of these arrays are allocated by default. The general principle the existing code employs is that expensive separable calculations should store precomputed values (e.g. of log(x1) or sin(x2)) in 1D arrays, and then return sums/products of these values when needed. This avoids repeated expensive calculations (by storing values), and it avoids using too much memory (by storing 3D results). The manner in which scratch arrays are used is a question strictly of performance, not accuracy.
• The file name without extension (my_coord) should be added as an option in 3 places in configure.py: the list of all valid choices, and twice in lists of choices that are not compatible with GR. Search the file for coord to find these.

Symmetry preservation in Athena++

Hydrodynamics and MHD test problems with physical symmetries are discriminating tests of computational astrophysics codes. In particular, linear instabilities may amplify any small floating-point differences (as small as 1 unit in the last place (ULP)), resulting in visible asymmetries in late-time solutions. The issues typically become more pronounced when using schemes with high-order accuracy; we have observed solutions with symmetry errors many orders of magnitude larger when using time/xorder=3 PPM than when using time/xorder=2 PLM reconstruction, for example.

While the underlying integration algorithm of Athena++ is dimensionally unsplit and is well-suited to preserving symmetries in theory, there are several (related) challenges when attempting to preserve symmetry with floating-point operations in practice:

• Both absolute and relative round-off errors are asymmetric about any value != 0.0
• Floating-point arithmetic is non-associative
• Compilers may enable value-unsafe optimizations of floating-point operations

Although symmetry-preservation is a desirable solver behavior, it may not be worth the developer effort (or even feasible) to maintain certain complex symmetry relationships for all possible test problems and physics. The design choices can be difficult, especially if the required changes would degrade the clarity, generality, or performance of the overall solver (see Style Guide).

Athena++ has been designed to preserve certain test problem symmetries exactly to double precision accuracy. This capability is automatic; no additional options when Configuring or Running the Code are required to preserve reflective symmetry in a coordinate direction (e.g. transformation in x1 about x2-x3 plane), as long as the following constraints are satisfied:

• Uniform mesh using Cartesian coordinates
  • Confirmed working with AMR/SMR without MPI
• Symmetric initial condition
• Hydrodynamics problems only
  • Constrained transport (CT) magnetic field updates currently break symmetry due to non-associativity of operations in field/calculate_corner_e.cpp
• Plane of symmetry located at x1=0.0; this restriction is sufficient for any number of MeshBlocks and nx1 (even or odd valued ---> reflection either about cell faces x1f=0.0 or volume-centers x1v=0.0, respectively)
• Known bug: Configured without MPI. Solvers using unrefined grids with MPI might preserve symmetry; AMR with MPI likely won't.
• Compiled with a value-safe floating-point arithmetic mode

This has been tested up to 64-bit floating-point precision for a variety of resolutions, additional physics, Riemann solvers, and reconstruction methods.

Miscellaneous notes for users:

• When performing symmetry tests and using HDF5 output, the code should be configured with -h5double for maximum precision.
• The Analysis Tools provided in the repository includes the athdf() class in athena_read.py for HDF5 data. This reader was designed to analyze both AMR/SMR and unrefined results. When reading AMR/SMR results into Python, the script interpolates quantities to a common grid and may introduce asymmetries. However, when provided with unrefined, uniform grid HDF5 data, the reader is guaranteed to preserve any symmetries present in the raw floating-point data. Alternatively, the h5py module can be used to directly load the .athdf files.

Compiler-dependent floating-point math modes

The Intel compiler is more aggressive with math optimizations than either the Clang or GNU Compiler Collection (GCC). Both clang++ and g++ disable -ffast-math at all optimization levels, by default. In contrast, the default icc math mode is -fp-model fast=1, which is value-unsafe. Studies of the above symmetric problems in Athena++ have shown that compiling with -fp-model strict is required to preserve symmetry with the Intel compiler. However, this flag prevents the use of OpenMP 4.0 SIMD vectorization extensions; thus, the performance of the code drops dramatically. The less restrictive -fp-model precise option is value-safe, but is insufficient to prevent the creation of asymmetries in solutions to the 2D Rayleigh-Taylor problem, for example.

According to Intel fp-model documentation, the only differences between strict and precise is that the former 1) also enables floating-point exception semantics and 2) disables floating-point contractions, including fused multiply-add (FMA) instructions. Both modes ensure reproducibility by prohibiting any reordering/change of associativity of floating-point operations:

disables optimizations that can change the result of floating-point calculations, which is required for strict ANSI conformance.

Future work will ensure symmetry-preservation with the Intel compiler while enabling SIMD vectorization. The ability of OpenMP threading to reorder floating-point operations of inside SIMD loops can be toggled via the environment variable KMP_DETERMINISTIC_REDUCTION; setting it to 0 may be sufficient to maintain problem symmetry with -fp-model precise:

Enables (1) or disables (0) the use of a specific ordering of the reduction operations for implementing the reduction clause for an OpenMP* parallel region. This has the effect that, for a given number of threads, in a given parallel region, for a given data set and reduction operation, a floating point reduction done for an OpenMP reduction clause will have a consistent floating point result from run to run, since round-off errors will be identical.

Tips for Athena++ developers

When adding new physics or changing existing Athena++ source code, it is important to keep the following concerns for preserving symmetry in mind:

• User-specified initial conditions, source terms, refinement conditions, etc. that depend on a function of position should refer to the cell-interface positions, if possible. For example, f(x1f) is preferred to f(x1v).
• Any operation that is expressed in separate per-direction calculations, such as functions OperatorX1(), OperatorX2(), and OperatorX3(), should use identical sequences of floating-point operations. For example, earlier versions of the PiecewiseLinearX*() functions directly divided a quantity for reconstructions along x1, whereas an inverse was precomputed and stored for reconstructions along x2:

// PLMx1()
Real& dx_im2 = pco->dx1v(i-2);
dql = (q(nin,k,j,i-1) - q(nin,k,j,i-2))/dx_im2;
// PLMx2()
Real dx2jm2i = 1.0/pco->dx2v(j-2);
dql = (q(nin,k,j-1,i) - q(nin,k,j-2,i))*dx2jm2i;

The latter sequence incurs an intermediate rounding error that does not occur in the first formulation. This difference created visible artifacts in the Liska-Wendroff implosion problem at high-resolutions.

• The Regression Testing suite should be periodically extended to check symmetry maintenance of new problems. Unlike issues with error convergence, symmetry violation bugs can be introduced by very minor changes and may require significant developer effort to debug.
• Developers should be sure that stencils/expressions that reference mesh quantities indexed across a range of cells in x1, x2, or x3 greater than +/- 1 are written without biased floating-point associativity. For example, in reconstruct/ppm.cpp, the original calculation of approximations to the second-derivative of input AthenaArray q along the x1 direction (indexed by i in variable names) was expressed as:

d2qc_im1(i) = q_im2(n,i) - 2.0*q_im1(n,i) + q    (n,i);
d2qc    (i) = q_im1(n,i) - 2.0*q    (n,i) + q_ip1(n,i); //(CD eq 85a) (no 1/2)
d2qc_ip1(i) = q    (n,i) - 2.0*q_ip1(n,i) + q_ip2(n,i);

Even though these stencils are "centered" (d2qc expression equally weights im1 and ip1 quantities) in infinite precision arithmetic, they are not centered in finite-precision arithmetic. The solver always adds/subtracts the lower indexed quantities first due to implicit left-to-right associativity in floating-point arithmetic. Hence, the round-off error from subtracting the first two terms (lower i) "dominates" the round-off error in the final result and introduces asymmetries in the solution. Simply replacing the above lines with:

d2qc_im1(i) = q_im2(n,i) + q    (n,i) - 2.0*q_im1(n,i) ;
d2qc    (i) = q_im1(n,i) + q_ip1(n,i) - 2.0*q(n,i); //(CD eq 85a) (no 1/2)
d2qc_ip1(i) = q    (n,i) + q_ip2(n,i) - 2.0*q_ip1(n,i) ;

guarantees that the ip1 and im1 quantities can be swapped in the middle expression and return the identical floating-point value, for example.

Implementation in Coordinates class

The class is currently designed such that the calculation of cell face positions in the Coordinates() class constructor uses a special inlined function UniformMeshGeneratorX1() in mesh.hpp. This function parameterizes the real cell faces in the calling MeshBlock by [-0.5, 0.5] instead of the [0,1] range used by the DefaultMeshGeneratorX1() function for user-defined mesh generators. The same function is used to compute the boundary positions of each MeshBlock. It was designed to preserve exact floating point symmetry of x1f if x1min=-x1max and work for nx1 even (central face at x1f=0) or odd (central cell has x1v=0). Furthermore, dx1f=dx1v are guaranteed to be identical and constant when x1rat=1.0 for a Cartesian mesh.