Core/Propagator/GenericRungeKuttaStepperTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Definitions/Direction.hpp"
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/Definitions/Units.hpp"
0016 #include "Acts/EventData/BoundTrackParameters.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/MagneticField/ConstantBField.hpp"
0019 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0020 #include "Acts/MagneticField/MagneticFieldError.hpp"
0021 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0022 #include "Acts/Propagator/ButcherTableau.hpp"
0023 #include "Acts/Propagator/ConstrainedStep.hpp"
0024 #include "Acts/Propagator/GenericRungeKuttaStepper.hpp"
0025 #include "Acts/Propagator/HelixStepper.hpp"
0026 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0027
0028 #include <cmath>
0029 #include <functional>
0030 #include <limits>
0031 #include <memory>
0032 #include <optional>
0033 #include <ostream>
0034 #include <stdexcept>
0035 #include <vector>
0036
0037 using namespace Acts;
0038 using namespace Acts::UnitLiterals;
0039 using Acts::VectorHelpers::makeVector4;
0040
0041 // The test contexts print the name of the tableau
0042 BOOST_TEST_DONT_PRINT_LOG_VALUE(std::shared_ptr<const Acts::ButcherTableau>)
0043 BOOST_TEST_DONT_PRINT_LOG_VALUE(Acts::GenericRungeKuttaStepper::ErrorEstimation)
0044
0045 namespace ActsTests {
0046
0047 namespace {
0048
0049 const auto tgContext = GeometryContext::dangerouslyDefaultConstruct();
0050 const MagneticFieldContext mfContext;
0051
0052 constexpr double infinity = std::numeric_limits<double>::infinity();
0053
0054 /// A smooth, non-uniform field with an optional analytic gradient
0055 ///
0056 ///   B = (s a y z, -s a x z, B0 + s a (x^2 + y^2)), s = B0 / L^2, a = 1/2
0057 class SmoothField final : public MagneticFieldProvider {
0058  public:
0059   struct Cache {
0060     explicit Cache(const MagneticFieldContext& /*mctx*/) {}
0061   };
0062
0063   SmoothField(double b0, double length, bool withGradient = false)
0064       : m_b0(b0), m_length(length), m_withGradient(withGradient) {}
0065
0066   MagneticFieldProvider::Cache makeCache(
0067       const MagneticFieldContext& mctx) const override {
0068     return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
0069   }
0070
0071   Result<Vector3> getField(
0072       const Vector3& position,
0073       MagneticFieldProvider::Cache& /*cache*/) const override {
0074     return Result<Vector3>::success(evaluate(position).field);
0075   }
0076
0077   bool providesFieldGradient() const override { return m_withGradient; }
0078
0079   Result<FieldAndGradient> getFieldAndGradient(
0080       const Vector3& position,
0081       MagneticFieldProvider::Cache& /*cache*/) const override {
0082     ++m_gradientLookups;
0083     if (!m_withGradient) {
0084       return Result<FieldAndGradient>::failure(
0085           MagneticFieldError::NotImplemented);
0086     }
0087     return Result<FieldAndGradient>::success(evaluate(position));
0088   }
0089
0090   /// Number of calls of getFieldAndGradient
0091   std::size_t gradientLookups() const { return m_gradientLookups; }
0092
0093  private:
0094   FieldAndGradient evaluate(const Vector3& p) const {
0095     constexpr double a = 0.5;
0096     const double s = m_b0 / (m_length * m_length);
0097     const double x = p.x();
0098     const double y = p.y();
0099     const double z = p.z();
0100     FieldAndGradient result;
0101     result.field =
0102         Vector3(s * a * y * z, -s * a * x * z, m_b0 + s * a * (x * x + y * y));
0103     // clang-format off
0104     result.gradient <<          0.,  s * a * z,  s * a * y,
0105                         -s * a * z,         0., -s * a * x,
0106                      2 * s * a * x, 2 * s * a * y,        0.;
0107     // clang-format on
0108     return result;
0109   }
0110
0111   double m_b0;
0112   double m_length;
0113   bool m_withGradient;
0114   mutable std::size_t m_gradientLookups = 0;
0115 };
0116
0117 /// A constant field that fails outside of |x| < xMax
0118 class BoundedField final : public MagneticFieldProvider {
0119  public:
0120   struct Cache {
0121     explicit Cache(const MagneticFieldContext& /*mctx*/) {}
0122   };
0123
0124   BoundedField(const Vector3& field, double xMax)
0125       : m_field(field), m_xMax(xMax) {}
0126
0127   MagneticFieldProvider::Cache makeCache(
0128       const MagneticFieldContext& mctx) const override {
0129     return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
0130   }
0131
0132   Result<Vector3> getField(
0133       const Vector3& p,
0134       MagneticFieldProvider::Cache& /*cache*/) const override {
0135     if (std::abs(p.x()) >= m_xMax) {
0136       return Result<Vector3>::failure(MagneticFieldError::OutOfBounds);
0137     }
0138     return Result<Vector3>::success(m_field);
0139   }
0140
0141  private:
0142   Vector3 m_field;
0143   double m_xMax;
0144 };
0145
0146 /// A constant field that counts its lookups
0147 class CountingField final : public MagneticFieldProvider {
0148  public:
0149   struct Cache {
0150     explicit Cache(const MagneticFieldContext& /*mctx*/) {}
0151   };
0152
0153   explicit CountingField(const Vector3& field) : m_field(field) {}
0154
0155   MagneticFieldProvider::Cache makeCache(
0156       const MagneticFieldContext& mctx) const override {
0157     return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
0158   }
0159
0160   Result<Vector3> getField(
0161       const Vector3& /*p*/,
0162       MagneticFieldProvider::Cache& /*cache*/) const override {
0163     ++m_lookups;
0164     return Result<Vector3>::success(m_field);
0165   }
0166
0167   std::size_t lookups() const { return m_lookups; }
0168
0169  private:
0170   Vector3 m_field;
0171   mutable std::size_t m_lookups = 0;
0172 };
0173
0174 FreeVector makeStart(const Vector3& pos, const Vector3& dir, double qop) {
0175   FreeVector start = FreeVector::Zero();
0176   start.segment<3>(eFreePos0) = pos;
0177   start[eFreeTime] = 1.;
0178   start.segment<3>(eFreeDir0) = dir.normalized();
0179   start[eFreeQOverP] = qop;
0180   return start;
0181 }
0182
0183 /// Propagate the free parameters over the path length and return the state.
0184 template <typename stepper_t>
0185 typename stepper_t::State propagateFree(
0186     const stepper_t& stepper, typename stepper_t::Options options,
0187     const FreeVector& start, double pathLength, bool covTransport,
0188     const ParticleHypothesis& particleHypothesis = ParticleHypothesis::pion()) {
0189   typename stepper_t::State state = stepper.makeState(options);
0190   // Initialise from the start parameters, so that the direction columns of
0191   // the bound-to-free jacobian are orthogonal to the start direction.
0192   stepper.initialize(
0193       state,
0194       BoundTrackParameters::createCurvilinear(
0195           makeVector4(start.segment<3>(eFreePos0), start[eFreeTime]),
0196           start.segment<3>(eFreeDir0), start[eFreeQOverP],
0197           covTransport ? std::optional<BoundMatrix>(BoundMatrix::Identity())
0198                        : std::nullopt,
0199           particleHypothesis));
0200   // Keep the exact start values, the round trip through the bound
0201   // parameters changes them at the level of the round-off.
0202   state.pars = start;
0203   while (std::abs(pathLength - state.pathAccumulated) > 1e-12) {
0204     state.stepSize.release(ConstrainedStep::Type::Navigator);
0205     state.stepSize.update(pathLength - state.pathAccumulated,
0206                           ConstrainedStep::Type::Navigator);
0207     auto res = stepper.step(state, Direction::fromScalar(pathLength), nullptr);
0208     BOOST_REQUIRE(res.ok());
0209   }
0210   return state;
0211 }
0212
0213 GenericRungeKuttaStepper::Options fixedStepOptions(double stepSize) {
0214   GenericRungeKuttaStepper::Options options(tgContext, mfContext);
0215   options.errorEstimation = GenericRungeKuttaStepper::ErrorEstimation::None;
0216   options.initialStepSize = stepSize;
0217   options.maxStepSize = stepSize;
0218   return options;
0219 }
0220
0221 }  // namespace
0222
0223 BOOST_AUTO_TEST_SUITE(PropagatorSuite)
0224
0225 BOOST_AUTO_TEST_CASE(butcher_tableau_validation) {
0226   BOOST_CHECK_THROW(ButcherTableau("sizes", 1, 0, {0., 1.}, {{}}, {1.}, {}),
0227                     std::invalid_argument);
0228   BOOST_CHECK_THROW(
0229       ButcherTableau("node", 1, 0, {0., 1.}, {{}, {0.5}}, {0.5, 0.5}, {}),
0230       std::invalid_argument);
0231   BOOST_CHECK_THROW(ButcherTableau("first", 1, 0, {1.}, {{}}, {1.}, {}),
0232                     std::invalid_argument);
0233   BOOST_CHECK_NO_THROW(
0234       ButcherTableau("heun", 2, 1, {0., 1.}, {{}, {1.}}, {0.5, 0.5}, {1., 0.}));
0235 }
0236
0237 BOOST_AUTO_TEST_CASE(generic_runge_kutta_stepper_without_tableau) {
0238   auto field = std::make_shared<ConstantBField>(Vector3(0., 0., 2_T));
0239   BOOST_CHECK_THROW(GenericRungeKuttaStepper(
0240                         GenericRungeKuttaStepper::Config{field, nullptr}),
0241                     std::invalid_argument);
0242 }
0243
0244 /// Check the order conditions of the built-in tableaus up to order 5. The
0245 /// Verner 9(8) coefficients were checked against all conditions up to order 9
0246 /// with 60-digit arithmetic when they were added.
0247 BOOST_AUTO_TEST_CASE(butcher_tableau_order_conditions) {
0248   using Vec = Eigen::VectorXd;
0249   using Mat = Eigen::MatrixXd;
0250
0251   // Value of each condition for the weights w, with its order and its
0252   // expected value 1 / gamma(tree)
0253   struct Condition {
0254     unsigned order;
0255     double expected;
0256     std::function<double(const Vec&, const Vec&, const Mat&)> value;
0257   };
0258   const std::vector<Condition> conditions = {
0259       {1, 1., [](auto& w, auto&, auto&) { return w.sum(); }},
0260       {2, 1. / 2, [](auto& w, auto& c, auto&) { return w.dot(c); }},
0261       {3, 1. / 3,
0262        [](auto& w, auto& c, auto&) { return w.dot(Vec(c.cwiseProduct(c))); }},
0263       {3, 1. / 6, [](auto& w, auto& c, auto& a) { return w.dot(Vec(a * c)); }},
0264       {4, 1. / 4,
0265        [](auto& w, auto& c, auto&) {
0266          return w.dot(Vec(c.array().pow(3).matrix()));
0267        }},
0268       {4, 1. / 8,
0269        [](auto& w, auto& c, auto& a) {
0270          return w.dot(Vec(c.cwiseProduct(a * c)));
0271        }},
0272       {4, 1. / 12,
0273        [](auto& w, auto& c, auto& a) {
0274          return w.dot(Vec(a * c.cwiseProduct(c)));
0275        }},
0276       {4, 1. / 24,
0277        [](auto& w, auto& c, auto& a) { return w.dot(Vec(a * a * c)); }},
0278       {5, 1. / 5,
0279        [](auto& w, auto& c, auto&) {
0280          return w.dot(Vec(c.array().pow(4).matrix()));
0281        }},
0282       {5, 1. / 10,
0283        [](auto& w, auto& c, auto& a) {
0284          return w.dot(Vec(c.cwiseProduct(c).cwiseProduct(a * c)));
0285        }},
0286       {5, 1. / 20,
0287        [](auto& w, auto& c, auto& a) {
0288          return w.dot(Vec((a * c).cwiseProduct(a * c)));
0289        }},
0290       {5, 1. / 15,
0291        [](auto& w, auto& c, auto& a) {
0292          return w.dot(Vec(c.cwiseProduct(a * c.cwiseProduct(c))));
0293        }},
0294       {5, 1. / 30,
0295        [](auto& w, auto& c, auto& a) {
0296          return w.dot(Vec(c.cwiseProduct(a * a * c)));
0297        }},
0298       {5, 1. / 20,
0299        [](auto& w, auto& c, auto& a) {
0300          return w.dot(Vec(a * c.array().pow(3).matrix()));
0301        }},
0302       {5, 1. / 40,
0303        [](auto& w, auto& c, auto& a) {
0304          return w.dot(Vec(a * c.cwiseProduct(a * c)));
0305        }},
0306       {5, 1. / 60,
0307        [](auto& w, auto& c, auto& a) {
0308          return w.dot(Vec(a * a * c.cwiseProduct(c)));
0309        }},
0310       {5, 1. / 120,
0311        [](auto& w, auto& c, auto& a) { return w.dot(Vec(a * a * a * c)); }},
0312   };
0313
0314   for (const auto& tableau :
0315        {ButcherTableau::classicalRk4(), ButcherTableau::dormandPrince54(),
0316         ButcherTableau::verner98()}) {
0317     const std::size_t s = tableau->stages();
0318     Vec c(s);
0319     Vec b(s);
0320     Vec bEmbedded = Vec::Zero(s);
0321     Mat a = Mat::Zero(s, s);
0322     for (std::size_t i = 0; i < s; ++i) {
0323       c[i] = tableau->c(i);
0324       b[i] = tableau->b(i);
0325       if (tableau->hasEmbedded()) {
0326         bEmbedded[i] = tableau->bEmbedded(i);
0327       }
0328       for (std::size_t j = 0; j < i; ++j) {
0329         a(i, j) = tableau->a(i, j);
0330       }
0331     }
0332
0333     BOOST_TEST_CONTEXT(tableau->name()) {
0334       for (const auto& condition : conditions) {
0335         if (condition.order <= tableau->order()) {
0336           CHECK_CLOSE_ABS(condition.value(b, c, a), condition.expected, 1e-14);
0337         }
0338         if (tableau->hasEmbedded() &&
0339             condition.order <= tableau->embeddedOrder()) {
0340           CHECK_CLOSE_ABS(condition.value(bEmbedded, c, a), condition.expected,
0341                           1e-14);
0342         }
0343       }
0344       // The solution does not have a higher order than stated
0345       if (tableau->order() == 4) {
0346         BOOST_CHECK_GT(std::abs(conditions[8].value(b, c, a) - 1. / 5), 1e-6);
0347       }
0348       // The embedded solution differs from the solution at the next order
0349       if (tableau->hasEmbedded() && tableau->embeddedOrder() < 5) {
0350         bool differs = false;
0351         for (const auto& condition : conditions) {
0352           if (condition.order == tableau->embeddedOrder() + 1) {
0353             differs |= std::abs(condition.value(bEmbedded, c, a) -
0354                                 condition.expected) > 1e-6;
0355           }
0356         }
0357         BOOST_CHECK(differs);
0358       }
0359     }
0360   }
0361 }
0362
0363 /// In a constant field the reference stepper reproduces the exact helix.
0364 BOOST_AUTO_TEST_CASE(generic_runge_kutta_stepper_constant_field) {
0365   auto field = std::make_shared<ConstantBField>(Vector3(0.1_T, -0.2_T, 2_T));
0366   const FreeVector start =
0367       makeStart(Vector3(1., 2., 3.), Vector3(4., -5., 6.), -1. / 0.7_GeV);
0368   const double pathLength = 1_m;
0369
0370   HelixStepper helix(field);
0371   HelixStepper::Options helixOptions(tgContext, mfContext);
0372   auto helixState = propagateFree(helix, helixOptions, start, pathLength, true);
0373
0374   GenericRungeKuttaStepper reference(field);
0375   GenericRungeKuttaStepper::Options options(tgContext, mfContext);
0376   options.stepTolerance = 1e-12;
0377   options.initialStepSize = 1_cm;
0378   auto state = propagateFree(reference, options, start, pathLength, true);
0379
0380   CHECK_CLOSE_ABS(state.pars, helixState.pars, 1e-9);
0381   CHECK_CLOSE_ABS(state.derivative, helixState.derivative, 1e-9);
0382
0383   // The reference projects the jacobian of the direction onto the unit sphere
0384   // after each step, so compare the curvilinear jacobians.
0385   const BoundMatrix refJac = reference.transportToCurvilinear(state);
0386   const BoundMatrix helixJac = helix.transportToCurvilinear(helixState);
0387   CHECK_CLOSE_OR_SMALL(refJac, helixJac, 1e-8, 1e-8);
0388 }
0389
0390 namespace {
0391
0392 struct ConvergenceCase {
0393   std::shared_ptr<const ButcherTableau> tableau;
0394   double qop{};
0395   std::vector<int> stepCounts;
0396
0397   friend std::ostream& operator<<(std::ostream& os,
0398                                   const ConvergenceCase& testCase) {
0399     return os << testCase.tableau->name();
0400   }
0401 };
0402
0403 // The ninth-order error needs a strongly bent track to stay above the
0404 // round-off while the step size halves
0405 const std::vector<ConvergenceCase> convergenceCases = {
0406     {ButcherTableau::classicalRk4(), 1. / 1_GeV, {8, 16, 32}},
0407     {ButcherTableau::dormandPrince54(), 1. / 1_GeV, {8, 16, 32}},
0408     {ButcherTableau::verner98(), 1. / 0.1_GeV, {8, 12, 18}},
0409 };
0410
0411 }  // namespace
0412
0413 /// The fixed-step global error falls with the order of the tableau.
0414 BOOST_DATA_TEST_CASE(generic_runge_kutta_stepper_convergence,
0415                      boost::unit_test::data::make(convergenceCases), testCase) {
0416   auto field = std::make_shared<SmoothField>(2_T, 1_m);
0417   const FreeVector start =
0418       makeStart(Vector3(100., -50., 20.), Vector3(1., 0.3, 0.2), testCase.qop);
0419   const double pathLength = 1_m;
0420
0421   // Every consistent tableau converges to the same solution, so a wrong order
0422   // still shows in the ratios.
0423   GenericRungeKuttaStepper::Config config{field, ButcherTableau::verner98()};
0424   const GenericRungeKuttaStepper reference(config);
0425   const FreeVector referencePars =
0426       propagateFree(reference, fixedStepOptions(pathLength / 1024), start,
0427                     pathLength, false)
0428           .pars;
0429
0430   config.tableau = testCase.tableau;
0431   const GenericRungeKuttaStepper stepper(config);
0432   std::vector<double> errors;
0433   for (int nSteps : testCase.stepCounts) {
0434     const FreeVector end =
0435         propagateFree(stepper, fixedStepOptions(pathLength / nSteps), start,
0436                       pathLength, false)
0437             .pars;
0438     errors.push_back(
0439         (end - referencePars).segment<3>(eFreePos0).cwiseAbs().maxCoeff());
0440   }
0441
0442   BOOST_TEST_CONTEXT(testCase.tableau->name()
0443                      << " errors " << errors[0] << " " << errors[1] << " "
0444                      << errors[2]) {
0445     for (std::size_t i = 0; i + 1 < errors.size(); ++i) {
0446       const double ratio = errors[i] / errors[i + 1];
0447       const double expectedRatio =
0448           std::pow(static_cast<double>(testCase.stepCounts[i + 1]) /
0449                        testCase.stepCounts[i],
0450                    testCase.tableau->order());
0451       BOOST_CHECK_GT(ratio, 0.6 * expectedRatio);
0452       BOOST_CHECK_LT(ratio, 1.6 * expectedRatio);
0453     }
0454   }
0455 }
0456
0457 /// The chain-rule jacobian matches finite differences of the full
0458 /// propagation, and it needs the field gradient to do so.
0459 BOOST_AUTO_TEST_CASE(generic_runge_kutta_stepper_jacobian) {
0460   const FreeVector start =
0461       makeStart(Vector3(100., -50., 20.), Vector3(1., 0.3, 0.2), -1. / 1_GeV);
0462   const double pathLength = 1_m;
0463   const auto options = fixedStepOptions(pathLength / 16);
0464
0465   auto field = std::make_shared<SmoothField>(2_T, 1_m);
0466   const GenericRungeKuttaStepper stepper(field);
0467
0468   FreeMatrix numeric;
0469   for (std::size_t j = 0; j < eFreeSize; ++j) {
0470     // A larger delta keeps the round-off of positions of order 1 m small
0471     const double delta =
0472         j == eFreeQOverP ? 1e-4 * std::abs(start[eFreeQOverP]) : 1e-4;
0473     FreeVector plus = start;
0474     FreeVector minus = start;
0475     plus[j] += delta;
0476     minus[j] -= delta;
0477     numeric.col(j) =
0478         (propagateFree(stepper, options, plus, pathLength, false).pars -
0479          propagateFree(stepper, options, minus, pathLength, false).pars) /
0480         (2. * delta);
0481   }
0482
0483   const auto analytic =
0484       propagateFree(stepper, options, start, pathLength, true);
0485   // Both normalise the end direction, so the projected jacobian compares
0486   // directly.
0487   CHECK_CLOSE_OR_SMALL(analytic.jacTransport, numeric, 1e-5, 1e-5);
0488
0489   // Without the gradient term the jacobian is wrong
0490   auto noGradientOptions = options;
0491   noGradientOptions.includeFieldGradient = false;
0492   const auto noGradient =
0493       propagateFree(stepper, noGradientOptions, start, pathLength, true);
0494   BOOST_CHECK_GT((noGradient.jacTransport - numeric).cwiseAbs().maxCoeff(),
0495                  1e-3);
0496
0497   // The analytic gradient of the provider gives the same jacobian as the
0498   // finite differences
0499   auto analyticField = std::make_shared<SmoothField>(2_T, 1_m, true);
0500   const GenericRungeKuttaStepper analyticStepper(analyticField);
0501   const auto withProviderGradient =
0502       propagateFree(analyticStepper, options, start, pathLength, true);
0503   CHECK_CLOSE_OR_SMALL(withProviderGradient.jacTransport, analytic.jacTransport,
0504                        1e-6, 1e-8);
0505   // Only the stepper with the analytic field uses the provider gradient
0506   BOOST_CHECK_GT(analyticField->gradientLookups(), 0u);
0507   BOOST_CHECK_EQUAL(field->gradientLookups(), 0u);
0508 }
0509
0510 namespace {
0511
0512 using ErrorEstimation = GenericRungeKuttaStepper::ErrorEstimation;
0513
0514 struct AdaptiveCase {
0515   std::shared_ptr<const ButcherTableau> tableau;
0516   ErrorEstimation errorEstimation{};
0517
0518   friend std::ostream& operator<<(std::ostream& os,
0519                                   const AdaptiveCase& testCase) {
0520     return os << testCase.tableau->name() << " "
0521               << (testCase.errorEstimation == ErrorEstimation::Embedded
0522                       ? "embedded"
0523                       : "step doubling");
0524   }
0525 };
0526
0527 const std::vector<AdaptiveCase> adaptiveCases = {
0528     {ButcherTableau::dormandPrince54(), ErrorEstimation::Embedded},
0529     {ButcherTableau::verner98(), ErrorEstimation::Embedded},
0530     {ButcherTableau::classicalRk4(), ErrorEstimation::StepDoubling},
0531     {ButcherTableau::dormandPrince54(), ErrorEstimation::StepDoubling},
0532 };
0533
0534 }  // namespace
0535
0536 /// The adaptive step size meets the tolerance.
0537 BOOST_DATA_TEST_CASE(generic_runge_kutta_stepper_adaptive,
0538                      boost::unit_test::data::make(adaptiveCases), testCase) {
0539   auto field = std::make_shared<SmoothField>(2_T, 1_m);
0540   const FreeVector start =
0541       makeStart(Vector3(100., -50., 20.), Vector3(1., 0.3, 0.2), 1. / 1_GeV);
0542   const double pathLength = 1_m;
0543
0544   GenericRungeKuttaStepper::Config config{field, ButcherTableau::verner98()};
0545   // There is no closed-form solution in this field. The reference is the
0546   // fixed-step Verner 9(8) solution with 1024 steps. Its error is many orders
0547   // of magnitude below the tolerances, also for the adaptive Verner 9(8)
0548   // case, which takes far fewer steps.
0549   const FreeVector referencePars =
0550       propagateFree(GenericRungeKuttaStepper(config),
0551                     fixedStepOptions(pathLength / 1024), start, pathLength,
0552                     false)
0553           .pars;
0554
0555   config.tableau = testCase.tableau;
0556   const GenericRungeKuttaStepper stepper(config);
0557   GenericRungeKuttaStepper::Options options(tgContext, mfContext);
0558   options.errorEstimation = testCase.errorEstimation;
0559   options.initialStepSize = 10_m;
0560   double previousError = infinity;
0561   for (double tolerance : {1e-4, 1e-7, 1e-10}) {
0562     options.stepTolerance = tolerance;
0563     const auto state =
0564         propagateFree(stepper, options, start, pathLength, false);
0565     const double error = (state.pars - referencePars)
0566                              .segment<3>(eFreePos0)
0567                              .cwiseAbs()
0568                              .maxCoeff();
0569     BOOST_TEST_CONTEXT(testCase << " tolerance " << tolerance << " error "
0570                                 << error << " steps " << state.nSteps) {
0571       BOOST_CHECK_GT(state.statistics.nRejectedSteps, 0u);
0572       // The global error is the sum of the local errors
0573       BOOST_CHECK_LT(error, 10. * tolerance * state.nSteps);
0574       BOOST_CHECK_LT(error, previousError);
0575     }
0576     previousError = error;
0577   }
0578 }
0579
0580 /// The embedded error estimate needs a tableau with embedded weights.
0581 BOOST_AUTO_TEST_CASE(generic_runge_kutta_stepper_embedded_without_weights) {
0582   auto field = std::make_shared<SmoothField>(2_T, 1_m);
0583   const GenericRungeKuttaStepper stepper(
0584       GenericRungeKuttaStepper::Config{field, ButcherTableau::classicalRk4()});
0585   GenericRungeKuttaStepper::Options options(tgContext, mfContext);
0586   BOOST_REQUIRE(options.errorEstimation == ErrorEstimation::Embedded);
0587   BOOST_CHECK_THROW(stepper.makeState(options), std::invalid_argument);
0588
0589   // The options in the state can change after makeState
0590   options.errorEstimation = ErrorEstimation::StepDoubling;
0591   auto state = stepper.makeState(options);
0592   stepper.initialize(state, BoundTrackParameters::createCurvilinear(
0593                                 Vector4::Zero(), Vector3::UnitX(), 1. / 1_GeV,
0594                                 std::nullopt, ParticleHypothesis::pion()));
0595   state.options.errorEstimation = ErrorEstimation::Embedded;
0596   BOOST_CHECK_THROW(
0597       static_cast<void>(stepper.step(state, Direction::Forward(), nullptr)),
0598       std::invalid_argument);
0599 }
0600
0601 /// Without rejections, step doubling equals fixed steps of h/2.
0602 BOOST_AUTO_TEST_CASE(generic_runge_kutta_stepper_step_doubling_halves) {
0603   auto field = std::make_shared<SmoothField>(2_T, 1_m);
0604   const GenericRungeKuttaStepper stepper(
0605       GenericRungeKuttaStepper::Config{field, ButcherTableau::classicalRk4()});
0606   const FreeVector start =
0607       makeStart(Vector3(100., -50., 20.), Vector3(1., 0.3, 0.2), 1. / 1_GeV);
0608
0609   GenericRungeKuttaStepper::Options options(tgContext, mfContext);
0610   options.errorEstimation = ErrorEstimation::StepDoubling;
0611   options.stepTolerance = infinity;
0612   options.initialStepSize = 10_cm;
0613   options.maxStepSize = 10_cm;
0614   const auto doubling = propagateFree(stepper, options, start, 1_m, true);
0615   const auto halves =
0616       propagateFree(stepper, fixedStepOptions(5_cm), start, 1_m, true);
0617
0618   BOOST_CHECK_EQUAL(doubling.nSteps, 10u);
0619   BOOST_CHECK_EQUAL(halves.nSteps, 20u);
0620   BOOST_CHECK_EQUAL(doubling.statistics.nRejectedSteps, 0u);
0621   CHECK_CLOSE_ABS(doubling.pars, halves.pars, 1e-12);
0622   CHECK_CLOSE_OR_SMALL(doubling.jacTransport, halves.jacTransport, 1e-12,
0623                        1e-12);
0624 }
0625
0626 /// An FSAL tableau reuses the field of its last stage at the end.
0627 BOOST_AUTO_TEST_CASE(generic_runge_kutta_stepper_first_same_as_last) {
0628   BOOST_CHECK(ButcherTableau::dormandPrince54()->firstSameAsLast());
0629   BOOST_CHECK(!ButcherTableau::classicalRk4()->firstSameAsLast());
0630   BOOST_CHECK(!ButcherTableau::verner98()->firstSameAsLast());
0631
0632   const FreeVector start =
0633       makeStart(Vector3::Zero(), Vector3(1., 0.3, 0.2), 1. / 1_GeV);
0634   for (const auto& tableau :
0635        {ButcherTableau::classicalRk4(), ButcherTableau::dormandPrince54()}) {
0636     auto field = std::make_shared<CountingField>(Vector3(0., 0., 2_T));
0637     const GenericRungeKuttaStepper stepper(
0638         GenericRungeKuttaStepper::Config{field, tableau});
0639     propagateFree(stepper, fixedStepOptions(10_cm), start, 1_m, false);
0640
0641     // The start, then per step the stages after the first and the end
0642     const std::size_t s = tableau->stages();
0643     const std::size_t expected =
0644         1 + 10 * (tableau->firstSameAsLast() ? s - 1 : s);
0645     BOOST_TEST_CONTEXT(tableau->name()) {
0646       BOOST_CHECK_EQUAL(field->lookups(), expected);
0647     }
0648   }
0649 }
0650
0651 /// The jacobian in the error estimate makes the adaptive steps smaller.
0652 BOOST_DATA_TEST_CASE(generic_runge_kutta_stepper_jacobian_in_error_estimate,
0653                      boost::unit_test::data::make(std::vector<ErrorEstimation>{
0654                          ErrorEstimation::Embedded,
0655                          ErrorEstimation::StepDoubling}),
0656                      errorEstimation) {
0657   auto field = std::make_shared<SmoothField>(2_T, 1_m);
0658   const GenericRungeKuttaStepper stepper(field);
0659   GenericRungeKuttaStepper::Options options(tgContext, mfContext);
0660   options.errorEstimation = errorEstimation;
0661   options.stepTolerance = 1e-7;
0662   options.initialStepSize = 10_m;
0663   const FreeVector start =
0664       makeStart(Vector3(100., -50., 20.), Vector3(1., 0.3, 0.2), 1. / 1_GeV);
0665
0666   const auto withoutJacobian =
0667       propagateFree(stepper, options, start, 1_m, true);
0668   options.jacobianInErrorEstimate = true;
0669   const auto withJacobian = propagateFree(stepper, options, start, 1_m, true);
0670
0671   BOOST_TEST_CONTEXT("steps " << withoutJacobian.nSteps << " "
0672                               << withJacobian.nSteps) {
0673     BOOST_CHECK_GT(withJacobian.nSteps, withoutJacobian.nSteps);
0674   }
0675   CHECK_CLOSE_ABS(withJacobian.pars, withoutJacobian.pars, 1e-6);
0676 }
0677
0678 /// The time of a neutral particle does not depend on q/p, because its charge
0679 /// and so its momentum hypothesis are fixed.
0680 BOOST_AUTO_TEST_CASE(generic_runge_kutta_stepper_neutral_time) {
0681   auto field = std::make_shared<ConstantBField>(Vector3::Zero());
0682   const GenericRungeKuttaStepper stepper(field);
0683   const auto options = fixedStepOptions(10_cm);
0684   const FreeVector start =
0685       makeStart(Vector3(1., 2., 3.), Vector3(1., 0.3, 0.2), 1. / 0.5_GeV);
0686
0687   const auto neutral = propagateFree(stepper, options, start, 1_m, true,
0688                                      ParticleHypothesis::pion0());
0689   BOOST_CHECK_EQUAL(neutral.jacTransport(eFreeTime, eFreeQOverP), 0.);
0690
0691   const auto charged = propagateFree(stepper, options, start, 1_m, true,
0692                                      ParticleHypothesis::pion());
0693   BOOST_CHECK_NE(charged.jacTransport(eFreeTime, eFreeQOverP), 0.);
0694 }
0695
0696 /// A failed field lookup returns the error and leaves the state unchanged,
0697 /// both for a stage and for the end of the step.
0698 ///
0699 /// The explicit Euler method has a single stage at the start, so only the
0700 /// lookup at the end of the step can fail.
0701 BOOST_DATA_TEST_CASE(generic_runge_kutta_stepper_field_failure,
0702                      boost::unit_test::data::make(
0703                          std::vector<std::shared_ptr<const ButcherTableau>>{
0704                              ButcherTableau::classicalRk4(),
0705                              std::make_shared<const ButcherTableau>(
0706                                  "Euler", 1, 0, std::vector<double>{0.},
0707                                  std::vector<std::vector<double>>{{}},
0708                                  std::vector<double>{1.},
0709                                  std::vector<double>{})}),
0710                      tableau) {
0711   auto field = std::make_shared<BoundedField>(Vector3(0., 0., 2_T), 50_cm);
0712   const GenericRungeKuttaStepper stepper(
0713       GenericRungeKuttaStepper::Config{field, tableau});
0714   const auto options = fixedStepOptions(1_m);
0715   const FreeVector start =
0716       makeStart(Vector3::Zero(), Vector3(1., 0., 0.), 1. / 1_GeV);
0717
0718   auto state = stepper.makeState(options);
0719   stepper.initialize(
0720       state, BoundTrackParameters::createCurvilinear(
0721                  makeVector4(start.segment<3>(eFreePos0), start[eFreeTime]),
0722                  start.segment<3>(eFreeDir0), start[eFreeQOverP],
0723                  BoundMatrix::Identity(), ParticleHypothesis::pion()));
0724   const FreeVector parsBefore = state.pars;
0725   const FreeMatrix jacBefore = state.jacTransport;
0726
0727   BOOST_TEST_CONTEXT(tableau->name()) {
0728     auto res = stepper.step(state, Direction::Forward(), nullptr);
0729     BOOST_REQUIRE(!res.ok());
0730     BOOST_CHECK(res.error() == MagneticFieldError::OutOfBounds);
0731     BOOST_CHECK_EQUAL(state.pars, parsBefore);
0732     BOOST_CHECK_EQUAL(state.jacTransport, jacBefore);
0733     BOOST_CHECK_EQUAL(state.pathAccumulated, 0.);
0734     BOOST_CHECK_EQUAL(state.nSteps, 0u);
0735   }
0736 }
0737
0738 BOOST_AUTO_TEST_SUITE_END()
0739
0740 }  // namespace ActsTests