step/stepcore/gslsolver.cc

0001 /*
0002     SPDX-FileCopyrightText: 2007 Vladimir Kuznetsov <ks.vladimir@gmail.com>
0003
0004     SPDX-License-Identifier: GPL-2.0-or-later
0005 */
0006
0007 #include "gslsolver.h"
0008
0009 #ifdef STEPCORE_WITH_GSL
0010
0011 #include "util.h"
0012 #include <cstring>
0013 #include <cmath>
0014
0015 namespace StepCore
0016 {
0017
0018 STEPCORE_META_OBJECT(GslGenericSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslGenericSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "GSL generic solver"), MetaObject::ABSTRACT, STEPCORE_SUPER_CLASS(Solver),)
0019
0020 STEPCORE_META_OBJECT(GslSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "GSL non-adaptive solver"), MetaObject::ABSTRACT,
0021     STEPCORE_SUPER_CLASS(GslGenericSolver),)
0022
0023 STEPCORE_META_OBJECT(GslAdaptiveSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslAdaptiveSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "GSL adaptive solver"), MetaObject::ABSTRACT,
0024     STEPCORE_SUPER_CLASS(GslGenericSolver),)
0025
0026 STEPCORE_META_OBJECT(GslRK2Solver, QT_TRANSLATE_NOOP("ObjectClass", "GslRK2Solver"), QT_TRANSLATE_NOOP("ObjectDescription", "Runge-Kutta second-order solver from GSL library"),
0027                         0, STEPCORE_SUPER_CLASS(GslSolver),)
0028 STEPCORE_META_OBJECT(GslAdaptiveRK2Solver, QT_TRANSLATE_NOOP("ObjectClass", "GslAdaptiveRK2Solver"), QT_TRANSLATE_NOOP("ObjectDescription", "Adaptive Runge-Kutta second-order solver from GSL library"),
0029                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
0030
0031 STEPCORE_META_OBJECT(GslRK4Solver, QT_TRANSLATE_NOOP("ObjectClass", "GslRK4Solver"), QT_TRANSLATE_NOOP("ObjectDescription", "Runge-Kutta classical fourth-order solver from GSL library"),
0032                         0, STEPCORE_SUPER_CLASS(GslSolver),)
0033 STEPCORE_META_OBJECT(GslAdaptiveRK4Solver, QT_TRANSLATE_NOOP("ObjectClass", "GslAdaptiveRK4Solver"), QT_TRANSLATE_NOOP("ObjectDescription", "Adaptive Runge-Kutta classical fourth-order solver from GSL library"),
0034                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
0035
0036 STEPCORE_META_OBJECT(GslRKF45Solver, QT_TRANSLATE_NOOP("ObjectClass", "GslRKF45Solver"), QT_TRANSLATE_NOOP("ObjectDescription", "Runge-Kutta-Fehlberg (4,5) solver from GSL library"),
0037                         0, STEPCORE_SUPER_CLASS(GslSolver),)
0038 STEPCORE_META_OBJECT(GslAdaptiveRKF45Solver, QT_TRANSLATE_NOOP("ObjectClass", "GslAdaptiveRKF45Solver"), QT_TRANSLATE_NOOP("ObjectDescription", "Adaptive Runge-Kutta-Fehlberg (4,5) solver from GSL library"),
0039                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
0040
0041 STEPCORE_META_OBJECT(GslRKCKSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslRKCKSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "Runge-Kutta Cash-Karp (4,5) solver from GSL library"),
0042                         0, STEPCORE_SUPER_CLASS(GslSolver),)
0043 STEPCORE_META_OBJECT(GslAdaptiveRKCKSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslAdaptiveRKCKSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "Adaptive Runge-Kutta Cash-Karp (4,5) solver from GSL library"),
0044                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
0045
0046 STEPCORE_META_OBJECT(GslRK8PDSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslRK8PDSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "Runge-Kutta Prince-Dormand (8,9) solver from GSL library"),
0047                         0, STEPCORE_SUPER_CLASS(GslSolver),)
0048 STEPCORE_META_OBJECT(GslAdaptiveRK8PDSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslAdaptiveRK8PDSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "Adaptive Runge-Kutta Prince-Dormand (8,9) solver from GSL library"),
0049                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
0050
0051 STEPCORE_META_OBJECT(GslRK2IMPSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslRK2IMPSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "Runge-Kutta implicit second-order solver from GSL library"),
0052                         0, STEPCORE_SUPER_CLASS(GslSolver),)
0053 STEPCORE_META_OBJECT(GslAdaptiveRK2IMPSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslAdaptiveRK2IMPSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "Adaptive Runge-Kutta implicit second-order solver from GSL library"),
0054                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
0055
0056 STEPCORE_META_OBJECT(GslRK4IMPSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslRK4IMPSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "Runge-Kutta implicit fourth-order solver from GSL library"),
0057                         0, STEPCORE_SUPER_CLASS(GslSolver),)
0058 STEPCORE_META_OBJECT(GslAdaptiveRK4IMPSolver, QT_TRANSLATE_NOOP("ObjectClass", "GslAdaptiveRK4IMPSolver"), QT_TRANSLATE_NOOP("ObjectDescription", "Adaptive Runge-Kutta implicit fourth-order solver from GSL library"),
0059                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
0060
0061
0062 void GslGenericSolver::init()
0063 {
0064     _yerr.resize(_dimension);
0065     _ytemp.resize(_dimension);
0066     _ydiff.resize(_dimension);
0067     _dydt_in.resize(_dimension);
0068     _dydt_out.resize(_dimension);
0069
0070     _gslStep = gsl_odeiv_step_alloc(_gslStepType, _dimension);
0071     STEPCORE_ASSERT_NOABORT(nullptr != _gslStep);
0072
0073     _gslSystem.function = gslFunction;
0074     _gslSystem.jacobian = nullptr;
0075     _gslSystem.dimension = _dimension;
0076     _gslSystem.params = this;
0077
0078     if(_adaptive) {
0079         _gslControl = gsl_odeiv_control_y_new(_toleranceAbs, _toleranceRel);
0080         STEPCORE_ASSERT_NOABORT(nullptr != _gslControl);
0081         _gslEvolve = gsl_odeiv_evolve_alloc(_dimension);
0082         STEPCORE_ASSERT_NOABORT(nullptr != _gslEvolve);
0083     } else {
0084         _gslControl = nullptr;
0085         _gslEvolve = nullptr;
0086     }
0087 }
0088
0089 void GslGenericSolver::fini()
0090 {
0091     if(_gslStep != nullptr) gsl_odeiv_step_free(_gslStep);
0092     if(_gslControl != nullptr) gsl_odeiv_control_free(_gslControl);
0093     if(_gslEvolve != nullptr) gsl_odeiv_evolve_free(_gslEvolve);
0094 }
0095
0096 int GslGenericSolver::gslFunction(double t, const double* y, double* f, void* params)
0097 {
0098     GslGenericSolver* s = static_cast<GslGenericSolver*>(params);
0099     return s->_function(t, y, nullptr, f, nullptr, s->_params);
0100 }
0101
0102 int GslGenericSolver::doCalcFn(double* t, const VectorXd* y,
0103             const VectorXd* yvar, VectorXd* f, VectorXd* fvar)
0104 {
0105     //int ret = GSL_ODEIV_FN_EVAL(&_gslSystem, *t, y, _ydiff);
0106     int ret = _function(*t, y->data(), yvar?yvar->data():nullptr, f ? f->data() : _ydiff.data(),
0107                         fvar?fvar->data():nullptr, _params);
0108     //if(f != NULL) std::memcpy(f, _ydiff, _dimension*sizeof(*f));
0109     return ret;
0110     //_hasSavedState = true;
0111 }
0112
0113 int GslGenericSolver::doEvolve(double* t, double t1, VectorXd* y, VectorXd* yvar)
0114 {
0115     STEPCORE_ASSERT_NOABORT(*t + _stepSize != *t);
0116     STEPCORE_ASSERT_NOABORT(*t != t1);
0117     //STEPCORE_ASSERT_NOABORT(_dimension != 0);
0118
0119     /*
0120     if(_hasSavedState) {
0121         std::memcpy(_ydiff_in, _ydiff_out, _dimension*sizeof(*_ydiff_in));
0122     } else {
0123         GSL_ODEIV_FN_EVAL(&_gslSystem, *t, y, _ydiff_in);
0124         _hasSavedState = true;
0125     }
0126     */
0127
0128     int gsl_result;
0129     _ytemp = *y;
0130
0131     if(!_adaptive) {
0132         gsl_result = GSL_ODEIV_FN_EVAL(&_gslSystem, *t, y->data(), _dydt_in.data());
0133         if(gsl_result != 0) return gsl_result;
0134     }
0135
0136     while(*t < t1) {
0137         double tt = *t;
0138         if(_adaptive) {
0139             gsl_odeiv_evolve_reset(_gslEvolve); // XXX
0140             gsl_result = gsl_odeiv_evolve_apply(_gslEvolve, _gslControl, _gslStep, &_gslSystem,
0141                                             &tt, t1, &_stepSize, _ytemp.data());
0142             _yerr = VectorXd::Map(_gslEvolve->yerr, _dimension);
0143         } else {
0144             STEPCORE_ASSERT_NOABORT(t1-tt > _stepSize/100);
0145             gsl_result = gsl_odeiv_step_apply(_gslStep, tt, (_stepSize < t1-tt ? _stepSize : t1-tt),
0146                                               _ytemp.data(), _yerr.data(), _dydt_in.data(),
0147                                               _dydt_out.data(), &_gslSystem);
0148             tt = _stepSize < t1-tt ? tt + _stepSize : t1;
0149         }
0150         if(gsl_result != 0) return gsl_result;
0151
0152         // XXX: Do we need to check it ?
0153         _localError = 0;
0154         _localErrorRatio = 0;
0155         for(int i=0; i<_dimension; ++i) {
0156             double error = fabs(_yerr[i]);
0157             if(error > _localError) _localError = error;
0158             double errorRatio = error / (_toleranceAbs + _toleranceRel * fabs(_ytemp[i]));
0159             if(errorRatio > _localErrorRatio) _localErrorRatio = errorRatio;
0160         }
0161         if(_localErrorRatio > 1.1) return ToleranceError;
0162
0163         *t = tt;
0164         *y = _ytemp;
0165         if(!_adaptive) _dydt_in = _dydt_out;
0166     }
0167
0168     return OK;
0169 }
0170
0171 } // namespace StepCore
0172
0173 #endif // STEPCORE_WITH_GSL
0174