File indexing completed on 2024-04-14 03:49:26

 /*
     SPDX-FileCopyrightText: 2007 Vladimir Kuznetsov <ks.vladimir@gmail.com>
 
     SPDX-License-Identifier: GPL-2.0-or-later
 */
 
 /** \file gslsolver.h
  *  \brief GslGenericSolver, GslSolver and GslAdaptiveSolver classes
  */
 
 #ifndef STEPCORE_GSLSOLVER_H
 #define STEPCORE_GSLSOLVER_H
 
 // HACK: CMake does not passes definitions to moc
 #if defined(STEPCORE_WITH_GSL) || defined(Q_MOC_RUN)
 
 #include "solver.h"
 #include "object.h"
 
 #include <gsl/gsl_odeiv.h>
 
 namespace StepCore
 {
 
 /** \ingroup solvers
  *  \brief Adaptive and non-adaptive solvers from GSL library
  *  
  *  See https://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html#Ordinary-Differential-Equations
  *  and https://en.wikipedia.org/wiki/Numerical_ordinary_differential_equations
  *
  */
 class GslGenericSolver: public Solver
 {
     STEPCORE_OBJECT(GslGenericSolver)
 
 public:
     /** Constructs GslSolver */
     GslGenericSolver(double stepSize, bool adaptive, const gsl_odeiv_step_type* gslStepType)
         : Solver(stepSize), _adaptive(adaptive), _gslStepType(gslStepType) { init(); }
     /** Constructs GslSolver */
     GslGenericSolver(int dimension, Function function, void* params, double stepSize,
                 bool adaptive, const gsl_odeiv_step_type* gslStepType)
         : Solver(dimension, function, params, stepSize),
             _adaptive(adaptive), _gslStepType(gslStepType) { init(); }
     /** Copy constructor */
     GslGenericSolver(const GslGenericSolver& gslSolver)
         : Solver(gslSolver), _adaptive(gslSolver._adaptive),
             _gslStepType(gslSolver._gslStepType) { init(); }
 
     ~GslGenericSolver() { fini(); }
 
     void setDimension(int dimension) { fini(); _dimension = dimension; init(); }
     void setToleranceAbs(double toleranceAbs) { fini(); _toleranceAbs = toleranceAbs; init(); }
     void setToleranceRel(double toleranceRel) { fini(); _toleranceRel = toleranceRel; init(); }
 
     int doCalcFn(double* t, const VectorXd* y, const VectorXd* yvar,
                             VectorXd* f = nullptr, VectorXd* fvar = nullptr);
     int doEvolve(double* t, double t1, VectorXd* y, VectorXd* yvar);
 
 protected:
     static int gslFunction(double t, const double* y, double* f, void* params);
     void init();
     void fini();
 
     bool   _adaptive;
 
     //gsl_odeiv_control*  _gslControl;
     //gsl_odeiv_evolve*   _gslEvolve;
     VectorXd _yerr;
     VectorXd _ytemp;
     VectorXd _ydiff;
     VectorXd _dydt_in;
     VectorXd _dydt_out;
 
     const gsl_odeiv_step_type* _gslStepType;
     gsl_odeiv_system   _gslSystem;
     gsl_odeiv_step*    _gslStep;
     gsl_odeiv_control* _gslControl;
     gsl_odeiv_evolve*  _gslEvolve;
 };
 
 /** \ingroup solvers
  *  \brief Non-adaptive solvers from GSL library
  */
 class GslSolver: public GslGenericSolver
 {
     STEPCORE_OBJECT(GslSolver)
 public:
     GslSolver(double stepSize, const gsl_odeiv_step_type* gslStepType):
                             GslGenericSolver(stepSize, false, gslStepType) {}
     GslSolver(int dimension, Function function, void* params, double stepSize,
                             const gsl_odeiv_step_type* gslStepType)
                     : GslGenericSolver(dimension, function, params, stepSize, false, gslStepType) {}
     GslSolver(const GslSolver& gslSolver): GslGenericSolver(gslSolver) {}
 };
 
 /** \ingroup solvers
  *  \brief Adaptive solvers from GSL library
  */
 class GslAdaptiveSolver: public GslGenericSolver
 {
     STEPCORE_OBJECT(GslAdaptiveSolver)
 public:
     explicit GslAdaptiveSolver(const gsl_odeiv_step_type* gslStepType):
                             GslGenericSolver(1, true, gslStepType) {}
     GslAdaptiveSolver(int dimension, Function function, void* params,
                             const gsl_odeiv_step_type* gslStepType)
                     : GslGenericSolver(dimension, function, params, 1, true, gslStepType) {}
     GslAdaptiveSolver(const GslAdaptiveSolver& gslSolver): GslGenericSolver(gslSolver) {}
 };
 
 #define STEPCORE_DECLARE_GSLSOLVER(Class, type) \
 class Gsl##Class##Solver: public GslSolver { \
     STEPCORE_OBJECT(Gsl##Class##Solver) \
 public: \
     Gsl##Class##Solver(double stepSize = 0.01): GslSolver(stepSize, gsl_odeiv_step_##type) {} \
     Gsl##Class##Solver(int dimension, Function function, void* params, double stepSize) \
                  : GslSolver(dimension, function, params, stepSize, gsl_odeiv_step_##type) {} \
     Gsl##Class##Solver(const Gsl##Class##Solver& gslSolver): GslSolver(gslSolver) {} \
 };
 
 #define STEPCORE_DECLARE_GSLASOLVER(Class, type) \
 class GslAdaptive##Class##Solver: public GslAdaptiveSolver { \
     STEPCORE_OBJECT(GslAdaptive##Class##Solver) \
 public: \
     GslAdaptive##Class##Solver(): GslAdaptiveSolver(gsl_odeiv_step_##type) {} \
     GslAdaptive##Class##Solver(int dimension, Function function, void* params) \
                     : GslAdaptiveSolver(dimension, function, params, gsl_odeiv_step_##type) {} \
     GslAdaptive##Class##Solver(const GslAdaptive##Class##Solver& gslSolver): GslAdaptiveSolver(gslSolver) {} \
 };
 
 /** \ingroup solvers
  *  \class GslRK2Solver
  *  \brief Runge-Kutta second-order solver from GSL library
  */
 STEPCORE_DECLARE_GSLSOLVER(RK2, rk2)
 
 /** \ingroup solvers
  *  \class GslAdaptiveRK2Solver
  *  \brief Adaptive Runge-Kutta second-order solver from GSL library
  */
 STEPCORE_DECLARE_GSLASOLVER(RK2, rk2)
 
 /** \ingroup solvers
  *  \class GslRK4Solver
  *  \brief Runge-Kutta classical fourth-order solver from GSL library
  */
 STEPCORE_DECLARE_GSLSOLVER(RK4, rk4)
 
 /** \ingroup solvers
  *  \class GslAdaptiveRK4Solver
  *  \brief Adaptive Runge-Kutta classical fourth-order solver from GSL library 
  */
 STEPCORE_DECLARE_GSLASOLVER(RK4, rk4)
 
 /** \ingroup solvers
  *  \class GslRKF45Solver
  *  \brief Runge-Kutta-Fehlberg (4,5) solver from GSL library
  */
 STEPCORE_DECLARE_GSLSOLVER(RKF45, rkf45)
 
 /** \ingroup solvers
  *  \class AdaptiveGslRKF45Solver
  *  \brief Adaptive Runge-Kutta-Fehlberg (4,5) solver from GSL library
  */
 STEPCORE_DECLARE_GSLASOLVER(RKF45, rkf45)
 
 /** \ingroup solvers
  *  \class GslRKCKSolver
  *  \brief Runge-Kutta Cash-Karp (4,5) solver from GSL library
  */
 STEPCORE_DECLARE_GSLSOLVER(RKCK, rkck)
 
 /** \ingroup solvers
  *  \class AdaptiveGslRKCKSolver
  *  \brief Adaptive Runge-Kutta Cash-Karp (4,5) solver from GSL library
  */
 STEPCORE_DECLARE_GSLASOLVER(RKCK, rkck)
 
 /** \ingroup solvers
  *  \class GslRK8PDSolver
  *  \brief Runge-Kutta Prince-Dormand (8,9) solver from GSL library
  */
 STEPCORE_DECLARE_GSLSOLVER(RK8PD, rk8pd)
 
 /** \ingroup solvers
  *  \class GslAdaptiveRK8PDSolver
  *  \brief Adaptive Runge-Kutta Prince-Dormand (8,9) solver from GSL library
  */
 STEPCORE_DECLARE_GSLASOLVER(RK8PD, rk8pd)
 
 /** \ingroup solvers
  *  \class GslRK2IMPSolver
  *  \brief Runge-Kutta implicit second-order solver from GSL library
  */
 STEPCORE_DECLARE_GSLSOLVER(RK2IMP, rk2imp)
 
 /** \ingroup solvers
  *  \class GslAdaptiveRK2IMPSolver
  *  \brief Adaptive Runge-Kutta Prince-Dormand (8,9) solver from GSL library
  */
 STEPCORE_DECLARE_GSLASOLVER(RK2IMP, rk2imp)
 
 /** \ingroup solvers
  *  \class GslRK4IMPSolver
  *  \brief Runge-Kutta implicit fourth-order solver from GSL library
  */
 STEPCORE_DECLARE_GSLSOLVER(RK4IMP, rk4imp)
 
 /** \ingroup solvers
  *  \class GslAdaptiveRK4IMPSolver
  *  \brief Runge-Kutta implicit fourth-order solver from GSL library
  */
 STEPCORE_DECLARE_GSLASOLVER(RK4IMP, rk4imp)
 
 } // namespace StepCore
 
 #endif // defined(STEPCORE_WITH_GSL) || defined(Q_MOC_RUN)
 
 #endif // STEPCORE_GSLSOLVER_H