OMCompiler/SimulationRuntime/cpp/Solver/CVode/CVode.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * This file belongs to the OpenModelica Run-Time System | ||
| 3 | * | ||
| 4 | * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings | ||
| 5 | * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights | ||
| 6 | * reserved. | ||
| 7 | * | ||
| 8 | * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE | ||
| 9 | * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY | ||
| 10 | * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S | ||
| 11 | * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL | ||
| 12 | * VERSION 3, ACCORDING TO RECIPIENTS CHOICE. | ||
| 13 | * | ||
| 14 | * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License | ||
| 15 | * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs: | ||
| 16 | * http://www.openmodelica.org or https://github.com/OpenModelica/ or | ||
| 17 | * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU | ||
| 18 | * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW | ||
| 19 | * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause. | ||
| 20 | * | ||
| 21 | * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 22 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY | ||
| 23 | * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF | ||
| 24 | * OSMC-PL. | ||
| 25 | * | ||
| 26 | */ | ||
| 27 | |||
| 28 | /** @addtogroup solverCvode | ||
| 29 | * | ||
| 30 | * @{ | ||
| 31 | */ | ||
| 32 | #include <Core/ModelicaDefine.h> | ||
| 33 | #include <Core/Modelica.h> | ||
| 34 | #include <Solver/CVode/CVode.h> | ||
| 35 | #include <Core/Math/Functions.h> | ||
| 36 | |||
| 37 | |||
| 38 | ✗ | Cvode::Cvode(IMixedSystem* system, ISolverSettings* settings) | |
| 39 | : SolverDefaultImplementation(system, settings), | ||
| 40 | ✗ | _cvodesettings(dynamic_cast<ISolverSettings*>(_settings)), | |
| 41 | ✗ | _cvodeMem(NULL), | |
| 42 | ✗ | _sunctx(NULL), | |
| 43 | ✗ | _z(NULL), | |
| 44 | ✗ | _zInit(NULL), | |
| 45 | ✗ | _zWrite(NULL), | |
| 46 | ✗ | _dimSys(0), | |
| 47 | ✗ | _cv_rt(0), | |
| 48 | ✗ | _outStps(0), | |
| 49 | ✗ | _locStps(0), | |
| 50 | ✗ | _idid(0), | |
| 51 | ✗ | _hOut(0.0), | |
| 52 | ✗ | _tOut(0.0), | |
| 53 | ✗ | _tZero(0.0), | |
| 54 | ✗ | _zeroSign(NULL), | |
| 55 | ✗ | _absTol(NULL), | |
| 56 | ✗ | _cvode_initialized(false), | |
| 57 | ✗ | _tLastEvent(0.0), | |
| 58 | ✗ | _event_n(0), | |
| 59 | ✗ | _properties(NULL), | |
| 60 | ✗ | _continuous_system(NULL), | |
| 61 | ✗ | _event_system(NULL), | |
| 62 | ✗ | _mixed_system(NULL), | |
| 63 | ✗ | _time_system(NULL), | |
| 64 | ✗ | _numberOfOdeEvaluations(0), | |
| 65 | ✗ | _delta(NULL), | |
| 66 | ✗ | _deltaInv(NULL), | |
| 67 | ✗ | _ysave(NULL), | |
| 68 | ✗ | _colorOfColumn(NULL), | |
| 69 | ✗ | _jacobianAIndex(NULL), | |
| 70 | ✗ | _jacobianALeadindex(NULL), | |
| 71 | ✗ | _CV_absTol(), | |
| 72 | ✗ | _tLastWrite(-1.0), | |
| 73 | ✗ | _bWritten(false), | |
| 74 | ✗ | _zeroFound(false), | |
| 75 | ✗ | _CV_y0(), | |
| 76 | ✗ | _CV_y(), | |
| 77 | ✗ | _CV_yWrite(), | |
| 78 | ✗ | _CV_ySolver(NULL), | |
| 79 | ✗ | _CV_linSol(NULL), | |
| 80 | ✗ | _CV_J(NULL), | |
| 81 | ✗ | _maxColors(0), | |
| 82 | ✗ | _jacobianANonzeros(0) | |
| 83 | { | ||
| 84 | ✗ | _data = ((void*) this); | |
| 85 | |||
| 86 | #ifdef RUNTIME_PROFILING | ||
| 87 | if (MeasureTime::getInstance() != NULL) | ||
| 88 | { | ||
| 89 | measureTimeFunctionsArray = new std::vector<MeasureTimeData*>(7, NULL); //0 calcFunction //1 solve ... //6 solver statistics | ||
| 90 | (*measureTimeFunctionsArray)[0] = new MeasureTimeData("calcFunction"); | ||
| 91 | (*measureTimeFunctionsArray)[1] = new MeasureTimeData("solve"); | ||
| 92 | (*measureTimeFunctionsArray)[2] = new MeasureTimeData("writeOutput"); | ||
| 93 | (*measureTimeFunctionsArray)[3] = new MeasureTimeData("evaluateZeroFuncs"); | ||
| 94 | (*measureTimeFunctionsArray)[4] = new MeasureTimeData("initialize"); | ||
| 95 | (*measureTimeFunctionsArray)[5] = new MeasureTimeData("stepCompleted"); | ||
| 96 | (*measureTimeFunctionsArray)[6] = new MeasureTimeData("solverStatistics"); | ||
| 97 | |||
| 98 | MeasureTime::addResultContentBlock(system->getModelName(), "cvode", measureTimeFunctionsArray); | ||
| 99 | measuredFunctionStartValues = MeasureTime::getZeroValues(); | ||
| 100 | measuredFunctionEndValues = MeasureTime::getZeroValues(); | ||
| 101 | solveFunctionStartValues = MeasureTime::getZeroValues(); | ||
| 102 | solveFunctionEndValues = MeasureTime::getZeroValues(); | ||
| 103 | solverValues = new MeasureTimeValuesSolver(); | ||
| 104 | |||
| 105 | delete (*measureTimeFunctionsArray)[6]->_sumMeasuredValues; | ||
| 106 | (*measureTimeFunctionsArray)[6]->_sumMeasuredValues = solverValues; | ||
| 107 | } | ||
| 108 | else | ||
| 109 | { | ||
| 110 | measureTimeFunctionsArray = new std::vector<MeasureTimeData*>(); | ||
| 111 | measuredFunctionStartValues = NULL; | ||
| 112 | measuredFunctionEndValues = NULL; | ||
| 113 | solveFunctionStartValues = NULL; | ||
| 114 | solveFunctionEndValues = NULL; | ||
| 115 | solverValues = NULL; | ||
| 116 | } | ||
| 117 | #endif | ||
| 118 | ✗ | } | |
| 119 | |||
| 120 | ✗ | Cvode::~Cvode() | |
| 121 | { | ||
| 122 | ✗ | if (_z) | |
| 123 | ✗ | delete[] _z; | |
| 124 | ✗ | if (_zInit) | |
| 125 | ✗ | delete[] _zInit; | |
| 126 | ✗ | if (_zeroSign) | |
| 127 | ✗ | delete[] _zeroSign; | |
| 128 | ✗ | if (_absTol) | |
| 129 | ✗ | delete[] _absTol; | |
| 130 | ✗ | if (_zWrite) | |
| 131 | ✗ | delete[] _zWrite; | |
| 132 | ✗ | if (_cvode_initialized) | |
| 133 | { | ||
| 134 | ✗ | N_VDestroy_Serial(_CV_y0); | |
| 135 | ✗ | N_VDestroy_Serial(_CV_y); | |
| 136 | ✗ | N_VDestroy_Serial(_CV_yWrite); | |
| 137 | ✗ | N_VDestroy_Serial(_CV_absTol); | |
| 138 | ✗ | N_VDestroy_Serial(_CV_ySolver); | |
| 139 | ✗ | SUNMatDestroy(_CV_J); | |
| 140 | ✗ | SUNLinSolFree(_CV_linSol); | |
| 141 | ✗ | CVodeFree(&_cvodeMem); | |
| 142 | ✗ | SUNContext_Free(&_sunctx); | |
| 143 | } | ||
| 144 | |||
| 145 | ✗ | if (_colorOfColumn) | |
| 146 | ✗ | delete[] _colorOfColumn; | |
| 147 | ✗ | if (_delta) | |
| 148 | ✗ | delete[] _delta; | |
| 149 | ✗ | if (_deltaInv) | |
| 150 | ✗ | delete[] _deltaInv; | |
| 151 | ✗ | if (_ysave) | |
| 152 | ✗ | delete[] _ysave; | |
| 153 | |||
| 154 | #ifdef RUNTIME_PROFILING | ||
| 155 | delete measuredFunctionStartValues; | ||
| 156 | delete measuredFunctionEndValues; | ||
| 157 | delete solveFunctionStartValues; | ||
| 158 | delete solveFunctionEndValues; | ||
| 159 | /* solverValues is owned by (*measureTimeFunctionsArray)[6] and freed by ~MeasureTimeData() */ | ||
| 160 | #endif | ||
| 161 | ✗ | } | |
| 162 | |||
| 163 | ✗ | void Cvode::initialize() | |
| 164 | { | ||
| 165 | ✗ | _properties = dynamic_cast<ISystemProperties*>(_system); | |
| 166 | ✗ | _continuous_system = dynamic_cast<IContinuous*>(_system); | |
| 167 | ✗ | _event_system = dynamic_cast<IEvent*>(_system); | |
| 168 | ✗ | _mixed_system = dynamic_cast<IMixedSystem*>(_system); | |
| 169 | ✗ | _time_system = dynamic_cast<ITime*>(_system); | |
| 170 | ✗ | IGlobalSettings* global_settings = dynamic_cast<ISolverSettings*>(_cvodesettings)->getGlobalSettings(); | |
| 171 | // Kennzeichnung, dass initialize()() (vor der Integration) aufgerufen wurde | ||
| 172 | ✗ | _idid = 5000; | |
| 173 | ✗ | _tLastEvent = 0.0; | |
| 174 | ✗ | _event_n = 0; | |
| 175 | ✗ | SolverDefaultImplementation::initialize(); | |
| 176 | ✗ | _dimSys = _continuous_system->getDimContinuousStates(); | |
| 177 | ✗ | _dimZeroFunc = _event_system->getDimZeroFunc(); | |
| 178 | |||
| 179 | ✗ | if (_dimSys == 0) | |
| 180 | ✗ | _dimSys = 1; // introduce dummy state | |
| 181 | |||
| 182 | ✗ | if (_dimSys <= 0) | |
| 183 | { | ||
| 184 | ✗ | _idid = -1; | |
| 185 | ✗ | throw ModelicaSimulationError(SOLVER, "Cvode::initialize()"); | |
| 186 | } | ||
| 187 | else | ||
| 188 | { | ||
| 189 | // Allocate state vectors, stages and temporary arrays | ||
| 190 | ✗ | if (_z) | |
| 191 | ✗ | delete[] _z; | |
| 192 | ✗ | if (_zInit) | |
| 193 | ✗ | delete[] _zInit; | |
| 194 | ✗ | if (_zWrite) | |
| 195 | ✗ | delete[] _zWrite; | |
| 196 | ✗ | if (_zeroSign) | |
| 197 | ✗ | delete[] _zeroSign; | |
| 198 | ✗ | if (_absTol) | |
| 199 | ✗ | delete[] _absTol; | |
| 200 | ✗ | if (_delta) | |
| 201 | ✗ | delete[] _delta; | |
| 202 | ✗ | if (_deltaInv) | |
| 203 | ✗ | delete[] _deltaInv; | |
| 204 | ✗ | if (_ysave) | |
| 205 | ✗ | delete[] _ysave; | |
| 206 | |||
| 207 | ✗ | _z = new double[_dimSys]; | |
| 208 | ✗ | _zInit = new double[_dimSys]; | |
| 209 | ✗ | _zWrite = new double[_dimSys]; | |
| 210 | ✗ | _zeroSign = new int[_dimZeroFunc]; | |
| 211 | ✗ | _absTol = new double[_dimSys]; | |
| 212 | ✗ | _delta = new double[_dimSys]; | |
| 213 | ✗ | _deltaInv = new double[_dimSys]; | |
| 214 | ✗ | _ysave = new double[_dimSys]; | |
| 215 | |||
| 216 | ✗ | memset(_z, 0, _dimSys * sizeof(double)); | |
| 217 | ✗ | memset(_zInit, 0, _dimSys * sizeof(double)); | |
| 218 | ✗ | memset(_ysave, 0, _dimSys * sizeof(double)); | |
| 219 | |||
| 220 | // Counter initialisieren | ||
| 221 | ✗ | _outStps = 0; | |
| 222 | |||
| 223 | ✗ | if (_cvodesettings->getDenseOutput()) | |
| 224 | { | ||
| 225 | // Ausgabeschrittweite | ||
| 226 | ✗ | _hOut = global_settings->gethOutput(); | |
| 227 | |||
| 228 | } | ||
| 229 | |||
| 230 | // Allocate memory for the solver | ||
| 231 | ✗ | if (SUNContext_Create(SUN_COMM_NULL, &_sunctx) != SUN_SUCCESS) | |
| 232 | ✗ | throw ModelicaSimulationError(SOLVER,"SUNDIALS_ERROR: SUNContext_Create failed"); | |
| 233 | /* Mute SUNDIALS' own logger: package level messages go to stderr/stdout by default, and we report solver failures ourselves. */ | ||
| 234 | { | ||
| 235 | ✗ | SUNLogger _sunlogger = NULL; | |
| 236 | ✗ | if (SUNContext_GetLogger(_sunctx, &_sunlogger) == SUN_SUCCESS && _sunlogger != NULL) { | |
| 237 | ✗ | SUNLogger_SetErrorFilename(_sunlogger, ""); | |
| 238 | ✗ | SUNLogger_SetWarningFilename(_sunlogger, ""); | |
| 239 | ✗ | SUNLogger_SetInfoFilename(_sunlogger, ""); | |
| 240 | ✗ | SUNLogger_SetDebugFilename(_sunlogger, ""); | |
| 241 | } | ||
| 242 | } | ||
| 243 | |||
| 244 | ✗ | _cvodeMem = CVodeCreate(CV_BDF, _sunctx); | |
| 245 | ✗ | if (check_flag((void*)_cvodeMem, "CVodeCreate", 0)) | |
| 246 | { | ||
| 247 | ✗ | _idid = -5; | |
| 248 | ✗ | throw ModelicaSimulationError(SOLVER,/*_idid,_tCurrent,*/"Cvode::initialize()"); | |
| 249 | } | ||
| 250 | |||
| 251 | // | ||
| 252 | // Make Cvode ready for integration | ||
| 253 | // | ||
| 254 | |||
| 255 | // Set initial values for CVODE | ||
| 256 | ✗ | _continuous_system->evaluateAll(IContinuous::CONTINUOUS); | |
| 257 | ✗ | _continuous_system->getContinuousStates(_zInit); | |
| 258 | ✗ | memcpy(_z, _zInit, _dimSys * sizeof(double)); | |
| 259 | |||
| 260 | // Get nominal values | ||
| 261 | ✗ | _absTol[0] = 1.0; // in case of dummy state | |
| 262 | ✗ | _continuous_system->getNominalStates(_absTol); | |
| 263 | ✗ | for (int i = 0; i < _dimSys; i++) | |
| 264 | ✗ | _absTol[i] *= dynamic_cast<ISolverSettings*>(_cvodesettings)->getATol(); | |
| 265 | |||
| 266 | ✗ | _CV_y0 = N_VMake_Serial(_dimSys, _zInit, _sunctx); | |
| 267 | ✗ | _CV_y = N_VMake_Serial(_dimSys, _z, _sunctx); | |
| 268 | ✗ | _CV_yWrite = N_VMake_Serial(_dimSys, _zWrite, _sunctx); | |
| 269 | ✗ | _CV_absTol = N_VMake_Serial(_dimSys, _absTol, _sunctx); | |
| 270 | ✗ | _CV_ySolver = N_VNew_Serial(_dimSys, _sunctx); | |
| 271 | |||
| 272 | ✗ | if (check_flag((void*)_CV_y0, "N_VMake_Serial", 0)) | |
| 273 | { | ||
| 274 | ✗ | _idid = -5; | |
| 275 | ✗ | throw ModelicaSimulationError(SOLVER, "Cvode::initialize()"); | |
| 276 | } | ||
| 277 | |||
| 278 | // Initialize Cvode (Initial values are required) | ||
| 279 | ✗ | _idid = CVodeInit(_cvodeMem, CV_fCallback, _tCurrent, _CV_y0); | |
| 280 | ✗ | if (_idid < 0) | |
| 281 | { | ||
| 282 | ✗ | _idid = -5; | |
| 283 | ✗ | throw ModelicaSimulationError(SOLVER, "Cvode::initialize()"); | |
| 284 | } | ||
| 285 | |||
| 286 | // Set Tolerances | ||
| 287 | ✗ | _idid = CVodeSVtolerances(_cvodeMem, dynamic_cast<ISolverSettings*>(_cvodesettings)->getRTol(), _CV_absTol); // RTOL and ATOL | |
| 288 | ✗ | if (_idid < 0) | |
| 289 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::initialize()"); | |
| 290 | |||
| 291 | // Set the pointer to user-defined data | ||
| 292 | ✗ | _idid = CVodeSetUserData(_cvodeMem, _data); | |
| 293 | ✗ | if (_idid < 0) | |
| 294 | ✗ | throw ModelicaSimulationError(SOLVER, "Cvode::initialize()"); | |
| 295 | |||
| 296 | ✗ | _idid = CVodeSetInitStep(_cvodeMem, 1e-6); // INITIAL STEPSIZE | |
| 297 | ✗ | if (_idid < 0) | |
| 298 | ✗ | throw ModelicaSimulationError(SOLVER, "Cvode::initialize()"); | |
| 299 | |||
| 300 | ✗ | _idid = CVodeSetMaxOrd(_cvodeMem, 5); // Max Order | |
| 301 | ✗ | if (_idid < 0) | |
| 302 | ✗ | throw ModelicaSimulationError(SOLVER, "CVoder::initialize()"); | |
| 303 | |||
| 304 | ✗ | _idid = CVodeSetMaxConvFails(_cvodeMem, 100); // Maximale Fehler im Konvergenztest | |
| 305 | ✗ | if (_idid < 0) | |
| 306 | ✗ | throw ModelicaSimulationError(SOLVER, "CVoder::initialize()"); | |
| 307 | |||
| 308 | ✗ | _idid = CVodeSetStabLimDet(_cvodeMem, SUNTRUE); // Stability Detection | |
| 309 | ✗ | if (_idid < 0) | |
| 310 | ✗ | throw ModelicaSimulationError(SOLVER, "CVoder::initialize()"); | |
| 311 | |||
| 312 | ✗ | _idid = CVodeSetMinStep(_cvodeMem, dynamic_cast<ISolverSettings*>(_cvodesettings)->getLowerLimit()); // MINIMUM STEPSIZE | |
| 313 | ✗ | if (_idid < 0) | |
| 314 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::initialize()"); | |
| 315 | |||
| 316 | ✗ | _idid = CVodeSetMaxStep(_cvodeMem, global_settings->getEndTime() / 10.0); // MAXIMUM STEPSIZE | |
| 317 | ✗ | if (_idid < 0) | |
| 318 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::initialize()"); | |
| 319 | |||
| 320 | ✗ | _idid = CVodeSetMaxNonlinIters(_cvodeMem, 5); // Max number of iterations | |
| 321 | ✗ | if (_idid < 0) | |
| 322 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::initialize()"); | |
| 323 | ✗ | _idid = CVodeSetMaxErrTestFails(_cvodeMem, 100); | |
| 324 | ✗ | if (_idid < 0) | |
| 325 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::initialize()"); | |
| 326 | |||
| 327 | ✗ | _idid = CVodeSetMaxNumSteps(_cvodeMem, 1e3); // Max Number of steps | |
| 328 | ✗ | if (_idid < 0) | |
| 329 | ✗ | throw ModelicaSimulationError(SOLVER,/*_idid,_tCurrent,*/"Cvode::initialize()"); | |
| 330 | |||
| 331 | // Initialize dense linear solver | ||
| 332 | ✗ | _CV_J = SUNDenseMatrix(_dimSys, _dimSys, _sunctx); | |
| 333 | #ifdef USE_SUNDIALS_LAPACK | ||
| 334 | _CV_linSol = SUNLinSol_LapackDense(_CV_ySolver, _CV_J, _sunctx); | ||
| 335 | #else | ||
| 336 | ✗ | _CV_linSol = SUNLinSol_Dense(_CV_ySolver, _CV_J, _sunctx); | |
| 337 | #endif | ||
| 338 | ✗ | _idid = CVodeSetLinearSolver(_cvodeMem, _CV_linSol, _CV_J); | |
| 339 | ✗ | if (_idid < 0) | |
| 340 | ✗ | throw ModelicaSimulationError(SOLVER, "Cvode::initialize()"); | |
| 341 | |||
| 342 | // Use own jacobian matrix | ||
| 343 | // Check if Colored Jacobians are worth to use | ||
| 344 | ✗ | _maxColors = _system->getAMaxColors(); | |
| 345 | ✗ | if (_maxColors < _dimSys && _continuous_system->getDimContinuousStates() > 0) | |
| 346 | { | ||
| 347 | // _idid = CVodeSetJacFn(_cvodeMem, &CV_JCallback); | ||
| 348 | // initializeColoredJac(); | ||
| 349 | } | ||
| 350 | |||
| 351 | ✗ | if (_idid < 0) | |
| 352 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::initialize()"); | |
| 353 | |||
| 354 | ✗ | if (_dimZeroFunc) | |
| 355 | { | ||
| 356 | ✗ | _idid = CVodeRootInit(_cvodeMem, _dimZeroFunc, &CV_ZerofCallback); | |
| 357 | |||
| 358 | ✗ | memset(_zeroSign, 0, _dimZeroFunc * sizeof(int)); | |
| 359 | ✗ | _idid = CVodeSetRootDirection(_cvodeMem, _zeroSign); | |
| 360 | ✗ | if (_idid < 0) | |
| 361 | ✗ | throw ModelicaSimulationError(SOLVER,/*_idid,_tCurrent,*/"CVode::initialize()"); | |
| 362 | ✗ | memset(_zeroSign, -1, _dimZeroFunc * sizeof(int)); | |
| 363 | ✗ | memset(_zeroVal, -1, _dimZeroFunc * sizeof(int)); | |
| 364 | |||
| 365 | } | ||
| 366 | |||
| 367 | |||
| 368 | ✗ | _cvode_initialized = true; | |
| 369 | |||
| 370 | ✗ | LOGGER_WRITE("Cvode: initialized", LC_SOLVER, LL_DEBUG); | |
| 371 | } | ||
| 372 | ✗ | } | |
| 373 | |||
| 374 | ✗ | void Cvode::solve(const SOLVERCALL action) | |
| 375 | { | ||
| 376 | ✗ | bool writeEventOutput = (_settings->getGlobalSettings()->getOutputPointType() == OPT_ALL); | |
| 377 | ✗ | bool writeOutput = !(_settings->getGlobalSettings()->getOutputPointType() == OPT_NONE); | |
| 378 | |||
| 379 | #ifdef RUNTIME_PROFILING | ||
| 380 | MEASURETIME_REGION_DEFINE(cvodeSolveFunctionHandler, "solve"); | ||
| 381 | if (MeasureTime::getInstance() != NULL) | ||
| 382 | { | ||
| 383 | MEASURETIME_START(solveFunctionStartValues, cvodeSolveFunctionHandler, "solve"); | ||
| 384 | } | ||
| 385 | #endif | ||
| 386 | |||
| 387 | ✗ | if (_cvodesettings && _system) | |
| 388 | { | ||
| 389 | // Solver und System für Integration vorbereiten | ||
| 390 | ✗ | if ((action & RECORDCALL) && (action & FIRST_CALL)) | |
| 391 | { | ||
| 392 | #ifdef RUNTIME_PROFILING | ||
| 393 | MEASURETIME_REGION_DEFINE(cvodeInitializeHandler, "CVodeInitialize"); | ||
| 394 | if (MeasureTime::getInstance() != NULL) | ||
| 395 | { | ||
| 396 | MEASURETIME_START(measuredFunctionStartValues, cvodeInitializeHandler, "CVodeInitialize"); | ||
| 397 | } | ||
| 398 | #endif | ||
| 399 | |||
| 400 | ✗ | initialize(); | |
| 401 | |||
| 402 | #ifdef RUNTIME_PROFILING | ||
| 403 | if (MeasureTime::getInstance() != NULL) | ||
| 404 | { | ||
| 405 | MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[4], cvodeInitializeHandler); | ||
| 406 | } | ||
| 407 | #endif | ||
| 408 | |||
| 409 | ✗ | if (writeOutput) | |
| 410 | ✗ | writeToFile(0, _tCurrent, _h); | |
| 411 | ✗ | _tLastWrite = 0; | |
| 412 | |||
| 413 | ✗ | return; | |
| 414 | } | ||
| 415 | |||
| 416 | ✗ | if ((action & RECORDCALL) && !(action & FIRST_CALL)) | |
| 417 | { | ||
| 418 | ✗ | writeToFile(_accStps, _tCurrent, _h); | |
| 419 | ✗ | return; | |
| 420 | } | ||
| 421 | |||
| 422 | // Nach einem TimeEvent wird der neue Zustand recorded | ||
| 423 | ✗ | if (action & RECALL) | |
| 424 | { | ||
| 425 | ✗ | _firstStep = true; | |
| 426 | ✗ | if (writeEventOutput) | |
| 427 | ✗ | writeToFile(0, _tCurrent, _h); | |
| 428 | ✗ | if (writeOutput) | |
| 429 | ✗ | writeCVodeOutput(_tCurrent, _h, _locStps); | |
| 430 | ✗ | _continuous_system->getContinuousStates(_z); | |
| 431 | } | ||
| 432 | |||
| 433 | // Solver soll fortfahren | ||
| 434 | ✗ | _solverStatus = ISolver::CONTINUE; | |
| 435 | |||
| 436 | ✗ | while ((_solverStatus & ISolver::CONTINUE) && !_interrupt) | |
| 437 | { | ||
| 438 | // Zuvor wurde initialize aufgerufen und hat funktioniert => RESET IDID | ||
| 439 | ✗ | if (_idid == 5000) | |
| 440 | ✗ | _idid = 0; | |
| 441 | |||
| 442 | // Solveraufruf | ||
| 443 | ✗ | if (_idid == 0) | |
| 444 | { | ||
| 445 | // Zähler zurücksetzen | ||
| 446 | ✗ | _accStps = 0; | |
| 447 | ✗ | _locStps = 0; | |
| 448 | |||
| 449 | // Solverstart | ||
| 450 | ✗ | CVodeCore(); | |
| 451 | |||
| 452 | } | ||
| 453 | |||
| 454 | // Integration war nicht erfolgreich und wurde auch nicht vom User unterbrochen | ||
| 455 | ✗ | if (_idid != 0 && _idid != 1) | |
| 456 | { | ||
| 457 | ✗ | _solverStatus = ISolver::SOLVERERROR; | |
| 458 | //throw ModelicaSimulationError(SOLVER,_idid,_tCurrent,"CVode::solve()"); | ||
| 459 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::solve()"); | |
| 460 | } | ||
| 461 | |||
| 462 | // Abbruchkriterium (erreichen der Endzeit) | ||
| 463 | ✗ | else if ((_tEnd - _tCurrent) <= dynamic_cast<ISolverSettings*>(_cvodesettings)->getEndTimeTol()) | |
| 464 | ✗ | _solverStatus = DONE; | |
| 465 | } | ||
| 466 | |||
| 467 | ✗ | _firstCall = false; | |
| 468 | |||
| 469 | } | ||
| 470 | else | ||
| 471 | { | ||
| 472 | |||
| 473 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::solve()"); | |
| 474 | } | ||
| 475 | |||
| 476 | #ifdef RUNTIME_PROFILING | ||
| 477 | if (MeasureTime::getInstance() != NULL) | ||
| 478 | { | ||
| 479 | MEASURETIME_END(solveFunctionStartValues, solveFunctionEndValues, (*measureTimeFunctionsArray)[1], cvodeSolveFunctionHandler); | ||
| 480 | |||
| 481 | long int nst, nfe, nsetups, netf, nni, ncfn; | ||
| 482 | int qlast, qcur; | ||
| 483 | sunrealtype h0u, hlast, hcur, tcur; | ||
| 484 | |||
| 485 | int flag; | ||
| 486 | |||
| 487 | flag = CVodeGetIntegratorStats(_cvodeMem, &nst, &nfe, &nsetups, &netf, &qlast, &qcur, &h0u, &hlast, &hcur, &tcur); | ||
| 488 | flag = CVodeGetNonlinSolvStats(_cvodeMem, &nni, &ncfn); | ||
| 489 | |||
| 490 | MeasureTimeValuesSolver solverVals = MeasureTimeValuesSolver(nfe, netf); | ||
| 491 | (*measureTimeFunctionsArray)[6]->_sumMeasuredValues->_numCalcs += nst; | ||
| 492 | (*measureTimeFunctionsArray)[6]->_sumMeasuredValues->add(&solverVals); | ||
| 493 | } | ||
| 494 | #endif | ||
| 495 | } | ||
| 496 | ✗ | bool Cvode::isInterrupted() | |
| 497 | { | ||
| 498 | ✗ | if (_interrupt) | |
| 499 | { | ||
| 500 | ✗ | _solverStatus = DONE; | |
| 501 | ✗ | return true; | |
| 502 | } | ||
| 503 | else | ||
| 504 | { | ||
| 505 | return false; | ||
| 506 | } | ||
| 507 | } | ||
| 508 | ✗ | void Cvode::CVodeCore() | |
| 509 | { | ||
| 510 | ✗ | _idid = CVodeReInit(_cvodeMem, _tCurrent, _CV_y); | |
| 511 | ✗ | _idid = CVodeSetStopTime(_cvodeMem, _tEnd); | |
| 512 | ✗ | _idid = CVodeSetInitStep(_cvodeMem, 1e-12); | |
| 513 | ✗ | if (_idid < 0) | |
| 514 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::ReInit"); | |
| 515 | |||
| 516 | ✗ | bool writeEventOutput = (_settings->getGlobalSettings()->getOutputPointType() == OPT_ALL); | |
| 517 | ✗ | bool writeOutput = !(_settings->getGlobalSettings()->getOutputPointType() == OPT_NONE); | |
| 518 | |||
| 519 | ✗ | while ((_solverStatus & ISolver::CONTINUE) && !_interrupt) | |
| 520 | { | ||
| 521 | ✗ | _cv_rt = CVode(_cvodeMem, _tEnd, _CV_y, &_tCurrent, CV_ONE_STEP); | |
| 522 | |||
| 523 | ✗ | _idid = CVodeGetNumSteps(_cvodeMem, &_locStps); | |
| 524 | ✗ | if (_idid != CV_SUCCESS) | |
| 525 | ✗ | throw ModelicaSimulationError(SOLVER, "CVodeGetNumSteps failed. The cvode mem pointer is NULL"); | |
| 526 | |||
| 527 | ✗ | _idid = CVodeGetLastStep(_cvodeMem, &_h); | |
| 528 | ✗ | if (_idid != CV_SUCCESS) | |
| 529 | ✗ | throw ModelicaSimulationError(SOLVER, "CVodeGetLastStep failed. The cvode mem pointer is NULL"); | |
| 530 | |||
| 531 | //set completed step to system and check if terminate was called | ||
| 532 | ✗ | if (_continuous_system->stepCompleted(_tCurrent)) | |
| 533 | ✗ | _solverStatus = DONE; | |
| 534 | |||
| 535 | //Check if there was at least one output-point within the last solver interval | ||
| 536 | // -> Write output if true | ||
| 537 | ✗ | if (writeOutput) | |
| 538 | { | ||
| 539 | |||
| 540 | try | ||
| 541 | { | ||
| 542 | ✗ | writeCVodeOutput(_tCurrent, _h, _locStps); | |
| 543 | } | ||
| 544 | ✗ | catch (std::exception& ex) | |
| 545 | { | ||
| 546 | |||
| 547 | ✗ | if (!(_cv_rt == CV_ROOT_RETURN)) // if a zero crossing was dected before the event iteration was called and writeoutput throws and error | |
| 548 | // for this time step the event iteration evaluates the system with corrected values. | ||
| 549 | { | ||
| 550 | ✗ | string error_msg = string("CVode write output failed") + ex.what(); | |
| 551 | ✗ | throw ModelicaSimulationError(SOLVER, error_msg); | |
| 552 | } | ||
| 553 | ✗ | } | |
| 554 | } | ||
| 555 | |||
| 556 | #ifdef RUNTIME_PROFILING | ||
| 557 | MEASURETIME_REGION_DEFINE(cvodeStepCompletedHandler, "CVodeStepCompleted"); | ||
| 558 | if (MeasureTime::getInstance() != NULL) | ||
| 559 | { | ||
| 560 | MEASURETIME_START(measuredFunctionStartValues, cvodeStepCompletedHandler, "CVodeStepCompleted"); | ||
| 561 | } | ||
| 562 | #endif | ||
| 563 | |||
| 564 | |||
| 565 | |||
| 566 | #ifdef RUNTIME_PROFILING | ||
| 567 | if (MeasureTime::getInstance() != NULL) | ||
| 568 | { | ||
| 569 | MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[5], cvodeStepCompletedHandler); | ||
| 570 | } | ||
| 571 | #endif | ||
| 572 | |||
| 573 | // Perform state selection | ||
| 574 | ✗ | bool state_selection = stateSelection(); | |
| 575 | ✗ | if (state_selection) | |
| 576 | ✗ | _continuous_system->getContinuousStates(_z); | |
| 577 | |||
| 578 | ✗ | _zeroFound = false; | |
| 579 | |||
| 580 | // Check if step was successful | ||
| 581 | ✗ | if (check_flag(&_cv_rt, "CVode", 1)) | |
| 582 | { | ||
| 583 | ✗ | _solverStatus = ISolver::SOLVERERROR; | |
| 584 | ✗ | break; | |
| 585 | } | ||
| 586 | |||
| 587 | // A root was found | ||
| 588 | ✗ | if ((_cv_rt == CV_ROOT_RETURN) && !isInterrupted()) | |
| 589 | { | ||
| 590 | // CVode is setting _tCurrent to the time where the first event occurred | ||
| 591 | ✗ | double _abs = fabs(_tLastEvent - _tCurrent); | |
| 592 | ✗ | _zeroFound = true; | |
| 593 | |||
| 594 | ✗ | if ((_abs < 1e-3) && _event_n == 0) | |
| 595 | { | ||
| 596 | ✗ | _tLastEvent = _tCurrent; | |
| 597 | ✗ | _event_n++; | |
| 598 | } | ||
| 599 | ✗ | else if ((_abs < 1e-3) && (_event_n >= 1 && _event_n < 500)) | |
| 600 | { | ||
| 601 | ✗ | _event_n++; | |
| 602 | } | ||
| 603 | ✗ | else if ((_abs >= 1e-3)) | |
| 604 | { | ||
| 605 | //restart event counter | ||
| 606 | ✗ | _tLastEvent = _tCurrent; | |
| 607 | ✗ | _event_n = 0; | |
| 608 | } | ||
| 609 | else{ | ||
| 610 | ✗ | std::stringstream zeros; | |
| 611 | ✗ | _idid = CVodeGetRootInfo(_cvodeMem, _zeroSign); | |
| 612 | ✗ | for (int i = 0; i < _dimZeroFunc; i++){ | |
| 613 | ✗ | if (_zeroSign[i] != 0) | |
| 614 | ✗ | zeros << i << " "; | |
| 615 | } | ||
| 616 | ✗ | throw ModelicaSimulationError(EVENT_HANDLING, "Number of events of zero function(s) " + zeros.str() + "exceeded in time interval " + to_string(_abs) + " at time " + to_string(_tCurrent)); | |
| 617 | ✗ | } | |
| 618 | |||
| 619 | // CVode has interpolated the states at time 'tCurrent' | ||
| 620 | ✗ | _time_system->setTime(_tCurrent); | |
| 621 | |||
| 622 | // To get steep steps in the result file, two value points (P1 and P2) must be added | ||
| 623 | // | ||
| 624 | // Y | (P2) X........... | ||
| 625 | // | : | ||
| 626 | // | : | ||
| 627 | // |........X (P1) | ||
| 628 | // |----------------------------------> | ||
| 629 | // | ^ t | ||
| 630 | // _tCurrent | ||
| 631 | |||
| 632 | // Write the values of (P1) | ||
| 633 | ✗ | if (writeEventOutput) | |
| 634 | { | ||
| 635 | |||
| 636 | try | ||
| 637 | { | ||
| 638 | ✗ | _continuous_system->evaluateAll(IContinuous::CONTINUOUS); | |
| 639 | } | ||
| 640 | ✗ | catch (std::exception& ex) | |
| 641 | { | ||
| 642 | |||
| 643 | // if a zero crossing was dected before the event iteration was called and evalutateAll throws and error | ||
| 644 | // for this time step the event iteration evaluates the system with corrected values. | ||
| 645 | |||
| 646 | ✗ | } | |
| 647 | |||
| 648 | ✗ | writeToFile(0, _tCurrent, _h); | |
| 649 | } | ||
| 650 | |||
| 651 | ✗ | _idid = CVodeGetRootInfo(_cvodeMem, _zeroSign); | |
| 652 | |||
| 653 | ✗ | for (int i = 0; i < _dimZeroFunc; i++) | |
| 654 | ✗ | _events[i] = bool(_zeroSign[i]); | |
| 655 | |||
| 656 | ✗ | if (_mixed_system->handleSystemEvents(_events)) | |
| 657 | { | ||
| 658 | // State variables were reinitialized, thus we have to give these values to the cvode-solver | ||
| 659 | // Take care about the memory regions, _z is the same like _CV_y | ||
| 660 | ✗ | _continuous_system->getContinuousStates(_z); | |
| 661 | } | ||
| 662 | } | ||
| 663 | |||
| 664 | ✗ | if ((_zeroFound || state_selection) && !isInterrupted()) | |
| 665 | { | ||
| 666 | // Write the values of (P2) | ||
| 667 | ✗ | if (writeEventOutput) | |
| 668 | { | ||
| 669 | // If we want to write the event-results, we should evaluate the whole system again | ||
| 670 | ✗ | _continuous_system->evaluateAll(IContinuous::CONTINUOUS); | |
| 671 | ✗ | writeToFile(0, _tCurrent, _h); | |
| 672 | } | ||
| 673 | |||
| 674 | ✗ | _idid = CVodeReInit(_cvodeMem, _tCurrent, _CV_y); | |
| 675 | ✗ | if (_idid < 0) | |
| 676 | ✗ | throw ModelicaSimulationError(SOLVER, "CVode::ReInit()"); | |
| 677 | |||
| 678 | // Der Eventzeitpunkt kann auf der Endzeit liegen (Time-Events). In diesem Fall wird der Solver beendet, da CVode sonst eine interne Warnung schmeißt | ||
| 679 | ✗ | if (_tCurrent == _tEnd) | |
| 680 | ✗ | _cv_rt = CV_TSTOP_RETURN; | |
| 681 | ✗ | if (_continuous_system->stepCompleted(_tCurrent)) | |
| 682 | ✗ | _solverStatus = DONE; | |
| 683 | } | ||
| 684 | |||
| 685 | // Zähler für die Anzahl der ausgegebenen Schritte erhöhen | ||
| 686 | ✗ | ++_outStps; | |
| 687 | ✗ | _tLastSuccess = _tCurrent; | |
| 688 | |||
| 689 | ✗ | if (_cv_rt == CV_TSTOP_RETURN) | |
| 690 | { | ||
| 691 | ✗ | _time_system->setTime(_tEnd); | |
| 692 | //Solver has finished calculation - calculate the final values | ||
| 693 | ✗ | _continuous_system->setContinuousStates(NV_DATA_S(_CV_y)); | |
| 694 | ✗ | _continuous_system->evaluateAll(IContinuous::CONTINUOUS); | |
| 695 | ✗ | if (writeOutput) | |
| 696 | ✗ | writeToFile(0, _tEnd, _h); | |
| 697 | |||
| 698 | ✗ | _accStps += _locStps; | |
| 699 | ✗ | _solverStatus = DONE; | |
| 700 | } | ||
| 701 | } | ||
| 702 | ✗ | } | |
| 703 | |||
| 704 | ✗ | void Cvode::setTimeOut(unsigned int time_out) | |
| 705 | { | ||
| 706 | ✗ | SimulationMonitor::setTimeOut(time_out); | |
| 707 | ✗ | } | |
| 708 | ✗ | void Cvode::stop() | |
| 709 | { | ||
| 710 | ✗ | SimulationMonitor::stop(); | |
| 711 | ✗ | } | |
| 712 | |||
| 713 | ✗ | void Cvode::writeCVodeOutput(const double &time, const double &h, const int &stp) | |
| 714 | { | ||
| 715 | #ifdef RUNTIME_PROFILING | ||
| 716 | MEASURETIME_REGION_DEFINE(cvodeWriteOutputHandler, "CVodeWriteOutput"); | ||
| 717 | if (MeasureTime::getInstance() != NULL) | ||
| 718 | { | ||
| 719 | MEASURETIME_START(measuredFunctionStartValues, cvodeWriteOutputHandler, "CVodeWriteOutput"); | ||
| 720 | } | ||
| 721 | #endif | ||
| 722 | |||
| 723 | ✗ | if (stp > 0) | |
| 724 | { | ||
| 725 | ✗ | if (_cvodesettings->getDenseOutput()) | |
| 726 | { | ||
| 727 | ✗ | _bWritten = false; | |
| 728 | /* double *oldValues = NULL;*/ | ||
| 729 | |||
| 730 | //We have to find all output-points within the last solver step | ||
| 731 | ✗ | while (_tLastWrite + dynamic_cast<ISolverSettings*>(_cvodesettings)->getGlobalSettings()->gethOutput() <= time) | |
| 732 | { | ||
| 733 | ✗ | if (!_bWritten) | |
| 734 | { | ||
| 735 | ✗ | _continuous_system->restoreOldValues(); | |
| 736 | ////Rescue the calculated derivatives | ||
| 737 | // oldValues = new double[_continuous_system->getDimRHS()]; | ||
| 738 | // _continuous_system->getRHS(oldValues); | ||
| 739 | |||
| 740 | } | ||
| 741 | ✗ | _bWritten = true; | |
| 742 | ✗ | _tLastWrite = _tLastWrite + dynamic_cast<ISolverSettings*>(_cvodesettings)->getGlobalSettings()->gethOutput(); | |
| 743 | //Get the state vars at the output-point (interpolated) | ||
| 744 | ✗ | _idid = CVodeGetDky(_cvodeMem, _tLastWrite, 0, _CV_yWrite); | |
| 745 | ✗ | _time_system->setTime(_tLastWrite); | |
| 746 | ✗ | _continuous_system->setContinuousStates(NV_DATA_S(_CV_yWrite)); | |
| 747 | ✗ | _continuous_system->evaluateAll(IContinuous::CONTINUOUS); | |
| 748 | ✗ | SolverDefaultImplementation::writeToFile(stp, _tLastWrite, h); | |
| 749 | } //end if time -_tLastWritten | ||
| 750 | ✗ | if (_bWritten) | |
| 751 | { | ||
| 752 | ✗ | _time_system->setTime(time); | |
| 753 | ✗ | _continuous_system->setContinuousStates(_z); | |
| 754 | ✗ | _continuous_system->restoreNewValues(); | |
| 755 | /* _continuous_system->setStateDerivatives(oldValues); | ||
| 756 | delete[] oldValues;*/ | ||
| 757 | } | ||
| 758 | ✗ | else if (time == _tEnd && _tLastWrite != time) | |
| 759 | { | ||
| 760 | ✗ | _idid = CVodeGetDky(_cvodeMem, time, 0, _CV_y); | |
| 761 | ✗ | _time_system->setTime(time); | |
| 762 | ✗ | _continuous_system->setContinuousStates(NV_DATA_S(_CV_y)); | |
| 763 | ✗ | _continuous_system->evaluateAll(IContinuous::CONTINUOUS); | |
| 764 | ✗ | SolverDefaultImplementation::writeToFile(stp, _tEnd, h); | |
| 765 | } | ||
| 766 | } | ||
| 767 | else | ||
| 768 | { | ||
| 769 | ✗ | SolverDefaultImplementation::writeToFile(stp, time, h); | |
| 770 | } | ||
| 771 | } | ||
| 772 | #ifdef RUNTIME_PROFILING | ||
| 773 | if (MeasureTime::getInstance() != NULL) | ||
| 774 | { | ||
| 775 | MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[2], cvodeWriteOutputHandler); | ||
| 776 | } | ||
| 777 | #endif | ||
| 778 | ✗ | } | |
| 779 | |||
| 780 | ✗ | bool Cvode::stateSelection() | |
| 781 | { | ||
| 782 | ✗ | return SolverDefaultImplementation::stateSelection(); | |
| 783 | } | ||
| 784 | ✗ | int Cvode::calcFunction(const double& time, const double* y, double* f) | |
| 785 | { | ||
| 786 | #ifdef RUNTIME_PROFILING | ||
| 787 | MEASURETIME_REGION_DEFINE(cvodeCalcFunctionHandler, "CVodeCalcFunction"); | ||
| 788 | if (MeasureTime::getInstance() != NULL) | ||
| 789 | { | ||
| 790 | MEASURETIME_START(measuredFunctionStartValues, cvodeCalcFunctionHandler, "CVodeCalcFunction"); | ||
| 791 | } | ||
| 792 | #endif | ||
| 793 | |||
| 794 | int returnValue = 0; | ||
| 795 | try | ||
| 796 | { | ||
| 797 | ✗ | f[0] = 0.0; // in case of dummy state | |
| 798 | ✗ | _time_system->setTime(time); | |
| 799 | ✗ | _continuous_system->setContinuousStates(y); | |
| 800 | ✗ | _continuous_system->evaluateODE(IContinuous::CONTINUOUS); | |
| 801 | ✗ | _continuous_system->getRHS(f); | |
| 802 | ✗ | _numberOfOdeEvaluations++; | |
| 803 | } //workaround until exception can be catch from c- libraries | ||
| 804 | ✗ | catch (std::exception & ex) | |
| 805 | { | ||
| 806 | |||
| 807 | //cerr << "CVode integration error: " << diagnostic_information(ex); | ||
| 808 | returnValue = 1; | ||
| 809 | ✗ | } | |
| 810 | |||
| 811 | #ifdef RUNTIME_PROFILING | ||
| 812 | if (MeasureTime::getInstance() != NULL) | ||
| 813 | { | ||
| 814 | MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[0], cvodeCalcFunctionHandler); | ||
| 815 | } | ||
| 816 | #endif | ||
| 817 | |||
| 818 | ✗ | return returnValue; | |
| 819 | } | ||
| 820 | |||
| 821 | ✗ | int Cvode::CV_fCallback(double t, N_Vector y, N_Vector ydot, void *user_data) | |
| 822 | { | ||
| 823 | ✗ | return ((Cvode*)user_data)->calcFunction(t, NV_DATA_S(y), NV_DATA_S(ydot)); | |
| 824 | |||
| 825 | } | ||
| 826 | |||
| 827 | ✗ | void Cvode::giveZeroVal(const double &t, const double *y, double *zeroValue) | |
| 828 | { | ||
| 829 | #ifdef RUNTIME_PROFILING | ||
| 830 | MEASURETIME_REGION_DEFINE(cvodeEvalZeroHandler, "evaluateZeroFuncs"); | ||
| 831 | if (MeasureTime::getInstance() != NULL) | ||
| 832 | { | ||
| 833 | MEASURETIME_START(measuredFunctionStartValues, cvodeEvalZeroHandler, "evaluateZeroFuncs"); | ||
| 834 | } | ||
| 835 | #endif | ||
| 836 | |||
| 837 | ✗ | _time_system->setTime(t); | |
| 838 | ✗ | _continuous_system->setContinuousStates(y); | |
| 839 | |||
| 840 | // System aktualisieren | ||
| 841 | ✗ | _continuous_system->evaluateZeroFuncs(IContinuous::DISCRETE); | |
| 842 | |||
| 843 | ✗ | _event_system->getZeroFunc(zeroValue); | |
| 844 | |||
| 845 | #ifdef RUNTIME_PROFILING | ||
| 846 | if (MeasureTime::getInstance() != NULL) | ||
| 847 | { | ||
| 848 | MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[3], cvodeEvalZeroHandler); | ||
| 849 | } | ||
| 850 | #endif | ||
| 851 | ✗ | } | |
| 852 | |||
| 853 | ✗ | int Cvode::CV_ZerofCallback(double t, N_Vector y, double *zeroval, void *user_data) | |
| 854 | { | ||
| 855 | ✗ | ((Cvode*)user_data)->giveZeroVal(t, NV_DATA_S(y), zeroval); | |
| 856 | |||
| 857 | ✗ | return (0); | |
| 858 | } | ||
| 859 | |||
| 860 | ✗ | int Cvode::CV_JCallback(long int N, double t, N_Vector y, N_Vector fy, SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) | |
| 861 | { | ||
| 862 | ✗ | return ((Cvode*)user_data)->calcJacobian(t, N, tmp1, tmp2, tmp3, NV_DATA_S(y), fy, Jac); | |
| 863 | |||
| 864 | } | ||
| 865 | |||
| 866 | ✗ | int Cvode::calcJacobian(double t, long int N, N_Vector fHelp, N_Vector errorWeight, N_Vector jthCol, double* y, N_Vector fy, SUNMatrix Jac) | |
| 867 | { | ||
| 868 | try | ||
| 869 | { | ||
| 870 | int l, g; | ||
| 871 | double fnorm, minInc, *f_data, *fHelp_data, *errorWeight_data, h, srur, delta_inv; | ||
| 872 | |||
| 873 | ✗ | f_data = NV_DATA_S(fy); | |
| 874 | ✗ | errorWeight_data = NV_DATA_S(errorWeight); | |
| 875 | ✗ | fHelp_data = NV_DATA_S(fHelp); | |
| 876 | |||
| 877 | |||
| 878 | //Get relevant info | ||
| 879 | ✗ | _idid = CVodeGetErrWeights(_cvodeMem, errorWeight); | |
| 880 | ✗ | if (_idid < 0) | |
| 881 | { | ||
| 882 | ✗ | _idid = -5; | |
| 883 | ✗ | throw ModelicaSimulationError(SOLVER, "Cvode::calcJacobian()"); | |
| 884 | } | ||
| 885 | ✗ | _idid = CVodeGetCurrentStep(_cvodeMem, &h); | |
| 886 | ✗ | if (_idid < 0) | |
| 887 | { | ||
| 888 | ✗ | _idid = -5; | |
| 889 | ✗ | throw ModelicaSimulationError(SOLVER, "Cvode::calcJacobian()"); | |
| 890 | } | ||
| 891 | |||
| 892 | srur = sqrt(UROUND); | ||
| 893 | |||
| 894 | ✗ | fnorm = N_VWrmsNorm(fy, errorWeight); | |
| 895 | ✗ | minInc = (fnorm != 0.0) ? | |
| 896 | ✗ | (1000.0 * abs(h) * UROUND * N * fnorm) : 1.0; | |
| 897 | |||
| 898 | ✗ | for (int j = 0; j < N; j++) | |
| 899 | { | ||
| 900 | ✗ | _delta[j] = max(srur*abs(y[j]), minInc / errorWeight_data[j]); | |
| 901 | } | ||
| 902 | |||
| 903 | ✗ | for (int j = 0; j < N; j++) | |
| 904 | { | ||
| 905 | ✗ | _deltaInv[j] = 1 / _delta[j]; | |
| 906 | } | ||
| 907 | |||
| 908 | ✗ | if (_jacobianANonzeros != 0) | |
| 909 | { | ||
| 910 | ✗ | for (int color = 1; color <= _maxColors; color++) | |
| 911 | { | ||
| 912 | ✗ | for (int k = 0; k < _dimSys; k++) | |
| 913 | { | ||
| 914 | ✗ | if ((_colorOfColumn[k]) == color) | |
| 915 | { | ||
| 916 | ✗ | _ysave[k] = y[k]; | |
| 917 | ✗ | y[k] += _delta[k]; | |
| 918 | } | ||
| 919 | } | ||
| 920 | |||
| 921 | ✗ | calcFunction(t, y, fHelp_data); | |
| 922 | |||
| 923 | ✗ | for (int k = 0; k < _dimSys; k++) | |
| 924 | { | ||
| 925 | ✗ | if ((_colorOfColumn[k]) == color) | |
| 926 | { | ||
| 927 | ✗ | y[k] = _ysave[k]; | |
| 928 | |||
| 929 | ✗ | int startOfColumn = k * _dimSys; | |
| 930 | ✗ | for (int j = _jacobianALeadindex[k]; j < _jacobianALeadindex[k + 1]; j++) | |
| 931 | { | ||
| 932 | ✗ | l = _jacobianAIndex[j]; | |
| 933 | ✗ | g = l + startOfColumn; | |
| 934 | ✗ | SM_DATA_D(Jac)[g] = (fHelp_data[l] - f_data[l]) * _deltaInv[k]; | |
| 935 | } | ||
| 936 | } | ||
| 937 | } | ||
| 938 | } | ||
| 939 | } | ||
| 940 | |||
| 941 | |||
| 942 | |||
| 943 | |||
| 944 | |||
| 945 | } | ||
| 946 | |||
| 947 | |||
| 948 | |||
| 949 | |||
| 950 | |||
| 951 | |||
| 952 | /* | ||
| 953 | //Calculation of J without colouring | ||
| 954 | for (j = 0; j < N; j++) | ||
| 955 | { | ||
| 956 | |||
| 957 | |||
| 958 | //N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol); | ||
| 959 | |||
| 960 | _ysave[j] = y[j]; | ||
| 961 | |||
| 962 | y[j] += _delta[j]; | ||
| 963 | |||
| 964 | calcFunction(t, y, fHelp_data); | ||
| 965 | |||
| 966 | y[j] = _ysave[j]; | ||
| 967 | |||
| 968 | delta_inv = 1.0/_delta[j]; | ||
| 969 | N_VLinearSum(delta_inv, fHelp, -delta_inv, fy, jthCol); | ||
| 970 | |||
| 971 | for(int i=0; i<_dimSys; ++i) | ||
| 972 | { | ||
| 973 | SM_DATA_D(Jac)[i+j*_dimSys] = NV_Ith_S(jthCol,i); | ||
| 974 | } | ||
| 975 | |||
| 976 | //DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol); | ||
| 977 | } | ||
| 978 | */ | ||
| 979 | |||
| 980 | //workaround until exception can be catch from c- libraries | ||
| 981 | ✗ | catch (std::exception & ex) | |
| 982 | { | ||
| 983 | |||
| 984 | ✗ | cerr << "CVode integration error: " << ex.what(); | |
| 985 | return 1; | ||
| 986 | ✗ | } | |
| 987 | |||
| 988 | |||
| 989 | ✗ | return 0; | |
| 990 | } | ||
| 991 | |||
| 992 | ✗ | void Cvode::initializeColoredJac() | |
| 993 | { | ||
| 994 | |||
| 995 | ✗ | if (_colorOfColumn) | |
| 996 | ✗ | delete[] _colorOfColumn; | |
| 997 | ✗ | _colorOfColumn = new int[_dimSys]; | |
| 998 | ✗ | _system->getAColorOfColumn(_colorOfColumn, _dimSys); | |
| 999 | |||
| 1000 | // _system->getJacobian(_jacobianA); | ||
| 1001 | //_jacobianANonzeros = boost::numeric::bindings::traits::spmatrix_num_nonzeros (_jacobianA); | ||
| 1002 | // _jacobianAIndex = bindings::begin_index_minor(_jacobianA); | ||
| 1003 | //_jacobianALeadindex = bindings::begin_index_major(_jacobianA); | ||
| 1004 | |||
| 1005 | ✗ | } | |
| 1006 | |||
| 1007 | ✗ | int Cvode::reportErrorMessage(ostream& messageStream) | |
| 1008 | { | ||
| 1009 | ✗ | if (_solverStatus == ISolver::SOLVERERROR) | |
| 1010 | { | ||
| 1011 | ✗ | if (_idid == -1) | |
| 1012 | messageStream << "Invalid system dimension." << std::endl; | ||
| 1013 | ✗ | if (_idid == -2) | |
| 1014 | messageStream << "Method not implemented." << std::endl; | ||
| 1015 | ✗ | if (_idid == -3) | |
| 1016 | messageStream << "No valid system/settings available." << std::endl; | ||
| 1017 | ✗ | if (_idid == -11) | |
| 1018 | messageStream << "Step size too small." << std::endl; | ||
| 1019 | } | ||
| 1020 | |||
| 1021 | ✗ | else if (_solverStatus == ISolver::USER_STOP) | |
| 1022 | { | ||
| 1023 | ✗ | messageStream << "Simulation terminated by user at t: " << _tCurrent << std::endl; | |
| 1024 | } | ||
| 1025 | |||
| 1026 | ✗ | return _idid; | |
| 1027 | } | ||
| 1028 | |||
| 1029 | ✗ | void Cvode::writeSimulationInfo() | |
| 1030 | { | ||
| 1031 | // don't write before memory has been initialized | ||
| 1032 | ✗ | if (_cvodeMem == NULL) | |
| 1033 | ✗ | return; | |
| 1034 | |||
| 1035 | long int nst, nfe, nsetups, nni, ncfn, netf; | ||
| 1036 | long int nfQe, netfQ; | ||
| 1037 | long int nfSe, nfeS, nsetupsS, nniS, ncfnS; | ||
| 1038 | long int nfQSe, netfQS; | ||
| 1039 | |||
| 1040 | int qlast, qcur; | ||
| 1041 | sunrealtype h0u, hlast, hcur, tcur; | ||
| 1042 | |||
| 1043 | int flag; | ||
| 1044 | |||
| 1045 | ✗ | flag = CVodeGetIntegratorStats(_cvodeMem, &nst, &nfe, &nsetups, &netf, &qlast, &qcur, &h0u, &hlast, &hcur, &tcur); | |
| 1046 | |||
| 1047 | ✗ | flag = CVodeGetNonlinSolvStats(_cvodeMem, &nni, &ncfn); | |
| 1048 | |||
| 1049 | ✗ | LOGGER_WRITE("Cvode: number steps = " + to_string(nst), LC_SOLVER, LL_INFO); | |
| 1050 | ✗ | LOGGER_WRITE("Cvode: function evaluations 'f' = " + to_string(nfe), LC_SOLVER, LL_INFO); | |
| 1051 | ✗ | LOGGER_WRITE("Cvode: linear solver setups 'nsetups' = " + to_string(nsetups), LC_SOLVER, LL_INFO); | |
| 1052 | ✗ | LOGGER_WRITE("Cvode: nonlinear iterations 'nni' = " + to_string(nni), LC_SOLVER, LL_INFO); | |
| 1053 | ✗ | LOGGER_WRITE("Cvode: convergence failures 'ncfn' = " + to_string(ncfn), LC_SOLVER, LL_INFO); | |
| 1054 | ✗ | LOGGER_WRITE("Cvode: number of evaluateODE calls 'eODE' = " + to_string(_numberOfOdeEvaluations), LC_SOLVER, LL_INFO); | |
| 1055 | |||
| 1056 | //// Solver | ||
| 1057 | //outputStream << "\nSolver: " << getName() | ||
| 1058 | // << "\nVerfahren: "; | ||
| 1059 | //if(_cvodesettings->iMethod == EulerSettings::EULERFORWARD) | ||
| 1060 | // outputStream << " Expliziter Cvode\n\n"; | ||
| 1061 | //else if(_cvodesettings->iMethod == EulerSettings::EULERBACKWARD) | ||
| 1062 | // outputStream << " Impliziter Cvode\n\n"; | ||
| 1063 | |||
| 1064 | //// System | ||
| 1065 | //outputStream | ||
| 1066 | // << "Dimension des Systems (ODE): " << (int)_dimSys << "\n"; | ||
| 1067 | //// Status, Anzahl Schritte, Nullstellenzeugs | ||
| 1068 | //SolverDefaultImplementation::writeSimulationInfo(outputStream); | ||
| 1069 | |||
| 1070 | //// Nullstellensuche | ||
| 1071 | //if (_cvodesettings->iZeroSearchMethod == SolverSettings::NO_ZERO_SEARCH) | ||
| 1072 | //{ | ||
| 1073 | // outputStream << "Nullstellensuche: Keine\n\n" << endl; | ||
| 1074 | //} | ||
| 1075 | //else | ||
| 1076 | //{ | ||
| 1077 | // /*if (_cvodesettings->iZeroSearchMethod == SolverSettings::BISECTION) | ||
| 1078 | // { | ||
| 1079 | // outputStream << "Nullstellensuche: Bisektion\n" << endl; | ||
| 1080 | // } | ||
| 1081 | // else | ||
| 1082 | // {*/ | ||
| 1083 | // outputStream << "Nullstellensuche: Lineare Interpolation\n" << endl; | ||
| 1084 | // /*}*/ | ||
| 1085 | |||
| 1086 | //} | ||
| 1087 | |||
| 1088 | //// Schritteweite | ||
| 1089 | //outputStream | ||
| 1090 | // << "ausgegebene Schritte: " << _outStps << "\n" | ||
| 1091 | // << "Anfangsschrittweite: " << _cvodesettings->dH_init << "\n" | ||
| 1092 | // << "Ausgabeschrittweite: " << dynamic_cast<ISolverSettings*>(_cvodesettings)->getGlobalSettings()->gethOutput() << "\n" | ||
| 1093 | // << "Obere Grenze für Schrittweite: " << _hUpLim << "\n\n"; | ||
| 1094 | //// Status | ||
| 1095 | //outputStream | ||
| 1096 | // << "Solver-Status: " << _idid << "\n\n"; | ||
| 1097 | } | ||
| 1098 | |||
| 1099 | ✗ | int Cvode::check_flag(void *flagvalue, const char *funcname, int opt) | |
| 1100 | { | ||
| 1101 | int *errflag; | ||
| 1102 | |||
| 1103 | /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ | ||
| 1104 | |||
| 1105 | ✗ | if (opt == 0 && flagvalue == NULL) | |
| 1106 | { | ||
| 1107 | ✗ | fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", funcname); | |
| 1108 | ✗ | return (1); | |
| 1109 | } | ||
| 1110 | |||
| 1111 | /* Check if flag < 0 */ | ||
| 1112 | |||
| 1113 | ✗ | else if (opt == 1) | |
| 1114 | { | ||
| 1115 | errflag = (int *)flagvalue; | ||
| 1116 | ✗ | if (*errflag < 0) | |
| 1117 | { | ||
| 1118 | ✗ | fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", funcname, *errflag); | |
| 1119 | ✗ | return (1); | |
| 1120 | } | ||
| 1121 | } | ||
| 1122 | |||
| 1123 | /* Check if function returned NULL pointer - no memory allocated */ | ||
| 1124 | |||
| 1125 | ✗ | else if (opt == 2 && flagvalue == NULL) | |
| 1126 | { | ||
| 1127 | ✗ | fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", funcname); | |
| 1128 | ✗ | return (1); | |
| 1129 | } | ||
| 1130 | |||
| 1131 | return (0); | ||
| 1132 | } | ||
| 1133 | /** @} */ // end of solverCvode | ||
| 1134 |