celeritas/field/FieldSubstepper.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file celeritas/field/FieldSubstepper.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/math/NumericLimits.hh"
0015 #include "corecel/math/SoftEqual.hh"
0016
0017 #include "FieldDriverOptions.hh"
0018 #include "Types.hh"
0019
0020 #include "detail/FieldUtils.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Advance the field state by a single substep based on user tolerances.
0027  *
0028  * The substep length is based on the radius of curvature for the step,
0029  * ensuring that the "miss distance" (sagitta, the distance between the
0030  * straight-line arc and the furthest point) is less than the \c delta_chord
0031  * option. This target length is reduced into sub-substeps if necessary to meet
0032  * a targeted relative error `epsilon_rel_max` based on the position and
0033  * momentum update.
0034  *
0035  * This iteratively reduces the given step length until the sagitta is no more
0036  * than \c delta_chord . The sagitta is calculated as the projection of the
0037  * mid-step point onto the line between the start and end-step points.
0038  *
0039  * Each iteration reduces the step length by a factor of no more than \c
0040  * min_chord_shrink , but is based on an approximate "exact" correction factor
0041  * if the chord length is very small and the curve is circular.
0042  * The sagitta \em h is related to the chord length \em s and radius of
0043  * curvature \em r with the trig expression: \f[
0044    r - h = r \cos \frac{s}{2r}
0045   \f]
0046  * For small chord lengths or a large radius, we expand
0047  * \f$ \cos \theta \sim 1 - \frac{\theta^2}{2} \f$, giving a radius of
0048  * curvature \f[
0049    r = \frac{s^2}{8h} \; .
0050    \f]
0051  * Given a trial step (chord length) \em s and resulting sagitta of \em h,
0052  * the exact step needed to give a chord length of \f$ \epsilon = {} \f$ \c
0053  * delta_chord is \f[
0054    s' = s \sqrt{\frac{\epsilon}{h}} \,.
0055  * \f]
0056  *
0057  * \note This class is based on G4ChordFinder and G4MagIntegratorDriver.
0058  */
0059 template<class IntegratorT>
0060 class FieldSubstepper
0061 {
0062   public:
0063     // Construct with options data and the integrator
0064     inline CELER_FUNCTION FieldSubstepper(FieldDriverOptions const& options,
0065                                           IntegratorT&& integrate);
0066
0067     // For a given trial step, advance by a sub_step within a tolerance error
0068     inline CELER_FUNCTION Substep operator()(real_type step,
0069                                              OdeState const& state);
0070
0071     //// TESTABLE HELPER FUNCTIONS ////
0072
0073     // An adaptive step size control from G4MagIntegratorDriver
0074     // Move this to private after all tests with non-uniform field are done
0075     inline CELER_FUNCTION Substep accurate_advance(real_type step,
0076                                                    OdeState const& state,
0077                                                    real_type hinitial) const;
0078
0079     //// ACCESSORS (TODO: refactor) ////
0080
0081     CELER_FUNCTION short int max_substeps() const
0082     {
0083         return options_.max_substeps;
0084     }
0085
0086     CELER_FUNCTION real_type minimum_step() const
0087     {
0088         return options_.minimum_step;
0089     }
0090
0091     // TODO: this should be field propagator data
0092     CELER_FUNCTION real_type delta_intersection() const
0093     {
0094         return options_.delta_intersection;
0095     }
0096
0097   private:
0098     //// DATA ////
0099
0100     // Driver configuration
0101     FieldDriverOptions const& options_;
0102
0103     // Integrator for this field driver
0104     IntegratorT integrate_;
0105
0106     // Maximum chord length based on a previous estimate
0107     real_type max_chord_{numeric_limits<real_type>::infinity()};
0108
0109     //// TYPES ////
0110
0111     //! A helper output for private member functions
0112     struct ChordSearch
0113     {
0114         Substep end;  //!< Step taken and post-step state
0115         real_type err_sq;  //!< Square of the truncation error
0116     };
0117
0118     struct Integration
0119     {
0120         Substep end;  //!< Step taken and post-step state
0121         real_type proposed_length;  //!< Proposed next step size
0122     };
0123
0124     //// HEPER FUNCTIONS ////
0125
0126     // Find the next acceptable chord whose sagitta is less than delta_chord
0127     inline CELER_FUNCTION ChordSearch
0128     find_next_chord(real_type step, OdeState const& state) const;
0129
0130     // Advance for a given step and evaluate the next predicted step.
0131     inline CELER_FUNCTION Integration
0132     integrate_step(real_type step, OdeState const& state) const;
0133
0134     // Advance within the truncated error and estimate a good next step size
0135     inline CELER_FUNCTION Integration one_good_step(real_type step,
0136                                                     OdeState const& state) const;
0137
0138     // Propose a next step size from a given step size and associated error
0139     inline CELER_FUNCTION real_type new_step_scale(real_type error_sq) const;
0140
0141     //// COMMON PROPERTIES ////
0142
0143     static CELER_CONSTEXPR_FUNCTION real_type half() { return 0.5; }
0144 };
0145
0146 //---------------------------------------------------------------------------//
0147 // DEDUCTION GUIDES
0148 //---------------------------------------------------------------------------//
0149 template<class IntegratorT>
0150 CELER_FUNCTION FieldSubstepper(FieldDriverOptions const&,
0151                                IntegratorT&&) -> FieldSubstepper<IntegratorT>;
0152
0153 //---------------------------------------------------------------------------//
0154 // INLINE DEFINITIONS
0155 //---------------------------------------------------------------------------//
0156 /*!
0157  * Construct with options and the step advancement functor.
0158  */
0159 template<class IntegratorT>
0160 CELER_FUNCTION
0161 FieldSubstepper<IntegratorT>::FieldSubstepper(FieldDriverOptions const& options,
0162                                               IntegratorT&& integrate)
0163     : options_(options)
0164     , integrate_(::celeritas::forward<IntegratorT>(integrate))
0165 {
0166     CELER_EXPECT(options_);
0167 }
0168
0169 //---------------------------------------------------------------------------//
0170 /*!
0171  * Adaptive step control based on G4ChordFinder and G4MagIntegratorDriver.
0172  *
0173  * \param step maximum step length
0174  * \param state starting state
0175  * \return substep and updated state
0176  *
0177  * For a given trial step, advance by a sub-step within a required tolerance
0178  * and update the current state (position and momentum).  For an efficient
0179  * adaptive integration, the proposed chord of which the sagitta (the
0180  * largest distance from the curved trajectory to the chord) is smaller than
0181  * a reference distance (dist_chord) will be accepted if its stepping error is
0182  * within a reference accuracy. Otherwise, the more accurate step integration
0183  * (advance_accurate) will be performed.
0184  */
0185 template<class IntegratorT>
0186 CELER_FUNCTION Substep
0187 FieldSubstepper<IntegratorT>::operator()(real_type step, OdeState const& state)
0188 {
0189     if (step <= options_.minimum_step)
0190     {
0191         // If the input is a very tiny step, do a "quick advance".
0192         Substep result;
0193         result.state = integrate_(step, state).end_state;
0194         result.length = step;
0195         return result;
0196     }
0197
0198     // Calculate the next chord length (and get an end state "for free") based
0199     // on delta_chord, reusing previous estimates
0200     ChordSearch next
0201         = this->find_next_chord(celeritas::min(step, max_chord_), state);
0202     CELER_ASSERT(next.end.length <= step);
0203     if (next.end.length < step)
0204     {
0205         // Chord length was reduced due to constraints: save the estimate for
0206         // the next potential field advance inside the propagation loop
0207         max_chord_ = next.end.length * (1 / options_.min_chord_shrink);
0208     }
0209
0210     if (next.err_sq > 1)
0211     {
0212         // Discard the original end state and advance more accurately with the
0213         // newly proposed (reduced) step
0214         real_type next_step = step * this->new_step_scale(next.err_sq);
0215         next.end = this->accurate_advance(next.end.length, state, next_step);
0216     }
0217
0218     CELER_ENSURE(next.end.length > 0 && next.end.length <= step);
0219     return next.end;
0220 }
0221
0222 //---------------------------------------------------------------------------//
0223 /*!
0224  * Find the maximum step length that satisfies a maximum "miss distance".
0225  */
0226 template<class IntegratorT>
0227 CELER_FUNCTION auto FieldSubstepper<IntegratorT>::find_next_chord(
0228     real_type step, OdeState const& state) const -> ChordSearch
0229 {
0230     bool succeeded = false;
0231     auto remaining_steps = options_.max_nsteps;
0232     FieldIntegration integrated;
0233
0234     do
0235     {
0236         // Try with the proposed step
0237         integrated = integrate_(step, state);
0238
0239         // Check whether the distance to the chord is smaller than the
0240         // reference
0241         real_type dchord = detail::distance_chord(
0242             state.pos, integrated.mid_state.pos, integrated.end_state.pos);
0243
0244         if (dchord > options_.delta_chord + options_.dchord_tol)
0245         {
0246             // Estimate a new trial chord with a relative scale
0247             real_type scale_step = max(std::sqrt(options_.delta_chord / dchord),
0248                                        options_.min_chord_shrink);
0249             step *= scale_step;
0250         }
0251         else
0252         {
0253             succeeded = true;
0254         }
0255     } while (!succeeded && --remaining_steps > 0);
0256
0257     // Update step, position and momentum
0258     ChordSearch result;
0259     result.end.length = step;
0260     result.end.state = integrated.end_state;
0261     result.err_sq = detail::rel_err_sq(integrated.err_state, step, state.mom)
0262                     / ipow<2>(options_.epsilon_rel_max);
0263
0264     return result;
0265 }
0266
0267 //---------------------------------------------------------------------------//
0268 /*!
0269  * Accurate advance for an adaptive step control.
0270  *
0271  * Perform an adaptive step integration for a proposed step or a series of
0272  * sub-steps within a required tolerance until the the accumulated curved path
0273  * is equal to the input step length.
0274  *
0275  * TODO: maybe this should be moved out of the substepper into the propagation
0276  * loop?
0277  */
0278 template<class IntegratorT>
0279 CELER_FUNCTION Substep FieldSubstepper<IntegratorT>::accurate_advance(
0280     real_type step, OdeState const& state, real_type hinitial) const
0281 {
0282     CELER_ASSERT(step > 0);
0283
0284     // Set an initial proposed step and evaluate the minimum threshold
0285     real_type end_curve_length = step;
0286
0287     // Use a pre-defined initial step size if it is smaller than the input
0288     // step length and larger than the per-million fraction of the step length.
0289     // Otherwise, use the input step length for the first trial.
0290     // TODO: review whether this approach is an efficient bootstrapping.
0291     real_type h
0292         = ((hinitial > options_.initial_step_tol * step) && (hinitial < step))
0293               ? hinitial
0294               : step;
0295     real_type h_threshold = options_.epsilon_step * step;
0296
0297     // Output with the next good step
0298     Integration result;
0299     result.end.state = state;
0300
0301     // Perform integration
0302     bool succeeded = false;
0303     real_type curve_length = 0;
0304     auto remaining_steps = options_.max_nsteps;
0305
0306     do
0307     {
0308         CELER_ASSERT(h > 0);
0309         result = this->integrate_step(h, result.end.state);
0310
0311         curve_length += result.end.length;
0312
0313         if (h < h_threshold || curve_length >= end_curve_length)
0314         {
0315             succeeded = true;
0316         }
0317         else
0318         {
0319             h = celeritas::min(
0320                 celeritas::max(result.proposed_length, options_.minimum_step),
0321                 end_curve_length - curve_length);
0322         }
0323     } while (!succeeded && --remaining_steps > 0);
0324
0325     // Curve length may be slightly longer than step due to roundoff in
0326     // accumulation
0327     CELER_ENSURE(curve_length > 0
0328                  && (curve_length <= step || soft_equal(curve_length, step)));
0329     result.end.length = min(curve_length, step);
0330     return result.end;
0331 }
0332
0333 //---------------------------------------------------------------------------//
0334 /*!
0335  * Advance for a given step and evaluate the next predicted step.
0336  *
0337  * Helper function for accurate_advance.
0338  */
0339 template<class IntegratorT>
0340 CELER_FUNCTION auto FieldSubstepper<IntegratorT>::integrate_step(
0341     real_type step, OdeState const& state) const -> Integration
0342 {
0343     CELER_EXPECT(step > 0);
0344
0345     // Output with a next proposed step
0346     Integration result;
0347
0348     if (step > options_.minimum_step)
0349     {
0350         result = this->one_good_step(step, state);
0351     }
0352     else
0353     {
0354         // Do an integration step for a small step (a.k.a quick advance)
0355         FieldIntegration integrated = integrate_(step, state);
0356
0357         // Update position and momentum
0358         result.end.state = integrated.end_state;
0359         result.end.length = step;
0360
0361         // Compute a proposed new step
0362         real_type err_sq
0363             = detail::rel_err_sq(integrated.err_state, step, state.mom)
0364               / ipow<2>(options_.epsilon_rel_max);
0365         result.proposed_length = step * this->new_step_scale(err_sq);
0366     }
0367
0368     return result;
0369 }
0370
0371 //---------------------------------------------------------------------------//
0372 /*!
0373  * Advance within a relative truncation error and estimate a good step size
0374  * for the next integration.
0375  */
0376 template<class IntegratorT>
0377 CELER_FUNCTION auto FieldSubstepper<IntegratorT>::one_good_step(
0378     real_type step, OdeState const& state) const -> Integration
0379 {
0380     // Perform integration for adaptive step control with the truncation error
0381     bool succeeded = false;
0382     size_type remaining_steps = options_.max_nsteps;
0383     real_type err_sq;
0384     FieldIntegration integrated;
0385
0386     do
0387     {
0388         integrated = integrate_(step, state);
0389
0390         err_sq = detail::rel_err_sq(integrated.err_state, step, state.mom)
0391                  / ipow<2>(options_.epsilon_rel_max);
0392
0393         if (err_sq > 1)
0394         {
0395             // Truncation error too large, reduce stepsize with a low bound
0396             step *= max(this->new_step_scale(err_sq),
0397                         options_.max_stepping_decrease);
0398         }
0399         else
0400         {
0401             // Success or possibly nan!
0402             succeeded = true;
0403         }
0404     } while (!succeeded && --remaining_steps > 0);
0405
0406     // Update state, step taken by this trial and the next predicted step
0407     Integration result;
0408     result.end.state = integrated.end_state;
0409     result.end.length = step;
0410     result.proposed_length
0411         = step
0412           * min(this->new_step_scale(err_sq), options_.max_stepping_increase);
0413
0414     return result;
0415 }
0416
0417 //---------------------------------------------------------------------------//
0418 /*!
0419  * Estimate the new predicted step size based on the error estimate.
0420  */
0421 template<class IntegratorT>
0422 CELER_FUNCTION real_type
0423 FieldSubstepper<IntegratorT>::new_step_scale(real_type err_sq) const
0424 {
0425     CELER_ASSERT(err_sq >= 0);
0426     return options_.safety
0427            * fastpow(err_sq,
0428                      half() * (err_sq > 1 ? options_.pshrink : options_.pgrow));
0429 }
0430
0431 //---------------------------------------------------------------------------//
0432 }  // namespace celeritas