include/Geant4/G4BulirschStoer.hh

0001 // ********************************************************************
0002 // * License and Disclaimer                                           *
0003 // *                                                                  *
0004 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0005 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0006 // * conditions of the Geant4 Software License,  included in the file *
0007 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0008 // * include a list of copyright holders.                             *
0009 // *                                                                  *
0010 // * Neither the authors of this software system, nor their employing *
0011 // * institutes,nor the agencies providing financial support for this *
0012 // * work  make  any representation or  warranty, express or implied, *
0013 // * regarding  this  software system or assume any liability for its *
0014 // * use.  Please see the license in the file  LICENSE  and URL above *
0015 // * for the full disclaimer and the limitation of liability.         *
0016 // *                                                                  *
0017 // * This  code  implementation is the result of  the  scientific and *
0018 // * technical work of the GEANT4 collaboration.                      *
0019 // * By using,  copying,  modifying or  distributing the software (or *
0020 // * any work based  on the software)  you  agree  to acknowledge its *
0021 // * use  in  resulting  scientific  publications,  and indicate your *
0022 // * acceptance of all terms of the Geant4 Software license.          *
0023 // ********************************************************************
0024 //
0025 // G4BulirschStoer
0026 //
0027 // Class description:
0028 //
0029 // The Bulirsch-Stoer is a controlled driver that adjusts both step size
0030 // and order of the method. The algorithm uses the modified midpoint and
0031 // a polynomial extrapolation computes the solution.
0032
0033 // Author: Dmitry Sorokin (CERN, Google Summer of Code 2016), 13.02.2018
0034 // Supervision: John Apostolakis (CERN)
0035 // --------------------------------------------------------------------
0036 #ifndef G4BULIRSCH_STOER_HH
0037 #define G4BULIRSCH_STOER_HH
0038
0039 #include "G4ModifiedMidpoint.hh"
0040
0041 #include "G4FieldTrack.hh"
0042
0043 /**
0044  * @brief G4BulirschStoer is a controlled driver that adjusts both step size
0045  * and order of the method. The algorithm uses the modified midpoint and
0046  * a polynomial extrapolation computes the solution.
0047  */
0048
0049 class G4BulirschStoer
0050 {
0051   public:
0052
0053     enum class step_result { success, fail };
0054
0055     /**
0056      * Constructor for G4BulirschStoer.
0057      *  @param[in] equation Pointer to the provided equation of motion.
0058      *  @param[in] nvar The number of integration variables.
0059      *  @param[in] eps_rel Relative tolerance.
0060      *  @param[in] max_dt Maximum allowed time step.
0061      */
0062     G4BulirschStoer(G4EquationOfMotion* equation, G4int nvar,
0063                     G4double eps_rel, G4double max_dt = DBL_MAX);
0064
0065     /**
0066      * Default Destructor.
0067      */
0068     ~G4BulirschStoer() = default;
0069
0070     /**
0071      * Modifiers.
0072      */
0073     inline void set_max_dt(G4double max_dt);
0074     inline void set_max_relative_error(G4double eps_rel);
0075
0076     /**
0077      * Stepper method.
0078      *  @param[in] in Initial position.
0079      *  @param[in] dxdt dxdt for mid-point calculation.
0080      *  @param[out] t The updated step.
0081      *  @param[out] out Updated position.
0082      *  @param[in,out] dt Step size.
0083      *  @returns success if step is not rejected.
0084      */
0085     step_result try_step(const G4double in[], const G4double dxdt[],
0086                          G4double& t, G4double out[], G4double& dt);
0087
0088     /**
0089      * Resets the internal state of the stepper.
0090      */
0091     void reset();
0092
0093     /**
0094      * Setter and getter for the equation of motion.
0095      */
0096     inline void SetEquationOfMotion(G4EquationOfMotion* equation);
0097     inline G4EquationOfMotion* GetEquationOfMotion() const;
0098
0099     /**
0100      * Returns the number of integration variables.
0101      */
0102     inline G4int GetNumberOfVariables() const;
0103
0104   private:
0105
0106     /**
0107      * Polynomial extrapolation.
0108      */
0109     void extrapolate(std::size_t k, G4double xest[]);
0110
0111     /**
0112      * Calculates the optimal step size for a given error and stage number.
0113      */
0114     G4double calc_h_opt(G4double h, G4double error, std::size_t k) const;
0115
0116     /**
0117      * Calculates the optimal stage number.
0118      */
0119     G4bool set_k_opt(std::size_t k, G4double& dt);
0120
0121     /**
0122      * Utilities.
0123      */
0124     G4bool in_convergence_window(G4int k) const;
0125     G4bool should_reject(G4double error, G4int k) const;
0126
0127   private:
0128
0129     /** Maximum number of stages. */
0130     const static G4int m_k_max = 8;
0131
0132     /** Number of vars to be integrated. */
0133     G4int fnvar;
0134
0135     /** Relative tolerance. */
0136     G4double m_eps_rel;
0137
0138     /** Modified midpoint algorithm. */
0139     G4ModifiedMidpoint m_midpoint;
0140
0141     /** Flags for step. */
0142     G4bool m_last_step_rejected{false};
0143     G4bool m_first{true};
0144
0145     /** Last step size. */
0146     G4double m_dt_last{0.0};
0147
0148     /** Max allowed time step. */
0149     G4double m_max_dt;
0150
0151     /** Crude estimate of optimal order. */
0152     G4int m_current_k_opt;
0153
0154     /** Error estimate. */
0155     G4double m_err[G4FieldTrack::ncompSVEC];
0156
0157     /** Stores the successive interval counts. */
0158     G4int m_interval_sequence[m_k_max+1];
0159
0160     /** Extrapolation coeffs (Neville's algorithm). */
0161     G4double m_coeff[m_k_max+1][m_k_max];
0162
0163     /** Costs for interval count. */
0164     G4int m_cost[m_k_max+1];
0165
0166     /** Sequence of states for extrapolation. */
0167     G4double m_table[m_k_max][G4FieldTrack::ncompSVEC];
0168
0169     /** Optimal step size. */
0170     G4double h_opt[m_k_max+1];
0171
0172     /** Work per unit step. */
0173     G4double work[m_k_max+1];
0174 };
0175
0176 #include "G4BulirschStoer.icc"
0177
0178 #endif