File indexing completed on 2024-03-03 03:59:18

0001 /*
0002     SPDX-FileCopyrightText: 2007 Vladimir Kuznetsov <ks.vladimir@gmail.com>
0003 
0004     SPDX-License-Identifier: GPL-2.0-or-later
0005 */
0006 
0007 /** \file gslsolver.h
0008  *  \brief GslGenericSolver, GslSolver and GslAdaptiveSolver classes
0009  */
0010 
0011 #ifndef STEPCORE_GSLSOLVER_H
0012 #define STEPCORE_GSLSOLVER_H
0013 
0014 // HACK: CMake does not passes definitions to moc
0015 #if defined(STEPCORE_WITH_GSL) || defined(Q_MOC_RUN)
0016 
0017 #include "solver.h"
0018 #include "object.h"
0019 
0020 #include <gsl/gsl_odeiv.h>
0021 
0022 namespace StepCore
0023 {
0024 
0025 /** \ingroup solvers
0026  *  \brief Adaptive and non-adaptive solvers from GSL library
0027  *  
0028  *  See https://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html#Ordinary-Differential-Equations
0029  *  and https://en.wikipedia.org/wiki/Numerical_ordinary_differential_equations
0030  *
0031  */
0032 class GslGenericSolver: public Solver
0033 {
0034     STEPCORE_OBJECT(GslGenericSolver)
0035 
0036 public:
0037     /** Constructs GslSolver */
0038     GslGenericSolver(double stepSize, bool adaptive, const gsl_odeiv_step_type* gslStepType)
0039         : Solver(stepSize), _adaptive(adaptive), _gslStepType(gslStepType) { init(); }
0040     /** Constructs GslSolver */
0041     GslGenericSolver(int dimension, Function function, void* params, double stepSize,
0042                 bool adaptive, const gsl_odeiv_step_type* gslStepType)
0043         : Solver(dimension, function, params, stepSize),
0044             _adaptive(adaptive), _gslStepType(gslStepType) { init(); }
0045     /** Copy constructor */
0046     GslGenericSolver(const GslGenericSolver& gslSolver)
0047         : Solver(gslSolver), _adaptive(gslSolver._adaptive),
0048             _gslStepType(gslSolver._gslStepType) { init(); }
0049 
0050     ~GslGenericSolver() { fini(); }
0051 
0052     void setDimension(int dimension) { fini(); _dimension = dimension; init(); }
0053     void setToleranceAbs(double toleranceAbs) { fini(); _toleranceAbs = toleranceAbs; init(); }
0054     void setToleranceRel(double toleranceRel) { fini(); _toleranceRel = toleranceRel; init(); }
0055 
0056     int doCalcFn(double* t, const VectorXd* y, const VectorXd* yvar,
0057                             VectorXd* f = nullptr, VectorXd* fvar = nullptr);
0058     int doEvolve(double* t, double t1, VectorXd* y, VectorXd* yvar);
0059 
0060 protected:
0061     static int gslFunction(double t, const double* y, double* f, void* params);
0062     void init();
0063     void fini();
0064 
0065     bool   _adaptive;
0066 
0067     //gsl_odeiv_control*  _gslControl;
0068     //gsl_odeiv_evolve*   _gslEvolve;
0069     VectorXd _yerr;
0070     VectorXd _ytemp;
0071     VectorXd _ydiff;
0072     VectorXd _dydt_in;
0073     VectorXd _dydt_out;
0074 
0075     const gsl_odeiv_step_type* _gslStepType;
0076     gsl_odeiv_system   _gslSystem;
0077     gsl_odeiv_step*    _gslStep;
0078     gsl_odeiv_control* _gslControl;
0079     gsl_odeiv_evolve*  _gslEvolve;
0080 };
0081 
0082 /** \ingroup solvers
0083  *  \brief Non-adaptive solvers from GSL library
0084  */
0085 class GslSolver: public GslGenericSolver
0086 {
0087     STEPCORE_OBJECT(GslSolver)
0088 public:
0089     GslSolver(double stepSize, const gsl_odeiv_step_type* gslStepType):
0090                             GslGenericSolver(stepSize, false, gslStepType) {}
0091     GslSolver(int dimension, Function function, void* params, double stepSize,
0092                             const gsl_odeiv_step_type* gslStepType)
0093                     : GslGenericSolver(dimension, function, params, stepSize, false, gslStepType) {}
0094     GslSolver(const GslSolver& gslSolver): GslGenericSolver(gslSolver) {}
0095 };
0096 
0097 /** \ingroup solvers
0098  *  \brief Adaptive solvers from GSL library
0099  */
0100 class GslAdaptiveSolver: public GslGenericSolver
0101 {
0102     STEPCORE_OBJECT(GslAdaptiveSolver)
0103 public:
0104     explicit GslAdaptiveSolver(const gsl_odeiv_step_type* gslStepType):
0105                             GslGenericSolver(1, true, gslStepType) {}
0106     GslAdaptiveSolver(int dimension, Function function, void* params,
0107                             const gsl_odeiv_step_type* gslStepType)
0108                     : GslGenericSolver(dimension, function, params, 1, true, gslStepType) {}
0109     GslAdaptiveSolver(const GslAdaptiveSolver& gslSolver): GslGenericSolver(gslSolver) {}
0110 };
0111 
0112 #define STEPCORE_DECLARE_GSLSOLVER(Class, type) \
0113 class Gsl##Class##Solver: public GslSolver { \
0114     STEPCORE_OBJECT(Gsl##Class##Solver) \
0115 public: \
0116     Gsl##Class##Solver(double stepSize = 0.01): GslSolver(stepSize, gsl_odeiv_step_##type) {} \
0117     Gsl##Class##Solver(int dimension, Function function, void* params, double stepSize) \
0118                  : GslSolver(dimension, function, params, stepSize, gsl_odeiv_step_##type) {} \
0119     Gsl##Class##Solver(const Gsl##Class##Solver& gslSolver): GslSolver(gslSolver) {} \
0120 };
0121 
0122 #define STEPCORE_DECLARE_GSLASOLVER(Class, type) \
0123 class GslAdaptive##Class##Solver: public GslAdaptiveSolver { \
0124     STEPCORE_OBJECT(GslAdaptive##Class##Solver) \
0125 public: \
0126     GslAdaptive##Class##Solver(): GslAdaptiveSolver(gsl_odeiv_step_##type) {} \
0127     GslAdaptive##Class##Solver(int dimension, Function function, void* params) \
0128                     : GslAdaptiveSolver(dimension, function, params, gsl_odeiv_step_##type) {} \
0129     GslAdaptive##Class##Solver(const GslAdaptive##Class##Solver& gslSolver): GslAdaptiveSolver(gslSolver) {} \
0130 };
0131 
0132 /** \ingroup solvers
0133  *  \class GslRK2Solver
0134  *  \brief Runge-Kutta second-order solver from GSL library
0135  */
0136 STEPCORE_DECLARE_GSLSOLVER(RK2, rk2)
0137 
0138 /** \ingroup solvers
0139  *  \class GslAdaptiveRK2Solver
0140  *  \brief Adaptive Runge-Kutta second-order solver from GSL library
0141  */
0142 STEPCORE_DECLARE_GSLASOLVER(RK2, rk2)
0143 
0144 /** \ingroup solvers
0145  *  \class GslRK4Solver
0146  *  \brief Runge-Kutta classical fourth-order solver from GSL library
0147  */
0148 STEPCORE_DECLARE_GSLSOLVER(RK4, rk4)
0149 
0150 /** \ingroup solvers
0151  *  \class GslAdaptiveRK4Solver
0152  *  \brief Adaptive Runge-Kutta classical fourth-order solver from GSL library 
0153  */
0154 STEPCORE_DECLARE_GSLASOLVER(RK4, rk4)
0155 
0156 /** \ingroup solvers
0157  *  \class GslRKF45Solver
0158  *  \brief Runge-Kutta-Fehlberg (4,5) solver from GSL library
0159  */
0160 STEPCORE_DECLARE_GSLSOLVER(RKF45, rkf45)
0161 
0162 /** \ingroup solvers
0163  *  \class AdaptiveGslRKF45Solver
0164  *  \brief Adaptive Runge-Kutta-Fehlberg (4,5) solver from GSL library
0165  */
0166 STEPCORE_DECLARE_GSLASOLVER(RKF45, rkf45)
0167 
0168 /** \ingroup solvers
0169  *  \class GslRKCKSolver
0170  *  \brief Runge-Kutta Cash-Karp (4,5) solver from GSL library
0171  */
0172 STEPCORE_DECLARE_GSLSOLVER(RKCK, rkck)
0173 
0174 /** \ingroup solvers
0175  *  \class AdaptiveGslRKCKSolver
0176  *  \brief Adaptive Runge-Kutta Cash-Karp (4,5) solver from GSL library
0177  */
0178 STEPCORE_DECLARE_GSLASOLVER(RKCK, rkck)
0179 
0180 /** \ingroup solvers
0181  *  \class GslRK8PDSolver
0182  *  \brief Runge-Kutta Prince-Dormand (8,9) solver from GSL library
0183  */
0184 STEPCORE_DECLARE_GSLSOLVER(RK8PD, rk8pd)
0185 
0186 /** \ingroup solvers
0187  *  \class GslAdaptiveRK8PDSolver
0188  *  \brief Adaptive Runge-Kutta Prince-Dormand (8,9) solver from GSL library
0189  */
0190 STEPCORE_DECLARE_GSLASOLVER(RK8PD, rk8pd)
0191 
0192 /** \ingroup solvers
0193  *  \class GslRK2IMPSolver
0194  *  \brief Runge-Kutta implicit second-order solver from GSL library
0195  */
0196 STEPCORE_DECLARE_GSLSOLVER(RK2IMP, rk2imp)
0197 
0198 /** \ingroup solvers
0199  *  \class GslAdaptiveRK2IMPSolver
0200  *  \brief Adaptive Runge-Kutta Prince-Dormand (8,9) solver from GSL library
0201  */
0202 STEPCORE_DECLARE_GSLASOLVER(RK2IMP, rk2imp)
0203 
0204 /** \ingroup solvers
0205  *  \class GslRK4IMPSolver
0206  *  \brief Runge-Kutta implicit fourth-order solver from GSL library
0207  */
0208 STEPCORE_DECLARE_GSLSOLVER(RK4IMP, rk4imp)
0209 
0210 /** \ingroup solvers
0211  *  \class GslAdaptiveRK4IMPSolver
0212  *  \brief Runge-Kutta implicit fourth-order solver from GSL library
0213  */
0214 STEPCORE_DECLARE_GSLASOLVER(RK4IMP, rk4imp)
0215 
0216 } // namespace StepCore
0217 
0218 #endif // defined(STEPCORE_WITH_GSL) || defined(Q_MOC_RUN)
0219 
0220 #endif // STEPCORE_GSLSOLVER_H
0221