From 8535376efccaf09a375ed1e556a6fa174c6d28ea Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 10:29:23 -0500 Subject: [PATCH 01/35] Reapply "Prohibit ratracing from outside" This reverts commit b4ff3bbb01eacc53ca8326d152f1f5177df73454. --- src/geocel/vg/VecgeomTrackView.hh | 39 ++++++---------------------- src/geocel/vg/detail/BVHNavigator.hh | 13 +++++++--- 2 files changed, 17 insertions(+), 35 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 10a30e3c0c..19120024a7 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -425,32 +425,12 @@ CELER_FUNCTION Real3 VecgeomTrackView::normal() const /*! * Find the distance to the next geometric boundary. * - * This function is allowed to be called from the exterior for ray tracing. + * \todo To support ray tracing from unknown distances from the world, we + * could add another special function \c find_first_step . */ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step() { - if (this->is_outside()) - { - // SPECIAL CASE: find distance to interior from outside world volume - auto const& world_pv = this->world(); - next_step_ = world_pv.DistanceToIn(detail::to_vector(pos_), - detail::to_vector(dir_), - vecgeom::kInfLength); - - next_step_ = max(next_step_, this->extra_push()); - vgnext_.Clear(); - Propagation result; - result.distance = next_step_; - result.boundary = next_step_ < vecgeom::kInfLength; - - if (result.boundary) - { - vgnext_.Push(&world_pv); - vgnext_.SetBoundaryState(true); - } - - return result; - } + CELER_EXPECT(!this->is_outside()); return this->find_next_step(vecgeom::kInfLength); } @@ -570,22 +550,19 @@ CELER_FUNCTION void VecgeomTrackView::move_to_boundary() */ CELER_FUNCTION void VecgeomTrackView::cross_boundary() { + CELER_EXPECT(!this->is_outside()); CELER_EXPECT(this->is_on_boundary()); CELER_EXPECT(this->is_next_boundary()); // Relocate to next tracking volume (maybe across multiple boundaries) if (vgnext_.Top() != nullptr) { - if (!CELERITAS_VECGEOM_SURFACE || !vgstate_.IsOutside()) - { - // In surf model, relocation does not work from [OUTSIDE] - Navigator::RelocateToNextVolume(detail::to_vector(this->pos_), - detail::to_vector(this->dir_), + Navigator::RelocateToNextVolume(detail::to_vector(this->pos_), + detail::to_vector(this->dir_), #if CELERITAS_VECGEOM_SURFACE - *next_surf_, + *next_surf_, #endif - vgnext_); - } + vgnext_); } vgstate_ = vgnext_; diff --git a/src/geocel/vg/detail/BVHNavigator.hh b/src/geocel/vg/detail/BVHNavigator.hh index f393896db2..a62658db06 100644 --- a/src/geocel/vg/detail/BVHNavigator.hh +++ b/src/geocel/vg/detail/BVHNavigator.hh @@ -19,7 +19,9 @@ #include #include +#include "corecel/Macros.hh" #include "corecel/math/Algorithms.hh" +#include "geocel/vg/VecgeomTypes.hh" #ifdef VECGEOM_ENABLE_CUDA # include @@ -31,9 +33,6 @@ # include "VgNavStateWrapper.hh" #endif -#include "corecel/Macros.hh" -#include "geocel/vg/VecgeomTypes.hh" - namespace celeritas { namespace detail @@ -69,6 +68,7 @@ class BVHNavigator } } + CELER_ASSERT(vol); path.Push(vol); VgReal3 currentpoint(point); @@ -89,8 +89,12 @@ class BVHNavigator break; } - currentpoint = daughterlocalpoint; + CELER_ASSERT(vol); path.Push(vol); + + // Transform to daughter + currentpoint = daughterlocalpoint; + // Only exclude the placed volume once since we could enter it // again via a different volume history. exclude = nullptr; @@ -100,6 +104,7 @@ class BVHNavigator CELER_FUNCTION static void RelocatePoint(VgReal3 const& localpoint, NavState& path) { + CELER_EXPECT(!path.IsOutside()); VgPlacedVol const* currentmother = path.Top(); VgReal3 transformed = localpoint; do From 70b0ffe5ca86314c4ca28f316ec7b2ec959944fc Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 09:49:07 -0500 Subject: [PATCH 02/35] Add safety with lmit --- src/geocel/vg/VecgeomTrackView.hh | 9 ++------- src/geocel/vg/detail/BVHNavigator.hh | 8 +++----- src/geocel/vg/detail/SurfNavigator.hh | 7 ++++--- 3 files changed, 9 insertions(+), 15 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 19120024a7..0677664726 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -511,13 +511,8 @@ CELER_FUNCTION real_type VecgeomTrackView::find_safety(real_type max_radius) CELER_EXPECT(!this->is_on_boundary()); CELER_EXPECT(max_radius > 0); - real_type safety = Navigator::ComputeSafety(detail::to_vector(this->pos()), - vgstate_ -#if VECGEOM_VERSION >= 0x200000 - , - max_radius -#endif - ); + real_type safety = Navigator::ComputeSafety( + detail::to_vector(this->pos()), vgstate_, max_radius); safety = min(safety, max_radius); // Since the reported "safety" is negative if we've moved slightly beyond diff --git a/src/geocel/vg/detail/BVHNavigator.hh b/src/geocel/vg/detail/BVHNavigator.hh index a62658db06..340d7ca218 100644 --- a/src/geocel/vg/detail/BVHNavigator.hh +++ b/src/geocel/vg/detail/BVHNavigator.hh @@ -236,11 +236,9 @@ class BVHNavigator public: // Computes the isotropic safety from the globalpoint. - CELER_FUNCTION static double - ComputeSafety(VgReal3 const& globalpoint, - NavState const& state, - vg_real_type safety - = std::numeric_limits::infinity()) + CELER_FUNCTION static double ComputeSafety(VgReal3 const& globalpoint, + NavState const& state, + vg_real_type safety) { VgPlacedVol const* pvol = state.Top(); vecgeom::Transformation3D m; diff --git a/src/geocel/vg/detail/SurfNavigator.hh b/src/geocel/vg/detail/SurfNavigator.hh index d8a3a07f3f..aa02f31ab6 100644 --- a/src/geocel/vg/detail/SurfNavigator.hh +++ b/src/geocel/vg/detail/SurfNavigator.hh @@ -70,12 +70,13 @@ class SurfNavigator /// @param globalpoint Point in global coordinates /// @param state Path where to compute safety /// @return Isotropic safe distance - CELER_FUNCTION static vg_real_type - ComputeSafety(VgReal3 const& globalpoint, NavState const& state) + CELER_FUNCTION static vg_real_type ComputeSafety(VgReal3 const& globalpoint, + NavState const& state, + Precision limit) { auto safety = vgbrep::protonav::BVHSurfNavigator::ComputeSafety( - globalpoint, state); + globalpoint, state, limit); return safety; } From eacba71a45ef59fe3fd0ed460a7eda311e3a8f6e Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 11:05:20 -0500 Subject: [PATCH 03/35] Remove world accessor --- src/geocel/vg/VecgeomTrackView.hh | 20 ++++---------------- 1 file changed, 4 insertions(+), 16 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 0677664726..cbf4e28722 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -193,9 +193,6 @@ class VecgeomTrackView // Whether the next distance-to-boundary is to a surface inline CELER_FUNCTION bool is_next_boundary() const; - // Get a reference to the world volume - inline CELER_FUNCTION VgPlacedVol const& world() const; - // Get a reference to the current volume instance inline CELER_FUNCTION VgPlacedVol const& physical_volume() const; @@ -285,9 +282,11 @@ VecgeomTrackView::operator=(Initializer_t const& init) // Set up current state and locate daughter volume vgstate_.Clear(); #if CELERITAS_VECGEOM_SURFACE - auto world = vecgeom::NavigationState::WorldId(); + // VecGeom's BVHSurfNav takes `int pvol_id ` but vecgeom's navtuple/index + // return NavIndex_t via `NavInd` :( + VgPlacedVolumeInt world = vecgeom::NavigationState::WorldId(); #else - auto const* world = &this->world(); + auto const* world = params_.scalars.world(); #endif // LocatePointIn sets `vgstate_` constexpr bool contains_point = true; @@ -644,17 +643,6 @@ CELER_FUNCTION bool VecgeomTrackView::is_next_boundary() const } } -//---------------------------------------------------------------------------// -/*! - * Get a reference to the world volume instance. - */ -CELER_FUNCTION auto VecgeomTrackView::world() const -> VgPlacedVol const& -{ - auto* pv = params_.scalars.world(); - CELER_ENSURE(pv); - return *pv; -} - //---------------------------------------------------------------------------// /*! * Get a reference to the current volume. From 94ba06416ca957603b72d62887566ed416f6a834 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 11:19:27 -0500 Subject: [PATCH 04/35] Update type based on BVH::LevelInside --- src/geocel/vg/VecgeomTypes.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geocel/vg/VecgeomTypes.hh b/src/geocel/vg/VecgeomTypes.hh index 405a618c31..d23143898c 100644 --- a/src/geocel/vg/VecgeomTypes.hh +++ b/src/geocel/vg/VecgeomTypes.hh @@ -57,7 +57,7 @@ namespace celeritas //---------------------------------------------------------------------------// using VgSurfaceInt = long; -using VgPlacedVolumeInt = int; +using VgPlacedVolumeInt = long; using vg_real_type = VECGEOM_PRECISION_NAMESPACE::Precision; #if defined(VECGEOM_BVH_SINGLE) || defined(__DOXYGEN__) From 53526f07bb1f21210ff7abe6a0c2c675d94c50f3 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 11:37:02 -0500 Subject: [PATCH 05/35] WIP: use vg bvh for 2.x solids --- src/geocel/vg/detail/SolidsNavigator.hh | 70 ++++++++----------------- src/geocel/vg/detail/SurfNavigator.hh | 7 +-- 2 files changed, 24 insertions(+), 53 deletions(-) diff --git a/src/geocel/vg/detail/SolidsNavigator.hh b/src/geocel/vg/detail/SolidsNavigator.hh index 0bd456db6b..10a321550e 100644 --- a/src/geocel/vg/detail/SolidsNavigator.hh +++ b/src/geocel/vg/detail/SolidsNavigator.hh @@ -11,8 +11,7 @@ #include #include #include -#include -#include +#include #include "corecel/Macros.hh" #include "corecel/Types.hh" @@ -38,79 +37,56 @@ class SolidsNavigator public: using VgPlacedVol = VgPlacedVolume; -#if CELER_VGNAV == CELER_VGNAV_PATH - using NavState = vecgeom::NavStatePath; -#else using NavState = detail::VgNavStateWrapper; -#endif + using NavImpl = vecgeom::BVHNavigator; //-----------------------------------------------------------------------// // Locate a point in the geometry hierarchy CELER_FUNCTION static void LocatePointIn(VgPlacedVol const* vol, - VgReal3 const& point, - NavState& nav, + VgReal3 const& localpos, + NavState& state, bool top, VgPlacedVol const* exclude = nullptr) { - ScopedVgNavState temp_nav{nav}; - if (exclude) - { - // Exclude the volume from the search - vecgeom::GlobalLocator::LocateGlobalPointExclVolume( - vol, exclude, point, temp_nav, top); - } - else - { - // TODO: eliminate this branch by always using Excl - // Locate the point in the volume hierarchy - vecgeom::GlobalLocator::LocateGlobalPoint( - vol, point, temp_nav, top); - } + ScopedVgNavState temp_state{state}; + NavImpl::LocatePointIn(vol, localpos, temp_state, top, exclude); } //-----------------------------------------------------------------------// - // FIXME: this *crosses* the volume CELER_FUNCTION static vg_real_type - ComputeStepAndNextVolume(VgReal3 const& glpos, - VgReal3 const& gldir, + ComputeStepAndNextVolume(VgReal3 const& pos, + VgReal3 const& dir, vg_real_type step_limit, NavState const& in_state, NavState& out_state) { - auto* curr_volume = in_state.Top()->GetLogicalVolume(); - - // simple dispatch implementation - ScopedVgNavState temp_out_state{out_state}; - auto* navigator = curr_volume->GetNavigator(); - real_type step = navigator->ComputeStepAndPropagatedState( - glpos, gldir, step_limit, in_state, temp_out_state); - - return step; + ScopedVgNavState temp_state{out_state}; + // TODO: maybe use ComputeStepToApproachNextVolume and eliminate next + // state? + return NavImpl::ComputeStepAndNextVolume( + pos, dir, step_limit, in_state, temp_state); } //-----------------------------------------------------------------------// // Computes the isotropic safety from the globalpoint - CELER_FUNCTION static double - ComputeSafety(VgReal3 const& glpos, - NavState const& curr, - vg_real_type safety + CELER_FUNCTION static vg_real_type + ComputeSafety(VgReal3 const& pos, + NavState const& state, + vg_real_type limit = std::numeric_limits::infinity()) { - auto* navigator = curr.Top()->GetLogicalVolume()->GetNavigator(); - real_type result - = navigator->GetSafetyEstimator()->ComputeSafety(glpos, curr); - result = vecCore::math::Min(result, safety); - - return result; + return NavImpl::ComputeSafety(pos, state, limit); } //-----------------------------------------------------------------------// // Relocate a state that was returned from ComputeStepAndNextVolume - CELER_FUNCTION static void - RelocateToNextVolume(VgReal3 const&, VgReal3 const&, NavState&) + CELER_FUNCTION static void RelocateToNextVolume(VgReal3 const& pos, + VgReal3 const& dir, + NavState& state) { - // Relocation is done previously :( + ScopedVgNavState temp_state{state}; + return NavImpl::RelocateToNextVolume(pos, dir, temp_state); } }; diff --git a/src/geocel/vg/detail/SurfNavigator.hh b/src/geocel/vg/detail/SurfNavigator.hh index aa02f31ab6..97cbec47c9 100644 --- a/src/geocel/vg/detail/SurfNavigator.hh +++ b/src/geocel/vg/detail/SurfNavigator.hh @@ -95,12 +95,7 @@ class SurfNavigator NavState& out_state, VgSurfaceInt& hitsurf) { - if (step_limit <= 0) - { - in_state.CopyTo(&out_state); - out_state.SetBoundaryState(false); - return step_limit; - } + CELER_EXPECT(step_limit > 0); ScopedVgNavState temp_out_state{out_state}; auto step = vgbrep::protonav::BVHSurfNavigator< From 95054ceadc9fc5472f1f01f8d010aac33439ffff Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 13:26:43 -0500 Subject: [PATCH 06/35] Delete unused methods --- src/geocel/vg/detail/BVHNavigator.hh | 314 ++++++-------------------- src/geocel/vg/detail/SurfNavigator.hh | 20 +- 2 files changed, 74 insertions(+), 260 deletions(-) diff --git a/src/geocel/vg/detail/BVHNavigator.hh b/src/geocel/vg/detail/BVHNavigator.hh index 340d7ca218..b9843c9349 100644 --- a/src/geocel/vg/detail/BVHNavigator.hh +++ b/src/geocel/vg/detail/BVHNavigator.hh @@ -101,140 +101,6 @@ class BVHNavigator } } - CELER_FUNCTION static void - RelocatePoint(VgReal3 const& localpoint, NavState& path) - { - CELER_EXPECT(!path.IsOutside()); - VgPlacedVol const* currentmother = path.Top(); - VgReal3 transformed = localpoint; - do - { - path.Pop(); - transformed = currentmother->GetTransformation()->InverseTransform( - transformed); - currentmother = path.Top(); - } while (currentmother - && (currentmother->IsAssembly() - || !currentmother->UnplacedContains(transformed))); - - if (currentmother) - { - path.Pop(); - return LocatePointIn(currentmother, transformed, path, false); - } - } - - private: - // Computes a step in the current volume from the localpoint into localdir, - // taking step_limit into account. If a volume is hit, the function calls - // out_state.SetBoundaryState(true) and hitcandidate is set to the hit - // daughter volume, or kept unchanged if the current volume is left. - CELER_FUNCTION static double - ComputeStepAndHit(VgReal3 const& localpoint, - VgReal3 const& localdir, - vg_real_type step_limit, - NavState const& in_state, - NavState& out_state, - VgPlacedVol const*& hitcandidate) - { - if (step_limit <= 0) - { - // We don't need to ask any solid, step not limited by geometry. - in_state.CopyTo(&out_state); - out_state.SetBoundaryState(false); - return 0; - } - - vg_real_type step = step_limit; - VgPlacedVol const* pvol = in_state.Top(); - - // need to calc DistanceToOut first - step = pvol->DistanceToOut(localpoint, localdir, step_limit); - - if (step < 0) - step = 0; - - if (pvol->GetDaughters().size() > 0) - { - auto bvh - = vecgeom::BVHManager::GetBVH(pvol->GetLogicalVolume()->id()); - bvh->CheckDaughterIntersections( - localpoint, localdir, step, pvol, hitcandidate); - } - - // now we have the candidates and we prepare the out_state - in_state.CopyTo(&out_state); - if (step == vecgeom::kInfLength && step_limit > 0) - { - out_state.SetBoundaryState(true); - do - { - out_state.Pop(); - } while (out_state.Top()->IsAssembly()); - - return vecgeom::kTolerance; - } - - // Is geometry further away than physics step? - if (step > step_limit) - { - // Then this is a phyics step and we don't need to do anything. - out_state.SetBoundaryState(false); - return step_limit; - } - - // Otherwise it is a geometry step and we push the point to the - // boundary. - out_state.SetBoundaryState(true); - - if (step < 0) - { - step = 0; - } - - return step; - } - - // Computes a step in the current volume from the localpoint into localdir, - // until the next daughter bounding box, taking step_limit into account. - CELER_FUNCTION static double ApproachNextVolume(VgReal3 const& localpoint, - VgReal3 const& localdir, - vg_real_type step_limit, - NavState const& in_state) - { - vg_real_type step = step_limit; - VgPlacedVol const* pvol = in_state.Top(); - - if (pvol->GetDaughters().size() > 0) - { - auto bvh - = vecgeom::BVHManager::GetBVH(pvol->GetLogicalVolume()->id()); - // bvh->CheckDaughterIntersections(localpoint, localdir, step, - // pvol, hitcandidate); - bvh->ApproachNextDaughter(localpoint, localdir, step, pvol); - // Make sure we don't "step" on next boundary - step -= 10 * vecgeom::kTolerance; - } - - if (step == vecgeom::kInfLength && step_limit > 0) - return 0; - - // Is geometry further away than physics step? - if (step > step_limit) - { - // Then this is a phyics step and we don't need to do anything. - return step_limit; - } - - if (step < 0) - { - step = 0; - } - - return step; - } - - public: // Computes the isotropic safety from the globalpoint. CELER_FUNCTION static double ComputeSafety(VgReal3 const& globalpoint, NavState const& state, @@ -258,88 +124,6 @@ class BVHNavigator return safety; } - // Computes a step from the globalpoint (which must be in the current - // volume) into globaldir, taking step_limit into account. If a volume is - // hit, the function calls out_state.SetBoundaryState(true) and relocates - // the state to the next volume. - CELER_FUNCTION static double - ComputeStepAndPropagatedState(VgReal3 const& globalpoint, - VgReal3 const& globaldir, - vg_real_type step_limit, - NavState const& in_state, - NavState& out_state, - vg_real_type push = 0) - { - // If we are on the boundary, push a bit more - if (in_state.IsOnBoundary()) - { - push += kBoundaryPush; - } - if (step_limit < push) - { - // Go as far as the step limit says, assuming there is no boundary. - // TODO: Does this make sense? - in_state.CopyTo(&out_state); - out_state.SetBoundaryState(false); - return step_limit; - } - step_limit -= push; - - // calculate local point/dir from global point/dir - VgReal3 localpoint; - VgReal3 localdir; - // Impl::DoGlobalToLocalTransformation(in_state, globalpoint, - // globaldir, localpoint, localdir); - vecgeom::Transformation3D m; - in_state.TopMatrix(m); - localpoint = m.Transform(globalpoint); - localdir = m.TransformDirection(globaldir); - // The user may want to move point from boundary before computing the - // step - localpoint += push * localdir; - - VgPlacedVol const* hitcandidate = nullptr; - vg_real_type step = ComputeStepAndHit( - localpoint, localdir, step_limit, in_state, out_state, hitcandidate); - step += push; - - if (out_state.IsOnBoundary()) - { - // Relocate the point after the step to refine out_state. - localpoint += (step + kBoundaryPush) * localdir; - - if (!hitcandidate) - { - // We didn't hit a daughter but instead we're exiting the - // current volume. - RelocatePoint(localpoint, out_state); - } - else - { - // Otherwise check if we're directly entering other daughters - // transitively. - localpoint - = hitcandidate->GetTransformation()->Transform(localpoint); - LocatePointIn(hitcandidate, localpoint, out_state, false); - } - - if (out_state.Top() != nullptr) - { - while (out_state.Top()->IsAssembly() - || out_state.HasSamePathAsOther(in_state)) - { - out_state.Pop(); - } - CELER_ASSERT(!out_state.Top() - ->GetLogicalVolume() - ->GetUnplacedVolume() - ->IsAssembly()); - } - } - - return step; - } - // Computes a step from the globalpoint (which must be in the current // volume) into globaldir, taking step_limit into account. If a volume is // hit, the function calls out_state.SetBoundaryState(true) and @@ -353,7 +137,7 @@ class BVHNavigator vg_real_type step_limit, NavState const& in_state, NavState& out_state, - vg_real_type push = 0) + vg_real_type push) { // If we are on the boundary, push a bit more if (in_state.IsOnBoundary()) @@ -418,30 +202,6 @@ class BVHNavigator return step; } - // Computes a step from the globalpoint (which must be in the current - // volume) into globaldir, taking step_limit into account. - CELER_FUNCTION static vg_real_type - ComputeStepToApproachNextVolume(VgReal3 const& globalpoint, - VgReal3 const& globaldir, - vg_real_type step_limit, - NavState const& in_state) - { - // calculate local point/dir from global point/dir - VgReal3 localpoint; - VgReal3 localdir; - // Impl::DoGlobalToLocalTransformation(in_state, globalpoint, - // globaldir, localpoint, localdir); - vecgeom::Transformation3D m; - in_state.TopMatrix(m); - localpoint = m.Transform(globalpoint); - localdir = m.TransformDirection(globaldir); - - vg_real_type step - = ApproachNextVolume(localpoint, localdir, step_limit, in_state); - - return step; - } - // Relocate a state that was returned from ComputeStepAndNextVolume: It // recursively locates the pushed point in the containing volume. CELER_FUNCTION static void RelocateToNextVolume(VgReal3 const& globalpoint, @@ -473,6 +233,78 @@ class BVHNavigator ->IsAssembly()); } } + + private: + // USED ONLY BY ComputeStepAndNextVolume + // Computes a step in the current volume from the localpoint into localdir, + // taking step_limit into account. If a volume is hit, the function calls + // out_state.SetBoundaryState(true) and hitcandidate is set to the hit + // daughter volume, or kept unchanged if the current volume is left. + CELER_FUNCTION static double + ComputeStepAndHit(VgReal3 const& localpoint, + VgReal3 const& localdir, + vg_real_type step_limit, + NavState const& in_state, + NavState& out_state, + VgPlacedVol const*& hitcandidate) + { + if (step_limit <= 0) + { + // We don't need to ask any solid, step not limited by geometry. + in_state.CopyTo(&out_state); + out_state.SetBoundaryState(false); + return 0; + } + + vg_real_type step = step_limit; + VgPlacedVol const* pvol = in_state.Top(); + + // need to calc DistanceToOut first + step = pvol->DistanceToOut(localpoint, localdir, step_limit); + + if (step < 0) + step = 0; + + if (pvol->GetDaughters().size() > 0) + { + auto bvh + = vecgeom::BVHManager::GetBVH(pvol->GetLogicalVolume()->id()); + bvh->CheckDaughterIntersections( + localpoint, localdir, step, pvol, hitcandidate); + } + + // now we have the candidates and we prepare the out_state + in_state.CopyTo(&out_state); + if (step == vecgeom::kInfLength && step_limit > 0) + { + out_state.SetBoundaryState(true); + do + { + out_state.Pop(); + } while (out_state.Top()->IsAssembly()); + + return vecgeom::kTolerance; + } + + // Is geometry further away than physics step? + if (step > step_limit) + { + // Then this is a phyics step and we don't need to do anything. + out_state.SetBoundaryState(false); + return step_limit; + } + + // Otherwise it is a geometry step and we push the point to the + // boundary. + out_state.SetBoundaryState(true); + + if (step < 0) + { + step = 0; + } + + return step; + } }; //---------------------------------------------------------------------------// diff --git a/src/geocel/vg/detail/SurfNavigator.hh b/src/geocel/vg/detail/SurfNavigator.hh index 97cbec47c9..02e8c67c12 100644 --- a/src/geocel/vg/detail/SurfNavigator.hh +++ b/src/geocel/vg/detail/SurfNavigator.hh @@ -108,25 +108,7 @@ class SurfNavigator return step; } - // Computes a step from the globalpoint (which must be in the current - // volume) into globaldir, taking step_limit into account. If a volume is - // hit, the function calls out_state.SetBoundaryState(true) and relocates - // the state to the next volume. - CELER_FUNCTION static vg_real_type - ComputeStepAndPropagatedState(VgReal3 const& globalpoint, - VgReal3 const& globaldir, - vg_real_type step_limit, - VgSurfaceInt& hit_surf, - NavState const& in_state, - NavState& out_state) - { - return ComputeStepAndNextVolume( - globalpoint, globaldir, step_limit, in_state, out_state, hit_surf); - } - - // Relocate a state that was returned from ComputeStepAndNextVolume: the - // surface model does this computation within ComputeStepAndNextVolume, so - // the relocation does nothing + // Relocate a state that was returned from ComputeStepAndNextVolume CELER_FUNCTION static void RelocateToNextVolume(VgReal3 const& globalpoint, VgReal3 const& globaldir, VgSurfaceInt hitsurf_index, From 4b723a8e122e6980db6ac012482d28134b748d18 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 13:55:23 -0500 Subject: [PATCH 07/35] Document fma parameters --- src/corecel/math/Algorithms.hh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/corecel/math/Algorithms.hh b/src/corecel/math/Algorithms.hh index 28fd2bfefd..6a66fbb3ff 100644 --- a/src/corecel/math/Algorithms.hh +++ b/src/corecel/math/Algorithms.hh @@ -537,6 +537,8 @@ inline CELER_FUNCTION T fastpow(T a, T b) /*! * Use fused multiply-add for generic calculations. * + * \f[ x \gets a \times b + y \f] + * * This provides a floating point specialization so that \c fma can be used in * code that is accelerated for floating point calculations but still works * correctly with integer arithmetic. From 13b40ad708b784f64493f9909ebe85d7340cc039 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 15:15:55 -0500 Subject: [PATCH 08/35] Move BVH bump into track view --- src/geocel/vg/VecgeomTrackView.hh | 50 +++++++++++++++----- src/geocel/vg/VecgeomTypes.hh | 6 +++ src/geocel/vg/detail/BVHNavigator.hh | 30 ++++++------ src/geocel/vg/detail/SolidsNavigator.hh | 2 - test/celeritas/field/FieldPropagator.test.cc | 7 ++- test/geocel/GeoTests.cc | 7 --- 6 files changed, 62 insertions(+), 40 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index cbf4e28722..92bf577179 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -54,8 +54,9 @@ namespace celeritas VecgeomTrackView geom(vg_params_ref, vg_state_ref, trackslot_id); \endcode * - * The "next distance" is cached as part of `find_next_step`, but it is only - * used when the immediate next call is `move_to_boundary`. + * \note VecGeom's solid model is imprcise and requires bumping. The bump + * values are arbitrary and may lead to geometry errors. The values of the bump + * constants do not account for the problem's length scale. */ class VecgeomTrackView { @@ -86,11 +87,6 @@ class VecgeomTrackView inline CELER_FUNCTION VecgeomTrackView& operator=(Initializer_t const& init); - //// STATIC ACCESSORS //// - - //! A tiny push to make sure tracks do not get stuck at boundaries - static CELER_CONSTEXPR_FUNCTION real_type extra_push() { return 1e-13; } - //// ACCESSORS //// //!@{ @@ -185,6 +181,22 @@ class VecgeomTrackView // Temporary data real_type next_step_{0}; + // Static data + + //! A tiny push to make sure tracks do not get stuck at boundaries (TODO: + //! DELETEME) + static constexpr real_type extra_push_ = 1e-13; + + //! Absolute tolerance used in internal vecgeom construction: + //! 1e-9 for double, 1e-3 for single + static constexpr vg_real_type abs_tolerance_ + = CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE ? 1e-9 : 1e-3; + + //! Tolerance used for solid model relocation bump + //! V1 solids by default; + static constexpr vg_real_type relocate_bump_ + = CELERITAS_VECGEOM_SURFACE ? 0 : 10 * abs_tolerance_; + //// HELPER FUNCTIONS //// // Whether any next distance-to-boundary has been found @@ -292,6 +304,7 @@ VecgeomTrackView::operator=(Initializer_t const& init) constexpr bool contains_point = true; Navigator::LocatePointIn( world, detail::to_vector(pos_), vgstate_, contains_point); + return *this; } @@ -465,7 +478,7 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) CELER_ASSERT((*next_surf_ != null_surface()) == vgnext_.IsOnBoundary()); } - next_step_ = max(next_step_, this->extra_push()); + next_step_ = max(next_step_, this->extra_push_); if (!this->is_next_boundary()) { @@ -482,9 +495,9 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) CELER_ENSURE(this->has_next_step()); CELER_ENSURE(result.distance > 0); - CELER_ENSURE(result.distance <= max(max_step, this->extra_push())); + CELER_ENSURE(result.distance <= max(max_step, this->extra_push_)); CELER_ENSURE(result.boundary || result.distance == max_step - || max_step < this->extra_push()); + || max_step < this->extra_push_); return result; } @@ -551,7 +564,22 @@ CELER_FUNCTION void VecgeomTrackView::cross_boundary() // Relocate to next tracking volume (maybe across multiple boundaries) if (vgnext_.Top() != nullptr) { - Navigator::RelocateToNextVolume(detail::to_vector(this->pos_), + VgReal3 bumped_pos; + if constexpr (CELERITAS_VECGEOM_SURFACE) + { + // Surface relocation is exact, assuming a well constructed + // geometry + bumped_pos = detail::to_vector(pos_); + } + else + { + for (auto i : range(3)) + { + bumped_pos[i] = fma(relocate_bump_, dir_[i], pos_[i]); + } + } + + Navigator::RelocateToNextVolume(bumped_pos, detail::to_vector(this->dir_), #if CELERITAS_VECGEOM_SURFACE *next_surf_, diff --git a/src/geocel/vg/VecgeomTypes.hh b/src/geocel/vg/VecgeomTypes.hh index d23143898c..34c4d471f8 100644 --- a/src/geocel/vg/VecgeomTypes.hh +++ b/src/geocel/vg/VecgeomTypes.hh @@ -38,6 +38,12 @@ # define CELER_VGNAV CELER_VGNAV_PATH #endif +#ifdef VECGEOM_SINGLE_PRECISION +static_assert(CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_FLOAT); +#else +static_assert(CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE); +#endif + namespace vecgeom { #if VECGEOM_VERSION >= 0x020000 diff --git a/src/geocel/vg/detail/BVHNavigator.hh b/src/geocel/vg/detail/BVHNavigator.hh index b9843c9349..6820d4c67f 100644 --- a/src/geocel/vg/detail/BVHNavigator.hh +++ b/src/geocel/vg/detail/BVHNavigator.hh @@ -48,8 +48,6 @@ class BVHNavigator using NavState = detail::VgNavStateWrapper; #endif - static constexpr vg_real_type kBoundaryPush = 10 * vecgeom::kTolerance; - //! Update path (which must be reset in advance) CELER_FUNCTION static void LocatePointIn(VgPlacedVol const* vol, @@ -136,18 +134,19 @@ class BVHNavigator VgReal3 const& globaldir, vg_real_type step_limit, NavState const& in_state, - NavState& out_state, - vg_real_type push) + NavState& out_state) { + // See VecgeomTrackView::relocate_bump_ ; this is from the VG 2.0 + // boundary, setting external push to zero + static constexpr vg_real_type kBoundaryPush = 10 * vecgeom::kTolerance; + // If we are on the boundary, push a bit more - if (in_state.IsOnBoundary()) - { - push += kBoundaryPush; - } + vg_real_type push = in_state.IsOnBoundary() ? kBoundaryPush : 0; + if (step_limit < push) { // Go as far as the step limit says, assuming there is no boundary. - // TODO: Does this make sense? + // TODO: DELETE in_state.CopyTo(&out_state); out_state.SetBoundaryState(false); return step_limit; @@ -176,6 +175,7 @@ class BVHNavigator { if (!hitcandidate) { + // Pop back up to the parent volume VgPlacedVol const* currentmother = out_state.Top(); VgReal3 transformed = localpoint; // Push the point inside the next volume. @@ -204,17 +204,15 @@ class BVHNavigator // Relocate a state that was returned from ComputeStepAndNextVolume: It // recursively locates the pushed point in the containing volume. - CELER_FUNCTION static void RelocateToNextVolume(VgReal3 const& globalpoint, - VgReal3 const& globaldir, - NavState& state) + CELER_FUNCTION static void + RelocateToNextVolume(VgReal3 const& globalpoint, + VgReal3 const& /* unused: globaldir */, + NavState& state) { - // Push the point inside the next volume. - VgReal3 pushed = globalpoint + kBoundaryPush * globaldir; - // Calculate local point from global point. vecgeom::Transformation3D m; state.TopMatrix(m); - VgReal3 localpoint = m.Transform(pushed); + VgReal3 localpoint = m.Transform(globalpoint); VgPlacedVol const* pvol = state.Top(); diff --git a/src/geocel/vg/detail/SolidsNavigator.hh b/src/geocel/vg/detail/SolidsNavigator.hh index 10a321550e..ab1fd31475 100644 --- a/src/geocel/vg/detail/SolidsNavigator.hh +++ b/src/geocel/vg/detail/SolidsNavigator.hh @@ -62,8 +62,6 @@ class SolidsNavigator NavState& out_state) { ScopedVgNavState temp_state{out_state}; - // TODO: maybe use ComputeStepToApproachNextVolume and eliminate next - // state? return NavImpl::ComputeStepAndNextVolume( pos, dir, step_limit, in_state, temp_state); } diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index e78a826869..138fcfd561 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -1222,13 +1222,12 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(electron_stuck)) = [&geo]() { return std::hypot(geo.pos()[0], geo.pos()[1]); }; EXPECT_SOFT_EQ(30.000000000000011, calc_radius()); - // NOTE: vecgeom 2.x-solids puts this position slightly *outside* the beam - // tube rather than *inside* if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) { - // TODO: VecGeom 2.x-solids starts to diverge here + // NOTE: vecgeom 2 solid model thinks r=30 + epsilon is *outside* the + // 30cm radius cyl EXPECT_EQ("vacuum_tube", this->volume_name(geo)); - GTEST_SKIP() << "FIXME: VecGeom 2.x-solid construction failure."; + GTEST_SKIP() << "VecGeom 2.x solid disagrees where the solid is"; } EXPECT_EQ("si_tracker", this->volume_name(geo)); { diff --git a/test/geocel/GeoTests.cc b/test/geocel/GeoTests.cc index 6499e67f94..6c6eae7e00 100644 --- a/test/geocel/GeoTests.cc +++ b/test/geocel/GeoTests.cc @@ -363,13 +363,6 @@ void FourLevelsGeoTest::test_detailed_tracking() const // Find the next step (to top edge of Shape1) but then scatter back // toward the sphere ASSERT_NO_THROW(next = geo.find_next_step(from_cm(10.0))); - if (test_->geometry_type() == "VecGeom" && using_solids_vg - && vecgeom_version >= Version{2, 0}) - { - // Solids VG navig issues here - both v1,v2 work the same though - EXPECT_GT(1e-12, to_cm(next.distance)); - GTEST_SKIP() << "FIXME: VG_solids navig issues: 1e-13 vs. 6"; - } EXPECT_SOFT_EQ(6, to_cm(next.distance)); geo.set_dir({-1, 0, 0}); EXPECT_VEC_SOFT_EQ((Real3{15, 10, 10}), to_cm(geo.pos())); From e0c42db976d994509167e4c62c635b8016f0b8c5 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 15:26:39 -0500 Subject: [PATCH 09/35] WIP: start comparing again --- test/geocel/GeoTests.cc | 26 ++++++-------------------- 1 file changed, 6 insertions(+), 20 deletions(-) diff --git a/test/geocel/GeoTests.cc b/test/geocel/GeoTests.cc index 6c6eae7e00..9a4d74802d 100644 --- a/test/geocel/GeoTests.cc +++ b/test/geocel/GeoTests.cc @@ -1548,13 +1548,9 @@ void SolidsGeoTest::test_trace() const ref.halfway_safeties[15] = 18.8833925371992; ref.halfway_safeties[16] = 42.8430141842906; } - if (test_->geometry_type() != "VecGeom" - || vecgeom_version < Version{2, 0} || CELERITAS_VECGEOM_SURFACE) - { - // TODO: VecGemo 2.x-solids still missing some shapes - auto tol = test_->tracking_tol(); - EXPECT_REF_NEAR(ref, result, tol); - } + + auto tol = test_->tracking_tol(); + EXPECT_REF_NEAR(ref, result, tol); } { SCOPED_TRACE("Lower +x"); @@ -1670,12 +1666,7 @@ void SolidsGeoTest::test_trace() const } auto tol = test_->tracking_tol(); - if (test_->geometry_type() != "VecGeom" - || vecgeom_version < Version{2, 0} || CELERITAS_VECGEOM_SURFACE) - { - // TODO: VecGemo 2.x-solids still missing some shapes - EXPECT_REF_NEAR(ref, result, tol); - } + EXPECT_REF_NEAR(ref, result, tol); } { SCOPED_TRACE("Middle +y"); @@ -1737,13 +1728,8 @@ void SolidsGeoTest::test_trace() const 74.5, }; - if (test_->geometry_type() != "VecGeom" - || vecgeom_version < Version{2, 0} || CELERITAS_VECGEOM_SURFACE) - { - // TODO: VecGemo 2.x-solids still missing some shapes - auto tol = test_->tracking_tol(); - EXPECT_REF_NEAR(ref, result, tol); - } + auto tol = test_->tracking_tol(); + EXPECT_REF_NEAR(ref, result, tol); } } From fce6f9fa4dccb57a934c18f0168c43bdb86b7ec6 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 15:37:34 -0500 Subject: [PATCH 10/35] Revert "Update type based on BVH::LevelInside" This reverts commit e89c24974f3b9b2a761374ae953b97210b6c5072. --- src/geocel/vg/VecgeomTypes.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geocel/vg/VecgeomTypes.hh b/src/geocel/vg/VecgeomTypes.hh index 34c4d471f8..bdf1c8bfaf 100644 --- a/src/geocel/vg/VecgeomTypes.hh +++ b/src/geocel/vg/VecgeomTypes.hh @@ -63,7 +63,7 @@ namespace celeritas //---------------------------------------------------------------------------// using VgSurfaceInt = long; -using VgPlacedVolumeInt = long; +using VgPlacedVolumeInt = int; using vg_real_type = VECGEOM_PRECISION_NAMESPACE::Precision; #if defined(VECGEOM_BVH_SINGLE) || defined(__DOXYGEN__) From bb33f4d6ce9025d97658fb0c7206f6ad824ba542 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 17:18:54 -0500 Subject: [PATCH 11/35] Remove extra push from vg track view --- src/geocel/vg/VecgeomTrackView.hh | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 92bf577179..980c0f9887 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -183,10 +183,6 @@ class VecgeomTrackView // Static data - //! A tiny push to make sure tracks do not get stuck at boundaries (TODO: - //! DELETEME) - static constexpr real_type extra_push_ = 1e-13; - //! Absolute tolerance used in internal vecgeom construction: //! 1e-9 for double, 1e-3 for single static constexpr vg_real_type abs_tolerance_ @@ -460,6 +456,7 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) { *next_surf_ = null_surface(); } + // TODO: vgnext is simply copied and the boundary flag optionally set next_step_ = Navigator::ComputeStepAndNextVolume(detail::to_vector(pos_), detail::to_vector(dir_), @@ -471,6 +468,7 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) *next_surf_ #endif ); + CELER_ASSERT(next_step_ > 0); if constexpr (CELERITAS_VECGEOM_SURFACE) { // Our accessor uses the next_surf_ state, but the temporary used for @@ -478,8 +476,6 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) CELER_ASSERT((*next_surf_ != null_surface()) == vgnext_.IsOnBoundary()); } - next_step_ = max(next_step_, this->extra_push_); - if (!this->is_next_boundary()) { // Soft equivalence between distance and max step is because the @@ -495,9 +491,8 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) CELER_ENSURE(this->has_next_step()); CELER_ENSURE(result.distance > 0); - CELER_ENSURE(result.distance <= max(max_step, this->extra_push_)); - CELER_ENSURE(result.boundary || result.distance == max_step - || max_step < this->extra_push_); + CELER_ENSURE(result.distance <= max_step); + CELER_ENSURE(result.boundary || result.distance == max_step); return result; } From eb2cd54ce0bdad1cec391a12f7772e24b30a88f9 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 17:28:26 -0500 Subject: [PATCH 12/35] Add error checking code to vecgeom track view --- src/geocel/vg/VecgeomTrackView.hh | 48 +++++++++++++++++++++++++++++-- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 980c0f9887..30433bdc2e 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -41,6 +41,11 @@ # include "detail/VgNavStateWrapper.hh" #endif +#if !CELER_DEVICE_COMPILE +# include "corecel/io/Logger.hh" +# include "corecel/io/Repr.hh" +#endif + namespace celeritas { //---------------------------------------------------------------------------// @@ -117,7 +122,7 @@ class VecgeomTrackView // Whether the track is exactly on a surface CELER_FORCEINLINE_FUNCTION bool is_on_boundary() const; //! Whether the last operation resulted in an error - CELER_FORCEINLINE_FUNCTION bool failed() const { return false; } + CELER_FORCEINLINE_FUNCTION bool failed() const { return failed_; } // Get the normal vector of the current surface inline CELER_FUNCTION Real3 normal() const; @@ -180,6 +185,7 @@ class VecgeomTrackView // Temporary data real_type next_step_{0}; + bool failed_{false}; // Static data @@ -290,7 +296,7 @@ VecgeomTrackView::operator=(Initializer_t const& init) // Set up current state and locate daughter volume vgstate_.Clear(); #if CELERITAS_VECGEOM_SURFACE - // VecGeom's BVHSurfNav takes `int pvol_id ` but vecgeom's navtuple/index + // VecGeom's BVHSurfNav takes `int pvol_id` but vecgeom's navtuple/index // return NavIndex_t via `NavInd` :( VgPlacedVolumeInt world = vecgeom::NavigationState::WorldId(); #else @@ -301,6 +307,15 @@ VecgeomTrackView::operator=(Initializer_t const& init) Navigator::LocatePointIn( world, detail::to_vector(pos_), vgstate_, contains_point); + if (CELER_UNLIKELY(vgstate_.IsOutside())) + { +#if !CELER_DEVICE_COMPILE + auto msg = CELER_LOG_LOCAL(error); + msg << "Failed to initialize geometry state at " << repr(pos_); +#endif + failed_ = true; + } + return *this; } @@ -468,7 +483,21 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) *next_surf_ #endif ); - CELER_ASSERT(next_step_ > 0); + + if (CELER_UNLIKELY(!(next_step_ > 0))) + { +#if !CELER_DEVICE_COMPILE + auto msg = CELER_LOG_LOCAL(error); + msg << "Failed to find next step at " << repr(pos_) << " along " + << repr(dir_) << ": computed step is " << repr(next_step_); +#endif + failed_ = true; + Propagation result; + result.distance = next_step_; + result.boundary = false; + return result; + } + if constexpr (CELERITAS_VECGEOM_SURFACE) { // Our accessor uses the next_surf_ state, but the temporary used for @@ -582,6 +611,19 @@ CELER_FUNCTION void VecgeomTrackView::cross_boundary() vgnext_); } + // Relocation should have changed volume + if (CELER_UNLIKELY(vgstate_.HasSamePathAsOther(vgnext_))) + { +#if !CELER_DEVICE_COMPILE + auto msg = CELER_LOG_LOCAL(error); + msg << "Failed to cross boundary: unique volume instance is the same " + "before and after at position " + << repr(pos_); +#endif + failed_ = true; + return; + } + vgstate_ = vgnext_; CELER_ENSURE(this->is_on_boundary()); From fdf143263b8e7f297c90bda61cc4d2d779d41ad5 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 25 Nov 2025 18:04:38 -0500 Subject: [PATCH 13/35] WIP: field propagation --- test/celeritas/field/FieldPropagator.test.cc | 84 +++++++++++++++----- 1 file changed, 65 insertions(+), 19 deletions(-) diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index 138fcfd561..f8d6e1c306 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -1031,7 +1031,8 @@ TEST_F(TwoBoxesTest, for (int i : range(2)) { SCOPED_TRACE(i); - auto result = propagate(radius * dtheta); + Propagation result; + EXPECT_NO_THROW(result = propagate(radius * dtheta)); if (result.boundary) { geo.cross_boundary(); @@ -1246,7 +1247,7 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(electron_stuck)) { // Surface geometry does not intersect the cylinder boundary, so // the track keeps going until the "looping" counter is hit - EXPECT_SOFT_EQ(1.0314309658010318e-13, result.distance); + EXPECT_LT(result.distance, 2e-13); EXPECT_FALSE(result.looping); } else @@ -1266,18 +1267,29 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(electron_stuck)) { auto integrate = make_mag_field_integrator( field, particle.charge()); + auto propagate = make_field_propagator(integrate, driver_options, particle, geo); - auto result = propagate(30); + Propagation result; + try + { + result = propagate(30); + } + catch (RuntimeError const& e) + { + if (using_surface_vg) + { + EXPECT_TRUE(geo.failed()); + GTEST_SKIP() + << "FIXME: VecGeom surface model fails: " << e.what(); + } + FAIL() << e.what(); + } + EXPECT_EQ(result.boundary, geo.is_on_boundary()); EXPECT_SOFT_NEAR( double{30}, static_cast(integrate.count()), 0.2); - if (using_surface_vg) - { - GTEST_SKIP() << "FIXME: VecGeom surface model fails a boundary " - "requirement."; - } ASSERT_TRUE(geo.is_on_boundary()); if (geo.check_normal()) @@ -1322,7 +1334,8 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) field, particle.charge()); auto propagate = make_field_propagator(integrate, driver_options, particle, geo); - auto result = propagate(1.39170198361108938e-05); + Propagation result; + EXPECT_NO_THROW(result = propagate(1.39170198361108938e-05)); EXPECT_EQ(result.boundary, geo.is_on_boundary()); EXPECT_EQ("em_calorimeter", this->volume_name(geo)); EXPECT_SOFT_EQ(125.00000000000001, calc_radius()); @@ -1336,8 +1349,26 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) geo.set_dir({-1.31178657592616127e-01, -8.29310561920304168e-01, -5.43172303859124073e-01}); - geo.cross_boundary(); // TODO: hack cross_boundary to handle this case - successful_reentry = (this->volume_name(geo) == "em_calorimeter"); + try + { + geo.cross_boundary(); + successful_reentry = (this->volume_name(geo) == "em_calorimeter"); + } + catch (RuntimeError const& e) + { + if (geo.failed()) + { + // VecGeom can set the 'failed' flag, which CheckedGeoTrackView + // turns into an error. Ignore that error for now and say it's + // a failure + successful_reentry = false; + } + else + { + // Something else happened... + throw; + } + } if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_ORANGE) { // ORANGE should successfully reenter. However, under certain @@ -1345,7 +1376,7 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) // so we don't test in all cases. EXPECT_EQ("em_calorimeter", this->volume_name(geo)); } - else + else if (!geo.failed()) { EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; } @@ -1359,11 +1390,14 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) CELER_LOG(warning) << "Reentry failed for " << cmake::core_geo << " geometry: post-propagation volume is " << this->volume_name(geo); + if (!scoped_log_.empty()) + { + CELER_LOG(warning) << scoped_log_; + } if (using_solids_vg && CELERITAS_VECGEOM_VERSION < 0x020000) { - // FIXME: VecGeom 1.x navigation failure (solids model) - EXPECT_EQ("world", this->volume_name(geo)); + EXPECT_EQ("em_calorimeter", this->volume_name(geo)); } else { @@ -1374,17 +1408,18 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) // slightly differently. Only surface model thinks the surface // was actually crossed, therefore the next step will find distinct // results - auto result = geo.find_next_step(1); + Propagation result; + EXPECT_NO_THROW(result = geo.find_next_step(1)); EXPECT_LT(result.distance, 2e-8); - if (result.distance < 1e-6) + if (result.distance < 1e-6 && result.boundary) { geo.move_to_boundary(); geo.cross_boundary(); // back into em_calorimeter } // then they are back to agreement - result = geo.find_next_step(1); + EXPECT_NO_THROW(result = geo.find_next_step(1)); EXPECT_EQ(result.distance, 1); EXPECT_FALSE(result.boundary); EXPECT_TRUE(geo.is_on_boundary()); @@ -1408,7 +1443,7 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) Propagation result; // This absurdly long step is because in the "failed" case the // track thinks it's in the world volume (nearly vacuum) - result = propagate(2.12621374950874703e+21); + EXPECT_NO_THROW(result = propagate(2.12621374950874703e+21)); if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_GEANT4 && result.boundary != geo.is_on_boundary()) @@ -1510,7 +1545,6 @@ TEST_F(CmseTest, coarse) } catch (RuntimeError const& e) { - // Failure during Geant4 propagation CELER_LOG(error) << e.what(); failed = true; break; @@ -1549,6 +1583,18 @@ TEST_F(CmseTest, coarse) expected_num_integration = {80659, 1670, 1956, 1092}; EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; } + else if (using_solids_vg && failed) + { + expected_num_boundary[1] = 37; + expected_num_step[1] = 179; + expected_num_intercept[1] = 614; + expected_num_integration[1] = 1666; + static char const* const expected_log_messages[] = { + R"(Moved internally from boundary but safety didn't increase: volume 18 from {10.32, -6.565, 796.9} to {10.32, -6.565, 796.9} (distance: 1.000e-4))", + R"(Failed to find next step at {10.32, -6.565, 796.9} along {0.6896, -0.1485, 0.7088})"}; + EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()) + << scoped_log_; + } else if (using_solids_vg) { // Bumped (platform-dependent!): counts change a bit From 24c0b52ea52e30b1555f669d6d6ee038e11008ba Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 27 Nov 2025 12:48:20 -0500 Subject: [PATCH 14/35] Add comments --- src/geocel/vg/VecgeomTrackView.hh | 3 +++ src/geocel/vg/detail/SurfNavigator.hh | 6 +----- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 30433bdc2e..094a1f527b 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -672,6 +672,9 @@ CELER_FUNCTION void VecgeomTrackView::move_internal(Real3 const& pos) * * This happens after a scattering event or movement inside a magnetic field. * It resets the calculated distance-to-boundary. + * + * \todo If on a boundary, determine as with ORANGE whether we should cancel + * the surface crossing. */ CELER_FUNCTION void VecgeomTrackView::set_dir(Real3 const& newdir) { diff --git a/src/geocel/vg/detail/SurfNavigator.hh b/src/geocel/vg/detail/SurfNavigator.hh index 02e8c67c12..b48e7b8c67 100644 --- a/src/geocel/vg/detail/SurfNavigator.hh +++ b/src/geocel/vg/detail/SurfNavigator.hh @@ -82,11 +82,7 @@ class SurfNavigator // Computes a step from the globalpoint (which must be in the current // volume) into globaldir, taking step_limit into account. If a volume is - // hit, the function calls out_state.SetBoundaryState(true) and relocates - // the state to the next volume. - // - // The surface model does automatic relocation, so this function does it as - // well. + // hit, the function calls out_state.SetBoundaryState(true) CELER_FUNCTION static vg_real_type ComputeStepAndNextVolume(VgReal3 const& globalpoint, VgReal3 const& globaldir, From 585a3007f9490bc33170d11207fb84ba40c98e85 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 11:01:34 -0500 Subject: [PATCH 15/35] Add length units to vecgeom debug output --- src/geocel/detail/LengthUnits.hh | 3 +++ src/geocel/vg/VecgeomTrackView.hh | 16 ++++++++++------ 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/src/geocel/detail/LengthUnits.hh b/src/geocel/detail/LengthUnits.hh index 2adbdc7f8f..20900cad0c 100644 --- a/src/geocel/detail/LengthUnits.hh +++ b/src/geocel/detail/LengthUnits.hh @@ -23,14 +23,17 @@ namespace lengthunits CELER_ICC meter{100}; CELER_ICC centimeter{1}; CELER_ICC millimeter{0.1}; +inline constexpr char const label[] = "cm"; #elif CELERITAS_UNITS == CELERITAS_UNITS_SI CELER_ICC meter{1}; CELER_ICC centimeter{0.01}; CELER_ICC millimeter{0.001}; +inline constexpr char const label[] = "m"; #elif CELERITAS_UNITS == CELERITAS_UNITS_CLHEP CELER_ICC meter{1000}; CELER_ICC centimeter{10}; CELER_ICC millimeter{1}; +inline constexpr char const label[] = "mm"; #else # error "CELERITAS_UNITS is undefined" #endif diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 094a1f527b..f021776273 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -44,6 +44,7 @@ #if !CELER_DEVICE_COMPILE # include "corecel/io/Logger.hh" # include "corecel/io/Repr.hh" +# include "geocel/detail/LengthUnits.hh" #endif namespace celeritas @@ -189,7 +190,7 @@ class VecgeomTrackView // Static data - //! Absolute tolerance used in internal vecgeom construction: + //! Absolute tolerance in vecgeom::kTolerance //! 1e-9 for double, 1e-3 for single static constexpr vg_real_type abs_tolerance_ = CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE ? 1e-9 : 1e-3; @@ -311,7 +312,8 @@ VecgeomTrackView::operator=(Initializer_t const& init) { #if !CELER_DEVICE_COMPILE auto msg = CELER_LOG_LOCAL(error); - msg << "Failed to initialize geometry state at " << repr(pos_); + msg << "Failed to initialize geometry state at " << repr(pos_) + << lengthunits::label; #endif failed_ = true; } @@ -488,8 +490,10 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) { #if !CELER_DEVICE_COMPILE auto msg = CELER_LOG_LOCAL(error); - msg << "Failed to find next step at " << repr(pos_) << " along " - << repr(dir_) << ": computed step is " << repr(next_step_); + msg << "Failed to find next step at " << repr(pos_) << ' ' + << lengthunits::label << " along " << repr(dir_) + << ": computed step is " << repr(next_step_) << ' ' + << lengthunits::label; #endif failed_ = true; Propagation result; @@ -617,8 +621,8 @@ CELER_FUNCTION void VecgeomTrackView::cross_boundary() #if !CELER_DEVICE_COMPILE auto msg = CELER_LOG_LOCAL(error); msg << "Failed to cross boundary: unique volume instance is the same " - "before and after at position " - << repr(pos_); + "before and after at " + << repr(pos_) << ' ' << lengthunits::label; #endif failed_ = true; return; From 0a573bc380266b1fa4fd6789f592b192d611a616 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 12:19:48 -0500 Subject: [PATCH 16/35] Clear failure flag on reinitialize --- src/geocel/vg/VecgeomTrackView.hh | 1 + 1 file changed, 1 insertion(+) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index f021776273..94d7012136 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -266,6 +266,7 @@ CELER_FUNCTION VecgeomTrackView& VecgeomTrackView::operator=(Initializer_t const& init) { CELER_EXPECT(is_soft_unit_vector(init.dir)); + failed_ = false; // Initialize direction dir_ = init.dir; From d6cbcbe12830dd893e52abd23a87d87805b0f87c Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 14:39:26 -0500 Subject: [PATCH 17/35] Vecgeom solid 1 and 2 pass vg tests --- test/geocel/GenericGeoResults.cc | 22 ++++ test/geocel/GenericGeoResults.hh | 6 + test/geocel/GenericGeoTestInterface.cc | 51 ++++---- test/geocel/GeoTests.cc | 169 ++++++++++++++++++------- 4 files changed, 175 insertions(+), 73 deletions(-) diff --git a/test/geocel/GenericGeoResults.cc b/test/geocel/GenericGeoResults.cc index af6ee20a93..8a58573933 100644 --- a/test/geocel/GenericGeoResults.cc +++ b/test/geocel/GenericGeoResults.cc @@ -28,6 +28,18 @@ namespace celeritas { namespace test { +namespace +{ + +template +void erase_after(std::vector& vec, std::size_t idx) +{ + auto iter = vec.begin() + std::min(idx, vec.size()); + vec.erase(iter, vec.end()); +} + +} // namespace + //---------------------------------------------------------------------------// ::testing::AssertionResult IsNormalEquiv(char const* expected_expr, char const* actual_expr, @@ -78,6 +90,16 @@ void GenericGeoTrackingResult::clear_boring_normals() } } +void GenericGeoTrackingResult::fail_at(std::size_t index) +{ + erase_after(volumes, index); + erase_after(volume_instances, index); + erase_after(distances, index); + erase_after(dot_normal, index); + erase_after(halfway_safeties, index); + volumes.push_back("[FAILURE]"); +} + void GenericGeoTrackingResult::print_expected() const { using std::cout; diff --git a/test/geocel/GenericGeoResults.hh b/test/geocel/GenericGeoResults.hh index 8ba79610e1..0663757dcd 100644 --- a/test/geocel/GenericGeoResults.hh +++ b/test/geocel/GenericGeoResults.hh @@ -74,6 +74,12 @@ struct GenericGeoTrackingResult // Whether surface normals are disabled bool disabled_surface_normal() const; + //! Add a failure sentinel at the end + void fail() { this->fail_at(volumes.size()); } + + // Add a failure sentinel at a certain index + void fail_at(std::size_t index); + // Print expected expression to cout void print_expected() const; }; diff --git a/test/geocel/GenericGeoTestInterface.cc b/test/geocel/GenericGeoTestInterface.cc index 623e57a523..aa81450270 100644 --- a/test/geocel/GenericGeoTestInterface.cc +++ b/test/geocel/GenericGeoTestInterface.cc @@ -49,30 +49,33 @@ auto GenericGeoTestInterface::track(Real3 const& pos, result.disable_surface_normal(); } -#define GGTI_EXPECT_NO_THROW(ACTION) \ - try \ - { \ - ACTION; \ - } \ - catch (CheckedGeoError const& e) \ - { \ - auto const& d = e.details(); \ - auto msg = CELER_LOG(debug); \ - msg << "Failed "; \ - if (!d.condition.empty()) \ - { \ - msg << '\'' << d.condition << "' "; \ - } \ - msg << "at " << d.file << ':' << d.line << " during '" << #ACTION \ - << "'"; \ - ADD_FAILURE() << d.what; \ - return result; \ - } \ - catch (std::exception const& e) \ - { \ - ADD_FAILURE() << "Caught exception during '" << #ACTION \ - << "': " << e.what() << ": " << geo; \ - return result; \ +#define GGTI_EXPECT_NO_THROW(ACTION) \ + try \ + { \ + ACTION; \ + } \ + catch (CheckedGeoError const& e) \ + { \ + auto const& d = e.details(); \ + { \ + auto msg = CELER_LOG(debug); \ + msg << "Failed "; \ + if (!d.condition.empty()) \ + { \ + msg << '\'' << d.condition << "' "; \ + } \ + msg << "at " << d.file << ':' << d.line << " during '" << #ACTION \ + << "'"; \ + } \ + CELER_LOG(error) << "Failed: " << d.what; \ + result.fail(); \ + return result; \ + } \ + catch (std::exception const& e) \ + { \ + ADD_FAILURE() << "Caught exception during '" << #ACTION \ + << "': " << e.what() << ": " << geo; \ + return result; \ } // Note: position is scaled according to test diff --git a/test/geocel/GeoTests.cc b/test/geocel/GeoTests.cc index 9a4d74802d..164aece042 100644 --- a/test/geocel/GeoTests.cc +++ b/test/geocel/GeoTests.cc @@ -24,6 +24,43 @@ #include "TestMacros.hh" #include "UnitUtils.hh" +/*! + * Allow a statement to skip the test when using a certain geometry. + * + * Run \c STATEMENT. When \c COND is true, the statement is expected to throw, + * and + */ +#define SHOULD_FAIL_WHEN(STATEMENT, COND) \ + do \ + { \ + bool succeeded_{false}; \ + try \ + { \ + STATEMENT; \ + succeeded_ = true; \ + } \ + catch (::celeritas::test::CheckedGeoError const& e) \ + { \ + if (COND) \ + { \ + GTEST_SKIP() << "Geometry fails when " << #COND << ": " \ + << e.details().what; \ + } \ + else \ + { \ + FAIL() << "'" << #STATEMENT \ + << "' failed unexpectedly: " << e.what(); \ + } \ + } \ + catch (std::exception const& e) \ + { \ + FAIL() << "'" << #STATEMENT \ + << "' failed unexpectedly: " << e.what(); \ + } \ + \ + EXPECT_EQ(!(COND), succeeded_); \ + } while (0) + namespace celeritas { namespace test @@ -198,6 +235,18 @@ void CmseGeoTest::test_trace() const 29.931406871193, 40.276406871193, 57.573593128807, 57.573593128807, 57.573593128807, 57.573593128807,}; // clang-format on + if (test_->geometry_type() == "VecGeom" + && vecgeom_version >= Version{2, 0} && using_solids_vg) + { + /* + * Failed to cross boundary: + * unique volume instance is the same before and after + * at {30, 30, 655.04101161124} [cm] along {0, 0, 1}, + * [FAILED] [ON BOUNDARY] in OCMS_PV/CMSE/MUON + */ + ref.fail_at(12); + } + auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } @@ -217,6 +266,18 @@ void CmseGeoTest::test_trace() const ref.halfway_safeties = {85, 267.5, 85.85, 60.4, 0.078366388350267, 2.343262600759, 0.078366388350244, 60.4, 85.85, 267.5, 460,}; // clang-format on + if (test_->geometry_type() == "VecGeom" + && vecgeom_version >= Version{2, 0} && using_solids_vg) + { + /* + * Failed to cross boundary: + * unique volume instance is the same before and after + * at {-295, 0, -48.5} [cm] along {1, 0, 0}, + * [FAILED] [ON BOUNDARY] in OCMS_PV/CMSE/MUON + */ + ref.fail_at(2); + } + auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } @@ -229,6 +290,18 @@ void CmseGeoTest::test_trace() const ref.distances = {12.495, 287.505, 530, 920}; ref.dot_normal = {}; // All normals are along track dir ref.halfway_safeties = {6.2475, 47.95, 242, 460}; + if (test_->geometry_type() == "VecGeom" + && vecgeom_version >= Version{2, 0} && using_solids_vg) + { + /* + * Failed to cross boundary: + * unique volume instance is the same before and after + * at {300, 0, 1328} [cm] along {1, 0, 0}, + * [FAILED] [ON BOUNDARY] in OCMS_PV/CMSE/OQUA@0 + */ + ref.fail_at(2); + } + auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } @@ -370,10 +443,11 @@ void FourLevelsGeoTest::test_detailed_tracking() const EXPECT_TRUE(geo.is_on_boundary()); // Check the distance to the sphere boundary again, then scatter - // into the sphere (this may be a "bump": 1e-13 for surface VG, Geant4; - // 1e-8 for volume VG; BUT exactly zero for ORANGE thanks to - // "reentrant" logic) - ASSERT_NO_THROW(next = geo.find_next_step(from_cm(20.0))); + // into the sphere (this may be a "bump": 1e-13 for Geant4; + // BUT exactly zero for ORANGE thanks to "reentrant" logic) + SHOULD_FAIL_WHEN(next = geo.find_next_step(from_cm(20.0)), + vecgeom_version >= Version{2} + && test_->geometry_type() == "VecGeom"); EXPECT_LE(next.distance, to_cm(1e-8)); ASSERT_TRUE(next.boundary); if (next.distance > 0) @@ -404,26 +478,9 @@ void FourLevelsGeoTest::test_detailed_tracking() const else { geo.set_dir(Real3{1, 0, 0}); - if (test_->geometry_type() == "VecGeom" - && CELERITAS_VECGEOM_SURFACE) - { - // Assertion failure in NavStateTuple::PushDaughterImpl: - // trying to push into a daughter but there are none - // (pv ID 1) - GTEST_SKIP() << "FIXME: vecgeom surface breaks"; - } - - ASSERT_NO_THROW(geo.cross_boundary()); - if (test_->geometry_type() == "VecGeom") - { - // FIXME: boundary crossing doesn't change volume like it - // should - EXPECT_EQ("Shape2", test_->volume_name(geo)); - } - else - { - EXPECT_EQ("Shape1", test_->volume_name(geo)); - } + SHOULD_FAIL_WHEN(geo.cross_boundary(), + test_->geometry_type() == "VecGeom"); + EXPECT_EQ("Shape1", test_->volume_name(geo)); geo.set_dir(Real3{-1, 0, 0}); ASSERT_NO_THROW(geo.cross_boundary()); } @@ -1282,12 +1339,6 @@ void ReplicaGeoTest::test_trace() const } delete_orange_safety(*test_, ref, result); - if (test_->geometry_type() != "VecGeom" - || vecgeom_version < Version{2, 0} || CELERITAS_VECGEOM_SURFACE) - { - // TODO: VecGemo 2.x-solids returns wrong distance values - EXPECT_REF_NEAR(ref, result, tol); - } } } @@ -1324,7 +1375,7 @@ void ReplicaGeoTest::test_volume_stack() const if (test_->geometry_type() != "VecGeom" || vecgeom_version < Version{2, 0} || CELERITAS_VECGEOM_SURFACE) { - // TODO: VecGemo 2.x-solids returns wrong volume instances + // TODO: VecGeom 2.x-solids returns wrong volume instances EXPECT_REF_EQ(ref, result); } } @@ -1547,6 +1598,17 @@ void SolidsGeoTest::test_trace() const ref.halfway_safeties[14] = 42.8397753718277; ref.halfway_safeties[15] = 18.8833925371992; ref.halfway_safeties[16] = 42.8430141842906; + + if (vecgeom_version >= Version{2}) + { + /* + * Unique volume instance is the same before and after + * Failed during cross_boundary: + * at {258.1, 0, 0.5} [cm] along {-1, 0, 0}, + * [FAILED] [ON BOUNDARY] in World_PV/polycone1_PV + */ + ref.fail_at(4); + } } auto tol = test_->tracking_tol(); @@ -1655,6 +1717,18 @@ void SolidsGeoTest::test_trace() const ref.halfway_safeties[13] = 19.0382940808067; ref.halfway_safeties[14] = 0.5; ref.halfway_safeties[17] = 28.6150602709819; + + if (vecgeom_version >= Version{2}) + { + /* + * Unique volume instance is the same before and after + * Failed during cross_boundary: + * at {-335.05, -125, 0.5} [cm] along {1, 0, 0}, + * [FAILED] [ON BOUNDARY] in World_PV/arb8b_PV + */ + ref.halfway_safeties[4] = 39.0470100365853; + ref.fail_at(5); + } } if (test_->geometry_type() == "Geant4" @@ -2382,7 +2456,6 @@ void TwoBoxesGeoTest::test_reentrant() const // Cross into the new volume, needed for optical physics (+; +,-) geo.cross_boundary(); - EXPECT_TRUE(geo.is_on_boundary()); if (check_normal) { EXPECT_NORMAL_EQUIV((Real3{1, 0, 0}), geo.normal()); @@ -2399,18 +2472,22 @@ void TwoBoxesGeoTest::test_reentrant() const EXPECT_EQ("world", test_->volume_name(geo)); // Cross back into previous volume (-; -,-) - if (CELERITAS_DEBUG && test_->geometry_type() == "Geant4") + // GeantGTV has an extra check because we know it can't do this :( + if (test_->geometry_type() == "Geant4") { - // GeantGTV has an extra check because we know it can't do this :( - EXPECT_THROW(geo.cross_boundary(), DebugError); - GTEST_SKIP() << "Consecutive boundary crossing fails for G4"; + // TODO: Geant4 does not allow crossing to new volume and returning + // to old + if (CELERITAS_DEBUG) + { + EXPECT_THROW(geo.cross_boundary(), DebugError); + } + GTEST_SKIP() << "Geant4 does not allow crossing to new volume and " + "returning to old"; } else { - // Typical case ASSERT_NO_THROW(geo.cross_boundary()); } - EXPECT_TRUE(geo.is_on_boundary()); if (check_normal) { @@ -2425,9 +2502,9 @@ void TwoBoxesGeoTest::test_reentrant() const { GTEST_SKIP() << "Unexpected failure to cross volume"; } - if ("[OUTSIDE]" == test_->volume_name(geo)) + if (geo.is_outside()) { - GTEST_SKIP() << "FIXME: Unexpected track location."; + GTEST_SKIP() << "FIXME: Unexpected track location"; } } @@ -2435,9 +2512,7 @@ void TwoBoxesGeoTest::test_reentrant() const { // VecGeom with surfaces seems to have issues here EXPECT_EQ("[OUTSIDE]", test_->volume_name(geo)); - { - GTEST_SKIP() << "FIXME: VecGeom v2.x-surface misses inner volume."; - } + GTEST_SKIP() << "FIXME: VecGeom v2.x-surface misses inner volume."; } EXPECT_EQ("inner", test_->volume_name(geo)); @@ -2544,7 +2619,8 @@ void TwoBoxesGeoTest::test_tangent() const // Crossing should *not* change volumes (-; -,-) { SCOPED_TRACE("trying to cross"); - ASSERT_NO_THROW(geo.cross_boundary()); + SHOULD_FAIL_WHEN(geo.cross_boundary(), + test_->geometry_type() == "VecGeom"); EXPECT_TRUE(geo.is_on_boundary()); if (test_->geometry_type() == "Geant4") { @@ -2552,11 +2628,6 @@ void TwoBoxesGeoTest::test_tangent() const EXPECT_EQ("world", test_->volume_name(geo)); GTEST_SKIP() << "Unexpected boundary crossing"; } - else if (test_->geometry_type() == "VecGeom" - && "world" == test_->volume_name(geo)) - { - GTEST_SKIP() << "Unexpected boundary crossing"; - } EXPECT_EQ("inner", test_->volume_name(geo)); } From 9b0eaa0f7d1f7af832eb190dbbce7faeaa202eb6 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 15:02:26 -0500 Subject: [PATCH 18/35] All vg, g4, orange passes geo tests --- test/geocel/GeoTests.cc | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/geocel/GeoTests.cc b/test/geocel/GeoTests.cc index 164aece042..d107f4d836 100644 --- a/test/geocel/GeoTests.cc +++ b/test/geocel/GeoTests.cc @@ -2620,11 +2620,13 @@ void TwoBoxesGeoTest::test_tangent() const { SCOPED_TRACE("trying to cross"); SHOULD_FAIL_WHEN(geo.cross_boundary(), - test_->geometry_type() == "VecGeom"); + test_->geometry_type() == "VecGeom" + && !CELERITAS_VECGEOM_SURFACE); EXPECT_TRUE(geo.is_on_boundary()); - if (test_->geometry_type() == "Geant4") + if (test_->geometry_type() == "Geant4" + || test_->geometry_type() == "VecGeom") { - // FIXME: Geant4 changes volumes :( + // FIXME: should reentrant volume changes be handled by geo track? EXPECT_EQ("world", test_->volume_name(geo)); GTEST_SKIP() << "Unexpected boundary crossing"; } From a0b5d54e5d9d45f47e5d602b9be0b949b5581ed3 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 15:03:55 -0500 Subject: [PATCH 19/35] VG solid 1 passes field propagator --- test/celeritas/field/FieldPropagator.test.cc | 139 ++----------------- 1 file changed, 11 insertions(+), 128 deletions(-) diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index f8d6e1c306..876559a689 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -1326,7 +1326,6 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) auto calc_radius = [&geo]() { return std::hypot(geo.pos()[0], geo.pos()[1]); }; - bool successful_reentry = false; { auto particle = this->make_particle_view( pdg::electron(), MevEnergy{3.27089632881079409e-02}); @@ -1352,23 +1351,19 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) try { geo.cross_boundary(); - successful_reentry = (this->volume_name(geo) == "em_calorimeter"); } - catch (RuntimeError const& e) + catch (CheckedGeoError const& e) { - if (geo.failed()) + if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_VECGEOM) { - // VecGeom can set the 'failed' flag, which CheckedGeoTrackView - // turns into an error. Ignore that error for now and say it's - // a failure - successful_reentry = false; + GTEST_SKIP() << e.what(); } else { - // Something else happened... - throw; + FAIL() << e.what(); } } + if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_ORANGE) { // ORANGE should successfully reenter. However, under certain @@ -1376,127 +1371,15 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(vecgeom_failure)) // so we don't test in all cases. EXPECT_EQ("em_calorimeter", this->volume_name(geo)); } - else if (!geo.failed()) - { - EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; - } - - if (!successful_reentry) + else { // Note that this is expected behavior in Geant4 and VecGeom, as // it is assumed that the track will actually change volumes at the // boundary (tangent tracks are not the norm and maybe not properly // handled). - CELER_LOG(warning) << "Reentry failed for " << cmake::core_geo - << " geometry: post-propagation volume is " - << this->volume_name(geo); - if (!scoped_log_.empty()) - { - CELER_LOG(warning) << scoped_log_; - } - - if (using_solids_vg && CELERITAS_VECGEOM_VERSION < 0x020000) - { - EXPECT_EQ("em_calorimeter", this->volume_name(geo)); - } - else - { - EXPECT_EQ("si_tracker", this->volume_name(geo)); - } - - // Interestingly, VecGeom surf and solid models see that surface - // slightly differently. Only surface model thinks the surface - // was actually crossed, therefore the next step will find distinct - // results - Propagation result; - EXPECT_NO_THROW(result = geo.find_next_step(1)); - EXPECT_LT(result.distance, 2e-8); - - if (result.distance < 1e-6 && result.boundary) - { - geo.move_to_boundary(); - geo.cross_boundary(); // back into em_calorimeter - } - - // then they are back to agreement - EXPECT_NO_THROW(result = geo.find_next_step(1)); - EXPECT_EQ(result.distance, 1); - EXPECT_FALSE(result.boundary); - EXPECT_TRUE(geo.is_on_boundary()); - EXPECT_EQ("em_calorimeter", this->volume_name(geo)); - EXPECT_SOFT_EQ(125.00000000000001, calc_radius()); - } - else - { - CELER_LOG(debug) << "Reentry succeeded: " << scoped_log_; - } - } - { - ScopedLogStorer scoped_log_{&celeritas::self_logger()}; - auto particle = this->make_particle_view( - pdg::electron(), MevEnergy{3.25917780979408864e-02}); - auto integrate = make_mag_field_integrator( - field, particle.charge()); - auto propagate - = make_field_propagator(integrate, driver_options, particle, geo); - - Propagation result; - // This absurdly long step is because in the "failed" case the - // track thinks it's in the world volume (nearly vacuum) - EXPECT_NO_THROW(result = propagate(2.12621374950874703e+21)); - - if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_GEANT4 - && result.boundary != geo.is_on_boundary()) - { - // FIXME: see #882 - GTEST_SKIP() << "The current fix fails with the Geant4 navigator"; - } - - EXPECT_EQ(result.boundary, geo.is_on_boundary()); - EXPECT_SOFT_NEAR(125, calc_radius(), 1e-2); - if (successful_reentry) - { - // ORANGE and *sometimes* vecgeom/geant4: extremely long - // propagation stopped by substep countdown - EXPECT_FALSE(result.boundary); - EXPECT_TRUE(result.looping); - EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; - - EXPECT_SOFT_EQ(12.02714054426572, result.distance); - EXPECT_EQ("em_calorimeter", this->volume_name(geo)); - EXPECT_EQ(573, integrate.count()); - EXPECT_TRUE(result.looping); - } - else - { - // Repeated substep bisection failed; particle is bumped - EXPECT_SOFT_NEAR(result.distance, 12.02714054426572, coarse_eps); - // Minor floating point differences could make this 98 or so - EXPECT_SOFT_NEAR( - real_type(573), real_type(integrate.count()), 0.05); - EXPECT_FALSE(result.boundary); // FIXME: should have reentered - EXPECT_TRUE(result.looping); // FIXME: looping?? - - if (scoped_log_.empty()) {} - else if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_GEANT4) - { - static char const* const expected_log_levels[] = {"error"}; - EXPECT_VEC_EQ(expected_log_levels, scoped_log_.levels()) - << scoped_log_; - } - else if (using_solids_vg) - { - static char const* const expected_log_messages[] = { - R"(Moved internally from boundary but safety didn't increase: volume 6 from {123.3, -20.82, -40.83} to {123.3, -20.82, -40.83} (distance: 1.000e-6))", - }; - EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()); - static char const* const expected_log_levels[] = {"warning"}; - EXPECT_VEC_EQ(expected_log_levels, scoped_log_.levels()); - } - else - { - ADD_FAILURE() << "Logged warning/error:" << scoped_log_; - } + // FIXME: see GeoTests: TwoBoxesGeoTest::test_tangent + EXPECT_EQ("si_tracker", this->volume_name(geo)); + GTEST_SKIP(); } } } @@ -1543,7 +1426,7 @@ TEST_F(CmseTest, coarse) { result = propagate(radius); } - catch (RuntimeError const& e) + catch (CheckedGeoError const& e) { CELER_LOG(error) << e.what(); failed = true; @@ -1591,7 +1474,7 @@ TEST_F(CmseTest, coarse) expected_num_integration[1] = 1666; static char const* const expected_log_messages[] = { R"(Moved internally from boundary but safety didn't increase: volume 18 from {10.32, -6.565, 796.9} to {10.32, -6.565, 796.9} (distance: 1.000e-4))", - R"(Failed to find next step at {10.32, -6.565, 796.9} along {0.6896, -0.1485, 0.7088})"}; + R"(Failed to find next step at {10.32, -6.565, 796.9} cm along {0.6896, -0.1485, 0.7088}: computed step is 0 cm)"}; EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()) << scoped_log_; } From 91d16a9d01e6148a596b55167397503daed527d0 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 15:59:27 -0500 Subject: [PATCH 20/35] OK now I guess it passes for real --- test/celeritas/CMakeLists.txt | 7 +- test/celeritas/field/FieldPropagator.test.cc | 131 +++++++++++++++---- 2 files changed, 106 insertions(+), 32 deletions(-) diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index 5e9a56ef45..392073ae9f 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -215,7 +215,8 @@ endif() #-----------------------------------------------------------------------------# # Field -if(CELERITAS_USE_VecGeom) +if(CELERITAS_CORE_GEO STREQUAL "VecGeom" + AND (CELERITAS_USE_CUDA OR VecGeom_VERSION VERSION_GREATER_EQUAL "2.0")) # VecGeom surface cannot rebuild the geometry set(_fieldprop_filter FILTER @@ -224,10 +225,6 @@ if(CELERITAS_USE_VecGeom) "SimpleCms*" "Cmse*" ) -else() - # Other geometries should run everything (and some tests e.g. CMSE are - # disabled) - set(_fieldprop_filter) endif() celeritas_add_device_test(field/Fields diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index 876559a689..33c6adbdde 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -813,8 +813,21 @@ TEST_F(TwoBoxesTest, electron_corner_hit) EXPECT_NORMAL_EQUIV((Real3{0, 1, 0}), geo.normal()); } - EXPECT_NO_THROW(geo.cross_boundary()); - EXPECT_EQ("world", this->volume_name(geo)); + if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) + { + /* + * unique volume instance is the same before and after + * at {-4.9930093111683, 5, 0} [cm] + * along {-0.99997563423471, 0.0069807547505593, 0}, + * [FAILED] [ON BOUNDARY] in world_PV/inner_PV + */ + EXPECT_THROW(geo.cross_boundary(), CheckedGeoError); + } + else + { + EXPECT_NO_THROW(geo.cross_boundary()); + EXPECT_EQ("world", this->volume_name(geo)); + } } { SCOPED_TRACE("Barely (correctly) misses y"); @@ -1032,16 +1045,34 @@ TEST_F(TwoBoxesTest, { SCOPED_TRACE(i); Propagation result; - EXPECT_NO_THROW(result = propagate(radius * dtheta)); - if (result.boundary) + result.distance = 0; + result.boundary = false; + + try { - geo.cross_boundary(); + if (!geo.failed()) + { + result = propagate(radius * dtheta); + } + } + catch (CheckedGeoError const& e) + { + CELER_LOG(error) << e.what(); } - boundary.push_back(result.boundary); + if (result.distance > 0) + { + boundary.push_back(result.boundary); + volumes.push_back(this->volume_name(geo)); + } + else + { + // Error sentinel + boundary.push_back(-1); + volumes.push_back("[FAILURE]"); + } distances.push_back(result.distance); substeps.push_back(integrate.count()); - volumes.push_back(this->volume_name(geo)); integrate.reset_count(); } } @@ -1073,10 +1104,54 @@ TEST_F(TwoBoxesTest, "world", }; - EXPECT_VEC_EQ(expected_boundary, boundary); - EXPECT_VEC_NEAR(expected_distances, distances, real_type{.1} * coarse_eps); - EXPECT_VEC_EQ(expected_substeps, substeps); - EXPECT_VEC_EQ(expected_volumes, volumes); + /* + * Failed to find next step: computed step is 0 cm + * at {0.0014111984735586, 5, 0} [cm] + * along {-0.96863764097147, -0.24847760561715, 0}, + * [FAILED] [ON BOUNDARY] in world_PV + */ + if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) + { + expected_boundary = {1, 0, 1, 0, 1, -1, 1, -1, -1, -1}; + expected_distances[1] = 0.00785398163397448; + expected_distances[3] = 0.00448798950512828; + expected_distances[5] = 0; + expected_distances[6] = 9.99996586224189e-09; + expected_distances[7] = 0; + expected_substeps[1] = 1; + expected_substeps[3] = 1; + expected_substeps[9] = 0; + expected_volumes = { + "inner", + "inner", + "inner", + "inner", + "inner", + "[FAILURE]", + "inner", + "[FAILURE]", + "[FAILURE]", + "[FAILURE]", + }; + } + else if (using_solids_vg) + { + expected_boundary[1] = 0; + expected_boundary[3] = 0; + expected_distances[1] = 0.00785398163397448; + expected_distances[3] = 0.00448798950512828; + expected_distances[5] = 1e-6; + expected_substeps[1] = 1; + expected_substeps[3] = 1; + expected_substeps[5] = 4; + expected_volumes.assign(expected_volumes.size(), "inner"); + } + + EXPECT_VEC_EQ(expected_boundary, boundary) << repr(boundary); + EXPECT_VEC_NEAR(expected_distances, distances, real_type{.1} * coarse_eps) + << repr(distances); + EXPECT_VEC_EQ(expected_substeps, substeps) << repr(substeps); + EXPECT_VEC_EQ(expected_volumes, volumes) << repr(volumes); } // Heuristic test: plotting points with finer propagation distance show a track @@ -1419,7 +1494,7 @@ TEST_F(CmseTest, coarse) int step_count = 0; int boundary_count = 0; int const max_steps = 10000; - while (!geo.is_outside() && step_count++ < max_steps) + while (!geo.is_outside() && !geo.failed() && step_count++ < max_steps) { Propagation result; try @@ -1429,12 +1504,26 @@ TEST_F(CmseTest, coarse) catch (CheckedGeoError const& e) { CELER_LOG(error) << e.what(); - failed = true; break; } if (result.boundary) { - geo.cross_boundary(); + try + { + geo.cross_boundary(); + } + catch (CheckedGeoError const& e) + { + if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_VECGEOM) + { + CELER_LOG(error) << e.details().what; + break; + } + else + { + FAIL() << e.what(); + } + } ++boundary_count; } } @@ -1466,7 +1555,7 @@ TEST_F(CmseTest, coarse) expected_num_integration = {80659, 1670, 1956, 1092}; EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; } - else if (using_solids_vg && failed) + else if (using_solids_vg) { expected_num_boundary[1] = 37; expected_num_step[1] = 179; @@ -1478,18 +1567,6 @@ TEST_F(CmseTest, coarse) EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()) << scoped_log_; } - else if (using_solids_vg) - { - // Bumped (platform-dependent!): counts change a bit - expected_num_boundary[1] = 101; - expected_num_step[1] = 6462; - expected_num_intercept[1] = 19551; - expected_num_integration[1] = 58282; - static char const* const expected_log_messages[] = { - R"(Moved internally from boundary but safety didn't increase: volume 18 from {10.32, -6.565, 796.9} to {10.32, -6.565, 796.9} (distance: 1.000e-4))"}; - EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()) - << scoped_log_; - } EXPECT_VEC_EQ(expected_num_boundary, num_boundary); EXPECT_VEC_EQ(expected_num_step, num_step); EXPECT_VEC_EQ(expected_num_intercept, num_intercept); From 0c0c3aa02fff8e1a0cc76103fb6e8d23d0f638c4 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 16:09:10 -0500 Subject: [PATCH 21/35] Field with vg2 solids passes --- test/celeritas/field/FieldPropagator.test.cc | 42 +++++++++++++------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index 33c6adbdde..01723f2972 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -1478,7 +1478,6 @@ TEST_F(CmseTest, coarse) std::vector num_integration; ScopedLogStorer scoped_log_{&celeritas::self_logger()}; - bool failed{false}; for (real_type radius : {5, 10, 20, 50}) { @@ -1538,8 +1537,9 @@ TEST_F(CmseTest, coarse) std::vector expected_num_step = {10001, 6450, 3236, 1303}; std::vector expected_num_intercept = {30419, 19521, 16170, 9956}; std::vector expected_num_integration = {80659, 58204, 41914, 26114}; + std::vector expected_log_messages; - if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_GEANT4 && failed) + if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_GEANT4) { // FIXME: this happens because of incorrect momentum update expected_num_boundary = {134, 37, 60, 40}; @@ -1547,13 +1547,17 @@ TEST_F(CmseTest, coarse) expected_num_intercept = {30419, 615, 16170, 9956}; expected_num_integration = {80659, 1670, 41914, 26114}; } - else if (using_surface_vg) + else if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) { - expected_num_boundary = {134, 37, 43, 16}; - expected_num_step = {10001, 179, 160, 63}; - expected_num_intercept = {30419, 615, 790, 414}; - expected_num_integration = {80659, 1670, 1956, 1092}; - EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; + expected_num_boundary = {1, 1, 2, 3}; + expected_num_step = {5, 3, 48, 29}; + expected_num_intercept = {17, 13, 236, 197}; + expected_num_integration = {40, 31, 612, 529}; + expected_log_messages = { + R"(Failed to cross boundary: unique volume instance is the same before and after at {4.033, -2.557, -286.0} cm)", + R"(Failed to cross boundary: unique volume instance is the same before and after at {6.349, 2.564, -280.4} cm)", + R"(Failed to cross boundary: unique volume instance is the same before and after at {19.03, 12.38, 351.4} cm)", + R"(Failed to cross boundary: unique volume instance is the same before and after at {51.18, 32.69, 657.9} cm)"}; } else if (using_solids_vg) { @@ -1561,16 +1565,24 @@ TEST_F(CmseTest, coarse) expected_num_step[1] = 179; expected_num_intercept[1] = 614; expected_num_integration[1] = 1666; - static char const* const expected_log_messages[] = { + expected_log_messages = { R"(Moved internally from boundary but safety didn't increase: volume 18 from {10.32, -6.565, 796.9} to {10.32, -6.565, 796.9} (distance: 1.000e-4))", R"(Failed to find next step at {10.32, -6.565, 796.9} cm along {0.6896, -0.1485, 0.7088}: computed step is 0 cm)"}; - EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()) - << scoped_log_; } - EXPECT_VEC_EQ(expected_num_boundary, num_boundary); - EXPECT_VEC_EQ(expected_num_step, num_step); - EXPECT_VEC_EQ(expected_num_intercept, num_intercept); - EXPECT_VEC_EQ(expected_num_integration, num_integration); + else if (using_surface_vg) + { + expected_num_boundary = {134, 37, 43, 16}; + expected_num_step = {10001, 179, 160, 63}; + expected_num_intercept = {30419, 615, 790, 414}; + expected_num_integration = {80659, 1670, 1956, 1092}; + EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; + } + EXPECT_VEC_EQ(expected_num_boundary, num_boundary) << repr(num_boundary); + EXPECT_VEC_EQ(expected_num_step, num_step) << repr(num_step); + EXPECT_VEC_EQ(expected_num_intercept, num_intercept) << repr(num_intercept); + EXPECT_VEC_EQ(expected_num_integration, num_integration) + << repr(num_integration); + EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()) << scoped_log_; } //---------------------------------------------------------------------------// From ccea451469e046226ac660cf9eb808b3f877ed84 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 16:14:28 -0500 Subject: [PATCH 22/35] Surface vg passes --- test/celeritas/field/FieldPropagator.test.cc | 24 ++++++++++++++------ 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index 01723f2972..2b8cec62d0 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -1110,7 +1110,7 @@ TEST_F(TwoBoxesTest, * along {-0.96863764097147, -0.24847760561715, 0}, * [FAILED] [ON BOUNDARY] in world_PV */ - if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) + if (CELERITAS_VECGEOM_VERSION >= 0x020000) { expected_boundary = {1, 0, 1, 0, 1, -1, 1, -1, -1, -1}; expected_distances[1] = 0.00785398163397448; @@ -1133,6 +1133,14 @@ TEST_F(TwoBoxesTest, "[FAILURE]", "[FAILURE]", }; + + if (using_surface_vg) + { + expected_boundary[8] = 1; + expected_substeps[8] = 1; + expected_substeps[9] = 1; + expected_volumes[8] = "inner"; + } } else if (using_solids_vg) { @@ -1567,15 +1575,17 @@ TEST_F(CmseTest, coarse) expected_num_integration[1] = 1666; expected_log_messages = { R"(Moved internally from boundary but safety didn't increase: volume 18 from {10.32, -6.565, 796.9} to {10.32, -6.565, 796.9} (distance: 1.000e-4))", - R"(Failed to find next step at {10.32, -6.565, 796.9} cm along {0.6896, -0.1485, 0.7088}: computed step is 0 cm)"}; + R"(Failed to find next step at {10.32, -6.565, 796.9} cm along {0.6896, -0.1485, 0.7088}: computed step is 0 cm)", + }; } else if (using_surface_vg) { - expected_num_boundary = {134, 37, 43, 16}; - expected_num_step = {10001, 179, 160, 63}; - expected_num_intercept = {30419, 615, 790, 414}; - expected_num_integration = {80659, 1670, 1956, 1092}; - EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; + expected_num_boundary = {106, 37, 43, 16}; + expected_num_step = {670, 179, 160, 63}; + expected_num_intercept = {2371, 615, 790, 414}; + expected_num_integration = {5931, 1670, 1956, 1092}; + expected_log_messages = { + R"(Failed to find next step at {9.956, 0.4182, 1833.} cm along {-0.2459, 0.6576, 0.7121}: computed step is -7.707e-13 cm)"}; } EXPECT_VEC_EQ(expected_num_boundary, num_boundary) << repr(num_boundary); EXPECT_VEC_EQ(expected_num_step, num_step) << repr(num_step); From c11e6c97da7d607496a5a2e5881c0c0642e27a3e Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 29 Nov 2025 18:35:57 -0500 Subject: [PATCH 23/35] Send 'failed' state back to propagation applier --- src/celeritas/alongstep/detail/PropagationApplier.hh | 5 +++++ src/celeritas/field/FieldPropagator.hh | 8 +++++++- src/celeritas/field/LinearPropagator.hh | 7 ++++++- src/geocel/Types.hh | 5 +++++ 4 files changed, 23 insertions(+), 2 deletions(-) diff --git a/src/celeritas/alongstep/detail/PropagationApplier.hh b/src/celeritas/alongstep/detail/PropagationApplier.hh index c0bccb4a93..3236a91a5e 100644 --- a/src/celeritas/alongstep/detail/PropagationApplier.hh +++ b/src/celeritas/alongstep/detail/PropagationApplier.hh @@ -114,6 +114,11 @@ PropagationApplierBaseImpl::operator()(CoreTrackView& track) #endif auto propagate = make_propagator(track); p = propagate(sim.step_length()); + if (CELER_UNLIKELY(p.failed)) + { + track.apply_errored(); + return; + } tracks_can_loop = propagate.tracks_can_loop(); CELER_ASSERT(p.distance > 0); #if CELERITAS_DEBUG diff --git a/src/celeritas/field/FieldPropagator.hh b/src/celeritas/field/FieldPropagator.hh index 6ccece3eaf..4d9ffc1b87 100644 --- a/src/celeritas/field/FieldPropagator.hh +++ b/src/celeritas/field/FieldPropagator.hh @@ -197,7 +197,13 @@ FieldPropagator::operator()(real_type step) -> result_type real_type const update_length = substep.length * linear_step.distance / chord.length; - if (!linear_step.boundary) + if (CELER_UNLIKELY(geo_.failed())) + { + // Time to die + result.failed = true; + return result; + } + else if (!linear_step.boundary) { // No boundary intersection along the chord: accept substep // movement inside the current volume and reset the remaining diff --git a/src/celeritas/field/LinearPropagator.hh b/src/celeritas/field/LinearPropagator.hh index 4285ad1b93..5199a77045 100644 --- a/src/celeritas/field/LinearPropagator.hh +++ b/src/celeritas/field/LinearPropagator.hh @@ -6,6 +6,7 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/Macros.hh" #include "corecel/math/Algorithms.hh" #include "geocel/Types.hh" @@ -78,7 +79,11 @@ CELER_FUNCTION auto LinearPropagator::operator()(real_type dist) result_type result = geo_.find_next_step(dist); - if (result.boundary) + if (CELER_UNLIKELY(geo_.failed())) + { + result.failed = true; + } + else if (result.boundary) { geo_.move_to_boundary(); } diff --git a/src/geocel/Types.hh b/src/geocel/Types.hh index 8a65d4ddc0..0761521e7b 100644 --- a/src/geocel/Types.hh +++ b/src/geocel/Types.hh @@ -109,12 +109,17 @@ struct GeoTrackInitializer * * The boundary flag means that the geometry is step limiting, but the surface * crossing must be called externally. + * + * \todo Change to a status type, or have separate "propagation" and "next + * step" result types. Currently this is copied from the geometry state to the + * propagation result because the track view failure flags are ephemeral. */ struct Propagation { real_type distance{0}; //!< Distance traveled bool boundary{false}; //!< True if hit a boundary before given distance bool looping{false}; //!< True if track is looping in the field propagator + bool failed{false}; //!< True if a failure state occurred }; //---------------------------------------------------------------------------// From 26257388484c827c2dcd0de52018b8d4aa80ae60 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 9 Dec 2025 18:33:45 -0600 Subject: [PATCH 24/35] Clean up a few things --- src/geocel/vg/detail/SurfNavigator.hh | 2 -- test/celeritas/field/FieldPropagator.test.cc | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/geocel/vg/detail/SurfNavigator.hh b/src/geocel/vg/detail/SurfNavigator.hh index b48e7b8c67..14187d5328 100644 --- a/src/geocel/vg/detail/SurfNavigator.hh +++ b/src/geocel/vg/detail/SurfNavigator.hh @@ -45,8 +45,6 @@ class SurfNavigator using SurfData = vgbrep::SurfData; using NavState = detail::VgNavStateWrapper; - static constexpr vg_real_type kBoundaryPush = 10 * vecgeom::kTolerance; - /// @brief Locates the point in the geometry volume tree /// @param pvol_id Placed volume id to be checked first /// @param point Point to be checked, in the local frame of pvol diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index 2b8cec62d0..ac619cca69 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -1308,7 +1308,7 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(electron_stuck)) if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) { - // NOTE: vecgeom 2 solid model thinks r=30 + epsilon is *outside* the + // NOTE: vecgeom 2 solid model thinks r=30 + epsilon is *inside* the // 30cm radius cyl EXPECT_EQ("vacuum_tube", this->volume_name(geo)); GTEST_SKIP() << "VecGeom 2.x solid disagrees where the solid is"; From 09f1eaa4e9c98ea361a76173aefd78ed21047a37 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 11 Dec 2025 11:24:25 -0500 Subject: [PATCH 25/35] Add edge case for starting near boundary at tangent --- test/geocel/GenericGeoTestInterface.cc | 12 ++- test/geocel/GeoTests.cc | 102 +++++++++++++++++++++++++ test/testdetail/TestMacrosImpl.hh | 30 ++++---- 3 files changed, 125 insertions(+), 19 deletions(-) diff --git a/test/geocel/GenericGeoTestInterface.cc b/test/geocel/GenericGeoTestInterface.cc index aa81450270..c2891df0a9 100644 --- a/test/geocel/GenericGeoTestInterface.cc +++ b/test/geocel/GenericGeoTestInterface.cc @@ -90,6 +90,9 @@ auto GenericGeoTestInterface::track(Real3 const& pos, auto from_native_length = [scale = unit_length.value](auto&& v) { return v / scale; }; auto const tol = this->tracking_tol(); + EqualOr> soft_eq_dist{ + tol.distance, tol.distance * unit_length.value}; + SoftZero is_bump(10 * soft_eq_dist.abs()); while (!geo.is_outside()) { @@ -97,7 +100,7 @@ auto GenericGeoTestInterface::track(Real3 const& pos, Propagation next; GGTI_EXPECT_NO_THROW(next = geo.find_next_step()); - if (SoftZero{tol.distance}(next.distance)) + if (is_bump(next.distance)) { // Add the point to the bump list for (auto p : geo.pos()) @@ -123,7 +126,8 @@ auto GenericGeoTestInterface::track(Real3 const& pos, real_type const half_distance = next.distance / 2; GGTI_EXPECT_NO_THROW(geo.move_internal(half_distance)); GGTI_EXPECT_NO_THROW(next = geo.find_next_step()); - EXPECT_SOFT_NEAR(next.distance, half_distance, tol.distance); + + EXPECT_SOFT_NEAR(next.distance, half_distance, soft_eq_dist); real_type safety{0}; GGTI_EXPECT_NO_THROW(safety = geo.find_safety()); @@ -146,7 +150,7 @@ auto GenericGeoTestInterface::track(Real3 const& pos, result.volumes.back() += "/" + this->volume_name(geo); } GGTI_EXPECT_NO_THROW(next = geo.find_next_step()); - EXPECT_SOFT_NEAR(next.distance, half_distance, tol.distance) + EXPECT_SOFT_NEAR(next.distance, half_distance, soft_eq_dist) << "reinitialized distance mismatch at index " << result.volumes.size() - 1 << ": " << geo; } @@ -253,7 +257,7 @@ size_type GenericGeoTestInterface::num_track_slots() const GenericGeoTrackingTolerance GenericGeoTestInterface::tracking_tol() const { GenericGeoTrackingTolerance result; - result.distance = SoftEqual{}.rel(); + result.distance = 5 * SoftEqual{}.rel(); result.normal = celeritas::sqrt_tol(); result.safety = result.distance; return result; diff --git a/test/geocel/GeoTests.cc b/test/geocel/GeoTests.cc index d107f4d836..f6d68b2a11 100644 --- a/test/geocel/GeoTests.cc +++ b/test/geocel/GeoTests.cc @@ -1922,6 +1922,108 @@ void SimpleCmsGeoTest::test_trace() const auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } + + // Edge case from FieldPropagator test + constexpr Real3 pos{ + -2.43293925496543e+01, -1.75522265870979e+01, 2.80918346435833e+02}; + constexpr Real3 dir{ + 7.01343313647855e-01, -6.43327996599957e-01, 3.06996164784077e-01}; + { + SCOPED_TRACE("edge case"); + auto result = test_->track(pos, dir); + result.volume_instances.clear(); + GenericGeoTrackingResult ref; + ref.volumes = { + "vacuum_tube", + "si_tracker", + "em_calorimeter", + "had_calorimeter", + "sc_solenoid", + "fe_muon_chambers", + "world", + }; + ref.distances = { + 12.743906479688, + 121.29083284073, + 53.606590194296, + 106.03009509735, + 105.51658247496, + 342.0533777001, + 719.28354304297, + }; + ref.dot_normal = { + 0.19238060078892, + 0.19238060078893, + 0.92504797421279, + 0.93820197534225, + 0.94626349101788, + 0.94878522291859, + 0.950872075671, + }; + ref.halfway_safeties = { + 0.61931256376742, + 40.222931079394, + 24.914281270933, + 49.898407971867, + 49.967525043337, + 162.41892241026, + 252.232351765063, + }; + ref.bumps = { + -24.329392549654, + -17.552226587098, + 280.91834643583, + 4.8849813083507e-14, + }; + delete_orange_safety(*test_, ref, result); + + if (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_FLOAT) + { + ref.bumps[3] = 5.245209e-06f; + } + + auto tol = test_->tracking_tol(); + EXPECT_REF_NEAR(ref, result, tol); + } + { + SCOPED_TRACE("edge case with flipped direction"); + auto result = test_->track(pos, -dir); + result.volume_instances.clear(); + GenericGeoTrackingResult ref; + ref.volumes = { + "si_tracker", + "em_calorimeter", + "had_calorimeter", + "sc_solenoid", + "fe_muon_chambers", + "world", + }; + ref.distances = { + 121.29083284073, + 53.606590194296, + 106.03009509735, + 105.51658247496, + 342.0533777001, + 662.64803982069, + }; + ref.dot_normal = { + 0.92504797421279, + 0.93820197534225, + 0.94626349101788, + 0.94878522291859, + 0.950872075671, + }; + ref.halfway_safeties = { + 40.222931079394, + 24.914281270933, + 49.898407971867, + 49.967525043337, + 162.41892241026, + 232.37188601505, + }; + auto tol = test_->tracking_tol(); + EXPECT_REF_NEAR(ref, result, tol); + } } //---------------------------------------------------------------------------// diff --git a/test/testdetail/TestMacrosImpl.hh b/test/testdetail/TestMacrosImpl.hh index 89200d33e6..ec8f0230bb 100644 --- a/test/testdetail/TestMacrosImpl.hh +++ b/test/testdetail/TestMacrosImpl.hh @@ -100,10 +100,12 @@ struct SoftPrecisionType /*! * Get a soft comparison function from a \c SOFT_NEAR argument. * - * If a floating point value, it defaults to + * If a floating point value, it defaults to EqualOrSoftEqual using explicit + * casts to the given value type to avoid errors with mixed-precision + * arithmetic. */ template -constexpr auto soft_comparator(CT&& cmp_or_tol) +constexpr auto make_soft_comparator(CT&& cmp_or_tol) { if constexpr (std::is_floating_point_v>) { @@ -133,8 +135,6 @@ IsSoftEquivImpl(typename BinaryOp::value_type expected, char const* actual_expr, BinaryOp comp) { - using value_type = typename BinaryOp::value_type; - if (comp(expected, actual)) { return ::testing::AssertionSuccess(); @@ -144,19 +144,19 @@ IsSoftEquivImpl(typename BinaryOp::value_type expected, ::testing::AssertionResult result = ::testing::AssertionFailure(); result << "Value of: " << actual_expr << "\n Actual: " << actual - << "\nExpected: " << expected_expr << "\nWhich is: " << expected - << '\n'; + << "\nExpected: " << expected_expr << "\nWhich is: " << expected; - if (SoftZero{comp.abs()}(expected)) + if (std::fabs(actual - expected) > comp.abs()) { - // Avoid divide by zero errors - result << "(Absolute error " << actual - expected - << " exceeds tolerance " << comp.abs() << ")"; + result << '\n' + << "- absolute error " << actual - expected + << " exceeds tolerance " << comp.abs(); } - else + if (expected != 0 + && std::fabs(actual - expected) > std::fabs(expected) * comp.rel()) { - result << "(Relative error " << (actual - expected) / expected - << " exceeds tolerance " << comp.rel() << ")"; + result << "\n- relative error " << actual / expected - 1 + << " exceeds tolerance " << comp.rel(); } return result; } @@ -214,7 +214,7 @@ template expected_expr, static_cast(actual), actual_expr, - soft_comparator(std::forward(cmp_or_tol))); + make_soft_comparator(std::forward(cmp_or_tol))); } //---------------------------------------------------------------------------// @@ -714,7 +714,7 @@ template expected_expr, actual, actual_expr, - soft_comparator(std::forward(cmp_or_tol))); + make_soft_comparator(std::forward(cmp_or_tol))); } //---------------------------------------------------------------------------// From b4f2a5ed151be27ff321f2efecaeda66837916a8 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 11 Dec 2025 11:35:07 -0500 Subject: [PATCH 26/35] Warn when starting on a boundary --- test/geocel/CheckedGeoTrackView.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/geocel/CheckedGeoTrackView.cc b/test/geocel/CheckedGeoTrackView.cc index e336a1d96c..12316403e3 100644 --- a/test/geocel/CheckedGeoTrackView.cc +++ b/test/geocel/CheckedGeoTrackView.cc @@ -134,6 +134,10 @@ CheckedGeoTrackView::operator=(GeoTrackInitializer const& init) *t_ = init; CGTV_VALIDATE_NOT_FAILED(*this, "initialization"); CGTV_VALIDATE(*this, !t_->is_outside(), << "initialized outside"); + if (t_->is_on_boundary()) + { + CELER_LOG_LOCAL(warning) << "Started on a boundary: " << *this; + } next_boundary_.reset(); return *this; } From bca596250c141a5f08c7871813315f1457e6a790 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 11 Dec 2025 11:36:25 -0500 Subject: [PATCH 27/35] Update geotest/field for vg2 solid init-on-boundary failure --- test/celeritas/field/FieldPropagator.test.cc | 15 ++++++++++----- test/geocel/GeoTests.cc | 14 ++++++++++++++ 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index ac619cca69..9953fdb87b 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -1306,14 +1306,13 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(electron_stuck)) = [&geo]() { return std::hypot(geo.pos()[0], geo.pos()[1]); }; EXPECT_SOFT_EQ(30.000000000000011, calc_radius()); - if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) + if (this->volume_name(geo) == "vacuum_tube") { // NOTE: vecgeom 2 solid model thinks r=30 + epsilon is *inside* the - // 30cm radius cyl - EXPECT_EQ("vacuum_tube", this->volume_name(geo)); - GTEST_SKIP() << "VecGeom 2.x solid disagrees where the solid is"; + // 30cm radius cyl because it's "on the surface" and on a boundary + EXPECT_TRUE(geo.is_on_boundary()); } - EXPECT_EQ("si_tracker", this->volume_name(geo)); + else { auto integrate = make_mag_field_integrator( field, particle.charge()); @@ -1369,6 +1368,12 @@ TEST_F(SimpleCmsTest, TEST_IF_CELERITAS_DOUBLE(electron_stuck)) FAIL() << e.what(); } + if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) + { + GTEST_SKIP() << "FIXME: VecGeom solid method didn't start in the " + "right place"; + } + EXPECT_EQ(result.boundary, geo.is_on_boundary()); EXPECT_SOFT_NEAR( double{30}, static_cast(integrate.count()), 0.2); diff --git a/test/geocel/GeoTests.cc b/test/geocel/GeoTests.cc index f6d68b2a11..1a80de6432 100644 --- a/test/geocel/GeoTests.cc +++ b/test/geocel/GeoTests.cc @@ -1982,6 +1982,12 @@ void SimpleCmsGeoTest::test_trace() const ref.bumps[3] = 5.245209e-06f; } + if (using_solids_vg && vecgeom_version >= Version{2, 0}) + { + // Track starts on boundary, avoiding bump + ref.bumps.clear(); + } + auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } @@ -2021,6 +2027,14 @@ void SimpleCmsGeoTest::test_trace() const 162.41892241026, 232.37188601505, }; + + if (using_solids_vg && vecgeom_version >= Version{2, 0}) + { + // Track fails immediately + ref = {}; + ref.fail(); + } + auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } From 67b24b4c6371e368cb054461df20c070e9754346 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 11 Dec 2025 11:51:38 -0500 Subject: [PATCH 28/35] VG 2 isn't "reference CI" build --- test/celeritas/GeantTestBase.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/celeritas/GeantTestBase.cc b/test/celeritas/GeantTestBase.cc index ff44f5b7ad..226feb981c 100644 --- a/test/celeritas/GeantTestBase.cc +++ b/test/celeritas/GeantTestBase.cc @@ -58,6 +58,12 @@ bool GeantTestBase::is_ci_build() // Config options are different return false; } + // Skip if VG 2.0 + if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_VECGEOM + && CELERITAS_VECGEOM_VERSION > 0x020000) + { + return false; + } // Check clhep/g4 versions auto clhep = Version::from_string(cmake::clhep_version); auto g4 = Version::from_string(cmake::geant4_version); From 59c436def25d77012db9326279499aeb317f5f4a Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 11 Dec 2025 11:57:08 -0500 Subject: [PATCH 29/35] Allow slight safety overestimates in geo track view --- test/geocel/CheckedGeoTrackView.cc | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/test/geocel/CheckedGeoTrackView.cc b/test/geocel/CheckedGeoTrackView.cc index 12316403e3..237eb1b5f8 100644 --- a/test/geocel/CheckedGeoTrackView.cc +++ b/test/geocel/CheckedGeoTrackView.cc @@ -271,12 +271,21 @@ Propagation CheckedGeoTrackView::find_next_step(real_type distance) && !started_on_boundary) { real_type safety = t_->find_safety(distance); - CGTV_VALIDATE(*this, - safety <= result.distance, - << "safety " << repr(safety) - << " exceeds actual distance " << repr(result.distance) - << " to boundary at " << t_->pos() << " in " - << t_->impl_volume_id().get()); + if (!(safety <= result.distance)) + { + auto const& units = this->unit_length(); + + CELER_LOG_LOCAL(warning) + << "Calculated safety " << repr(safety) + << " exceeds actual distance " << repr(result.distance) + << " to boundary at " << t_->pos() << " by " + << repr((safety - result.distance) / units.value) << " [" + << units.label << "]: " << *this; + CGTV_VALIDATE(*this, + safety <= 1.1 * result.distance, + << "calculated safety " << repr(safety) + << " is much too large"); + } } if (result.distance == 0) { From fe5a5397d575ebdabbecb83127b2d4ab5a2d7025 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 12 Dec 2025 10:45:22 -0500 Subject: [PATCH 30/35] Try bumping when initializing on boundary --- src/geocel/vg/VecgeomTrackView.hh | 46 +++++++++++++++++++++++-- src/geocel/vg/detail/SolidsNavigator.hh | 9 ++++- 2 files changed, 51 insertions(+), 4 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 94d7012136..6b1ffb72b3 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -17,6 +17,7 @@ #include "corecel/Macros.hh" #include "corecel/Types.hh" #include "corecel/cont/Span.hh" +#include "corecel/io/EnumStringMapper.hh" #include "corecel/math/ArraySoftUnit.hh" #include "corecel/math/ArrayUtils.hh" #include "corecel/sys/ThreadId.hh" @@ -309,12 +310,51 @@ VecgeomTrackView::operator=(Initializer_t const& init) Navigator::LocatePointIn( world, detail::to_vector(pos_), vgstate_, contains_point); + if constexpr (!CELERITAS_VECGEOM_SURFACE + && CELERITAS_VECGEOM_VERSION >= 0x020000) + { + // VecGeom 2 solids can start "on boundary" :( + // so it'll end up on the wrong side of an object + // Bump backward and try again + if (CELER_UNLIKELY(vgstate_.IsOnBoundary())) + { +#if !CELER_DEVICE_COMPILE + CELER_LOG_LOCAL(info) + << "Bumping geometry state on boundary backward by " + << relocate_bump_ << " at " << repr(pos_) + << lengthunits::label; +#endif + vgstate_.Clear(); + VgReal3 bumped_pos; + for (auto i : range(3)) + { + bumped_pos[i] = fma(-relocate_bump_, dir_[i], pos_[i]); + } + Navigator::LocatePointIn( + world, bumped_pos, vgstate_, contains_point); +#if !CELER_DEVICE_COMPILE + if (CELER_UNLIKELY(vgstate_.IsOnBoundary())) + { + CELER_LOG_LOCAL(warning) + << "Initialized geometry state on boundary at " + << repr(bumped_pos) << lengthunits::label; + } + else + { + CELER_LOG_LOCAL(info) + << "Reinitialized geometry state at " << repr(bumped_pos) + << lengthunits::label; + // vgstate_.SetBoundaryState(true); + } +#endif + } + } + if (CELER_UNLIKELY(vgstate_.IsOutside())) { #if !CELER_DEVICE_COMPILE - auto msg = CELER_LOG_LOCAL(error); - msg << "Failed to initialize geometry state at " << repr(pos_) - << lengthunits::label; + CELER_LOG_LOCAL(error) << "Failed to initialize geometry state at " + << repr(pos_) << lengthunits::label; #endif failed_ = true; } diff --git a/src/geocel/vg/detail/SolidsNavigator.hh b/src/geocel/vg/detail/SolidsNavigator.hh index ab1fd31475..7c6683ff8a 100644 --- a/src/geocel/vg/detail/SolidsNavigator.hh +++ b/src/geocel/vg/detail/SolidsNavigator.hh @@ -62,8 +62,15 @@ class SolidsNavigator NavState& out_state) { ScopedVgNavState temp_state{out_state}; + // FIXME: avoid pushing; and at least use a relative scale + static constexpr vg_real_type boundary_push = 10 * vecgeom::kTolerance; return NavImpl::ComputeStepAndNextVolume( - pos, dir, step_limit, in_state, temp_state); + pos, + dir, + step_limit, + in_state, + temp_state, + in_state.IsOnBoundary() ? boundary_push : 0); } //-----------------------------------------------------------------------// From 01b03d165abc8466f743185e821b85c898bb25e5 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 12 Dec 2025 10:45:28 -0500 Subject: [PATCH 31/35] Revert "Try bumping when initializing on boundary" This reverts commit fe5a5397d575ebdabbecb83127b2d4ab5a2d7025. --- src/geocel/vg/VecgeomTrackView.hh | 46 ++----------------------- src/geocel/vg/detail/SolidsNavigator.hh | 9 +---- 2 files changed, 4 insertions(+), 51 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 6b1ffb72b3..94d7012136 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -17,7 +17,6 @@ #include "corecel/Macros.hh" #include "corecel/Types.hh" #include "corecel/cont/Span.hh" -#include "corecel/io/EnumStringMapper.hh" #include "corecel/math/ArraySoftUnit.hh" #include "corecel/math/ArrayUtils.hh" #include "corecel/sys/ThreadId.hh" @@ -310,51 +309,12 @@ VecgeomTrackView::operator=(Initializer_t const& init) Navigator::LocatePointIn( world, detail::to_vector(pos_), vgstate_, contains_point); - if constexpr (!CELERITAS_VECGEOM_SURFACE - && CELERITAS_VECGEOM_VERSION >= 0x020000) - { - // VecGeom 2 solids can start "on boundary" :( - // so it'll end up on the wrong side of an object - // Bump backward and try again - if (CELER_UNLIKELY(vgstate_.IsOnBoundary())) - { -#if !CELER_DEVICE_COMPILE - CELER_LOG_LOCAL(info) - << "Bumping geometry state on boundary backward by " - << relocate_bump_ << " at " << repr(pos_) - << lengthunits::label; -#endif - vgstate_.Clear(); - VgReal3 bumped_pos; - for (auto i : range(3)) - { - bumped_pos[i] = fma(-relocate_bump_, dir_[i], pos_[i]); - } - Navigator::LocatePointIn( - world, bumped_pos, vgstate_, contains_point); -#if !CELER_DEVICE_COMPILE - if (CELER_UNLIKELY(vgstate_.IsOnBoundary())) - { - CELER_LOG_LOCAL(warning) - << "Initialized geometry state on boundary at " - << repr(bumped_pos) << lengthunits::label; - } - else - { - CELER_LOG_LOCAL(info) - << "Reinitialized geometry state at " << repr(bumped_pos) - << lengthunits::label; - // vgstate_.SetBoundaryState(true); - } -#endif - } - } - if (CELER_UNLIKELY(vgstate_.IsOutside())) { #if !CELER_DEVICE_COMPILE - CELER_LOG_LOCAL(error) << "Failed to initialize geometry state at " - << repr(pos_) << lengthunits::label; + auto msg = CELER_LOG_LOCAL(error); + msg << "Failed to initialize geometry state at " << repr(pos_) + << lengthunits::label; #endif failed_ = true; } diff --git a/src/geocel/vg/detail/SolidsNavigator.hh b/src/geocel/vg/detail/SolidsNavigator.hh index 7c6683ff8a..ab1fd31475 100644 --- a/src/geocel/vg/detail/SolidsNavigator.hh +++ b/src/geocel/vg/detail/SolidsNavigator.hh @@ -62,15 +62,8 @@ class SolidsNavigator NavState& out_state) { ScopedVgNavState temp_state{out_state}; - // FIXME: avoid pushing; and at least use a relative scale - static constexpr vg_real_type boundary_push = 10 * vecgeom::kTolerance; return NavImpl::ComputeStepAndNextVolume( - pos, - dir, - step_limit, - in_state, - temp_state, - in_state.IsOnBoundary() ? boundary_push : 0); + pos, dir, step_limit, in_state, temp_state); } //-----------------------------------------------------------------------// From cebffd1452abafcfc9fda8b78f3a098ddf4e9fd7 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 18 Dec 2025 12:31:46 -0500 Subject: [PATCH 32/35] Validate pointers from runtime in BVHNavigator --- src/geocel/vg/VecgeomParams.cc | 1 + src/geocel/vg/detail/VecgeomSetup.cu | 59 +++++++++++++++++++++++++++- src/geocel/vg/detail/VecgeomSetup.hh | 13 +++++- 3 files changed, 70 insertions(+), 3 deletions(-) diff --git a/src/geocel/vg/VecgeomParams.cc b/src/geocel/vg/VecgeomParams.cc index e436d36004..e9356d5239 100644 --- a/src/geocel/vg/VecgeomParams.cc +++ b/src/geocel/vg/VecgeomParams.cc @@ -837,6 +837,7 @@ void VecgeomParams::build_volume_tracking() check_bvh_device_pointers(); check_navindex_device_pointers(); + detail::check_other_device_pointers(); device_ownership_ = Ownership::value; } diff --git a/src/geocel/vg/detail/VecgeomSetup.cu b/src/geocel/vg/detail/VecgeomSetup.cu index 600019b8dc..471cb1e713 100644 --- a/src/geocel/vg/detail/VecgeomSetup.cu +++ b/src/geocel/vg/detail/VecgeomSetup.cu @@ -6,8 +6,13 @@ //---------------------------------------------------------------------------// #include "VecgeomSetup.hh" +#include #include +#if VECGEOM_VERSION >= 0x020000 +# include +#endif + #include "corecel/data/DeviceVector.hh" #if CELERITAS_VECGEOM_SURFACE @@ -39,7 +44,7 @@ struct BvhGetter pointer_type* dest{nullptr}; - CELER_FUNCTION void operator()(ThreadId tid) + __device__ void operator()(ThreadId tid) { CELER_EXPECT(tid == ThreadId{0}); *dest = vecgeom::cuda::BVHManager::GetBVH(0); @@ -62,6 +67,46 @@ struct NavIndexGetter } }; +//---------------------------------------------------------------------------// +//! Copy the logical volume pointer table +struct LogicalVolumesGetter +{ + using pointer_type = vecgeom::cuda::LogicalVolume const*; + static constexpr char const label[] = "logical-volumes"; + + pointer_type* dest{nullptr}; + + __device__ void operator()(ThreadId tid) + { + CELER_EXPECT(tid == ThreadId{0}); +#if VECGEOM_VERSION >= 0x020000 + *dest = vecgeom::globaldevicegeomdata::gDeviceLogicalVolumes; +#else + *dest = nullptr; +#endif + } +}; + +//---------------------------------------------------------------------------// +//! Copy the logical volume pointer table +struct PlacedVolumesGetter +{ + using pointer_type = vecgeom::cuda::VPlacedVolume const*; + static constexpr char const label[] = "placed-volumes"; + + pointer_type* dest{nullptr}; + + __device__ void operator()(ThreadId tid) + { + CELER_EXPECT(tid == ThreadId{0}); +#if VECGEOM_VERSION >= 0x020000 + *dest = vecgeom::globaldevicegeomdata::gCompactPlacedVolBuffer; +#else + *dest = nullptr; +#endif + } +}; + //---------------------------------------------------------------------------// //! Launch a kernel to copy a value from global memory template @@ -144,6 +189,18 @@ CudaPointers navindex_pointers_device() return result; } +void check_other_device_pointers() +{ + if constexpr (VECGEOM_VERSION < 0x020000) + return; + + CELER_VALIDATE(get_device_pointer() != nullptr, + << "failed to copy VG logical volumes to GPU"); + + CELER_VALIDATE(get_device_pointer() != nullptr, + << "failed to copy VG places volumes to GPU"); +} + #if CELER_VGNAV == CELER_VGNAV_TUPLE //---------------------------------------------------------------------------// /* diff --git a/src/geocel/vg/detail/VecgeomSetup.hh b/src/geocel/vg/detail/VecgeomSetup.hh index 60aed02bd2..cb278b3310 100644 --- a/src/geocel/vg/detail/VecgeomSetup.hh +++ b/src/geocel/vg/detail/VecgeomSetup.hh @@ -48,7 +48,11 @@ CudaPointers bvh_pointers_device(); //---------------------------------------------------------------------------// // Get pointers to the global nav index after setup, for consistency checking -CudaPointers navindex_pointers_device(); +CudaPointers navindex_pointers_device(); + +//---------------------------------------------------------------------------// +// Check consistency of BVH pointers for device data in vecgeom 2.x +void check_other_device_pointers(); //---------------------------------------------------------------------------// // Default-initialize navigation state because DeviceVector doesn't @@ -72,7 +76,12 @@ inline CudaPointers bvh_pointers_device() CELER_ASSERT_UNREACHABLE(); } -inline CudaPointers navindex_pointers_device() +inline CudaPointers navindex_pointers_device() +{ + CELER_ASSERT_UNREACHABLE(); +} + +inline void check_other_device_pointers() { CELER_ASSERT_UNREACHABLE(); } From c5e5463745458e44ce7efdb9dfb68cfe4039285c Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Wed, 4 Feb 2026 12:36:43 -0500 Subject: [PATCH 33/35] Update executor script --- scripts/cmake-presets/executor.json | 58 +++++++++++++++++------------ 1 file changed, 34 insertions(+), 24 deletions(-) diff --git a/scripts/cmake-presets/executor.json b/scripts/cmake-presets/executor.json index 0ef7e7c8b1..979172da86 100644 --- a/scripts/cmake-presets/executor.json +++ b/scripts/cmake-presets/executor.json @@ -14,27 +14,21 @@ "name": ".base", "hidden": true, "generator": "Ninja", - "inherits": [".spack-env"], "binaryDir": "${sourceDir}/build-${presetName}", "cacheVariables": { "BUILD_SHARED_LIBS": {"type": "BOOL", "value": "ON"}, "BUILD_TESTING": {"type": "BOOL", "value": "ON"}, - "CELERITAS_USE_covfie": {"type": "BOOL", "value": "ON"}, "CELERITAS_USE_CUDA": {"type": "BOOL", "value": "OFF"}, "CELERITAS_USE_DD4hep": {"type": "BOOL", "value": "OFF"}, "CELERITAS_USE_LArSoft": {"type": "BOOL", "value": "OFF"}, - "CELERITAS_USE_HepMC3": {"type": "BOOL", "value": "ON"}, "CELERITAS_USE_HIP": {"type": "BOOL", "value": "OFF"}, - "CELERITAS_USE_Geant4": {"type": "BOOL", "value": "ON"}, "CELERITAS_USE_MPI": {"type": "BOOL", "value": "OFF"}, "CELERITAS_USE_OpenMP": {"type": "BOOL", "value": "OFF"}, + "CELERITAS_USE_PNG": {"type": "BOOL", "value": "OFF"}, + "CELERITAS_USE_Python": {"type": "BOOL", "value": "ON"}, "CELERITAS_USE_ROOT": {"type": "BOOL", "value": "OFF"}, "CELERITAS_USE_VecGeom": {"type": "BOOL", "value": "OFF"}, "CELERITAS_BUILD_TESTS": {"type": "BOOL", "value": "ON"}, - "CELERITAS_BUILTIN_CLI11": {"type": "BOOL", "value": "OFF"}, - "CELERITAS_BUILTIN_G4VG": {"type": "BOOL", "value": "OFF"}, - "CELERITAS_BUILTIN_nlohmann_json": {"type": "BOOL", "value": "OFF"}, - "CELERITAS_BUILTIN_GTest": {"type": "BOOL", "value": "OFF"}, "CMAKE_C_COMPILER_LAUNCHER": null, "CMAKE_CXX_COMPILER_LAUNCHER": null, "CMAKE_EXPORT_COMPILE_COMMANDS": {"type": "BOOL", "value": "ON"}, @@ -47,10 +41,29 @@ "CMAKE_CXX_FLAGS": "-Wall -Wextra -Werror -Wno-error=deprecated -pedantic -fdiagnostics-color=always" } }, + { + "name": ".spackbase", + "hidden": true, + "generator": "Ninja", + "inherits": [".base", ".spack-env"], + "binaryDir": "${sourceDir}/build-${presetName}", + "cacheVariables": { + "CELERITAS_USE_covfie": {"type": "BOOL", "value": "ON"}, + "CELERITAS_USE_HepMC3": {"type": "BOOL", "value": "ON"}, + "CELERITAS_USE_Geant4": {"type": "BOOL", "value": "ON"}, + "CELERITAS_USE_PNG": {"type": "BOOL", "value": "ON"}, + "CELERITAS_USE_Python": {"type": "BOOL", "value": "ON"}, + "CELERITAS_BUILD_TESTS": {"type": "BOOL", "value": "ON"}, + "CELERITAS_BUILTIN_CLI11": {"type": "BOOL", "value": "OFF"}, + "CELERITAS_BUILTIN_G4VG": {"type": "BOOL", "value": "OFF"}, + "CELERITAS_BUILTIN_nlohmann_json": {"type": "BOOL", "value": "OFF"}, + "CELERITAS_BUILTIN_GTest": {"type": "BOOL", "value": "OFF"} + } + }, { "name": "iwyu", "displayName": "Include-what-you-use (see scripts/dev/run-iwyu.sh)", - "inherits": [".base"], + "inherits": [".spackbase"], "cacheVariables": { "CELERITAS_BUILD_TESTS": {"type": "BOOL", "value": "OFF"}, "CELERITAS_DEBUG": {"type": "BOOL", "value": "OFF"} @@ -58,7 +71,7 @@ }, { "name": "min", - "displayName": "Goldfinger minimal build", + "displayName": "Executor minimal build", "inherits": [".base", ".debug", "default"], "environment": { "CMAKE_PREFIX_PATH": "$env{SPACK_ROOT}/var/spack/environments/celeritas-min/.spack-env/view" @@ -72,7 +85,7 @@ { "name": "base", "displayName": "Goldfinger default options", - "inherits": [".ccache", ".base", ".debug", "default"], + "inherits": [".ccache", ".spackbase", ".debug", "default"], "binaryDir": "${sourceDir}/build", "cacheVariables": { "CELERITAS_BUILD_DOCS": {"type": "BOOL", "value": "ON"}, @@ -91,15 +104,10 @@ "CELERITAS_REAL_TYPE": "float" } }, - { - "name": "mini", - "displayName": "Minimal deps", - "inherits": ["minimal", ".base", ".debug"] - }, { "name": "si", "displayName": "With SI units", - "inherits": [".base", ".debug", "default"], + "inherits": [".spackbase", ".debug", "default"], "cacheVariables": { "CELERITAS_UNITS": "SI", "CELERITAS_CORE_GEO": "VecGeom", @@ -109,7 +117,7 @@ { "name": "vecgeom", "displayName": "With vecgeom", - "inherits": [".ccache", ".base", ".debug", "default"], + "inherits": [".ccache", ".spackbase", ".debug", "default"], "cacheVariables": { "CELERITAS_USE_VecGeom": {"type": "BOOL", "value": "ON"} } @@ -119,16 +127,18 @@ "displayName": "With local vecgeom", "inherits": [".ccache", ".base", ".debug", "default"], "environment": { - "CMAKE_PREFIX_PATH": "/opt/spack/opt/spack/tahoe/geant4/11.3.2/7oydgus:/opt/spack/opt/spack/tahoe/googletest/1.17.0/c332ssk:/opt/spack/opt/spack/tahoe/cli11/2.5.0/pmtt6hb:/opt/spack/opt/spack/tahoe/nlohmann-json/3.12.0/dkrlqn3:/opt/spack/opt/spack/tahoe/covfie/0.15.3/tlzukkw:/opt/spack/opt/spack/tahoe/zlib-ng/2.2.4/qbbzbte" + "CMAKE_PREFIX_PATH": "/opt/spack/opt/spack/tahoe/geant4/11.3.2/gykadvh:/opt/spack/opt/spack/tahoe/covfie/0.15.4/r5g2psb:/opt/spack/opt/spack/tahoe/nlohmann-json/3.12.0/7mlcviu:/opt/spack/opt/spack/tahoe/cli11/2.6.1/znuwv76:/opt/spack/opt/spack/tahoe/googletest/1.17.0/rmxgcak" }, "cacheVariables": { + "CELERITAS_BUILTIN_CLI11": {"type": "BOOL", "value": "OFF"}, "CELERITAS_BUILTIN_G4VG": {"type": "BOOL", "value": "ON"}, + "CELERITAS_BUILTIN_GTest": {"type": "BOOL", "value": "OFF"}, + "CELERITAS_BUILTIN_nlohmann_json": {"type": "BOOL", "value": "OFF"}, "CELERITAS_USE_HepMC3": {"type": "BOOL", "value": "OFF"}, "CELERITAS_USE_ROOT": {"type": "BOOL", "value": "OFF"}, "CELERITAS_USE_VecGeom": {"type": "BOOL", "value": "ON"}, "CELERITAS_VecGeom_SURFACE": {"type": "BOOL", "value": "ON"}, "VecGeom_DIR": {"type": "PATH", "value": "/Users/seth/Code/vecgeom/install/lib/cmake/VecGeom"} - } }, { @@ -150,7 +160,7 @@ { "name": "geant4", "displayName": "Using Geant4 navigation", - "inherits": [".base", ".debug", "default"], + "inherits": [".spackbase", ".debug", "default"], "cacheVariables": { "CELERITAS_CORE_GEO": "Geant4", "CELERITAS_UNITS": "CLHEP", @@ -162,7 +172,7 @@ { "name": "reldeb", "displayName": "With ORANGE in optimized+errcheck mode", - "inherits": [".reldeb", ".base"], + "inherits": [".reldeb", ".spackbase"], "cacheVariables": { "CELERITAS_USE_ROOT": {"type": "BOOL", "value": "OFF"} } @@ -178,7 +188,7 @@ { "name": "ndebug", "displayName": "With ORANGE in optimized mode", - "inherits": [".ccache", ".ndebug", ".base"], + "inherits": [".ccache", ".ndebug", ".spackbase"], "cacheVariables": { "CELERITAS_USE_ROOT": {"type": "BOOL", "value": "OFF"} } @@ -186,7 +196,7 @@ { "name": "ndebug-vecgeom", "displayName": "With vecgeom, clhep, dev g4 in optimized mode", - "inherits": [".ccache", ".base", ".ndebug", "vecgeom"], + "inherits": [".ccache", ".spackbase", ".ndebug", "vecgeom"], "cacheVariables": { "CELERITAS_USE_HepMC3": {"type": "BOOL", "value": "OFF"}, "CELERITAS_USE_ROOT": {"type": "BOOL", "value": "OFF"}, From 3500642e0ae24c35600be7cb6510e2d0abc56487 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Wed, 4 Feb 2026 13:15:25 -0500 Subject: [PATCH 34/35] Fix main geometry tests by increasing bump distance due to single-precision bvh All geo tests pass for vg 1.2.11, solids 2.0.0-rc7, surface 2.0.0-rc7, G4, ORANGE --- src/geocel/vg/VecgeomTrackView.hh | 13 +++--- test/geocel/GeoTests.cc | 67 +++---------------------------- 2 files changed, 11 insertions(+), 69 deletions(-) diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index ca15174ab1..dcca13098e 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -8,6 +8,7 @@ #include // NOTE: must include Global before most other vecgeom/veccore includes +#include #include #include #include @@ -19,6 +20,7 @@ #include "corecel/cont/Span.hh" #include "corecel/math/ArraySoftUnit.hh" #include "corecel/math/ArrayUtils.hh" +#include "corecel/math/NumericLimits.hh" #include "corecel/sys/ThreadId.hh" #include "geocel/Types.hh" @@ -190,15 +192,12 @@ class VecgeomTrackView // Static data - //! Absolute tolerance in vecgeom::kTolerance - //! 1e-9 for double, 1e-3 for single - static constexpr vg_real_type abs_tolerance_ - = CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE ? 1e-9 : 1e-3; - //! Tolerance used for solid model relocation bump - //! V1 solids by default; static constexpr vg_real_type relocate_bump_ - = CELERITAS_VECGEOM_SURFACE ? 0 : 10 * abs_tolerance_; + = CELERITAS_VECGEOM_SURFACE ? 0 + : std::is_same_v ? 1e-4f + : std::is_same_v ? 1e-5f + : 1e-8; //// HELPER FUNCTIONS //// diff --git a/test/geocel/GeoTests.cc b/test/geocel/GeoTests.cc index 569594eeb5..917d907707 100644 --- a/test/geocel/GeoTests.cc +++ b/test/geocel/GeoTests.cc @@ -393,7 +393,7 @@ void AtlasHgtdGeoTest::test_detailed_tracking() const // the updated point is still inside the original volume. SHOULD_FAIL_WHEN(geo.cross_boundary(), test_->geometry_type() == "VecGeom" - && vecgeom_version < Version{2}); + && using_solids_vg); EXPECT_EQ("HGTD", test_->volume_name(geo)); EXPECT_TRUE(geo.is_on_boundary()); @@ -522,18 +522,6 @@ void CmseGeoTest::test_trace() const 29.931406871193, 40.276406871193, 57.573593128807, 57.573593128807, 57.573593128807, 57.573593128807,}; // clang-format on - if (test_->geometry_type() == "VecGeom" - && vecgeom_version >= Version{2, 0} && using_solids_vg) - { - /* - * Failed to cross boundary: - * unique volume instance is the same before and after - * at {30, 30, 655.04101161124} [cm] along {0, 0, 1}, - * [FAILED] [ON BOUNDARY] in OCMS_PV/CMSE/MUON - */ - ref.fail_at(12); - } - auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } @@ -553,17 +541,6 @@ void CmseGeoTest::test_trace() const ref.halfway_safeties = {85, 267.5, 85.85, 60.4, 0.078366388350267, 2.343262600759, 0.078366388350244, 60.4, 85.85, 267.5, 460,}; // clang-format on - if (test_->geometry_type() == "VecGeom" - && vecgeom_version >= Version{2, 0} && using_solids_vg) - { - /* - * Failed to cross boundary: - * unique volume instance is the same before and after - * at {-295, 0, -48.5} [cm] along {1, 0, 0}, - * [FAILED] [ON BOUNDARY] in OCMS_PV/CMSE/MUON - */ - ref.fail_at(2); - } auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); @@ -577,17 +554,6 @@ void CmseGeoTest::test_trace() const ref.distances = {12.495, 287.505, 530, 920}; ref.dot_normal = {}; // All normals are along track dir ref.halfway_safeties = {6.2475, 47.95, 242, 460}; - if (test_->geometry_type() == "VecGeom" - && vecgeom_version >= Version{2, 0} && using_solids_vg) - { - /* - * Failed to cross boundary: - * unique volume instance is the same before and after - * at {300, 0, 1328} [cm] along {1, 0, 0}, - * [FAILED] [ON BOUNDARY] in OCMS_PV/CMSE/OQUA@0 - */ - ref.fail_at(2); - } auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); @@ -1943,17 +1909,6 @@ void SolidsGeoTest::test_trace() const ref.halfway_safeties[14] = 42.8397753718277; ref.halfway_safeties[15] = 18.8833925371992; ref.halfway_safeties[16] = 42.8430141842906; - - if (vecgeom_version >= Version{2}) - { - /* - * Unique volume instance is the same before and after - * Failed during cross_boundary: - * at {258.1, 0, 0.5} [cm] along {-1, 0, 0}, - * [FAILED] [ON BOUNDARY] in World_PV/polycone1_PV - */ - ref.fail_at(4); - } } auto tol = test_->tracking_tol(); @@ -2066,13 +2021,10 @@ void SolidsGeoTest::test_trace() const if (vecgeom_version >= Version{2}) { /* - * Unique volume instance is the same before and after - * Failed during cross_boundary: - * at {-335.05, -125, 0.5} [cm] along {1, 0, 0}, - * [FAILED] [ON BOUNDARY] in World_PV/arb8b_PV + * Issues near World_PV/arb8b_PV */ ref.halfway_safeties[4] = 39.0470100365853; - ref.fail_at(5); + ref.halfway_safeties[6] = 29.8360600858068; } } @@ -2336,6 +2288,8 @@ void SimpleCmsGeoTest::test_trace() const auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } + // Fails to find first step + if (!using_solids_vg || vecgeom_version < Version{2, 0}) { SCOPED_TRACE("edge case with flipped direction"); auto result = test_->track(pos, -dir); @@ -2373,13 +2327,6 @@ void SimpleCmsGeoTest::test_trace() const 232.37188601505, }; - if (using_solids_vg && vecgeom_version >= Version{2, 0}) - { - // Track fails immediately - ref = {}; - ref.fail(); - } - auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } @@ -3192,10 +3139,6 @@ void ZnenvGeoTest::test_trace() const auto tol = test_->tracking_tol(); fixup_orange(*test_, ref, result, "World"); - if (using_solids_vg && vecgeom_version >= Version{2, 0}) - { - GTEST_SKIP() << "FIXME: Znenv VecGeom model construction failure."; - } EXPECT_REF_NEAR(ref, result, tol); } } From bc5149e84b3070703c728ff6ce66aa6f63c9cd37 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Wed, 4 Feb 2026 13:24:49 -0500 Subject: [PATCH 35/35] Update field propagator for ORANGE This was broken before the last master merge --- test/celeritas/field/FieldPropagator.test.cc | 27 ++++++++++---------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index 88deb0b3a5..7d0e9733d3 100644 --- a/test/celeritas/field/FieldPropagator.test.cc +++ b/test/celeritas/field/FieldPropagator.test.cc @@ -1077,31 +1077,31 @@ TEST_F(TwoBoxesTest, } } - std::vector expected_boundary = {1, 1, 1, 1, 1, 0, 1, 0, 1, 0}; + std::vector expected_boundary = {1, 0, 1, 0, 1, 0, 1, 0, 1, 0}; std::vector expected_distances = { 0.0078534718906499, - 0.0028235332722979, + 0.0078539816339745, 0.0044879852658442, - 0.0028259738005751, - 1e-05, + 0.0044879895051283, 1e-05, + 1e-06, 9.9999658622419e-09, 1e-08, 9.9981633254417e-12, 1e-11, }; - std::vector expected_substeps = {1, 25, 1, 12, 1, 1, 1, 1, 1, 1}; + std::vector expected_substeps = {1, 1, 1, 1, 1, 4, 1, 1, 1, 1}; std::vector expected_volumes = { - "world", "inner", - "world", "inner", - "world", - "world", - "world", - "world", - "world", - "world", + "inner", + "inner", + "inner", + "inner", + "inner", + "inner", + "inner", + "inner", }; /* @@ -1120,6 +1120,7 @@ TEST_F(TwoBoxesTest, expected_distances[7] = 0; expected_substeps[1] = 1; expected_substeps[3] = 1; + expected_substeps[5] = 1; expected_substeps[9] = 0; expected_volumes = { "inner",