Motivation
Recent tests by @bwvdnbro and myself reveal that the current mechanism in SKIRT sometimes fails to build a Voronoi tessellation for site lists imported from hydro simulations, especially for a large number of sites/cells (and correspondingly high site/cell density). One can then still resample to an octree grid to perform the radiative transfer simulation (see pull request #95) but that is not always an acceptable solution. Furthermore, the behavior sometimes differs between compilers.
Description
This pull request attempts to address this situation with three changes:
Include the latest version of the Voro++ library written by Chris H. Rycroft (Harvard University/Lawrence Berkeley Laboratory) which can now be found at https://github.com/chr1shr/voro .
Rather than giving up when a problem arises, the Voronoi construction is allowed to continue. At the end, if there are invalid cells (empty volume, inconsistent neighbor lists), the sites corresponding to these cells are discarded and the Voronoi construction is retried (with a maximum of 5 iterations).
When compiled with the Intel compiler, the -fp-model precise flag is specified to block the compiler from using optimizations that cause non-standard floating point operations (e.g. reordering instructions resulting in different rounding). This brings the Intel compiler in line with standard practice.
Tests
We tested the updated code with Voronoi grids up to 16 million sites. The retry process needs 3 iterations in the worst case and the number of discarded sites is very small compared to the total number. Although these adjustments represent a risk of changing the physics of the input model, it seems that the effect is small in practice. Also, the behavior is now identical on all tested compilers (Intel, GCC, Clang).
The results of the functional tests are compatible with but not identical to previous results:
The order in which the cells and the neighbors within each cell are generated may differ. This does not change the SKIRT instrument output but may affect some of the probe output (for example, the grid plotting files list the edges to be plotted in a different order).
Some contrived test cases that have sites lined up on a regular Cartesian grid produce different output. Previously, a warning was issued about neighbors not being mutual; now these cells/sites are discarded.
Context
A distinguishing feature of the Voro++ library is that it carries out cell-based calculations, computing the Voronoi cell for each site individually, rather than computing the Voronoi tessellation as a global network of vertices and edges. It is therefore particularly well-suited for applications that require cell-based properties such as the cell volume, the centroid position or the number of faces or vertices. Equally important in the context of SKIRT, it is easy to distribute the work over parallel execution threads because, after setting up a common search data structure holding all sites, the calculations for the various cells are mutually independent.
This approach also has an important drawback. Because cells are handled independently of each other, floating pointing rounding errors sometimes cause inconsistencies where the results for one cell do not properly match the corresponding results for a neighboring cell. This most often happens when the generating sites are very close to each other and/or form certain hard-to-calculate geometries. These problems manifest themselves either as empty cells or as asymmetries in the neighbor lists. To handle these situations, as described above, the updated version of the SKIRT VoronoiMeshSnapshot class now removes the sites corresponding to these invalid cells from the input list.
Motivation Recent tests by @bwvdnbro and myself reveal that the current mechanism in SKIRT sometimes fails to build a Voronoi tessellation for site lists imported from hydro simulations, especially for a large number of sites/cells (and correspondingly high site/cell density). One can then still resample to an octree grid to perform the radiative transfer simulation (see pull request #95) but that is not always an acceptable solution. Furthermore, the behavior sometimes differs between compilers.
Description This pull request attempts to address this situation with three changes:
-fp-model preciseflag is specified to block the compiler from using optimizations that cause non-standard floating point operations (e.g. reordering instructions resulting in different rounding). This brings the Intel compiler in line with standard practice.Tests We tested the updated code with Voronoi grids up to 16 million sites. The retry process needs 3 iterations in the worst case and the number of discarded sites is very small compared to the total number. Although these adjustments represent a risk of changing the physics of the input model, it seems that the effect is small in practice. Also, the behavior is now identical on all tested compilers (Intel, GCC, Clang).
The results of the functional tests are compatible with but not identical to previous results:
Context A distinguishing feature of the Voro++ library is that it carries out cell-based calculations, computing the Voronoi cell for each site individually, rather than computing the Voronoi tessellation as a global network of vertices and edges. It is therefore particularly well-suited for applications that require cell-based properties such as the cell volume, the centroid position or the number of faces or vertices. Equally important in the context of SKIRT, it is easy to distribute the work over parallel execution threads because, after setting up a common search data structure holding all sites, the calculations for the various cells are mutually independent.
This approach also has an important drawback. Because cells are handled independently of each other, floating pointing rounding errors sometimes cause inconsistencies where the results for one cell do not properly match the corresponding results for a neighboring cell. This most often happens when the generating sites are very close to each other and/or form certain hard-to-calculate geometries. These problems manifest themselves either as empty cells or as asymmetries in the neighbor lists. To handle these situations, as described above, the updated version of the SKIRT VoronoiMeshSnapshot class now removes the sites corresponding to these invalid cells from the input list.