Fix thread-safety bug in `geom.c` by oskooi · Pull Request #75 · NanoComp/libctl

oskooi · 2026-04-26T15:53:38Z

Root cause: Two thread-safety bugs in utils/geom.c, exposed by Meep's recently added OMP parallization of set_chi1inv (NanoComp/meep#3166).

Fix 1 (primary — caused the crash of test_mode_coeffs.py): Remove the geom_fix_object_ptr(&o) call from geom_get_bounding_box(). For Prism objects, this call triggered reinit_prism(), which frees and re-mallocs all internal prism arrays. When multiple OMP threads call this concurrently on the same Prism, it causes double-free/use-after-free, corrupting heap metadata and producing the free(): unaligned chunk detected in tcache 2 error.

Fix 2 (secondary — would cause incorrect results): Change intersect_line_segment_with_prism() to use a stack-allocated variable-length array (VLA) (double slist[num_vertices + 2]) instead of the shared prsm->workspace.items buffer. This prevents data races during concurrent intersection calculations.

Verification: all 6 tests in Meep's test_mode_coeffs.py now pass with OMP_NUM_THREADS=1, 2, and 4.

@Luochenghuang

stevengj · 2026-04-26T21:20:52Z

Fix 1 (primary — caused the crash of test_mode_coeffs.py): Remove the geom_fix_object_ptr(&o) call from geom_get_bounding_box(). For Prism objects, this call triggered reinit_prism(), which frees and re-mallocs all internal prism arrays. When multiple OMP threads call this concurrently on the same Prism, it causes double-free/use-after-free, corrupting heap metadata and producing the free(): unaligned chunk detected in tcache 2 error.

This changes the semantics. But it's probably a good idea anyway since fix_object is expensive and should really only be done once at the beginning of all geometry calculations, and in fact we are already doing this.. (In Meep, it is called in the geom_epsilon constructor, and in MPB and pympb it is called in init_epsilon.

stevengj · 2026-04-26T21:21:32Z

Change intersect_line_segment_with_prism() to use a stack-allocated variable-length array (VLA) (double slist[num_vertices + 2]) instead of the shared prsm->workspace.items buffer.

This is potentially dangerous for a large prism because stack space is limited. In particular, it could crash for prisms more than about 100,000 elements.

stevengj · 2026-04-26T21:36:53Z

For the prism list, in practice the slist normally doesn't need to be equal to the number of vertices or anything close, since it is very unlikely you will have so many intersections. What we could do to make it both thread-safe and memory-safe would be something like:

Stack allocate a moderate number of intersections, e.g. double slist_stack[1024]; double *slist = slist_stack;, that will easily fit in the stack.
In intersect_line_with_prism, add an extra argument int slist_len for the length of slist (passing 1024)
Change the line from slist[num_intersections++] = s; to something like if (num_intersections++ < slist_len) slist[num_intersections-1] = s. Then, at the end of that loop, add a check if (num_intersections > slist_len) return num_intersections; to return immediately
In the caller intersect_line_segment_with_prism, check if (num_intersections > 1024) { slist = (double *) malloc(num_intersections * sizeof(double)); num_intersections = intersect_line_with_prism(prsm, pc, dc, slist, num_intersections); } to re-do the intersect_line_with_prism calculation with a heap-allocated array, as a slow fallback path. At the end of the function, also add if (num_intersections > 1024) free(slist); to deallocated it in this case.
Get rid of prsm->workspace, which is no longer used.

stevengj · 2026-05-01T19:39:29Z

I implemented my proposal and pushed a commit to this PR.

oskooi · 2026-05-08T06:30:35Z

The latest changes (c617512) are producing failures for the Meep unit test test_mode_coeffs.py when OMP_NUM_THREADS is set to 2 or 4:

OMP_NUM_THREADS=2:

malloc(): unaligned tcache chunk detected

OMP_NUM_THREADS=4:

degenerate plane in project_point_into_plane

