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"}, diff --git a/src/celeritas/field/FieldPropagator.hh b/src/celeritas/field/FieldPropagator.hh index 49d9f4f272..ed2d858884 100644 --- a/src/celeritas/field/FieldPropagator.hh +++ b/src/celeritas/field/FieldPropagator.hh @@ -181,7 +181,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 7463144570..1b86f4a443 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" @@ -56,7 +57,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/corecel/math/Algorithms.hh b/src/corecel/math/Algorithms.hh index cafba3accf..1922c4bbc3 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. 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 }; //---------------------------------------------------------------------------// 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/VecgeomData.hh b/src/geocel/vg/VecgeomData.hh index 4c6aadadf9..fd1e15fdd7 100644 --- a/src/geocel/vg/VecgeomData.hh +++ b/src/geocel/vg/VecgeomData.hh @@ -13,7 +13,6 @@ #include "corecel/Macros.hh" #include "corecel/Types.hh" #include "corecel/data/Collection.hh" -#include "corecel/data/CollectionBuilder.hh" #include "corecel/sys/ThreadId.hh" #include "geocel/Types.hh" 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/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 682694767f..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" @@ -41,6 +43,12 @@ # include "detail/VgNavStateWrapper.hh" #endif +#if !CELER_DEVICE_COMPILE +# include "corecel/io/Logger.hh" +# include "corecel/io/Repr.hh" +# include "geocel/detail/LengthUnits.hh" +#endif + namespace celeritas { //---------------------------------------------------------------------------// @@ -54,8 +62,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 +95,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 //// //!@{ @@ -121,7 +125,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; @@ -184,6 +188,16 @@ class VecgeomTrackView // Temporary data real_type next_step_{0}; + bool failed_{false}; + + // Static data + + //! Tolerance used for solid model relocation bump + static constexpr vg_real_type relocate_bump_ + = CELERITAS_VECGEOM_SURFACE ? 0 + : std::is_same_v ? 1e-4f + : std::is_same_v ? 1e-5f + : 1e-8; //// HELPER FUNCTIONS //// @@ -193,9 +207,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; @@ -254,6 +265,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; @@ -285,14 +297,27 @@ 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; 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_) + << lengthunits::label; +#endif + failed_ = true; + } + return *this; } @@ -426,32 +451,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); } @@ -469,6 +474,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_), @@ -480,6 +486,23 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) *next_surf_ #endif ); + + 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_) << ' ' + << lengthunits::label << " along " << repr(dir_) + << ": computed step is " << repr(next_step_) << ' ' + << lengthunits::label; +#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 @@ -487,14 +510,12 @@ 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 // BVH navigator subtracts and then re-adds a bump distance to the // step - CELER_ASSERT(soft_equal(next_step_, max(max_step, this->extra_push()))); + CELER_ASSERT(soft_equal(next_step_, max_step)); next_step_ = max_step; } @@ -504,9 +525,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; } @@ -532,13 +552,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 @@ -571,22 +586,47 @@ 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()) + VgReal3 bumped_pos; + if constexpr (CELERITAS_VECGEOM_SURFACE) + { + // Surface relocation is exact, assuming a well constructed + // geometry + bumped_pos = detail::to_vector(pos_); + } + else { - // In surf model, relocation does not work from [OUTSIDE] - Navigator::RelocateToNextVolume(detail::to_vector(this->pos_), - detail::to_vector(this->dir_), + 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_, + *next_surf_, #endif - vgnext_); - } + 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 " + << repr(pos_) << ' ' << lengthunits::label; +#endif + failed_ = true; + return; } vgstate_ = vgnext_; @@ -637,6 +677,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) { @@ -673,17 +716,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. diff --git a/src/geocel/vg/VecgeomTypes.hh b/src/geocel/vg/VecgeomTypes.hh index 405a618c31..bdf1c8bfaf 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 065d9e2285..dedf7b9495 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 @@ -49,12 +48,6 @@ class BVHNavigator using NavState = detail::VgNavStateWrapper; #endif -#ifdef VECGEOM_FLOAT_PRECISION - static constexpr vg_real_type kBoundaryPush = 10 * 1e-3f; -#else - static constexpr vg_real_type kBoundaryPush = 10 * 1e-9; -#endif - //! Update path (which must be reset in advance) CELER_FUNCTION static void LocatePointIn(VgPlacedVol const* vol, @@ -73,6 +66,7 @@ class BVHNavigator } } + CELER_ASSERT(vol); path.Push(vol); VgReal3 currentpoint(point); @@ -93,153 +87,22 @@ 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; } } - CELER_FUNCTION static void - RelocatePoint(VgReal3 const& localpoint, NavState& path) - { - 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, - 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; @@ -259,88 +122,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,18 +134,23 @@ class BVHNavigator VgReal3 const& globaldir, vg_real_type step_limit, NavState const& in_state, - NavState& out_state, - vg_real_type push = 0) + NavState& out_state) { + // See VecgeomTrackView::relocate_bump_ ; this is from the VG 2.0 + // boundary, setting external push to zero +#ifdef VECGEOM_FLOAT_PRECISION + static constexpr vg_real_type kBoundaryPush = 10 * 1e-3f; +#else + static constexpr vg_real_type kBoundaryPush = 10 * 1e-9; +#endif + // 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; @@ -393,6 +179,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. @@ -419,43 +206,17 @@ 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, - 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(); @@ -474,6 +235,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/SolidsNavigator.hh b/src/geocel/vg/detail/SolidsNavigator.hh index 0bd456db6b..ab1fd31475 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,54 @@ 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}; + 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 d1abe958bf..028074c6ad 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 @@ -70,22 +68,19 @@ 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; } // 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, @@ -94,12 +89,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< @@ -112,25 +102,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, 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(); } diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index 2288bfcb65..4080e83a9c 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -248,7 +248,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 @@ -257,10 +258,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/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); diff --git a/test/celeritas/field/FieldPropagator.test.cc b/test/celeritas/field/FieldPropagator.test.cc index 56a46a54a8..7d0e9733d3 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"); @@ -1031,51 +1044,123 @@ TEST_F(TwoBoxesTest, for (int i : range(2)) { SCOPED_TRACE(i); - auto result = propagate(radius * dtheta); - if (result.boundary) + Propagation result; + 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(); } } - 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", }; - 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 (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[5] = 1; + expected_substeps[9] = 0; + expected_volumes = { + "inner", + "inner", + "inner", + "inner", + "inner", + "[FAILURE]", + "inner", + "[FAILURE]", + "[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) + { + 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 @@ -1222,15 +1307,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()); - // NOTE: vecgeom 2.x-solids puts this position slightly *outside* the beam - // tube rather than *inside* - if (using_solids_vg && CELERITAS_VECGEOM_VERSION >= 0x020000) + if (this->volume_name(geo) == "vacuum_tube") { - // TODO: VecGeom 2.x-solids starts to diverge here - EXPECT_EQ("vacuum_tube", this->volume_name(geo)); - GTEST_SKIP() << "FIXME: VecGeom 2.x-solid construction failure."; + // NOTE: vecgeom 2 solid model thinks r=30 + epsilon is *inside* the + // 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()); @@ -1247,7 +1330,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 @@ -1267,18 +1350,35 @@ 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(); + } + + 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); - if (using_surface_vg) - { - GTEST_SKIP() << "FIXME: VecGeom surface model fails a boundary " - "requirement."; - } ASSERT_TRUE(geo.is_on_boundary()); if (geo.check_normal()) @@ -1315,7 +1415,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}); @@ -1323,7 +1422,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()); @@ -1337,132 +1437,38 @@ 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"); - if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_ORANGE) + try { - // ORANGE should successfully reenter. However, under certain - // system configurations, VecGeom will end up in the world volume, - // so we don't test in all cases. - EXPECT_EQ("em_calorimeter", this->volume_name(geo)); - } - else - { - EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; + geo.cross_boundary(); } - - if (!successful_reentry) + catch (CheckedGeoError const& e) { - // 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 (using_solids_vg && CELERITAS_VECGEOM_VERSION < 0x020000) + if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_VECGEOM) { - // FIXME: VecGeom 1.x navigation failure (solids model) - EXPECT_EQ("world", this->volume_name(geo)); + GTEST_SKIP() << e.what(); } else { - EXPECT_EQ("si_tracker", this->volume_name(geo)); + FAIL() << e.what(); } - - // 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 - auto result = geo.find_next_step(1); - EXPECT_LT(result.distance, 2e-8); - - if (result.distance < 1e-6) - { - geo.move_to_boundary(); - geo.cross_boundary(); // back into em_calorimeter - } - - // then they are back to agreement - 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) - 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) + if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_ORANGE) { - // 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); + // ORANGE should successfully reenter. However, under certain + // system configurations, VecGeom will end up in the world volume, + // so we don't test in all cases. 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_; - } + // 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). + // FIXME: see GeoTests: TwoBoxesGeoTest::test_tangent + EXPECT_EQ("si_tracker", this->volume_name(geo)); + GTEST_SKIP(); } } } @@ -1486,7 +1492,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}) { @@ -1502,23 +1507,36 @@ 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 { result = propagate(radius); } - catch (RuntimeError const& e) + catch (CheckedGeoError const& e) { - // Failure during Geant4 propagation 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; } } @@ -1533,8 +1551,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}; @@ -1542,30 +1561,44 @@ 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) { - // 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_; + expected_num_boundary[1] = 37; + expected_num_step[1] = 179; + expected_num_intercept[1] = 614; + 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)", + }; + } + else if (using_surface_vg) + { + 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); - EXPECT_VEC_EQ(expected_num_step, num_step); - EXPECT_VEC_EQ(expected_num_intercept, num_intercept); - EXPECT_VEC_EQ(expected_num_integration, num_integration); + 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_; } //---------------------------------------------------------------------------// diff --git a/test/geocel/CheckedGeoTrackView.cc b/test/geocel/CheckedGeoTrackView.cc index 7624830c4a..130df9b75c 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; } @@ -267,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) { 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 2b3432b033..d82d79ee26 100644 --- a/test/geocel/GenericGeoResults.hh +++ b/test/geocel/GenericGeoResults.hh @@ -69,6 +69,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/GeoTests.cc b/test/geocel/GeoTests.cc index 19adff0680..917d907707 100644 --- a/test/geocel/GeoTests.cc +++ b/test/geocel/GeoTests.cc @@ -25,6 +25,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 @@ -350,15 +387,13 @@ void AtlasHgtdGeoTest::test_detailed_tracking() const EXPECT_SOFT_EQ(344.45, to_cm(geo.pos()[2])); EXPECT_EQ("SPlate", test_->volume_name(geo)); EXPECT_TRUE(geo.is_on_boundary()); - geo.cross_boundary(); - if (test_->geometry_type() == "VecGeom" && vecgeom_version < Version{2}) - { - // VecGeom fails to cross the boundary! the internal bump along the - // path of travel doesn't change the Z coordinate, so it assumes - // the updated point is still inside the original volume. - EXPECT_EQ("SPlate", test_->volume_name(geo)); - return; - } + + // VecGeom fails to cross the boundary! the internal bump along the + // path of travel doesn't change the Z coordinate, so it assumes + // the updated point is still inside the original volume. + SHOULD_FAIL_WHEN(geo.cross_boundary(), + test_->geometry_type() == "VecGeom" + && using_solids_vg); EXPECT_EQ("HGTD", test_->volume_name(geo)); EXPECT_TRUE(geo.is_on_boundary()); @@ -506,6 +541,7 @@ 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 + auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } @@ -518,6 +554,7 @@ 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}; + auto tol = test_->tracking_tol(); EXPECT_REF_NEAR(ref, result, tol); } @@ -652,13 +689,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())); @@ -666,10 +696,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) @@ -700,26 +731,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()); } @@ -1604,12 +1618,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); - } } } @@ -1902,13 +1910,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"); @@ -2013,6 +2017,15 @@ 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}) + { + /* + * Issues near World_PV/arb8b_PV + */ + ref.halfway_safeties[4] = 39.0470100365853; + ref.halfway_safeties[6] = 29.8360600858068; + } } if (test_->geometry_type() == "Geant4" @@ -2024,12 +2037,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"); @@ -2091,13 +2099,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); } } @@ -2216,6 +2219,117 @@ 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; + } + + 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); + } + // 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); + 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); + } } //---------------------------------------------------------------------------// @@ -2750,7 +2864,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()); @@ -2767,18 +2880,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) { @@ -2793,9 +2910,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"; } } @@ -2803,9 +2920,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)); @@ -2912,19 +3027,17 @@ 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" + && !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"; } - else if (test_->geometry_type() == "VecGeom" - && "world" == test_->volume_name(geo)) - { - GTEST_SKIP() << "Unexpected boundary crossing"; - } EXPECT_EQ("inner", test_->volume_name(geo)); } @@ -3026,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); } } diff --git a/test/testdetail/TestMacrosImpl.hh b/test/testdetail/TestMacrosImpl.hh index 7390ffd8a4..f04f9aaf27 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))); } //---------------------------------------------------------------------------// @@ -729,7 +729,7 @@ template expected_expr, actual, actual_expr, - soft_comparator(std::forward(cmp_or_tol))); + make_soft_comparator(std::forward(cmp_or_tol))); } //---------------------------------------------------------------------------//