PointCloud Voronoi Construction Options

These options live under pointcloud_voronoi.construction: and control how the Voronoi tessellation is built once particles are loaded. All keys are optional; CGAL remains the production default when it is compiled into THOR.

dataset_type: 'pointcloud_voronoi'

pointcloud_voronoi:
  hilbert_bits: 8
  loader: "gadget"
  gadget:
    path: ./snapdir_099
    fields: ["PositionX", "PositionY", "PositionZ", "Density"]

  construction:
    voronoi_backend: cgal               # cgal, voro, or sycl
    fallback_backend: cgal              # cgal, voro, or none; SYCL only
    voro_threads: 1                     # host threads; serial is MPI-safe default
    bounded_tesselation: false           # debug-only; see below
    ghost_thickness: -1.0                # PBC/CUBE_OLD only; negative = auto (see below)
    ghost_safety_factor: 3.0             # C in the auto ghost_thickness estimate
    min_face_area: -1.0                  # -1 = no face-area filter
    post_cleanup_sanity_checks: true     # validate cell topology after build

Default backend: CGAL

With voronoi_backend: cgal (the default in a CGAL-enabled build), the Voronoi mesh is built by CGAL's exact-predicates Delaunay triangulation. This path is the long-standing production backend and is the safe choice for any new dataset. A build configured with THOR_ENABLE_CGAL=OFF and THOR_ENABLE_VORO=ON defaults to Voro++ instead.

Known operational caveat: launching several independent THOR runs that all construct CGAL Voronoi tessellations at the same time can deadlock inside CGAL, even though the runs are separate process instances. Avoid concurrent mesh construction on the same node; build once with mesh_cache_mode: "auto" or "save" and reuse the cached mesh for repeated runs. See Troubleshooting.

Construction options on the CGAL path:

Voro++ backend (opt in)

voronoi_backend: voro selects the optional host-side Voro++ backend. Build with THOR_ENABLE_VORO=ON; bundled Voro++ is used by default, or select THOR_VORO_PROVIDER=SYSTEM. Voro++ supports CUBE_OLD and PERIODIC, but not BOUNDED. Set voro_threads to the number of CPU cores available per MPI rank; it defaults to 1.

Mesh caching

Voronoi construction is the expensive part of loading a point cloud, so the neighbour connectivity can be written to a small binary file and reused. This also sidesteps the concurrent-CGAL deadlock described in Troubleshooting.

pointcloud_voronoi:
  mesh_cache_mode: "auto"
  mesh_save_path: /path/to/voronoi_mesh_cache.bin

When no explicit path is given, the cache name encodes the point count, interpolator, Hilbert bits (only when hilbert_sort is on, since otherwise they do not affect the mesh), boundary type, numeric precision, and a digest of the construction parameters (e.g. mesh_cache_pointcloud_voronoi_1000000_8_pbc_fp64_i32_p3f1c0a27b5d4e918.bin) so runs with different settings don't clobber each other — a parameter sweep keeps one cache per setting instead of rebuilding on every switch.

A cache is only reused when the point count, interpolator, boundary type and that same construction digest all match. The digest covers every knob that changes the tessellation — hilbert_bits, hilbert_sort, min_face_area, the resolved voronoi_backend, and the effective ghost shell — and deliberately ignores knobs that don't (post_cleanup_sanity_checks, voro_threads, cgal_marker_*), which therefore never force a rebuild.

Knobs are compared at the value the construction path will actually use, not as written — a setting this configuration ignores never forces a rebuild:

