Fix class name conflict, remove default initializers, and tweak field…

… driver (#504)

* Rename stepper result
* Update cmake presets
* Define rsqrt
* Minor consistency changes to driver and test
* Add MagFieldEquation test
* Delete default initializers
* More minor refactoring
* Use detail::ax for updating results
* Uniformly use ints for explicit 3-element looping

scripts/cmake-presets/goldfinger.json

@@ -10,6 +10,8 @@
       "cacheVariables": {
         "CELERITAS_USE_CUDA":    {"type": "BOOL",   "value": "OFF"},
         "CELERITAS_USE_HIP":     {"type": "BOOL",   "value": "OFF"},
+        "CELERITAS_USE_JSON":     {"type": "BOOL",   "value": "ON"},
+        "CELERITAS_USE_Geant4":   {"type": "BOOL",   "value": "ON"},
         "CMAKE_BUILD_TYPE":      {"type": "STRING", "value": "Debug"},
         "CMAKE_EXPORT_COMPILE_COMMANDS": {"type": "BOOL",   "value": "ON"},
         "CMAKE_OSX_DEPLOYMENT_TARGET": {"type": "STRING", "value": "12"},
@@ -27,7 +29,7 @@
         "CELERITAS_BUILD_DEMOS": {"type": "BOOL",   "value": "OFF"},
         "CELERITAS_BUILD_TESTS": {"type": "BOOL",   "value": "OFF"},
         "CELERITAS_USE_JSON":    {"type": "BOOL", "value": "ON"},
-        "CMAKE_CXX_INCLUDE_WHAT_YOU_USE": "$env{IWYU_ROOT}/bin/include-what-you-use;-Xiwyu;--no_fwd_decls;-Xiwyu;--no_comments;-Xiwyu;--transitive_includes_only"
+        "CMAKE_CXX_INCLUDE_WHAT_YOU_USE": "/opt/homebrew/bin/include-what-you-use;-Xiwyu;--no_fwd_decls;-Xiwyu;--no_comments;-Xiwyu;--transitive_includes_only"
       }
     },
     {

src/celeritas/em/distribution/UrbanMscScatter.hh

@@ -230,7 +230,7 @@ CELER_FUNCTION auto UrbanMscScatter::operator()(Engine& rng) -> MscInteraction
         if (length > 0)
         {
             result.displacement = this->sample_displacement_dir(rng, phi);
-            for (auto i : range(3))
+            for (int i = 0; i < 3; ++i)
             {
                 result.displacement[i] *= length;
             }

src/celeritas/em/interactor/MollerBhabhaInteractor.hh

@@ -11,7 +11,6 @@
 
 #include "corecel/Macros.hh"
 #include "corecel/Types.hh"
-#include "corecel/cont/Range.hh"
 #include "corecel/data/StackAllocator.hh"
 #include "corecel/math/Algorithms.hh"
 #include "corecel/math/ArrayUtils.hh"
@@ -182,7 +181,7 @@ CELER_FUNCTION Interaction MollerBhabhaInteractor::operator()(Engine& rng)
 
     // Calculate incident particle final direction
     Real3 inc_exiting_direction;
-    for (int i : range(3))
+    for (int i = 0; i < 3; ++i)
     {
         real_type inc_momentum_ijk       = inc_momentum_ * inc_direction_[i];
         real_type secondary_momentum_ijk = secondary_momentum

src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh

@@ -9,7 +9,6 @@
 
 #include "corecel/Macros.hh"
 #include "corecel/Types.hh"
-#include "corecel/cont/Range.hh"
 #include "corecel/data/StackAllocator.hh"
 #include "corecel/math/ArrayUtils.hh"
 #include "celeritas/Constants.hh"
@@ -145,7 +144,7 @@ CELER_FUNCTION Interaction MuBremsstrahlungInteractor::operator()(Engine& rng)
     Real3 gamma_dir = rotate(from_spherical(cost, phi(rng)), inc_direction_);
 
     Real3 inc_direction;
-    for (int i : range(3))
+    for (int i = 0; i < 3; ++i)
     {
         inc_direction[i] = value_as<Momentum>(particle_.momentum())
                                * inc_direction_[i]

src/celeritas/em/interactor/RayleighInteractor.hh

@@ -182,7 +182,7 @@ auto RayleighInteractor::evaluate_weight_and_prob() const -> SampleInput
     axpy(input.factor, b, &x);
 
     Real3 prob;
-    for (auto i : range(3))
+    for (int i = 0; i < 3; ++i)
     {
         input.weight[i] = (x[i] > this->fit_slice())
                               ? 1 - fastpow(1 + x[i], -n[i])

src/celeritas/em/interactor/detail/BremFinalStateHelper.hh

@@ -101,7 +101,7 @@ CELER_FUNCTION Interaction BremFinalStateHelper::operator()(
         = rotate(from_spherical(cost, sample_phi(rng)), inc_direction_);
 
     // Update parent particle direction
-    for (unsigned int i : range(3))
+    for (int i = 0; i < 3; ++i)
     {
         real_type inc_momentum_i   = inc_momentum_.value() * inc_direction_[i];
         real_type gamma_momentum_i = result.secondaries[0].energy.value()

src/celeritas/field/DormandPrinceStepper.hh

@@ -63,7 +63,7 @@ class DormandPrinceStepper
   public:
     //!@{
     //! Type aliases
-    using result_type = StepperResult;
+    using result_type = FieldStepperResult;
     //!@}
 
   public:
@@ -189,6 +189,7 @@ DormandPrinceStepper<E>::operator()(real_type       step,
     OdeState k7 = calc_rhs_(result.end_state);
 
     // The error estimate
+    result.err_state = {{0, 0, 0}, {0, 0, 0}};
     axpy(d71 * step, k1, &result.err_state);
     axpy(d73 * step, k3, &result.err_state);
     axpy(d74 * step, k4, &result.err_state);

src/celeritas/field/FieldDriver.hh

@@ -14,9 +14,8 @@
 #include "corecel/math/Algorithms.hh"
 
 #include "FieldDriverOptions.hh"
-#include "MagFieldEquation.hh"
-#include "RungeKuttaStepper.hh"
 #include "Types.hh"
+#include "detail/FieldUtils.hh"
 
 namespace celeritas
 {
@@ -173,9 +172,9 @@ FieldDriver<StepperT>::find_next_chord(real_type       step,
     // Output with a step control error
     ChordSearch output;
 
-    bool          succeeded       = false;
-    size_type     remaining_steps = options_.max_nsteps;
-    StepperResult result;
+    bool               succeeded       = false;
+    size_type          remaining_steps = options_.max_nsteps;
+    FieldStepperResult result;
 
     do
     {
@@ -197,7 +196,7 @@ FieldDriver<StepperT>::find_next_chord(real_type       step,
             // Estimate a new trial chord with a relative scale
             step *= max(std::sqrt(options_.delta_chord / dchord), half());
         }
-    } while (!succeeded && (--remaining_steps > 0));
+    } while (!succeeded && --remaining_steps > 0);
 
     // TODO: loop check and handle rare cases if happen
     CELER_ASSERT(succeeded);
@@ -291,7 +290,7 @@ FieldDriver<StepperT>::integrate_step(real_type       step,
     else
     {
         // Do an integration step for a small step (a.k.a quick advance)
-        StepperResult result = apply_step_(step, state);
+        FieldStepperResult result = apply_step_(step, state);
 
         // Update position and momentum
         output.end.state = result.end_state;
@@ -303,7 +302,7 @@ FieldDriver<StepperT>::integrate_step(real_type       step,
         // Compute a proposed new step
         CELER_ASSERT(output.end.step > 0);
         output.proposed_step = this->new_step_size(
-            step, dyerr / (step * options_.epsilon_step));
+            step, dyerr / (options_.epsilon_step * step));
     }
 
     return output;
@@ -326,7 +325,7 @@ FieldDriver<StepperT>::one_good_step(real_type step, const OdeState& state) cons
     bool          succeeded       = false;
     size_type     remaining_steps = options_.max_nsteps;
     real_type     errmax2;
-    StepperResult result;
+    FieldStepperResult result;
 
     do
     {
@@ -353,13 +352,12 @@ FieldDriver<StepperT>::one_good_step(real_type step, const OdeState& state) cons
     CELER_ASSERT(succeeded);
 
     // Update state, step taken by this trial and the next predicted step
-    output.end.state = result.end_state;
-    output.end.step  = step;
-    output.proposed_step
-        = (errmax2 > ipow<2>(options_.errcon))
-              ? options_.safety * step
-                    * fastpow(errmax2, half() * options_.pgrow)
-              : options_.max_stepping_increase * step;
+    output.end.state     = result.end_state;
+    output.end.step      = step;
+    output.proposed_step = (errmax2 > ipow<2>(options_.errcon))
+                               ? options_.safety * step
+                                     * fastpow(errmax2, half() * options_.pgrow)
+                               : options_.max_stepping_increase * step;
 
     return output;
 }

src/celeritas/field/MagFieldEquation.hh

@@ -15,6 +15,7 @@
 #include "celeritas/Quantities.hh"
 
 #include "Types.hh"
+#include "detail/FieldUtils.hh"
 
 namespace celeritas
 {
@@ -90,18 +91,18 @@ template<class FieldT>
 CELER_FUNCTION auto
 MagFieldEquation<FieldT>::operator()(const OdeState& y) const -> OdeState
 {
-    // Get a magnetic field value at a given position
-    auto&& mag_vec = calc_field_(y.pos);
+    CELER_EXPECT(y.mom[0] != 0 || y.mom[1] != 0 || y.mom[2] != 0);
+    real_type momentum_inv = celeritas::rsqrt(dot_product(y.mom, y.mom));
 
-    real_type momentum_mag2 = dot_product(y.mom, y.mom);
-    CELER_ASSERT(momentum_mag2 > 0);
-    real_type momentum_inv = 1 / std::sqrt(momentum_mag2);
-
-    // Evaluate the right-hand-side of the equation
+    // Evaluate the rate of change in particle's position per unit length: this
+    // is just the direction
     OdeState result;
+    result.pos = detail::ax(momentum_inv, y.mom);
 
-    axpy(momentum_inv, y.mom, &result.pos);
-    axpy(coeffi_ * momentum_inv, cross_product(y.mom, mag_vec), &result.mom);
+    // Calculate the magnetic field value at the current position
+    // to calculate the force on the particle
+    result.mom = detail::ax(coeffi_ * momentum_inv,
+                            cross_product(y.mom, calc_field_(y.pos)));
 
     return result;
 }

src/celeritas/field/RungeKuttaStepper.hh

@@ -31,7 +31,7 @@ class RungeKuttaStepper
   public:
     //!@{
     //! Type aliases
-    using result_type = StepperResult;
+    using result_type = FieldStepperResult;
     //!@}
 
   public:

src/celeritas/field/Types.hh

@@ -25,15 +25,15 @@ struct OdeState
     using MomentumUnits = units::MevMomentum;
     using Real3         = Array<real_type, 3>;
 
-    Real3 pos{0, 0, 0}; //!< Particle position
-    Real3 mom{0, 0, 0}; //!< Particle momentum
+    Real3 pos; //!< Particle position
+    Real3 mom; //!< Particle momentum
 };
 
 //---------------------------------------------------------------------------//
 /*!
  * The result of the integration stepper.
  */
-struct StepperResult
+struct FieldStepperResult
 {
     OdeState mid_state; //!< OdeState at the middle
     OdeState end_state; //!< OdeState at the end

src/celeritas/field/ZHelixStepper.hh

@@ -8,14 +8,13 @@
 #pragma once
 
 #include "corecel/Types.hh"
-#include "corecel/cont/Range.hh"
 #include "corecel/math/Algorithms.hh"
 
 #include "Types.hh"
-#include "UniformZField.hh"
 
 namespace celeritas
 {
+class UniformZField;
 //---------------------------------------------------------------------------//
 /*!
  * Analytically step along a helical path for a uniform Z magnetic field.
@@ -35,7 +34,7 @@ class ZHelixStepper
   public:
     //!@{
     //! Type aliases
-    using result_type = StepperResult;
+    using result_type = FieldStepperResult;
     //!@}
 
   public:
@@ -116,10 +115,10 @@ ZHelixStepper<E>::operator()(real_type step, const OdeState& beg_state) const
     result.end_state = this->move(step, radius, helicity, beg_state, rhs);
 
     // Solution are exact, but assign a tolerance for numerical treatments
-    for (auto i : range(3))
+    for (int i = 0; i < 3; ++i)
     {
-        result.err_state.pos[i] += ZHelixStepper::tolerance();
-        result.err_state.mom[i] += ZHelixStepper::tolerance();
+        result.err_state.pos[i] = ZHelixStepper::tolerance();
+        result.err_state.mom[i] = ZHelixStepper::tolerance();
     }
 
     return result;
@@ -156,14 +155,13 @@ CELER_FUNCTION OdeState ZHelixStepper<E>::move(real_type       step,
                                                const OdeState& beg_state,
                                                const OdeState& rhs) const
 {
-    OdeState end_state;
-
     // Solution for position and momentum after moving delta_phi on the helix
     real_type del_phi = (helicity == Helicity::positive) ? step / radius
                                                          : -step / radius;
     real_type sin_phi = std::sin(del_phi);
     real_type cos_phi = std::cos(del_phi);
 
+    OdeState end_state;
     end_state.pos = {(beg_state.pos[0] * cos_phi - beg_state.pos[1] * sin_phi),
                      (beg_state.pos[0] * sin_phi + beg_state.pos[1] * cos_phi),
                      beg_state.pos[2] + del_phi * radius * rhs.pos[2]};
@@ -173,7 +171,7 @@ CELER_FUNCTION OdeState ZHelixStepper<E>::move(real_type       step,
                      rhs.pos[2]};
 
     real_type momentum = norm(beg_state.mom);
-    for (auto i : range(3))
+    for (int i = 0; i < 3; ++i)
     {
         end_state.mom[i] *= momentum;
     }

src/celeritas/field/detail/FieldUtils.hh

@@ -43,7 +43,7 @@ template<class T>
 inline CELER_FUNCTION Array<T, 3> ax(T a, const Array<T, 3>& x)
 {
     Array<T, 3> result;
-    for (size_type i = 0; i != 3; ++i)
+    for (int i = 0; i < 3; ++i)
     {
         result[i] = a * x[i];
     }
@@ -57,13 +57,13 @@ inline CELER_FUNCTION Array<T, 3> ax(T a, const Array<T, 3>& x)
 inline CELER_FUNCTION Chord make_chord(const Real3& src, const Real3& dst)
 {
     Chord result;
-    for (size_type i = 0; i != 3; ++i)
+    for (int i = 0; i < 3; ++i)
     {
         result.dir[i] = dst[i] - src[i];
     }
     result.length        = norm(result.dir);
     const real_type norm = 1 / result.length;
-    for (size_type i = 0; i != 3; ++i)
+    for (int i = 0; i < 3; ++i)
     {
         result.dir[i] *= norm;
     }
@@ -88,7 +88,7 @@ inline CELER_FUNCTION real_type calc_miss_distance(const Real3& pos,
                                                    const Real3& target)
 {
     real_type delta_sq = 0;
-    for (size_type i = 0; i != 3; ++i)
+    for (int i = 0; i < 3; ++i)
     {
         delta_sq += ipow<2>(pos[i] - target[i] + distance * dir[i]);
     }
@@ -126,8 +126,10 @@ inline CELER_FUNCTION real_type truncation_error(real_type       step,
                                                  const OdeState& beg_state,
                                                  const OdeState& err_state)
 {
+    CELER_EXPECT(step > 0);
+    CELER_EXPECT(eps_rel_max > 0);
+
     // Evaluate tolerance and squre of the position and momentum accuracy
-    real_type eps_pos = eps_rel_max * step;
 
     real_type magvel2 = dot_product(beg_state.mom, beg_state.mom);
     real_type errpos2 = dot_product(err_state.pos, err_state.pos);
@@ -137,7 +139,7 @@ inline CELER_FUNCTION real_type truncation_error(real_type       step,
     CELER_ASSERT(errpos2 >= 0);
     CELER_ASSERT(magvel2 > 0);
 
-    errpos2 /= ipow<2>(eps_pos);
+    errpos2 /= ipow<2>(eps_rel_max * step);
     errvel2 /= (magvel2 * ipow<2>(eps_rel_max));
 
     // Return the square of the maximum truncation error
@@ -156,10 +158,10 @@ inline CELER_FUNCTION real_type distance_chord(const OdeState& beg_state,
                                                const OdeState& mid_state,
                                                const OdeState& end_state)
 {
-    Real3 beg_mid = mid_state.pos;
-    Real3 beg_end = end_state.pos;
+    Real3 beg_mid;
+    Real3 beg_end;
 
-    for (size_type i = 0; i != 3; ++i)
+    for (int i = 0; i < 3; ++i)
     {
         beg_mid[i] = mid_state.pos[i] - beg_state.pos[i];
         beg_end[i] = end_state.pos[i] - beg_state.pos[i];

src/celeritas/global/alongstep/detail/UrbanMsc.hh

@@ -170,7 +170,7 @@ CELER_FUNCTION void UrbanMsc::apply_step(CoreTrackView const& track,
         // Displacment during a boundary crossing is *not* OK
         CELER_ASSERT(!geo.is_on_boundary());
         Real3 new_pos;
-        for (int i : range(3))
+        for (int i = 0; i < 3; ++i)
         {
             new_pos[i] = geo.pos()[i] + msc_result.displacement[i];
         }

src/celeritas/random/distribution/GammaDistribution.hh

@@ -86,7 +86,7 @@ GammaDistribution<RealType>::GammaDistribution(real_type alpha, real_type beta)
     , beta_(beta)
     , alpha_p_(alpha < 1 ? alpha + 1 : alpha)
     , d_(alpha_p_ - real_type(1) / 3)
-    , c_(1 / std::sqrt(9 * d_))
+    , c_(celeritas::rsqrt(9 * d_))
 {
     CELER_EXPECT(alpha_ > 0);
     CELER_EXPECT(beta_ > 0);