From 28028ca748b7b453943f149ec5baf3b764df0cf3 Mon Sep 17 00:00:00 2001 From: "Seth R. Johnson" Date: Wed, 7 Sep 2022 14:06:58 -0400 Subject: [PATCH] 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 | 4 +- .../em/distribution/UrbanMscScatter.hh | 2 +- .../em/interactor/MollerBhabhaInteractor.hh | 3 +- .../interactor/MuBremsstrahlungInteractor.hh | 3 +- .../em/interactor/RayleighInteractor.hh | 2 +- .../interactor/detail/BremFinalStateHelper.hh | 2 +- src/celeritas/field/DormandPrinceStepper.hh | 3 +- src/celeritas/field/FieldDriver.hh | 30 +++---- src/celeritas/field/MagFieldEquation.hh | 19 +++-- src/celeritas/field/RungeKuttaStepper.hh | 2 +- src/celeritas/field/Types.hh | 6 +- src/celeritas/field/ZHelixStepper.hh | 16 ++-- src/celeritas/field/detail/FieldUtils.hh | 20 +++-- .../global/alongstep/detail/UrbanMsc.hh | 2 +- .../random/distribution/GammaDistribution.hh | 2 +- src/corecel/math/Algorithms.hh | 23 +++++ test/CMakeLists.txt | 1 + test/celeritas/field/FieldDriver.test.cc | 27 +++--- test/celeritas/field/MagFieldEquation.test.cc | 85 +++++++++++++++++++ test/celeritas/field/Steppers.test.cc | 2 +- test/corecel/math/Algorithms.test.cc | 17 ++++ 21 files changed, 202 insertions(+), 69 deletions(-) create mode 100644 test/celeritas/field/MagFieldEquation.test.cc diff --git a/scripts/cmake-presets/goldfinger.json b/scripts/cmake-presets/goldfinger.json index 3badf7ccb9..33dcc6cec3 100644 --- a/scripts/cmake-presets/goldfinger.json +++ b/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" } }, { diff --git a/src/celeritas/em/distribution/UrbanMscScatter.hh b/src/celeritas/em/distribution/UrbanMscScatter.hh index f9052ddca7..dfb964dd96 100644 --- a/src/celeritas/em/distribution/UrbanMscScatter.hh +++ b/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; } diff --git a/src/celeritas/em/interactor/MollerBhabhaInteractor.hh b/src/celeritas/em/interactor/MollerBhabhaInteractor.hh index cf0037b550..7f961b1c6d 100644 --- a/src/celeritas/em/interactor/MollerBhabhaInteractor.hh +++ b/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 diff --git a/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh b/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh index 29b1f8e1e6..1cfb40c1b1 100644 --- a/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh +++ b/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(particle_.momentum()) * inc_direction_[i] diff --git a/src/celeritas/em/interactor/RayleighInteractor.hh b/src/celeritas/em/interactor/RayleighInteractor.hh index 9fd3b2d19d..526483999e 100644 --- a/src/celeritas/em/interactor/RayleighInteractor.hh +++ b/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]) diff --git a/src/celeritas/em/interactor/detail/BremFinalStateHelper.hh b/src/celeritas/em/interactor/detail/BremFinalStateHelper.hh index 21ba50fc0c..4fb8baf3ee 100644 --- a/src/celeritas/em/interactor/detail/BremFinalStateHelper.hh +++ b/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() diff --git a/src/celeritas/field/DormandPrinceStepper.hh b/src/celeritas/field/DormandPrinceStepper.hh index f56de5adc7..9cfa19d480 100644 --- a/src/celeritas/field/DormandPrinceStepper.hh +++ b/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::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); diff --git a/src/celeritas/field/FieldDriver.hh b/src/celeritas/field/FieldDriver.hh index 4b5e283994..86c4de619d 100644 --- a/src/celeritas/field/FieldDriver.hh +++ b/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::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::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::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::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::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::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; } diff --git a/src/celeritas/field/MagFieldEquation.hh b/src/celeritas/field/MagFieldEquation.hh index 29df717f91..b4a611a551 100644 --- a/src/celeritas/field/MagFieldEquation.hh +++ b/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 CELER_FUNCTION auto MagFieldEquation::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; } diff --git a/src/celeritas/field/RungeKuttaStepper.hh b/src/celeritas/field/RungeKuttaStepper.hh index 5a1248449a..bc1716fe01 100644 --- a/src/celeritas/field/RungeKuttaStepper.hh +++ b/src/celeritas/field/RungeKuttaStepper.hh @@ -31,7 +31,7 @@ class RungeKuttaStepper public: //!@{ //! Type aliases - using result_type = StepperResult; + using result_type = FieldStepperResult; //!@} public: diff --git a/src/celeritas/field/Types.hh b/src/celeritas/field/Types.hh index 507cba9457..4dce7eeb94 100644 --- a/src/celeritas/field/Types.hh +++ b/src/celeritas/field/Types.hh @@ -25,15 +25,15 @@ struct OdeState using MomentumUnits = units::MevMomentum; using Real3 = Array; - 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 diff --git a/src/celeritas/field/ZHelixStepper.hh b/src/celeritas/field/ZHelixStepper.hh index 2835f0cfcd..374189b537 100644 --- a/src/celeritas/field/ZHelixStepper.hh +++ b/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::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::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::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; } diff --git a/src/celeritas/field/detail/FieldUtils.hh b/src/celeritas/field/detail/FieldUtils.hh index fddabf8212..c6ca923193 100644 --- a/src/celeritas/field/detail/FieldUtils.hh +++ b/src/celeritas/field/detail/FieldUtils.hh @@ -43,7 +43,7 @@ template inline CELER_FUNCTION Array ax(T a, const Array& x) { Array 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 ax(T a, const Array& 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]; diff --git a/src/celeritas/global/alongstep/detail/UrbanMsc.hh b/src/celeritas/global/alongstep/detail/UrbanMsc.hh index 0571faffe2..d6ddc52e17 100644 --- a/src/celeritas/global/alongstep/detail/UrbanMsc.hh +++ b/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]; } diff --git a/src/celeritas/random/distribution/GammaDistribution.hh b/src/celeritas/random/distribution/GammaDistribution.hh index 11a813164f..ab5f0ece49 100644 --- a/src/celeritas/random/distribution/GammaDistribution.hh +++ b/src/celeritas/random/distribution/GammaDistribution.hh @@ -86,7 +86,7 @@ GammaDistribution::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); diff --git a/src/corecel/math/Algorithms.hh b/src/corecel/math/Algorithms.hh index be0bfc949f..d5550a22f2 100644 --- a/src/corecel/math/Algorithms.hh +++ b/src/corecel/math/Algorithms.hh @@ -10,6 +10,7 @@ #include #include +#include "celeritas_config.h" #include "corecel/Assert.hh" #include "corecel/Macros.hh" @@ -316,5 +317,27 @@ inline CELER_FUNCTION T fastpow(T a, T b) return std::exp(b * std::log(a)); } +#ifdef __CUDACC__ +using ::rsqrt; +#else +//---------------------------------------------------------------------------// +/*! + * Calculate an inverse square root. + */ +inline CELER_FUNCTION double rsqrt(double value) +{ + return 1.0 / std::sqrt(value); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate an inverse square root. + */ +inline CELER_FUNCTION float rsqrt(float value) +{ + return 1.0f / std::sqrt(value); +} +#endif + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 526ee59f74..1003192efd 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -280,6 +280,7 @@ celeritas_add_test(celeritas/field/Steppers.test.cc) celeritas_add_test(celeritas/field/FieldDriver.test.cc) celeritas_add_test(celeritas/field/FieldPropagator.test.cc ${_needs_geo}) celeritas_add_test(celeritas/field/LinearPropagator.test.cc ${_needs_geo}) +celeritas_add_test(celeritas/field/MagFieldEquation.test.cc) #-------------------------------------# # Geo diff --git a/test/celeritas/field/FieldDriver.test.cc b/test/celeritas/field/FieldDriver.test.cc index 99ff4e1f20..44e74ff4c1 100644 --- a/test/celeritas/field/FieldDriver.test.cc +++ b/test/celeritas/field/FieldDriver.test.cc @@ -14,7 +14,6 @@ #include "celeritas/Constants.hh" #include "celeritas/Quantities.hh" #include "celeritas/field/DormandPrinceStepper.hh" -#include "celeritas/field/FieldDriver.hh" #include "celeritas/field/FieldDriverOptions.hh" #include "celeritas/field/MagFieldEquation.hh" #include "celeritas/field/Types.hh" @@ -50,7 +49,7 @@ class FieldDriverTest : public Test protected: // Field parameters - FieldDriverOptions field_params; + FieldDriverOptions driver_options; // Test parameters FieldTestParams test_params; @@ -71,16 +70,16 @@ make_mag_field_driver(FieldT&& field, options, Stepper_t{Equation_t{::celeritas::forward(field), charge}}}; } + //---------------------------------------------------------------------------// // TESTS //---------------------------------------------------------------------------// -TEST_F(FieldDriverTest, field_driver_host) +TEST_F(FieldDriverTest, types) { - // Construct FieldDriver auto driver = make_mag_field_driver( UniformField({0, 0, test_params.field_value}), - field_params, + driver_options, units::ElementaryCharge{-1}); // Make sure object is holding things by value @@ -91,6 +90,14 @@ TEST_F(FieldDriverTest, field_driver_host) // Size: field vector, q / c, reference to options EXPECT_EQ(sizeof(Real3) + sizeof(real_type) + sizeof(FieldDriverOptions*), sizeof(driver)); +} + +TEST_F(FieldDriverTest, revolutions) +{ + auto driver = make_mag_field_driver( + UniformField({0, 0, test_params.field_value}), + driver_options, + units::ElementaryCharge{-1}); // Test parameters and the sub-step size real_type circumference = 2 * constants::pi * test_params.radius; @@ -106,7 +113,7 @@ TEST_F(FieldDriverTest, field_driver_host) real_type total_step_length{0}; // Try the stepper by hstep for (num_revolutions * num_steps) times - real_type delta = field_params.errcon; + real_type delta = driver_options.errcon; for (int nr = 0; nr < test_params.revolutions; ++nr) { y_expected.pos @@ -129,11 +136,11 @@ TEST_F(FieldDriverTest, field_driver_host) total_step_length, circumference * test_params.revolutions, delta); } -TEST_F(FieldDriverTest, accurate_advance_host) +TEST_F(FieldDriverTest, accurate_advance) { auto driver = make_mag_field_driver( UniformField({0, 0, test_params.field_value}), - field_params, + driver_options, units::ElementaryCharge{-1}); // Test parameters and the sub-step size @@ -149,7 +156,7 @@ TEST_F(FieldDriverTest, accurate_advance_host) // Try the stepper by hstep for (num_revolutions * num_steps) times real_type total_curved_length{0}; - real_type delta = field_params.errcon; + real_type delta = driver_options.errcon; for (int nr = 0; nr < test_params.revolutions; ++nr) { @@ -171,4 +178,4 @@ TEST_F(FieldDriverTest, accurate_advance_host) // Check the total error, step/curve length EXPECT_LT(total_curved_length - circumference * test_params.revolutions, delta); -} \ No newline at end of file +} diff --git a/test/celeritas/field/MagFieldEquation.test.cc b/test/celeritas/field/MagFieldEquation.test.cc new file mode 100644 index 0000000000..d09812fafc --- /dev/null +++ b/test/celeritas/field/MagFieldEquation.test.cc @@ -0,0 +1,85 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2022 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/field/MagFieldEquation.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/field/MagFieldEquation.hh" + +#include "celeritas_test.hh" + +using celeritas::units::ElementaryCharge; + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// + +Real3 dummy_field(const Real3& pos) +{ + // Rotate and scale + Real3 result; + result[0] = real_type(0.5) * pos[1]; + result[1] = real_type(1.0) * pos[2]; + result[2] = real_type(2.0) * pos[0]; + return result; +} + +template +CELER_FUNCTION decltype(auto) +make_equation(FieldT&& field, ElementaryCharge charge) +{ + using Equation_t = celeritas::MagFieldEquation; + return Equation_t{::celeritas::forward(field), charge}; +} + +void print_expected(const OdeState& s) +{ + cout << "/*** BEGIN CODE ***/\n" + << "EXPECT_VEC_SOFT_EQ(Real3(" << repr(s.pos) << "), result.pos);\n" + << "EXPECT_VEC_SOFT_EQ(Real3(" << repr(s.mom) << "), result.mom);\n" + << "/*** END CODE ***/\n"; +} + +//---------------------------------------------------------------------------// + +TEST(MagFieldEquationTest, charged) +{ + auto eval = make_equation(dummy_field, ElementaryCharge{3}); + { + OdeState result = eval({{1, 2, 3}, {0, 0, 1}}); + EXPECT_VEC_SOFT_EQ(Real3({0, 0, 1}), result.pos); + EXPECT_VEC_SOFT_EQ(Real3({-0.02698132122, 0.00899377374, 0}), + result.mom); + } + { + OdeState result = eval({{0.5, -2, -1}, {1, 2, 3}}); + EXPECT_VEC_SOFT_EQ( + Real3({0.26726124191242, 0.53452248382485, 0.80178372573727}), + result.pos); + EXPECT_VEC_SOFT_EQ( + Real3({0.012018435696159, -0.009614748556927, 0.0024036871392318}), + result.mom); + } + if (CELERITAS_DEBUG) + { + EXPECT_THROW(eval({{0.5, -2, -1}, {0, 0, 0}}), DebugError); + } +} + +TEST(MagFieldEquationTest, neutral) +{ + auto eval = make_equation(dummy_field, ElementaryCharge{0}); + { + OdeState s{{0.5, -2, -1}, {0, 0, 1}}; + OdeState result = eval(s); + EXPECT_VEC_SOFT_EQ(Real3({0, 0, 1}), result.pos); + EXPECT_VEC_SOFT_EQ(Real3({0, 0, 0}), result.mom); + } +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/field/Steppers.test.cc b/test/celeritas/field/Steppers.test.cc index 6dd413dc9b..2f4aef9401 100644 --- a/test/celeritas/field/Steppers.test.cc +++ b/test/celeritas/field/Steppers.test.cc @@ -93,7 +93,7 @@ class SteppersTest : public Test expected_y.pos[2] = param.delta_z * (nr + 1) + i * 1.0e-6; for (CELER_MAYBE_UNUSED int j : celeritas::range(param.nsteps)) { - StepperResult result = stepper(hstep, y); + FieldStepperResult result = stepper(hstep, y); y = result.end_state; total_err2 diff --git a/test/corecel/math/Algorithms.test.cc b/test/corecel/math/Algorithms.test.cc index 274d0c8098..7a034a5134 100644 --- a/test/corecel/math/Algorithms.test.cc +++ b/test/corecel/math/Algorithms.test.cc @@ -215,3 +215,20 @@ TEST(MathTest, fastpow) EXPECT_TRUE((std::is_same::value)); } + +//---------------------------------------------------------------------------// + +TEST(MathTest, rsqrt) +{ + using celeritas::rsqrt; + + constexpr auto dblinf = std::numeric_limits::infinity(); + EXPECT_DOUBLE_EQ(0.5, rsqrt(4.0)); + EXPECT_DOUBLE_EQ(dblinf, rsqrt(0.0)); + EXPECT_DOUBLE_EQ(0.0, rsqrt(dblinf)); + + constexpr auto fltinf = std::numeric_limits::infinity(); + EXPECT_FLOAT_EQ(0.5f, rsqrt(4.0f)); + EXPECT_FLOAT_EQ(fltinf, rsqrt(0.0f)); + EXPECT_FLOAT_EQ(0.0f, rsqrt(fltinf)); +}