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