basic_rk_integrator.hpp

Go to the documentation of this file.

// ====================================================================
// This file is part of FlexibleSUSY.
//
// FlexibleSUSY is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// FlexibleSUSY is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FlexibleSUSY.  If not, see
// <http://www.gnu.org/licenses/>.
// ====================================================================
 
#ifndef BASIC_RK_INTEGRATOR_H
#define BASIC_RK_INTEGRATOR_H
 
#include <functional>
 
#include "rk.hpp"
 
namespace flexiblesusy {
 
namespace runge_kutta {
 
template <typename StateType, typename Derivs>
class Basic_rk_stepper {
public:
   double step(StateType&, const StateType&, double&, double,
               double, const StateType&, Derivs, int&) const;
private:
   void runge_kutta_step(const StateType&, const StateType&, double,
                         double, StateType&, StateType&, Derivs) const;
};
 
template <typename StateType, typename Derivs>
void Basic_rk_stepper<StateType, Derivs>::runge_kutta_step(
   const StateType& y, const StateType& dydx, double x,
   double h, StateType& yout, StateType& yerr, Derivs derivs) const
{
   rungeKuttaStep(y, dydx, x, h, yout, yerr, derivs);
}
 
template <typename StateType, typename Derivs>
double Basic_rk_stepper<StateType,Derivs>::step(
   StateType& y, const StateType& dydx, double& x, double htry,
   double eps, const StateType& yscal, Derivs derivs,
   int& max_step_dir) const
{
   return odeStepper(y, dydx, x, htry, eps, yscal, derivs, max_step_dir);
}
 
template <typename StateType,
          typename Derivs
          = std::function<StateType(double, const StateType&)>,
          typename Stepper = Basic_rk_stepper<StateType,Derivs> >
class Basic_rk_integrator {
public:
   void operator()(double start, double end, StateType& ystart,
                   Derivs derivs, double tolerance) const;
 
   void set_max_steps(int s) { max_steps = s; }
 
   int get_max_steps() const { return max_steps; }
 
private:
   int max_steps{400}; 
   Stepper stepper{};  
};
 
template <typename StateType, typename Derivs, typename Stepper>
void Basic_rk_integrator<StateType, Derivs, Stepper>::operator()(
   double start, double end, StateType& ystart, Derivs derivs,
   double tolerance) const
{
   const double guess = (start - end) * 0.1; // first step size
   const double hmin = (start - end) * tolerance * 1.0e-5;
   const auto rkqs = [this] (
      StateType& y, const StateType& dydx, double& x, double htry,
      double eps, const StateType& yscal, Derivs derivs,
      int& max_step_dir) -> double {
      return this->stepper.step(y, dydx, x, htry, eps,
                                yscal, derivs, max_step_dir);
   };
 
   integrateOdes(ystart, start, end, tolerance, guess, hmin,
                 derivs, rkqs, max_steps);
}
 
} // namespace runge_kutta
 
} // namespace flexiblesusy
 
#endif