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, Google Summer of Code 2016
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 class G4BulirschStoer
0044 {
0045   public:
0046
0047     enum class step_result
0048     {
0049       success,
0050       fail
0051     };
0052
0053     G4BulirschStoer( G4EquationOfMotion* equation, G4int nvar,
0054                      G4double eps_rel, G4double max_dt = DBL_MAX);
0055
0056     inline void set_max_dt(G4double max_dt);
0057     inline void set_max_relative_error(G4double eps_rel);
0058
0059     // Stepper method
0060     //
0061     step_result try_step(const G4double in[], const G4double dxdt[],
0062                          G4double& t, G4double out[], G4double& dt);
0063
0064     // Reset the internal state of the stepper
0065     //
0066     void reset();
0067
0068     inline void SetEquationOfMotion(G4EquationOfMotion* equation);
0069     inline G4EquationOfMotion* GetEquationOfMotion();
0070
0071     inline G4int GetNumberOfVariables() const;
0072
0073   private:
0074
0075     const static G4int m_k_max = 8;
0076
0077     void extrapolate(std::size_t k, G4double xest[]);
0078     G4double calc_h_opt(G4double h, G4double error, std::size_t k) const;
0079
0080     G4bool set_k_opt(std::size_t k, G4double& dt);
0081     G4bool in_convergence_window(G4int k) const;
0082     G4bool should_reject(G4double error, G4int k) const;
0083
0084     // Number of vars to be integrated
0085     G4int fnvar;
0086
0087     // Relative tolerance
0088     G4double m_eps_rel;
0089
0090     // Modified midpoint algorithm
0091     G4ModifiedMidpoint m_midpoint;
0092
0093     G4bool m_last_step_rejected{false};
0094     G4bool m_first{true};
0095
0096     G4double m_dt_last{0.0};
0097     // G4double m_t_last;
0098
0099     // Max allowed time step
0100     G4double m_max_dt;
0101
0102     G4int m_current_k_opt;
0103
0104     // G4double m_xnew[G4FieldTrack::ncompSVEC];
0105     G4double m_err[G4FieldTrack::ncompSVEC];
0106     // G4double m_dxdt[G4FieldTrack::ncompSVEC];
0107
0108     // Stores the successive interval counts
0109     G4int m_interval_sequence[m_k_max+1];
0110
0111     // Extrapolation coeffs (Neville’s algorithm)
0112     G4double m_coeff[m_k_max+1][m_k_max];
0113
0114     // Costs for interval count
0115     G4int m_cost[m_k_max+1];
0116
0117     // Sequence of states for extrapolation
0118     G4double m_table[m_k_max][G4FieldTrack::ncompSVEC];
0119
0120     // Optimal step size
0121     G4double h_opt[m_k_max+1];
0122
0123     // Work per unit step
0124     G4double work[m_k_max+1];
0125 };
0126
0127 #include "G4BulirschStoer.icc"
0128
0129 #endif