From 31ebcaf18401e900cd2deeb76842584a1e9216c7 Mon Sep 17 00:00:00 2001
From: Andreas Stefl <stefl.andreas@gmail.com>
Date: Thu, 20 Feb 2025 14:59:42 +0100
Subject: [PATCH] feat: Adaptive step size for the `AtlasStepper` (#4102)

Until now the `AtlasStepper` only reduced the step size by a factor of `1/2` if the accuracy target was not met. This can be slow and does not allow for increasing the step size which is important in inhomogeneous magnetic fields.

This PR adds our usual adaptive step size solution to the `AtlasStepper` which is identical in the `EigenStepper` and `SympyStepper`.

This also allows for apple to apple comparisons between the steppers as they will use the same number of steps to reach a certain path length.

<!-- This is an auto-generated comment: release notes by coderabbit.ai -->

## Summary by CodeRabbit

- **New Features**
  - Enhanced propagation performance with improved error estimation and dynamic step size adjustments, resulting in more robust and accurate operations.

- **Refactor**
  - Streamlined the underlying process by removing redundant calculations, ensuring more efficient and adaptive propagation behavior.

<!-- end of auto-generated comment: release notes by coderabbit.ai -->
---
 Core/include/Acts/Propagator/AtlasStepper.hpp | 51 ++++++++++++++++---
 1 file changed, 45 insertions(+), 6 deletions(-)

diff --git a/Core/include/Acts/Propagator/AtlasStepper.hpp b/Core/include/Acts/Propagator/AtlasStepper.hpp
index 16cafc17639..80eda451c13 100644
--- a/Core/include/Acts/Propagator/AtlasStepper.hpp
+++ b/Core/include/Acts/Propagator/AtlasStepper.hpp
@@ -1189,6 +1189,36 @@ class AtlasStepper {
     bool Helix = false;
     // if (std::abs(S) < m_cfg.helixStep) Helix = true;
 
+    const auto calcStepSizeScaling =
+        [&](const double errorEstimate_) -> double {
+      // For details about these values see ATL-SOFT-PUB-2009-001
+      constexpr double lower = 0.25;
+      constexpr double upper = 4.0;
+      // This is given by the order of the Runge-Kutta method
+      constexpr double exponent = 0.25;
+
+      double x = state.options.stepTolerance / errorEstimate_;
+
+      if constexpr (exponent == 0.25) {
+        // This is 3x faster than std::pow
+        x = std::sqrt(std::sqrt(x));
+      } else {
+        x = std::pow(x, exponent);
+      }
+
+      return std::clamp(x, lower, upper);
+    };
+
+    const auto isErrorTolerable = [&](const double errorEstimate_) {
+      // For details about these values see ATL-SOFT-PUB-2009-001
+      constexpr double marginFactor = 4.0;
+
+      return errorEstimate_ <= marginFactor * state.options.stepTolerance;
+    };
+
+    double EST = 0;
+    double initialH = h;
+
     std::size_t nStepTrials = 0;
     while (h != 0.) {
       nStepTrials++;
@@ -1269,12 +1299,13 @@ class AtlasStepper {
       //
       // This is (h2 * (sd.k1 - sd.k2 - sd.k3 + sd.k4).template lpNorm<1>())
       // in EigenStepper
-      double EST =
-          2. * h *
-          (std::abs((A1 + A6) - (A3 + A4)) + std::abs((B1 + B6) - (B3 + B4)) +
-           std::abs((C1 + C6) - (C3 + C4)));
-      if (std::abs(EST) > std::abs(state.options.stepTolerance)) {
-        h = h * .5;
+      EST = 2. * std::abs(h) *
+            (std::abs((A1 + A6) - (A3 + A4)) + std::abs((B1 + B6) - (B3 + B4)) +
+             std::abs((C1 + C6) - (C3 + C4)));
+      EST = std::max(1e-20, EST);
+      if (!isErrorTolerable(EST)) {
+        const double stepSizeScaling = calcStepSizeScaling(EST);
+        h *= stepSizeScaling;
         // neutralize the sign of h again
         state.stepSize.setAccuracy(h * propDir);
         //        dltm = 0.;
@@ -1414,6 +1445,14 @@ class AtlasStepper {
     ++state.nSteps;
     state.nStepTrials += nStepTrials;
 
+    const double stepSizeScaling = calcStepSizeScaling(EST);
+    const double nextAccuracy = std::abs(h * stepSizeScaling);
+    const double previousAccuracy = std::abs(state.stepSize.accuracy());
+    const double initialStepLength = std::abs(initialH);
+    if (nextAccuracy < initialStepLength || nextAccuracy > previousAccuracy) {
+      state.stepSize.setAccuracy(nextAccuracy);
+    }
+
     return h;
   }