intersect_line_segment_with_mesh added (b - last_s) when a hit fell at or past b, then the post-loop cleanup added the same length again because inside/last_s were not updated in the break branch. Drop the in-branch addition and let the post-loop check handle the inside-at-b case. Add test_cube_segments_vs_block: 1000 random unit-direction line segments through a unit cube mesh vs make_block, comparing interior lengths within 1e-4. This test surfaced the bug; max diff after the fix is ~4e-16. Also note that remove_duplicate_intersections should be unified with the prism slist dedup after NanoComp#75.

oskooi · 2026-05-20T21:43:20Z

This PR is ready to be merged. The errors I reported in my previous comment should be ignored as I discovered that the local build of Meep I was using for testing was linking to a different version of libctlgeom which did not include the changes in this PR.

I verified that the previously failing tests are now passing with this PR.

stevengj · 2026-05-21T10:38:21Z

We now need a PR enabling --enable-omp during CI.

…#81) * Fix init_prism centroid + remove racy per-query auto-fix from query paths Two related bugs surfaced by running the --enable-openmp test suite against test-prism, both pre-existing: 1) init_prism assigned prsm->centroid BEFORE the shift-to-o->center logic, then the shift only updated a local `centroid` variable. After construction, any prism whose center was modified ended up with prsm->centroid stale relative to the (shifted) vertices. 2) point_in_objectp, normal_to_object, and point_in_periodic_objectp each called geom_fix_object_ptr(&o) on every query. For prisms that path calls reinit_prism, which free()s + malloc()s shared per-prism arrays (vertices_p, vertices_top_p, ...). Under OpenMP the test was crashing with "malloc(): unaligned tcache chunk" / "free(): double free". This matches what PR #75 already did for geom_get_bounding_box -- callers are now responsible for calling geom_fix_object_ptr after mutating o.center. test-prism: after the random-shift branch, call geom_fix_object_ptr on the block and prism once, single-threaded, before the parallel test loops run. Verified at OMP_NUM_THREADS = 1, 2, 4, 8, 16, 32, 64, 128, 166: all clean (0 point failures, 0 segment failures, no heap corruption). make check passes. * Restore auto-fix in query wrappers; document them as not thread-safe Per review feedback, point_in_objectp / normal_to_object / point_in_periodic_objectp are restored to their original behavior (call geom_fix_object_ptr internally) rather than being silently turned into thin aliases for the *_fixed_* variants. Their doc comments now state they are NOT thread-safe and direct concurrent callers at the *_fixed_* entry points. test-prism's parallel loops switched from point_in_objectp / normal_to_object to point_in_fixed_objectp / normal_to_fixed_object, which are pure reads and safe to call from multiple threads against an object that was fixed up once beforehand. Combined with the init_prism centroid fix in the previous commit, make check passes cleanly across OMP_NUM_THREADS = 1, 2, 4, 8, 16, 64, 166 (21 trials each: 0 point-inclusion failures, 0 segment failures, no heap-corruption signatures).

