Fix general asymmetry of x1f, dx1v Coordinates terms; ensure exact float64 symmetry of PPM stencils

[felker]

Description

@zhuzh1983 first reported that PPM was breaking symmetry in the Rayleigh-Taylor test in November 2017. In December 2017, I discovered that PLM also broke symmetry due to asymmetric round-off error in the Coordinates class, but the resulting differences were very small. A non-general fix to Coordinates() was designed and fixed PLM, but the PPM solution was still asymmetric.

The issue persisted after redesigning the Reconstruction functions in Feb-Mar 2018. A general fix to Coordinates() for uniform, Cartesian meshes was designed, and the PPM-specific issue was diagnosed last week. Detailed report forthcoming.

Changes to x1 calculations, for example, include:

1. Replaced the calculation of x1f in Coordinates() class constructor to call a new inlined function UniformMeshGeneratorX1() in mesh.hpp. This parameterizes the real cell faces in the calling MeshBlock by [-0.5, 0.5] instead of the [0,1] used by the DefaultMeshGeneratorX1() function. 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 ).
2. Made dx1v=dx1f= constant for all cells in cartesian.cpp when x1rat=1.0.
3. PPM stencils were reorganized to eliminate biasing to lower-indexed quantities. Left-to-right associativity in floating point arithmetic is a subtle issue that can introduce asymmetric round-off error. For example,

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);

always added/subtracted the lower indexed quantities first, and hence this operation's round-off error dominated the final round-off error. I replaced such 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) ;

Explicit comments (signed "KGF") were added to those changes, since they will be hard to maintain. The stencil biasing was never was in issue in PLM, since those stencils are only 2 cells wide.

This also replaced all long int and long long int usages in the code with more portable precise-width int64_t 64-bit integers. See #95.

Testing and validation

Tested on Rayleigh-Taylor problem, uniform Cartesian mesh of 100x200 cells and one MeshBlock. I created an explicit, double-precision unit in the last place (ULP) comparison routine based on the 32-bit analysis from https://randomascii.wordpress.com/2012/01/11/tricks-with-the-floating-point-format/. I added symmetry traps in:

• calculate_fluxes.cpp: check wl, wr states
• add_flux_divergence.cpp: check u_out before and after update
• Reconstruction routines intermediate quantities.

rt-ppm-fixed-diffs.pdf shows the fixed solution in this test.

To-do

• Confirm that coarse_flag check is not needed for current use of UniformMeshGeneratorX1(), and the continued use of DefaultMeshGeneratorX1() in mesh.cpp is compatible. @tomidakn
• Consider adding the ULP symmetry checks to the repository for special DEBUG_SYMMETRY mode.
• Need to test with multiple MeshBlock configurations
• Need to test with SMR/AMR
• Add documentation to Wiki about Athena++ symmetry guarantees
• Add instructions to Wiki's Programmer Guide about floating point best practices to preserve symmetry
• Share report detailing the Coordinates redesign

[felker]

After discussion with @jmstone , I repeated the coarse_flag logic for noffset from the nonuniform Coordinates() calculations in the uniform calculation