step/stepcore

gslsolver.cc

00001 /* This file is part of StepCore library.
00002    Copyright (C) 2007 Vladimir Kuznetsov <ks.vladimir@gmail.com>
00003 
00004    StepCore library is free software; you can redistribute it and/or modify
00005    it under the terms of the GNU General Public License as published by
00006    the Free Software Foundation; either version 2 of the License, or
00007    (at your option) any later version.
00008 
00009    StepCore library is distributed in the hope that it will be useful,
00010    but WITHOUT ANY WARRANTY; without even the implied warranty of
00011    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012    GNU General Public License for more details.
00013 
00014    You should have received a copy of the GNU General Public License
00015    along with StepCore; if not, write to the Free Software
00016    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017 */
00018 
00019 #include "gslsolver.h"
00020 
00021 #ifdef STEPCORE_WITH_GSL
00022 
00023 #include "util.h"
00024 #include <cstring>
00025 #include <cmath>
00026 
00027 namespace StepCore
00028 {
00029 
00030 STEPCORE_META_OBJECT(GslGenericSolver, "GSL generic solver", MetaObject::ABSTRACT, STEPCORE_SUPER_CLASS(Solver),)
00031 
00032 STEPCORE_META_OBJECT(GslSolver, "GSL non-adaptive solver", MetaObject::ABSTRACT,
00033     STEPCORE_SUPER_CLASS(GslGenericSolver),)
00034 
00035 STEPCORE_META_OBJECT(GslAdaptiveSolver, "GSL adaptive solver", MetaObject::ABSTRACT,
00036     STEPCORE_SUPER_CLASS(GslGenericSolver),)
00037 
00038 STEPCORE_META_OBJECT(GslRK2Solver, "Runge-Kutta second-order solver from GSL library",
00039                         0, STEPCORE_SUPER_CLASS(GslSolver),)
00040 STEPCORE_META_OBJECT(GslAdaptiveRK2Solver, "Adaptive Runge-Kutta second-order solver from GSL library",
00041                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
00042 
00043 STEPCORE_META_OBJECT(GslRK4Solver, "Runge-Kutta classical fourth-order solver from GSL library",
00044                         0, STEPCORE_SUPER_CLASS(GslSolver),)
00045 STEPCORE_META_OBJECT(GslAdaptiveRK4Solver, "Adaptive Runge-Kutta classical fourth-order solver from GSL library",
00046                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
00047 
00048 STEPCORE_META_OBJECT(GslRKF45Solver, "Runge-Kutta-Fehlberg (4,5) solver from GSL library",
00049                         0, STEPCORE_SUPER_CLASS(GslSolver),)
00050 STEPCORE_META_OBJECT(GslAdaptiveRKF45Solver, "Adaptive Runge-Kutta-Fehlberg (4,5) solver from GSL library",
00051                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
00052 
00053 STEPCORE_META_OBJECT(GslRKCKSolver, "Runge-Kutta Cash-Karp (4,5) solver from GSL library",
00054                         0, STEPCORE_SUPER_CLASS(GslSolver),)
00055 STEPCORE_META_OBJECT(GslAdaptiveRKCKSolver, "Adaptive Runge-Kutta Cash-Karp (4,5) solver from GSL library",
00056                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
00057 
00058 STEPCORE_META_OBJECT(GslRK8PDSolver, "Runge-Kutta Prince-Dormand (8,9) solver from GSL library",
00059                         0, STEPCORE_SUPER_CLASS(GslSolver),)
00060 STEPCORE_META_OBJECT(GslAdaptiveRK8PDSolver, "Adaptive Runge-Kutta Prince-Dormand (8,9) solver from GSL library",
00061                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
00062 
00063 STEPCORE_META_OBJECT(GslRK2IMPSolver, "Runge-Kutta implicit second-order solver from GSL library",
00064                         0, STEPCORE_SUPER_CLASS(GslSolver),)
00065 STEPCORE_META_OBJECT(GslAdaptiveRK2IMPSolver, "Adaptive Runge-Kutta implicit second-order solver from GSL library",
00066                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
00067 
00068 STEPCORE_META_OBJECT(GslRK4IMPSolver, "Runge-Kutta implicit fourth-order solver from GSL library",
00069                         0, STEPCORE_SUPER_CLASS(GslSolver),)
00070 STEPCORE_META_OBJECT(GslAdaptiveRK4IMPSolver, "Adaptive Runge-Kutta implicit fource-order solver from GSL library",
00071                         0, STEPCORE_SUPER_CLASS(GslAdaptiveSolver),)
00072 
00073 
00074 void GslGenericSolver::init()
00075 {
00076     _yerr = new double[_dimension];
00077     _ytemp = new double[_dimension];
00078     _ydiff = new double[_dimension];
00079     _dydt_in  = new double[_dimension];
00080     _dydt_out = new double[_dimension];
00081 
00082     _gslStep = gsl_odeiv_step_alloc(_gslStepType, _dimension);
00083     STEPCORE_ASSERT_NOABORT(NULL != _gslStep);
00084 
00085     _gslSystem.function = gslFunction;
00086     _gslSystem.jacobian = NULL;
00087     _gslSystem.dimension = _dimension;
00088     _gslSystem.params = this;
00089 
00090     if(_adaptive) {
00091         _gslControl = gsl_odeiv_control_y_new(_toleranceAbs, _toleranceRel);
00092         STEPCORE_ASSERT_NOABORT(NULL != _gslControl);
00093         _gslEvolve = gsl_odeiv_evolve_alloc(_dimension);
00094         STEPCORE_ASSERT_NOABORT(NULL != _gslEvolve);
00095     } else {
00096         _gslControl = NULL;
00097         _gslEvolve = NULL;
00098     }
00099 }
00100 
00101 void GslGenericSolver::fini()
00102 {
00103     if(_gslStep != NULL) gsl_odeiv_step_free(_gslStep);
00104     if(_gslControl != NULL) gsl_odeiv_control_free(_gslControl);
00105     if(_gslEvolve != NULL) gsl_odeiv_evolve_free(_gslEvolve);
00106     delete[] _dydt_out;
00107     delete[] _dydt_in;
00108     delete[] _ydiff;
00109     delete[] _ytemp;
00110     delete[] _yerr;
00111 }
00112 
00113 int GslGenericSolver::gslFunction(double t, const double* y, double* f, void* params)
00114 {
00115     GslGenericSolver* s = static_cast<GslGenericSolver*>(params);
00116     return s->_function(t, y, 0, f, 0, s->_params);
00117 }
00118 
00119 int GslGenericSolver::doCalcFn(double* t, const double* y,
00120             const double* yvar, double* f, double* fvar)
00121 {
00122     //int ret = GSL_ODEIV_FN_EVAL(&_gslSystem, *t, y, _ydiff);
00123     int ret = _function(*t, y, yvar, f ? f : _ydiff, fvar, _params);
00124     //if(f != NULL) std::memcpy(f, _ydiff, _dimension*sizeof(*f));
00125     return ret;
00126     //_hasSavedState = true;
00127 }
00128 
00129 int GslGenericSolver::doEvolve(double* t, double t1, double* y, double* yvar)
00130 {
00131     STEPCORE_ASSERT_NOABORT(*t + _stepSize != *t);
00132     STEPCORE_ASSERT_NOABORT(*t != t1);
00133     //STEPCORE_ASSERT_NOABORT(_dimension != 0);
00134 
00135     /*
00136     if(_hasSavedState) {
00137         std::memcpy(_ydiff_in, _ydiff_out, _dimension*sizeof(*_ydiff_in));
00138     } else {
00139         GSL_ODEIV_FN_EVAL(&_gslSystem, *t, y, _ydiff_in);
00140         _hasSavedState = true;
00141     }
00142     */
00143 
00144     int gsl_result;
00145     std::memcpy(_ytemp, y, _dimension*sizeof(*_ytemp));
00146 
00147     if(!_adaptive) {
00148         gsl_result = GSL_ODEIV_FN_EVAL(&_gslSystem, *t, y, _dydt_in);
00149         if(gsl_result != 0) return gsl_result;
00150     }
00151 
00152     while(*t < t1) {
00153         double tt = *t;
00154         if(_adaptive) {
00155             gsl_odeiv_evolve_reset(_gslEvolve); // XXX
00156             gsl_result = gsl_odeiv_evolve_apply(_gslEvolve, _gslControl, _gslStep, &_gslSystem,
00157                                             &tt, t1, &_stepSize, _ytemp);
00158             std::memcpy(_yerr, _gslEvolve->yerr, _dimension*sizeof(*_yerr));
00159         } else {
00160             STEPCORE_ASSERT_NOABORT(t1-tt > _stepSize/100);
00161             gsl_result = gsl_odeiv_step_apply(_gslStep, tt, (_stepSize < t1-tt ? _stepSize : t1-tt),
00162                                                 _ytemp, _yerr, _dydt_in, _dydt_out, &_gslSystem);
00163             tt = _stepSize < t1-tt ? tt + _stepSize : t1;
00164         }
00165         if(gsl_result != 0) return gsl_result;
00166 
00167         // XXX: Do we need to check it ?
00168         _localError = 0;
00169         _localErrorRatio = 0;
00170         for(int i=0; i<_dimension; ++i) {
00171             double error = fabs(_yerr[i]);
00172             if(error > _localError) _localError = error;
00173             double errorRatio = error / (_toleranceAbs + _toleranceRel * fabs(_ytemp[i]));
00174             if(errorRatio > _localErrorRatio) _localErrorRatio = errorRatio;
00175         }
00176         if(_localErrorRatio > 1.1) return ToleranceError;
00177 
00178         *t = tt;
00179         std::memcpy(y, _ytemp, _dimension*sizeof(*y));
00180         if(!_adaptive) std::memcpy(_dydt_in, _dydt_out, _dimension*sizeof(*_dydt_in));
00181     }
00182 
00183     return OK;
00184 }
00185 
00186 } // namespace StepCore
00187 
00188 #endif // STEPCORE_WITH_GSL
00189
step/stepcore

gslsolver.cc

step/stepcore

kdeedu
