Ray direction sampling via rand(): overlap_check uses one constant direction, and all sites are confined to the positive octant

[cyb3ralbert]

Describe the Bug

Ray direction sampling for point-containment checks goes through libc rand() in a way that has three separate problems. The same idiom appears in four places — two in DAGMC, two in MOAB's GeomQueryTool:

Location	Notes
src/overlap_check/overlap.cpp:94	one direction for the entire run
src/geant4/DagSolid.cc:193-199	per Inside() call
MOAB src/GeomQueryTool.cpp:955 (point_in_volume)	fallback when uvw is null/zero
MOAB src/GeomQueryTool.cpp:1215 (find_volume)	same fallback

The DAGMC and MOAB copies are identical down to the const double magnitude line. MOAB line numbers are from 5.5.1 (DAGMC's default pull); the code is unchanged on MOAB master/develop today, at :944 and :1209.

The three problems, in the order I would rank them:

1. overlap_check fires one constant direction on every run. It links only dagmc, and there is no srand() anywhere in that binary, so rand() starts from the implementation default seed. The direction at :94 is therefore a compile-time constant in practice, identical for every geometry and every run.

2. All four sites sample only the positive octant. rand() returns [0, RAND_MAX], so all three components are non-negative and the normalized vector is confined to one eighth of the sphere.

3. DagSolid::Inside() duplicates MOAB's fallback. It computes a direction and passes it to point_in_volume, but MOAB does exactly the same thing when no direction is given, with the same code. Passing NULL would give a statistically identical result. The comment at DagSolid.cc:207 ("if uvw is not given, this function generate uvw ran") suggests this was known.

To be clear about severity: for a watertight closed volume, any non-degenerate direction gives the correct containment answer, so the octant restriction is a robustness issue rather than a wrong-answer issue. The point of choosing a direction at random is presumably to avoid a ray lying along an edge or passing exactly through a vertex, where the crossing count degenerates. Problem 1 is what I would call an actual defect: with a constant direction that protection does not exist, so a geometry badly aligned with this one vector produces the same wrong answer on every run, and re-running cannot help. The direction is used for both senses (overlap.cpp:39 negates it), so it is effectively one fixed line for the whole check.

To Reproduce

The constant direction can be shown without building DAGMC, since the expression is self-contained:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(void) {
  /* no srand() — same as overlap_check */
  double u = rand(), v = rand(), w = rand();
  double m = sqrt(u*u + v*v + w*w);
  printf("(%.6f, %.6f, %.6f)\n", u/m, v/m, w/m);
  return 0;
}

On glibc this prints (0.691871, 0.324763, 0.644860) every time — the first three values of the default stream, {1804289383, 846930886, 1681692777}, normalized.

Two wrinkles worth noting. In overlap.cpp:94 the three rand() calls are arguments to one constructor, and their relative order is unspecified in C++17 (indeterminately sequenced — well defined, just unordered), so which component receives which value is a compiler decision: gcc evaluates right-to-left and yields (0.644860, 0.324763, 0.691871), the mirror of the above. The values are also libc dependent, so the constant differs between glibc and macOS/musl. The result is a direction that is fixed for any given build but varies across toolchains for reasons unrelated to the geometry.

For the octant bias, sampling the same expression 1e6 times gives a mean direction of (0.5158, 0.5157, 0.5153), |mean| = 0.8931. Within the octant the distribution is not uniform either, since drawing uniformly in a cube and normalizing pushes density toward the cube corners.

Expected Behavior

A direction drawn uniformly over the whole sphere, from a stream that is explicitly seeded so runs are reproducible and the seed is under the caller's control.

Screenshots or Code Snippets

The shared idiom, as it appears in DagSolid.cc:193-199:

double u = rand();
double v = rand();
double w = rand();

const double magnitude = sqrt(u * u + v * v + w * w);
u /= magnitude;
v /= magnitude;
w /= magnitude;

The standard uniform-on-sphere method is already used elsewhere in this repository — src/dagmc/tools/ray_fire_test.cpp:48-56 (RNDVEC) samples z uniformly in [-1, 1] and the azimuth uniformly over 2π, and it seeds explicitly (randseed = 12345) and re-seeds when it needs to repeat a sequence:

inline void RNDVEC(CartVect& uvw, double& az) {
  double theta = az * denom * rand();
  double u = 2 * denom * rand() - 1;
  uvw[0] = sqrt(1 - u * u) * cos(theta);
  uvw[1] = sqrt(1 - u * u) * sin(theta);
  uvw[2] = u;
}

So this is an inconsistency within DAGMC rather than a missing technique, and no new dependency is needed: <random> is available, and DagSolid.cc:61 already includes Randomize.hh, so G4UniformRand() is in scope there.

Please complete the following information regarding your system:

• OS: Debian 12 (bookworm), gcc 12.2.0, glibc 2.36
• MOAB Version: source inspected at 5.5.1, and master/develop
• Physics codes versions installed with: n/a — this is a source-level report, not a runtime failure

Additional Context

A note on seeding elsewhere: src/tally/KDEMeshTally.cpp:561 also uses rand(), but for a scalar position along a sub-track, and that sampling is correct. Its issue is only the seed — srand(time(NULL)) at :68 unless the user passes seed= (documented in doc/usersguide/tally.rst:97), so sub-track points differ run to run by default. Because libc's stream is process-global, that seed also shifts directions drawn by geometry code in the same process, and vice versa — an unintended coupling between a tally and containment queries. src/make_watertight/tests/ also uses rand(), which looks fine as test code.

I am deliberately not claiming a multithreading bug: DAGMC ships a plain G4RunManager (src/geant4/app/exampleN01.cc:40, and the generator at src/geant4/generate_geant4:1391), and there is no G4MULTITHREADED anywhere under src/geant4/. The narrower current concern is that because DAGMC draws from libc's stream rather than Geant4's, this randomness is invisible to Geant4's own RNG management (/random/setSeeds, SetRandomNumberStore — the latter appears commented out at generate_geant4:1232), so a Geant4 run meant to be reproducible through those controls is not reproducible in the DagSolid part.

Questions before I write any patch. Every fix here moves numerical output, so I would rather agree scope first than open PRs you have to argue with.

1. This cannot be fully fixed inside DAGMC. point_in_volume called without a direction — as in src/dagmc/tests/dagmc_pointinvol_test.cpp:49 and dagmc_simple_test.cpp:154, and by downstream users — lands in MOAB's copy. Do you want this raised with MOAB as well, and would you rather open that conversation given you already track MOAB master in linux_upstream_test_moab.yml? I am happy to file it there if that is useful.
2. How much change in numerical output is acceptable? Are there regression tests with reference values that would need regenerating, and is now a reasonable time for that?
3. Is the MCNP5/6 KDEMeshTally path still maintained, or effectively frozen? It only builds under BUILD_MCNP5/BUILD_MCNP6/BUILD_CI_TESTS, and CI never sets the MCNP flags. If it is frozen I would rather leave it alone than churn it.
4. Does an overlap_check fix collide with Improvement on Overlap Check Tool #797? That issue proposes reworking the tool to fire rays along triangle edges instead of sampling points. If that rewrite is near, patching the current direction sampling may be wasted effort.

Suggested split, smallest blast radius first:

• overlap_check — <random> with an explicit fixed seed plus uniform-on-sphere sampling. This makes today's accidental determinism explicit and independent of libc and compiler, while actually covering the sphere; OverlappingVolumesTest already exercises this path. Optionally expose --seed, so a user who hits a degenerate direction has a way out, which today they do not.
• DagSolid — either drop the block and let MOAB's fallback handle it, or switch to G4UniformRand() plus uniform-on-sphere so the sampling participates in Geant4's RNG management. I would lean to the latter, but it depends on your answer to Q1.
• KDEMeshTally — only if the path is alive; minimally, drop the time(NULL) default. Feeding it from the transport code's own stream would be more correct but changes the tally API, so that is a separate conversation.

Happy to do all of this, one PR at a time, each with a doc/CHANGELOG.rst entry.