[yt-dev] round-off error in TotalQuantity

[Douglas Harvey Rudd]

Hi all,

I'm working on getting the ARTIO frontend working with a script used by the AGORA collaboration to measure density profiles.  This script uses the "TotalQuantity" quantity to sum all particle masses within a given sphere, then differences a nested set of spheres to measure the mean density in spherical shells.

I noticed a 2% difference at large radius from my own profiling routines, and tracked it down to the precision of the sum, which used the same precision as the input variable, float32 for ARTIO. Since there are a large number of particles, summing 1e5 x 10^-7 floats led to O(1%) error.

While this is clearly not an ideal way to measure density profiles, it's probably a good idea to help avoid such issues in future by changing TotalQuantity (and other places, like _TotalMass) to use a 64-bit accumulator, e.g.
        totals.append(data[field].sum(dtype=np.float64))

The alternative is to promote all float32 input variables to float64, which is wasteful of memory, or otherwise warn the user that they are likely to encounter loss of precision.

Douglas Rudd
Scientific Computing Consultant
Research Computing Center
drudd at uchicago.edu