Mirror the test-prism.c thread-safety verification (PR NanoComp#75 / commit 926527c): parallelize the two random-comparison loops in test-mesh.c that already use the pure _fixed_* query variants: - test_cube_vs_block (10000 random points, point_in_fixed_pobjectp) - test_cube_segments_vs_block (1000 random segments, intersect_line_segment_with_object — for MESH this dispatches directly to intersect_line_segment_with_mesh, which uses a stack-local mesh_hit_list with per-call heap fallback) Both objects (mesh + block) are constructed on the main thread before the parallel region, so init_mesh / BVH build happens once, single-threaded. After that, every query path through point_in_mesh, normal_to_mesh, find_closest_face, mesh_ray_all_intersections, ray_triangle_intersect, ray_bvh_node_intersect, and BVH traversal is a pure read against the mesh struct — no shared mutable state, no module-level statics (the only static module data is the const mesh_ray_dirs table). rand() is replaced with rand_r() seeded per iteration so results are deterministic across thread counts; this lets the assertions (0 mismatches, max_diff at machine precision) hold uniformly in serial and parallel runs. Verified with --enable-openmp at OMP_NUM_THREADS = 1, 2, 4, 8, 16, 64, 166: identical results (0 mismatches, max_diff = 4.44e-16) at every thread count, no heap corruption. 20 repeated trials at 64 threads: all clean.

@stevengj

* Add MESH geometric object type for arbitrary triangulated 3D surfaces Implement a new MESH primitive alongside SPHERE, CYLINDER, BLOCK, and PRISM. A mesh is defined by a vertex array and triangle index array, with a BVH (bounding volume hierarchy) for O(log N) queries. Core operations implemented: - point_in_mesh: ray casting with parity counting and deduplication - normal_to_mesh: BVH-accelerated closest-face search - get_mesh_bounding_box: from BVH root node - get_mesh_volume: divergence theorem (signed tetrahedron volumes) - intersect_line_segment_with_mesh: for subpixel smoothing integration - init_mesh: validates geometry, computes normals, fixes winding, builds BVH BVH uses SAH (surface area heuristic) binned construction with flat array layout. Ray-triangle intersection uses Moller-Trumbore algorithm. Closest-point queries use Eberly's algorithm with BVH pruning. Also adds make_mesh/make_mesh_with_center constructors, mesh copy/equal/destroy in geom-ctl-io.c, and 11 unit tests covering point-in, volume, bbox, normals, line-segment intersection, and comparison against the Block primitive. * Fix point_in_mesh robustness and reinit_mesh performance Two fixes: 1. point_in_mesh: use 3-ray majority vote instead of a single ray. A single ray can pass through a shared mesh edge, causing the deduplication to either miss or double-count a crossing and flip the inside/outside parity. Three rays in irrational directions make it extremely unlikely that more than one hits a degenerate case at the same query point. 2. reinit_mesh: skip BVH rebuild if already initialized. Unlike blocks/cylinders, mesh uses absolute coordinates and doesn't depend on the lattice basis. geom_fix_object_ptr is called hundreds of times during meep's init_sim via geometry copies; without the guard, each call rebuilt the BVH (~150ms for 9K triangles), causing ~40s of wasted time. * Lazy majority vote and mesh closure validation - point_in_mesh: cast 1 ray first, only fall back to 3-ray majority vote if deduplication removed near-duplicate t-values (indicating an edge/vertex hit). 3x faster for the common case (no edge hits). - init_mesh: validate mesh closure by building an edge-to-face-count hash table. Every edge must be shared by exactly 2 faces. Sets is_closed=false and prints a warning for open/non-manifold meshes. point_in_mesh returns false for all points on open meshes. - Add 2 new tests: open mesh detection, closed mesh detection (13 tests total, all passing). * Add closure validation edge case tests Test boundary edges (1 face per edge), non-manifold edges (4 faces sharing an edge), and isolated vertices (vertex unreferenced by any face). 16 tests total, all passing. * Eliminate redundant BVH build and point_in_mesh fallback 1. make_mesh_with_center: shift vertices before calling init_mesh so the BVH is only built once (was building twice when center was specified). 2. intersect_line_segment_with_mesh: determine inside/outside at the segment start from the full intersection list parity (count crossings before parameter 'a') instead of calling point_in_mesh as a fallback. Eliminates millions of extra ray casts during subpixel smoothing. ~20% faster set_epsilon on the Utah teapot (28.6s → 22.9s). * Minor cleanups: remove dead code, optimize BVH queries - Remove unused bvh_node_box helper function - find_closest_face: push nearer BVH child last so it's popped first, improving pruning and reducing nodes visited for normal queries - Extract mesh_triangle_vertices helper to avoid redundant vertex fetches. Add mesh_triangle_bbox_centroid to compute both in one pass during BVH binning (was fetching 3 vertices twice per face) * Add mesh API documentation to ctlgeom.h * Include generated ctlgeom-types.h and geom-ctl-io.c These files are normally auto-generated from geom.scm via gen-ctl-io (requires Guile). Including them allows building libctlgeom without Guile, which is the common case for meep's Python-only builds. * Address PR feedback: relative lengthscale tolerances and mesh input documentation * Fix 3 bugs found in code review: per-component winding, dynamic ray-hit buffer, vertex mutation rebuild * Drop vertex_hash mutation detection in reinit_mesh A non-NULL bvh is sufficient to skip rebuild: there is no documented way for callers to mutate mesh vertices post-construction, and none of our internal functions do so. Hash-based detection was also unsound on collisions. Removes the mesh_vertex_hash function, the vertex_hash field, and the now-irrelevant test_vertex_mutation case. * Fix double-counted interior length in mesh segment intersection intersect_line_segment_with_mesh added (b - last_s) when a hit fell at or past b, then the post-loop cleanup added the same length again because inside/last_s were not updated in the break branch. Drop the in-branch addition and let the post-loop check handle the inside-at-b case. Add test_cube_segments_vs_block: 1000 random unit-direction line segments through a unit cube mesh vs make_block, comparing interior lengths within 1e-4. This test surfaced the bug; max diff after the fix is ~4e-16. Also note that remove_duplicate_intersections should be unified with the prism slist dedup after #75. * Verify mesh thread-safety via parallel test loops Mirror the test-prism.c thread-safety verification (PR #75 / commit 926527c): parallelize the two random-comparison loops in test-mesh.c that already use the pure _fixed_* query variants: - test_cube_vs_block (10000 random points, point_in_fixed_pobjectp) - test_cube_segments_vs_block (1000 random segments, intersect_line_segment_with_object — for MESH this dispatches directly to intersect_line_segment_with_mesh, which uses a stack-local mesh_hit_list with per-call heap fallback) Both objects (mesh + block) are constructed on the main thread before the parallel region, so init_mesh / BVH build happens once, single-threaded. After that, every query path through point_in_mesh, normal_to_mesh, find_closest_face, mesh_ray_all_intersections, ray_triangle_intersect, ray_bvh_node_intersect, and BVH traversal is a pure read against the mesh struct — no shared mutable state, no module-level statics (the only static module data is the const mesh_ray_dirs table). rand() is replaced with rand_r() seeded per iteration so results are deterministic across thread counts; this lets the assertions (0 mismatches, max_diff at machine precision) hold uniformly in serial and parallel runs. Verified with --enable-openmp at OMP_NUM_THREADS = 1, 2, 4, 8, 16, 64, 166: identical results (0 mismatches, max_diff = 4.44e-16) at every thread count, no heap corruption. 20 repeated trials at 64 threads: all clean. * Apply suggestion from @stevengj --------- Co-authored-by: Steven G. Johnson <stevenj@mit.edu>

Drops the macos-latest matrix entry per maintainer preference (less useful: project is primarily Linux-targeted, macOS runners cost more minutes, and the existing Travis macOS coverage didn't exercise anything Linux didn't). Adds --enable-openmp to the configure step and matrices the run on OMP_NUM_THREADS = 1 and 8. Both values share the same build; only the runtime thread count differs. This catches regressions in the parallel test paths (test-prism's #pragma omp parallel for loops in test_point_inclusion, test_normal_to_object, and test_line_segment_intersection) that a serial-only run wouldn't flag — same matrix philosophy used to validate the prism thread-safety work in PR NanoComp#75 / NanoComp#81.

* Migrate CI from Travis to GitHub Actions Replaces the Travis configuration with an equivalent GitHub Actions workflow at .github/workflows/ci.yml. Same matrix (ubuntu-latest + macos-latest), same build pipeline (autogen.sh, make, make check, make install), same on-failure log dump for utils/*.log. Two minor modernizations: - Switches the Ubuntu guile package from guile-2.0-dev to guile-3.0-dev. guile-2.0 was dropped from Ubuntu after 20.04; libctl already builds and runs against guile 3.0. - macOS dependency install adds autoconf/automake/libtool explicitly rather than relying on the runner image to provide them. Workflow triggers on push to master, pull_request to master, and manual workflow_dispatch. fail-fast disabled so a failure on one OS doesn't cancel the other. Also removes the Travis CI badge from README.md and replaces it with a GitHub Actions badge pointing at the new workflow. * CI: drop macOS, add OMP_NUM_THREADS matrix to exercise parallel tests Drops the macos-latest matrix entry per maintainer preference (less useful: project is primarily Linux-targeted, macOS runners cost more minutes, and the existing Travis macOS coverage didn't exercise anything Linux didn't). Adds --enable-openmp to the configure step and matrices the run on OMP_NUM_THREADS = 1 and 8. Both values share the same build; only the runtime thread count differs. This catches regressions in the parallel test paths (test-prism's #pragma omp parallel for loops in test_point_inclusion, test_normal_to_object, and test_line_segment_intersection) that a serial-only run wouldn't flag — same matrix philosophy used to validate the prism thread-safety work in PR #75 / #81.

PR NanoComp#74 committed both ctlgeom-types.h and geom-ctl-io.c, but the versions it shipped were generated from a pre-PR-NanoComp#75 geom.scm that still had the prism workspace field. After NanoComp#75 removed workspace from geom.scm and from the C code in geom.c, the committed copies became stale. Without-Guile builds silently kept compiling because the two stale files were consistent with each other (prism_copy/destroy/equal operated on a workspace field that nothing else touched). With-Guile builds papered over it by regenerating ctlgeom-types.h to match the post-NanoComp#75 geom.scm, which silently differed from the committed copy. The previous commit (treat geom-ctl-io.c as checked-in) exposed the inconsistency by removing the regeneration safety net: with WITH_GUILE=true, ctlgeom-types.h gets regenerated fresh from geom.scm (no workspace, no number_list typedef), while the committed geom-ctl-io.c still tried to access prism->workspace, producing 'has no member named workspace' compile errors in prism_copy / prism_equal / prism_destroy. This commit aligns both committed files with current geom.scm: - ctlgeom-types.h: replaced with the gen-ctl-io output. Drops the unused number_list typedef and the workspace field from prism struct. Same mesh struct as before (vertices, face_indices, is_closed, internal). - geom-ctl-io.c: removed the stale workspace handling from prism_copy / prism_equal / prism_destroy. Kept the PIMPL hand-tweaks (extern mesh_init_internal / mesh_internal_free + their calls in mesh_copy / mesh_destroy) since the file is no longer auto-regenerated. Verified locally in both build modes: - WITHOUT Guile (WITH_GUILE=false, no regen): make + make check + OMP=1/8 clean. - WITH Guile (WITH_GUILE=true, ctlgeom-types.h regenerates and matches committed): make + make install + examples/ build (which exercises the cp utils/geom.c -> examples/geom.c MPB-style copy path) + make check + OMP=1/8 clean.

fix thread-safety bug in geom.c

0438715

oskooi added the bug label Apr 26, 2026

safer stack allocation

c617512

stevengj mentioned this pull request May 1, 2026

Mesh geometry import #74

Merged

2 tasks

Luochenghuang mentioned this pull request May 16, 2026

Add --enable-openmp to exercise ctlgeom thread-safety in test suite #79

Merged

stevengj merged commit 6de5034 into NanoComp:master May 21, 2026

oskooi deleted the thread_safety_fix branch May 21, 2026 19:29

Luochenghuang mentioned this pull request May 22, 2026

Fix init_prism centroid; make test-prism's parallel loops thread-safe #81

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Fix thread-safety bug in `geom.c` - #75

Fix thread-safety bug in `geom.c`#75
stevengj merged 2 commits into
NanoComp:masterfrom
oskooi:thread_safety_fix

oskooi commented Apr 26, 2026 •

edited

Loading

Uh oh!

stevengj commented Apr 26, 2026

Uh oh!

stevengj commented Apr 26, 2026 •

edited

Loading

Uh oh!

stevengj commented Apr 26, 2026 •

edited

Loading

Uh oh!

stevengj commented May 1, 2026

Uh oh!

oskooi commented May 8, 2026

Uh oh!

oskooi commented May 20, 2026

Uh oh!

stevengj commented May 21, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

oskooi commented Apr 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

stevengj commented Apr 26, 2026

Uh oh!

stevengj commented Apr 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

stevengj commented Apr 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

stevengj commented May 1, 2026

Uh oh!

oskooi commented May 8, 2026

Uh oh!

oskooi commented May 20, 2026

Uh oh!

stevengj commented May 21, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

oskooi commented Apr 26, 2026 •

edited

Loading

stevengj commented Apr 26, 2026 •

edited

Loading

stevengj commented Apr 26, 2026 •

edited

Loading