File indexing completed on 2024-04-28 07:39:46
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