• Backends are compared after resolution, so spelling out the value this build already defaults to is not a change. fallback_backend counts only under voronoi_backend: sycl, its only consumer.
• hilbert_bits counts only under hilbert_sort: true. It reaches the mesh through the Hilbert permutation, and with the sort off there is none.
• This rank's padded domain box counts always, whatever the sort is doing, because it is a construction input twice over. It is the rank's point cloud — under MPI each rank tessellates its own slab, and nothing else in the digest or the generated name tells two equal-sized ranks apart — and, with the sort on, it is also the Hilbert normalization box, so it fixes the row order the cached CSR is indexed against alongside the bit count. It moves with local_domain_pad (under PBC the padding is not clamped back to the global domain, so even a single rank sees it), with an explicit domain, and with the rank layout. The resolved box is compared rather than the keys behind it. For the sort the check is deliberately conservative: the digest is computed before the permutation, so it can only tell that the normalization changed, not whether the order did.
• min_face_area counts as the effective threshold, so every negative value is the same "disabled" setting. The face-area filter is a host-backend feature, so the key is ignored entirely under voronoi_backend: sycl with fallback_backend: none, where no host construction can run.
• ghost_thickness counts only where a ghost layer is actually built — a CGAL construction (voronoi_backend: cgal, or sycl with a CGAL fallback) on a CUBE or periodic domain. Voro++ tessellates periodic domains with native periodic containers, SYCL uses wall sentinels and native min-image wrapping, and bounded_tesselation has no ghost layer anywhere, so the key is ignored in all three. Where it does count, it is compared as the explicit value if one is set, as "auto" where auto sizing applies, and otherwise as the compiled default.
• ghost_safety_factor counts only when auto sizing is actually active — an explicit ghost_thickness, or a configuration that builds no ghosts, skips the estimate, and then the factor cannot affect the mesh.
• Whether the auto estimate will find radii counts too, wherever auto sizing is active — with no radius source the estimate warns and falls back to the compiled default shell, so "auto" is not a single shell. What is compared is that answer (together with the loader name), not the keys behind it: MaxDist reaching the loader through fields or the synthesized fields_required, the particle type, and, on gadget, smoothing_length suppressing the automatic request. So two gadget configs differing only in smoothing_length must not share a cache, while listing MaxDist explicitly on a gas run that would have requested it anyway is not a change. The debug-only gadget: {fields: all} is its own third answer rather than a yes or a no: it discovers whatever raw datasets the snapshot happens to carry, so whether it yields radii is not decidable from the config. It therefore never shares a cache with an equivalent explicit field list, in either direction. Whether gadget.aux_fields is set counts for the same reason: the sidecar is injected just before construction and contributes every dataset it holds, so one carrying MaxDist feeds the estimate radii the same config without it would not have. A warm field cache holding MaxDist is a radius source too, and the only one keyed on what the run actually did rather than on what the config says: a field-cache miss leaves MaxDist to the post-construction derivation, while a hit restores it before construction. So a cold first run and a warm second run of one unchanged config do not share a mesh cache, and a configured-but-cold cache does share with a config that has no field cache at all — in both cases matching what was actually built. Where MaxDist is requested outright, warming the cache is not a change and is not treated as one: the cold run materializes the recipe and writes those values, the warm run restores the same ones, so both reach construction with the same radii.
• VORONOI_PBC_GHOST is the one environment variable in the digest, folded in only when it is set and only under voronoi_backend: sycl on a periodic domain in a CGAL-enabled build — the one configuration that reads it. It swaps the native min-image candidate supply for the legacy ghost-padded one, so it builds a different neighbour graph; sharing a cache across the switch would make the diagnostic compare a mesh against itself (or, under auto, skip the diagnostic construction altogether). An unset variable — every normal run — hashes exactly as if the check were not there.

Positions are not hashed, so if you change particle positions but keep the same count you must delete the cache to force a rebuild. This also applies to the ghost shell in auto mode (ghost_thickness unset or negative): the thickness actually used is derived from the particle radii, so a different snapshot can resolve a different shell behind an identical digest. The same holds for an aux_fields sidecar: that it is configured is part of the identity, but which file it points at, and what that file contains, is not.

Caches written before the construction digest existed (payload version 1 or 2) are rejected with a warning and rebuilt — they record no construction parameters, so there is no way to tell what settings produced them.

