celeritas/field/DormandPrinceStepper.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2020-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/field/DormandPrinceStepper.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/math/Algorithms.hh"
0013
0014 #include "Types.hh"
0015
0016 namespace celeritas
0017 {
0018 //---------------------------------------------------------------------------//
0019 /*!
0020  * Integrate the field ODEs using Dormand-Prince RK5(4)7M.
0021  *
0022  * The algorithm, RK5(4)7M and the coefficients have been adapted from
0023  * J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
0024  * Journal Computational and Applied Mathematics, volume 6, no 1 (1980) and
0025  * the coefficients to locate the mid point are taken from L. F. Shampine,
0026  * "Some Practical Runge-Kutta Formulas", Mathematics of * Computation,
0027  * volume 46, number 17, pp 147 (1986).
0028  *
0029  * For a given ordinary differential equation, \f$dy/dx = f(x, y)\f$, the
0030  * fifth order solution, \f$ y_{n+1} \f$, an embedded fourth order solution,
0031  * \f$ y^{*}_{n+1} \f$, and the error estimate as difference between them are
0032  * as follows,
0033  * \f[
0034      y_{n+1}     = y_n + h \sum_{n=1}^{6} b_i  k_i + O(h^6)
0035      y^{*}_{n+1} = y_n + h \sum_{n=1}^{7} b*_i k_i + O(h^5)
0036      y_{error}   = y_{n+1} - y^{*}_{n+1} = \sum_{n=1}^{7} (b^{*}_i - b_i) k_i
0037    \f]
0038  * where \f$h\f$ is the step to advance and k_i is the right hand side of the
0039  * function at \f$x_n + h c_i\f$, and the coefficients (The Butcher table) for
0040  * Dormand-Prince RK5(4)7M are
0041  * \verbatim
0042    c_i | a_ij
0043    0   |
0044    1/5 | 1/5
0045    3/10| 3/40        9/40
0046    4/5 | 44/45      -56/15      32/9
0047    8/9 | 19372/6561 -25360/2187 64448/6561 -212/729
0048    1   | 9017/3168  -355/33     46732/5247  49/176  -5103/18656
0049    1   | 35/384      0          500/1113    125/192 -2187/6784    11/84
0050   ----------------------------------------------------------------------------
0051    b*_i| 35/384      0          500/1113    125/192 -2187/6784    11/84    0
0052    b_i | 5179/57600  0          7571/16695  393/640 -92097/339200 187/2100 1/40
0053   \endverbatim
0054  *
0055  * The result at the mid point is computed
0056  * \f[
0057      y_{n+1/2}   = y_n + (h/2) \Sigma_{n=1}^{7} c^{*}_i k_i
0058    \f]
0059  * with the coefficients \f$c^{*}\f$ taken from L. F. Shampine (1986).
0060  *
0061  * \todo Rename DormandPrinceIntegrator
0062  */
0063 template<class EquationT>
0064 class DormandPrinceStepper
0065 {
0066   public:
0067     //!@{
0068     //! \name Type aliases
0069     using result_type = FieldStepperResult;
0070     //!@}
0071
0072   public:
0073     //! Construct with the equation of motion
0074     explicit CELER_FUNCTION DormandPrinceStepper(EquationT&& eq)
0075         : calc_rhs_(::celeritas::forward<EquationT>(eq))
0076     {
0077     }
0078
0079     // Adaptive step size control
0080     CELER_FUNCTION result_type operator()(real_type step,
0081                                           OdeState const& beg_state) const;
0082
0083   private:
0084     // Functor to calculate the force applied to a particle
0085     EquationT calc_rhs_;
0086 };
0087
0088 //---------------------------------------------------------------------------//
0089 // DEDUCTION GUIDES
0090 //---------------------------------------------------------------------------//
0091 template<class EquationT>
0092 CELER_FUNCTION
0093 DormandPrinceStepper(EquationT&&) -> DormandPrinceStepper<EquationT>;
0094
0095 //---------------------------------------------------------------------------//
0096 // INLINE DEFINITIONS
0097 //---------------------------------------------------------------------------//
0098 /*!
0099  * Numerically integrate using the DormandPrince RK5(4)7M method.
0100  */
0101 template<class E>
0102 CELER_FUNCTION auto DormandPrinceStepper<E>::operator()(
0103     real_type step, OdeState const& beg_state) const -> result_type
0104 {
0105     using celeritas::axpy;
0106     using R = real_type;
0107
0108     // Coefficients for Dormand-Prince RK5(4)7M
0109     constexpr R a11 = 0.2;
0110
0111     constexpr R a21 = 0.075;
0112     constexpr R a22 = 0.225;
0113
0114     constexpr R a31 = 44 / R(45);
0115     constexpr R a32 = -56 / R(15);
0116     constexpr R a33 = 32 / R(9);
0117
0118     constexpr R a41 = 19372 / R(6561);
0119     constexpr R a42 = -25360 / R(2187);
0120     constexpr R a43 = 64448 / R(6561);
0121     constexpr R a44 = -212 / R(729);
0122
0123     constexpr R a51 = 9017 / R(3168);
0124     constexpr R a52 = -355 / R(33);
0125     constexpr R a53 = 46732 / R(5247);
0126     constexpr R a54 = 49 / R(176);
0127     constexpr R a55 = -5103 / R(18656);
0128
0129     constexpr R a61 = 35 / R(384);
0130     constexpr R a63 = 500 / R(1113);
0131     constexpr R a64 = 125 / R(192);
0132     constexpr R a65 = -2187 / R(6784);
0133     constexpr R a66 = 11 / R(84);
0134
0135     constexpr R d71 = a61 - 5179 / R(57600);
0136     constexpr R d73 = a63 - 7571 / R(16695);
0137     constexpr R d74 = a64 - 393 / R(640);
0138     constexpr R d75 = a65 + 92097 / R(339200);
0139     constexpr R d76 = a66 - 187 / R(2100);
0140     constexpr R d77 = -1 / R(40);
0141
0142     // Coefficients for the mid point calculation by Shampine
0143     constexpr R c71 = R(6025192743.) / R(30085553152.);
0144     constexpr R c73 = R(51252292925.) / R(65400821598.);
0145     constexpr R c74 = R(-2691868925.) / R(45128329728.);
0146     constexpr R c75 = R(187940372067.) / R(1594534317056.);
0147     constexpr R c76 = R(-1776094331.) / R(19743644256.);
0148     constexpr R c77 = R(11237099.) / R(235043384.);
0149
0150     result_type result;
0151
0152     // First step
0153     OdeState k1 = calc_rhs_(beg_state);
0154     OdeState state = beg_state;
0155     axpy(a11 * step, k1, &state);
0156
0157     // Second step
0158     OdeState k2 = calc_rhs_(state);
0159     state = beg_state;
0160     axpy(a21 * step, k1, &state);
0161     axpy(a22 * step, k2, &state);
0162
0163     // Third step
0164     OdeState k3 = calc_rhs_(state);
0165     state = beg_state;
0166     axpy(a31 * step, k1, &state);
0167     axpy(a32 * step, k2, &state);
0168     axpy(a33 * step, k3, &state);
0169
0170     // Fourth step
0171     OdeState k4 = calc_rhs_(state);
0172     state = beg_state;
0173     axpy(a41 * step, k1, &state);
0174     axpy(a42 * step, k2, &state);
0175     axpy(a43 * step, k3, &state);
0176     axpy(a44 * step, k4, &state);
0177
0178     // Fifth step
0179     OdeState k5 = calc_rhs_(state);
0180     state = beg_state;
0181     axpy(a51 * step, k1, &state);
0182     axpy(a52 * step, k2, &state);
0183     axpy(a53 * step, k3, &state);
0184     axpy(a54 * step, k4, &state);
0185     axpy(a55 * step, k5, &state);
0186
0187     // Sixth step
0188     OdeState k6 = calc_rhs_(state);
0189     result.end_state = beg_state;
0190     axpy(a61 * step, k1, &result.end_state);
0191     axpy(a63 * step, k3, &result.end_state);
0192     axpy(a64 * step, k4, &result.end_state);
0193     axpy(a65 * step, k5, &result.end_state);
0194     axpy(a66 * step, k6, &result.end_state);
0195
0196     // Seventh step: the final step
0197     OdeState k7 = calc_rhs_(result.end_state);
0198
0199     // The error estimate
0200     result.err_state = {{0, 0, 0}, {0, 0, 0}};
0201     axpy(d71 * step, k1, &result.err_state);
0202     axpy(d73 * step, k3, &result.err_state);
0203     axpy(d74 * step, k4, &result.err_state);
0204     axpy(d75 * step, k5, &result.err_state);
0205     axpy(d76 * step, k6, &result.err_state);
0206     axpy(d77 * step, k7, &result.err_state);
0207
0208     // The mid point
0209     real_type half_step = step / real_type(2);
0210     result.mid_state = beg_state;
0211     axpy(c71 * half_step, k1, &result.mid_state);
0212     axpy(c73 * half_step, k3, &result.mid_state);
0213     axpy(c74 * half_step, k4, &result.mid_state);
0214     axpy(c75 * half_step, k5, &result.mid_state);
0215     axpy(c76 * half_step, k6, &result.mid_state);
0216     axpy(c77 * half_step, k7, &result.mid_state);
0217
0218     return result;
0219 }
0220
0221 //---------------------------------------------------------------------------//
0222 }  // namespace celeritas