The file format is self-describing; the thor.common.mesh_cache Python reader and thor-tools common mesh-cache info can inspect one.

Row order and the identity column

Each cache carries a per-row ids column that records which particle of the original, unfiltered snapshot a row belongs to. It is on by default; mesh_cache_ids: false turns it off (recorded as id_kind = 0, no ids[] array in the file).

• id_kind = 2 (PARTICLE_IDS): ids[i] is the ParticleIDs value of row i, taken from the snapshot. Written whenever the (gadget) snapshot has a ParticleIDs dataset.
• id_kind = 1 (SNAPSHOT_INDEX): ids[i] is row i's position in the concatenated snapshot files, counted before any filtering, so it indexes the full snapshot arrays directly. Fallback when there is no ParticleIDs dataset.

Both kinds survive particle selection: the identity values pass through the same filter chain as the particle data, so they always point at the original snapshot. Only the gadget loader fills the column; sphere (no snapshot to map back to) and tipsy (not implemented) write id_kind = 0.

from thor.common.mesh_cache import MeshCache
import numpy as np

mesh = MeshCache("voronoi_mesh.bin")
print(mesh.id_provenance)          # 'PARTICLE_IDS' or 'SNAPSHOT_INDEX'

# Map mesh rows -> snapshot indices (PARTICLE_IDS case):
snap_ids = ...                     # snapshot ParticleIDs, concatenated file order
order = np.argsort(snap_ids)
row_to_snap = order[np.searchsorted(snap_ids, mesh.ids, sorter=order)]

Caches from older THOR versions have no ids column. Their rows only map back to the snapshot if the producing run used hilbert_sort: false and no filtering, in which case row order equals snapshot order.

Producing a mesh only

To use THOR purely as a Voronoi mesher — build the tessellation, write it, and stop before any radiative-transfer work — combine a save mode with the early exit:

pointcloud_voronoi:
  mesh_cache_mode: "save"
  mesh_save_path: /path/to/voronoi_mesh_cache.bin
  exit_after_mesh_write: true

THOR shuts down cleanly (MPI finalized) as soon as the cache is on disk.

SYCL Voronoi (experimental)

voronoi_backend: sycl selects a SYCL-native implementation of Ray et al. 2018 Meshless Voronoi on the GPU: each cell clips against its K nearest neighbors with a security-radius termination test, growing K via tiered retries. Cells that exhaust the tiers fall through to the configured host fallback so the final neighbor list stays geometrically correct. CGAL fallback remains global, while Voro++ can construct just the requested residual rows.

This backend will offer a consistent, SYCL-first parallelization within THOR.

Status: experimental

• Validated against CGAL byte-identically on the unit-test suite (thor_tests "[voronoi][sycl][periodic]") and on bench_voronoi at N = 1 M and N = 3 M for both CUBE_OLD and PERIODIC boundary modes.
• End-to-end on TNG50-4 (18.5 M particles, PBC) completes on both OMP and a 24 GB CUDA GPU, with the final neighbor list matching the CGAL reference on all converged cells. Residual cells (a few hundred) are routed to the per-cell CGAL fallback automatically.
• Not yet recommended as the default: we have not tested this on production, high-resolution cosmological simulations with highly inhomogeneous point distributions yet.

For a production rollout, leave voronoi_backend: cgal per-config and opt in case-by-case after spot-checking the SYCL output against CGAL on the specific dataset.

Construction options on the SYCL path

The legacy use_sycl_voronoi and use_cgal_fallback boolean keys remain accepted for existing configurations, with a deprecation warning. When both forms are present, voronoi_backend and fallback_backend take precedence.

The other construction keys (bounded_tesselation, ghost_thickness, min_face_area, post_cleanup_sanity_checks) retain their existing meaning on the SYCL path. Ghost shells are prepared only when CGAL construction may be needed; native Voro++ uses its container's CUBE walls or periodic dimensions instead. post_cleanup_sanity_checks validates the final merged CSR regardless of which backend produced each cell.