Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 68.5% 346 / 0 / 505
Functions: 70.0% 14 / 0 / 20
Branches: 30.7% 145 / 0 / 472
OMCompiler/SimulationRuntime/cpp/Solver/IDA/IDA.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 #include <Core/ModelicaDefine.h>
29 #include <Core/Modelica.h>
30 #include <Solver/IDA/IDA.h>
31 #include <Core/Math/Functions.h>
32
33
34 1 Ida::Ida(IMixedSystem* system, ISolverSettings* settings)
35 : SolverDefaultImplementation(system, settings),
36 1 _idasettings(dynamic_cast<ISolverSettings*>(_settings)),
37 1 _idaMem(NULL),
38 1 _sunctx(NULL),
39 /*_z(NULL),
40 _zInit(NULL),
41 _zWrite(NULL),
42 */
43 1 _y(NULL),
44 1 _yp(NULL),
45 1 _yInit(NULL),
46 1 _yWrite(NULL),
47 1 _ypWrite(NULL),
48 1 _dae_res(NULL),
49 1 _dimSys(0),
50 1 _dimAE(0),
51 1 _dimStates(0),
52 1 _cv_rt(0),
53 1 _outStps(0),
54 1 _locStps(0),
55 1 _idid(0),
56 1 _hOut(0.0),
57 1 _tOut(0.0),
58 1 _tZero(0.0),
59 1 _zeroSign(NULL),
60 1 _absTol(NULL),
61 1 _ida_initialized(false),
62 1 _tLastEvent(0.0),
63 1 _event_n(0),
64 1 _properties(NULL),
65 1 _continuous_system(NULL),
66 1 _event_system(NULL),
67 1 _mixed_system(NULL),
68 1 _time_system(NULL),
69 1 _delta(NULL),
70 1 _deltaInv(NULL),
71 1 _ysave(NULL),
72 1 _colorOfColumn (NULL),
73 1 _jacobianAIndex(NULL),
74 1 _jacobianALeadindex(NULL),
75 1 _CV_y0(),
76 1 _CV_y(),
77 1 _CV_yp(),
78 1 _CV_yWrite(),
79 1 _CV_ypWrite(),
80 1 _CV_absTol(),
81 1 _bWritten(false),
82 1 _zeroFound(false),
83 1 _maxColors(0),
84 1 _tLastWrite(-1.0),
85 1 _jacobianANonzeros(0)
86 {
87 1 _data = ((void*) this);
88 #ifdef RUNTIME_PROFILING
89 if(MeasureTime::getInstance() != NULL)
90 {
91 measureTimeFunctionsArray = new std::vector<MeasureTimeData*>(7, NULL); //0 calcFunction //1 solve ... //6 solver statistics
92 (*measureTimeFunctionsArray)[0] = new MeasureTimeData("calcFunction");
93 (*measureTimeFunctionsArray)[1] = new MeasureTimeData("solve");
94 (*measureTimeFunctionsArray)[2] = new MeasureTimeData("writeOutput");
95 (*measureTimeFunctionsArray)[3] = new MeasureTimeData("evaluateZeroFuncs");
96 (*measureTimeFunctionsArray)[4] = new MeasureTimeData("initialize");
97 (*measureTimeFunctionsArray)[5] = new MeasureTimeData("stepCompleted");
98 (*measureTimeFunctionsArray)[6] = new MeasureTimeData("solverStatistics");
99
100 MeasureTime::addResultContentBlock(system->getModelName(),"ida", measureTimeFunctionsArray);
101 measuredFunctionStartValues = MeasureTime::getZeroValues();
102 measuredFunctionEndValues = MeasureTime::getZeroValues();
103 solveFunctionStartValues = MeasureTime::getZeroValues();
104 solveFunctionEndValues = MeasureTime::getZeroValues();
105 solverValues = new MeasureTimeValuesSolver();
106
107 (*measureTimeFunctionsArray)[6]->_sumMeasuredValues = solverValues;
108 }
109 else
110 {
111 measureTimeFunctionsArray = new std::vector<MeasureTimeData*>();
112 measuredFunctionStartValues = NULL;
113 measuredFunctionEndValues = NULL;
114 solveFunctionStartValues = NULL;
115 solveFunctionEndValues = NULL;
116 solverValues = NULL;
117 }
118 #endif
119 1 }
120
121 2 Ida::~Ida()
122 {
123 /*
124 if (_z)
125 delete[] _z;
126 if (_zInit)
127 delete[] _zInit;
128 if (_zWrite)
129 delete[] _zWrite;
130 */
131
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_y)
132 1 delete[] _y;
133
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_yp)
134 1 delete[] _yp;
135
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_yInit)
136 1 delete[] _yInit;
137
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_yWrite)
138 1 delete[] _yWrite;
139
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_ypWrite)
140 1 delete[] _ypWrite;
141
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_dae_res)
142 1 delete[] _dae_res;
143
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_zeroSign)
144 1 delete[] _zeroSign;
145
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_absTol)
146 1 delete[] _absTol;
147
148
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_ida_initialized)
149 {
150 1 N_VDestroy_Serial(_CV_y0);
151 1 N_VDestroy_Serial(_CV_y);
152 1 N_VDestroy_Serial(_CV_yp);
153 1 N_VDestroy_Serial(_CV_yWrite);
154 1 N_VDestroy_Serial(_CV_absTol);
155 1 N_VDestroy_Serial(_ida_ySolver);
156 1 SUNMatDestroy(_ida_J);
157 1 SUNLinSolFree(_ida_linSol);
158 1 IDAFree(&_idaMem);
159 1 SUNContext_Free(&_sunctx);
160 }
161
162
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_colorOfColumn)
163 ✗ delete [] _colorOfColumn;
164
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(_delta)
165 1 delete [] _delta;
166
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(_deltaInv)
167 1 delete [] _deltaInv;
168
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(_ysave)
169 1 delete [] _ysave;
170
171 #ifdef RUNTIME_PROFILING
172 if(measuredFunctionStartValues)
173 delete measuredFunctionStartValues;
174 if(measuredFunctionEndValues)
175 delete measuredFunctionEndValues;
176 if(solveFunctionStartValues)
177 delete solveFunctionStartValues;
178 if(solveFunctionEndValues)
179 delete solveFunctionEndValues;
180 if(solverValues)
181 delete solverValues;
182 #endif
183 2 }
184
185 1 void Ida::initialize()
186 {
187
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _properties = dynamic_cast<ISystemProperties*>(_system);
188
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _continuous_system = dynamic_cast<IContinuous*>(_system);
189
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _event_system = dynamic_cast<IEvent*>(_system);
190 1 _mixed_system = dynamic_cast<IMixedSystem*>(_system);
191
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _time_system = dynamic_cast<ITime*>(_system);
192 1 IGlobalSettings* global_settings = dynamic_cast<ISolverSettings*>(_idasettings)->getGlobalSettings();
193 // Kennzeichnung, dass initialize()() (vor der Integration) aufgerufen wurde
194 1 _idid = 5000;
195 1 _tLastEvent = 0.0;
196 1 _event_n = 0;
197 1 SolverDefaultImplementation::initialize();
198
199 1 _dimStates = _continuous_system->getDimContinuousStates();
200 1 _dimZeroFunc = _event_system->getDimZeroFunc();
201 1 _dimAE = _continuous_system->getDimAE();
202
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(_dimAE>0)
203 1 _dimSys=_dimAE+ _dimStates;
204 else
205 ✗ _dimSys=_dimStates;
206
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_dimStates <= 0)
207
208 {
209 ✗ _idid = -1;
210 ✗ throw std::invalid_argument("Ida::initialize()");
211 }
212 else
213 {
214 // Allocate state vectors, stages and temporary arrays
215
216 /*if (_z)
217 delete[] _z;
218 if (_zInit)
219 delete[] _zInit;
220 if (_zWrite)
221 delete[] _zWrite;*/
222
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_y)
223 ✗ delete[] _y;
224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_yInit)
225 ✗ delete[] _yInit;
226
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_yWrite)
227 ✗ delete[] _yWrite;
228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_ypWrite)
229 ✗ delete[] _ypWrite;
230
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_yp)
231 ✗ delete[] _yp;
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_dae_res)
233 ✗ delete[] _dae_res;
234
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_zeroSign)
235 ✗ delete[] _zeroSign;
236
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_absTol)
237 ✗ delete[] _absTol;
238
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(_delta)
239 ✗ delete [] _delta;
240
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(_deltaInv)
241 ✗ delete [] _deltaInv;
242
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(_ysave)
243 ✗ delete [] _ysave;
244
245
246
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _y = new double[_dimSys];
247
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _yp = new double[_dimSys];
248
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _yInit = new double[_dimSys];
249
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _yWrite = new double[_dimSys];
250
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _ypWrite = new double[_dimSys];
251
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _dae_res = new double[_dimSys];
252 /*
253 _z = new double[_dimSys];
254 _zInit = new double[_dimSys];
255 _zWrite = new double[_dimSys];
256 */
257
258
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _zeroSign = new int[_dimZeroFunc];
259
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _absTol = new double[_dimSys];
260
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _delta =new double[_dimSys];
261
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _deltaInv =new double[_dimSys];
262
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 _ysave =new double[_dimSys];
263
264
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 memset(_y, 0, _dimSys * sizeof(double));
265 1 memset(_yp, 0, _dimSys * sizeof(double));
266 1 memset(_yInit, 0, _dimSys * sizeof(double));
267 1 memset(_ysave, 0, _dimSys * sizeof(double));
268
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 std::fill_n(_absTol, _dimSys, 1.0);
269 // Counter initialisieren
270 1 _outStps = 0;
271
272
1/2
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
 1 if (_idasettings->getDenseOutput())
273 {
274 // Ausgabeschrittweite
275 1 _hOut = global_settings->gethOutput();
276
277 }
278
279 // Allocate memory for the solver
280
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if (SUNContext_Create(SUN_COMM_NULL, &_sunctx) != SUN_SUCCESS)
281 ✗ throw ModelicaSimulationError(SOLVER,"SUNDIALS_ERROR: SUNContext_Create failed");
282 /* Mute SUNDIALS' own logger: package level messages go to stderr/stdout by default, and we report solver failures ourselves. */
283 {
284 1 SUNLogger _sunlogger = NULL;
285
2/4
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
 1 if (SUNContext_GetLogger(_sunctx, &_sunlogger) == SUN_SUCCESS && _sunlogger != NULL) {
286 1 SUNLogger_SetErrorFilename(_sunlogger, "");
287 1 SUNLogger_SetWarningFilename(_sunlogger, "");
288 1 SUNLogger_SetInfoFilename(_sunlogger, "");
289 1 SUNLogger_SetDebugFilename(_sunlogger, "");
290 }
291 }
292
293 1 _idaMem = IDACreate(_sunctx);
294
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if (check_flag((void*) _idaMem, "IDACreate", 0))
295 {
296 ✗ _idid = -5;
297 ✗ throw std::invalid_argument(/*_idid,_tCurrent,*/"Ida::initialize()");
298 }
299
300 //
301 // Make Ida ready for integration
302 //
303
304 // Set initial values for IDA
305 //_continuous_system->evaluateAll(IContinuous::CONTINUOUS);
306 1 _continuous_system->getContinuousStates(_yInit);
307
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 memcpy(_y, _yInit, _dimStates * sizeof(double));
308
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(_dimAE>0)
309 {
310 1 _mixed_system->getAlgebraicDAEVars(_yInit+_dimStates);
311 1 memcpy(_y+_dimStates, _yInit+_dimStates, _dimAE * sizeof(double));
312 1 _continuous_system->getContinuousStates(_yp);
313 }
314 // Get nominal values
315 1 _continuous_system->getNominalStates(_absTol);
316
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 1 time.
 10 for (int i = 0; i < _dimStates; i++)
317 9 _absTol[i] = dynamic_cast<ISolverSettings*>(_idasettings)->getATol();
318
319 1 _CV_y0 = N_VMake_Serial(_dimSys, _yInit, _sunctx);
320 1 _CV_y = N_VMake_Serial(_dimSys, _y, _sunctx);
321 1 _CV_yp = N_VMake_Serial(_dimSys, _yp, _sunctx);
322 1 _CV_yWrite = N_VMake_Serial(_dimSys, _yWrite, _sunctx);
323 1 _CV_ypWrite = N_VMake_Serial(_dimSys, _ypWrite, _sunctx);
324 1 _CV_absTol = N_VMake_Serial(_dimSys, _absTol, _sunctx);
325
326
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if (check_flag((void*) _CV_y0, "N_VMake_Serial", 0))
327 {
328 ✗ _idid = -5;
329 ✗ throw std::invalid_argument("Ida::initialize()");
330 }
331
332 //is already initialized: calcFunction(_tCurrent, NV_DATA_S(_CV_y0), NV_DATA_S(_CV_yp),NV_DATA_S(_CV_yp));
333
334 // Initialize Ida (Initial values are required)
335 1 _idid = IDAInit(_idaMem, rhsFunctionCB, _tCurrent, _CV_y0, _CV_yp);
336
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
337 {
338 ✗ _idid = -5;
339 ✗ throw std::invalid_argument("Ida::initialize()");
340 }
341 1 _idid = SUNContext_PushErrHandler(_sunctx, errOutputIDA, _data);
342
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
343 ✗ throw std::invalid_argument("IDA::initialize()");
344 // Set Tolerances
345 1 _idid = IDASVtolerances(_idaMem, dynamic_cast<ISolverSettings*>(_idasettings)->getRTol(), _CV_absTol); // RTOL and ATOL
346
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
347 ✗ throw std::invalid_argument("IDA::initialize()");
348
349 // Set the pointer to user-defined data
350 1 _idid = IDASetUserData(_idaMem, _data);
351
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
352 ✗ throw std::invalid_argument("IDA::initialize()");
353
354 1 _idid = IDASetInitStep(_idaMem, 1e-6); // INITIAL STEPSIZE
355
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
356 ✗ throw std::invalid_argument("Ida::initialize()");
357
358
359 1 _idid = IDASetMaxStep(_idaMem, global_settings->getEndTime() / 10.0); // MAXIMUM STEPSIZE
360
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
361 ✗ throw std::invalid_argument("IDA::initialize()");
362
363 1 _idid = IDASetMaxNonlinIters(_idaMem, 5); // Max number of iterations
364
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
365 ✗ throw std::invalid_argument("IDA::initialize()");
366 1 _idid = IDASetMaxErrTestFails(_idaMem, 100);
367
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
368 ✗ throw std::invalid_argument("IDA::initialize()");
369
370 1 _idid = IDASetMaxNumSteps(_idaMem, 1e3); // Max Number of steps
371
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
372 ✗ throw std::invalid_argument(/*_idid,_tCurrent,*/"IDA::initialize()");
373
374 // Initialize dense linear solver
375 1 _ida_ySolver = N_VNew_Serial(_dimSys, _sunctx);
376 1 _ida_J = SUNDenseMatrix(_dimSys, _dimSys, _sunctx);
377 1 _ida_linSol = SUNLinSol_Dense(_ida_ySolver, _ida_J, _sunctx);
378 1 _idid = IDASetLinearSolver(_idaMem, _ida_linSol, _ida_J);
379
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
380 ✗ throw std::invalid_argument("IDA::initialize()");
381
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(_dimAE>0)
382 {
383 1 _idid = IDASetSuppressAlg(_idaMem, SUNTRUE);
384
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 double* tmp = new double[_dimSys];
385
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 std::fill_n(tmp, _dimStates, 1.0);
386
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 std::fill_n(tmp+_dimStates, _dimAE, 0.0);
387 1 _idid = IDASetId(_idaMem, N_VMake_Serial(_dimSys, tmp, _sunctx));
388 1 delete [] tmp;
389
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
390 ✗ throw std::invalid_argument("IDA::initialize()");
391 }
392
393 // Use own jacobian matrix
394 //_idid = CVodeSetJacFn(_idaMem, &jacobianFunctionCB);
395 //if (_idid < 0)
396 // throw std::invalid_argument("IDA::initialize()");
397
398
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_dimZeroFunc)
399 {
400 1 _idid = IDARootInit(_idaMem, _dimZeroFunc, &zeroFunctionCB);
401
402 1 memset(_zeroSign, 0, _dimZeroFunc * sizeof(int));
403 1 _idid = IDASetRootDirection(_idaMem, _zeroSign);
404
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
405 ✗ throw std::invalid_argument(/*_idid,_tCurrent,*/"IDA::initialize()");
406 1 memset(_zeroSign, -1, _dimZeroFunc * sizeof(int));
407 1 memset(_zeroVal, -1, _dimZeroFunc * sizeof(int));
408
409 }
410
411
412 1 _ida_initialized = true;
413
414 //
415 // IDA is ready for integration
416 //
417 // BOOST_LOG_SEV(ida_lg::get(), ida_info) << "IDA initialized";
418 }
419 1 }
420
421 4 void Ida::solve(const SOLVERCALL action)
422 {
423 4 bool writeEventOutput = (_settings->getGlobalSettings()->getOutputPointType() == OPT_ALL);
424 4 bool writeOutput = !(_settings->getGlobalSettings()->getOutputPointType() == OPT_NONE);
425
426 #ifdef RUNTIME_PROFILING
427 MEASURETIME_REGION_DEFINE(idaSolveFunctionHandler, "solve");
428 if(MeasureTime::getInstance() != NULL)
429 {
430 MEASURETIME_START(solveFunctionStartValues, idaSolveFunctionHandler, "solve");
431 }
432 #endif
433
434
2/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 if (_idasettings && _system)
435 {
436 // Solver und System für Integration vorbereiten
437
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
 4 if ((action & RECORDCALL) && (action & FIRST_CALL))
438 {
439 #ifdef RUNTIME_PROFILING
440 MEASURETIME_REGION_DEFINE(idaInitializeHandler, "IDAInitialize");
441 if(MeasureTime::getInstance() != NULL)
442 {
443 MEASURETIME_START(measuredFunctionStartValues, idaInitializeHandler, "IDAInitialize");
444 }
445 #endif
446
447 1 initialize();
448
449 #ifdef RUNTIME_PROFILING
450 if(MeasureTime::getInstance() != NULL)
451 {
452 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[4], idaInitializeHandler);
453 }
454 #endif
455
456
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (writeOutput)
457 1 writeToFile(0, _tCurrent, _h);
458 1 _tLastWrite = 0;
459
460 1 return;
461 }
462
463
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
 3 if ((action & RECORDCALL) && !(action & FIRST_CALL))
464 {
465 2 writeToFile(_accStps, _tCurrent, _h);
466 2 return;
467 }
468
469 // Nach einem TimeEvent wird der neue Zustand recorded
470
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (action & RECALL)
471 {
472 ✗ _firstStep = true;
473 ✗ if (writeEventOutput)
474 ✗ writeToFile(0, _tCurrent, _h);
475 ✗ if (writeOutput)
476 ✗ writeIDAOutput(_tCurrent, _h, _locStps);
477 ✗ _continuous_system->getContinuousStates(_y);
478 }
479
480 // Solver soll fortfahren
481 1 _solverStatus = ISolver::CONTINUE;
482
483
3/4
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
 3 while ((_solverStatus & ISolver::CONTINUE) && !_interrupt )
484 {
485 // Zuvor wurde initialize aufgerufen und hat funktioniert => RESET IDID
486
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid == 5000)
487 ✗ _idid = 0;
488
489 // Solveraufruf
490
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (_idid == 0)
491 {
492 // Zähler zurücksetzen
493 1 _accStps = 0;
494 1 _locStps = 0;
495
496 // Solverstart
497 1 IDACore();
498
499 }
500
501 // Integration war nicht erfolgreich und wurde auch nicht vom User unterbrochen
502
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid != 0 && _idid != 1)
503 {
504 ✗ _solverStatus = ISolver::SOLVERERROR;
505 //throw std::invalid_argument(_idid,_tCurrent,"IDA::solve()");
506 ✗ throw std::invalid_argument("IDA::solve()");
507 }
508
509 // Abbruchkriterium (erreichen der Endzeit)
510
1/2
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
 1 else if ((_tEnd - _tCurrent) <= dynamic_cast<ISolverSettings*>(_idasettings)->getEndTimeTol())
511 1 _solverStatus = DONE;
512 }
513
514 1 _firstCall = false;
515
516 }
517 else
518 {
519
520 ✗ throw std::invalid_argument("IDA::solve()");
521 }
522
523 #ifdef RUNTIME_PROFILING
524 if(MeasureTime::getInstance() != NULL)
525 {
526 MEASURETIME_END(solveFunctionStartValues, solveFunctionEndValues, (*measureTimeFunctionsArray)[1], idaSolveFunctionHandler);
527
528 long int nst, nfe, nsetups, netf, nni, ncfn;
529 int qlast, qcur;
530 sunrealtype h0u, hlast, hcur, tcur;
531
532 int flag;
533
534 flag = IDAGetIntegratorStats(_idaMem, &nst, &nfe, &nsetups, &netf, &qlast, &qcur, &h0u, &hlast, &hcur, &tcur);
535 flag = IDAGetNonlinSolvStats(_idaMem, &nni, &ncfn);
536
537 MeasureTimeValuesSolver solverVals = MeasureTimeValuesSolver(nfe, netf);
538 (*measureTimeFunctionsArray)[6]->_sumMeasuredValues->_numCalcs += nst;
539 (*measureTimeFunctionsArray)[6]->_sumMeasuredValues->add(&solverVals);
540 }
541 #endif
542 }
543 4 bool Ida::isInterrupted()
544 {
545
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if(_interrupt)
546 {
547 ✗ _solverStatus = DONE;
548 ✗ return true;
549 }
550 else
551 {
552 return false;
553 }
554 }
555 1 void Ida::IDACore()
556 {
557 1 _idid = IDAReInit(_idaMem, _tCurrent, _CV_y,_CV_yp);
558 1 _idid = IDASetStopTime(_idaMem, _tEnd);
559 1 _idid = IDASetInitStep(_idaMem, 1e-12);
560
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (_idid < 0)
561 ✗ throw std::runtime_error("IDA::ReInit");
562
563 1 bool writeEventOutput = (_settings->getGlobalSettings()->getOutputPointType() == OPT_ALL);
564 1 bool writeOutput = !(_settings->getGlobalSettings()->getOutputPointType() == OPT_NONE);
565
566
3/4
✓ Branch 0 taken 13348 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 13348 times.
✗ Branch 3 not taken.
 13350 while ((_solverStatus & ISolver::CONTINUE) && !_interrupt )
567 {
568 13348 _cv_rt = IDASolve(_idaMem, _tEnd, &_tCurrent, _CV_y, _CV_yp, IDA_ONE_STEP);
569
570 13348 _idid = IDAGetNumSteps(_idaMem, &_locStps);
571
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13348 times.
 13348 if (_idid != IDA_SUCCESS)
572 ✗ throw std::runtime_error("IDAGetNumSteps failed. The ida mem pointer is NULL");
573
574 13348 _idid =IDAGetLastStep(_idaMem, &_h);
575
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13348 times.
 13348 if (_idid != IDA_SUCCESS)
576 ✗ throw std::runtime_error("IDAGetLastStep failed. The ida mem pointer is NULL");
577
578 //Check if there was at least one output-point within the last solver interval
579 // -> Write output if true
580
1/2
✓ Branch 0 taken 13348 times.
✗ Branch 1 not taken.
 13348 if (writeOutput)
581 {
582 13348 writeIDAOutput(_tCurrent, _h, _locStps);
583 }
584
585 #ifdef RUNTIME_PROFILING
586 MEASURETIME_REGION_DEFINE(idaStepCompletedHandler, "IDAStepCompleted");
587 if(MeasureTime::getInstance() != NULL)
588 {
589 MEASURETIME_START(measuredFunctionStartValues, idaStepCompletedHandler, "IDAStepCompleted");
590 }
591 #endif
592
593 //set completed step to system and check if terminate was called
594
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 13348 times.
 13348 if(_continuous_system->stepCompleted(_tCurrent))
595 ✗ _solverStatus = DONE;
596
597 #ifdef RUNTIME_PROFILING
598 if(MeasureTime::getInstance() != NULL)
599 {
600 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[5], idaStepCompletedHandler);
601 }
602 #endif
603
604 // Perform state selection
605 13348 bool state_selection = stateSelection();
606
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13348 times.
 13348 if (state_selection)
607 ✗ _continuous_system->getContinuousStates(_y);
608
609 13348 _zeroFound = false;
610
611 // Check if step was successful
612
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 13348 times.
 13348 if (check_flag(&_cv_rt, "IDA", 1))
613 {
614 ✗ _solverStatus = ISolver::SOLVERERROR;
615 ✗ break;
616 }
617
618 // A root was found
619
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 13346 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 13348 if ((_cv_rt == IDA_ROOT_RETURN) && !isInterrupted())
620 {
621 // IDA is setting _tCurrent to the time where the first event occurred
622 2 double _abs = fabs(_tLastEvent - _tCurrent);
623 2 _zeroFound = true;
624
625
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if ((_abs < 1e-3) && _event_n == 0)
626 {
627 ✗ _tLastEvent = _tCurrent;
628 ✗ _event_n++;
629 }
630
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 else if ((_abs < 1e-3) && (_event_n >= 1 && _event_n < 500))
631 {
632 ✗ _event_n++;
633 }
634
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 else if ((_abs >= 1e-3))
635 {
636 //restart event counter
637 2 _tLastEvent = _tCurrent;
638 2 _event_n = 0;
639 }
640 else
641 ✗ throw std::runtime_error("Number of events exceeded in time interval " + to_string(_abs) + " at time " + to_string(_tCurrent));
642
643 // IDA has interpolated the states at time 'tCurrent'
644 2 _time_system->setTime(_tCurrent);
645
646 // To get steep steps in the result file, two value points (P1 and P2) must be added
647 //
648 // Y | (P2) X...........
649 // | :
650 // | :
651 // |........X (P1)
652 // |---------------------------------->
653 // | ^ t
654 // _tCurrent
655
656 // Write the values of (P1)
657
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (writeEventOutput)
658 {
659
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(_dimAE>0)
660 {
661 2 _continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
662 }
663 else
664 {
665 ✗ _continuous_system->evaluateAll(IContinuous::CONTINUOUS);
666 }
667 2 writeToFile(0, _tCurrent, _h);
668 }
669
670 2 _idid = IDAGetRootInfo(_idaMem, _zeroSign);
671
672
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 2 times.
 20 for (int i = 0; i < _dimZeroFunc; i++)
673 18 _events[i] = bool(_zeroSign[i]);
674
675
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if (_mixed_system->handleSystemEvents(_events))
676 {
677 // State variables were reinitialized, thus we have to give these values to the ida-solver
678 // Take care about the memory regions, _z is the same like _CV_y
679 ✗ _continuous_system->getContinuousStates(_y);
680 ✗ if(_dimAE>0)
681 {
682 ✗ _mixed_system->getAlgebraicDAEVars(_y+_dimStates);
683 ✗ _continuous_system->getRHS(_yp);
684 }
685 ✗ calcFunction(_tCurrent, NV_DATA_S(_CV_y), NV_DATA_S(_CV_yp),_dae_res);
686
687 }
688 }
689
690
4/6
✓ Branch 0 taken 13346 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13346 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 13348 if ((_zeroFound || state_selection)&& !isInterrupted())
691 {
692 // Write the values of (P2)
693
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (writeEventOutput)
694 {
695 // If we want to write the event-results, we should evaluate the whole system again
696
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(_dimAE>0)
697 {
698 2 _continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
699 }
700 else
701 {
702 ✗ _continuous_system->evaluateAll(IContinuous::CONTINUOUS);
703 }
704 2 writeToFile(0, _tCurrent, _h);
705 }
706
707 2 _idid = IDAReInit(_idaMem, _tCurrent, _CV_y,_CV_yp);
708
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (_idid < 0)
709 ✗ throw std::runtime_error("IDA::ReInit()");
710
711 // Der Eventzeitpunkt kann auf der Endzeit liegen (Time-Events). In diesem Fall wird der Solver beendet, da IDA sonst eine interne Warnung schmeißt
712
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (_tCurrent == _tEnd)
713 ✗ _cv_rt = IDA_TSTOP_RETURN;
714 }
715
716 // Zähler für die Anzahl der ausgegebenen Schritte erhöhen
717 13348 ++_outStps;
718 13348 _tLastSuccess = _tCurrent;
719
720
2/2
✓ Branch 0 taken 13347 times.
✓ Branch 1 taken 1 time.
 13348 if (_cv_rt == IDA_TSTOP_RETURN)
721 {
722 1 _time_system->setTime(_tEnd);
723 1 _continuous_system->setContinuousStates(NV_DATA_S(_CV_y));
724
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(_dimAE>0)
725 {
726 1 _mixed_system->setAlgebraicDAEVars(NV_DATA_S(_CV_y)+_dimStates);
727 1 _continuous_system->setStateDerivatives(NV_DATA_S(_CV_yp));
728 1 _continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
729 }
730 else
731 {
732 ✗ _continuous_system->evaluateAll(IContinuous::CONTINUOUS);
733 }
734
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(writeOutput)
735 1 writeToFile(0, _tEnd, _h);
736
737 1 _accStps += _locStps;
738 1 _solverStatus = DONE;
739 }
740 }
741 1 }
742 ✗ void Ida::setTimeOut(unsigned int time_out)
743 {
744 ✗ SimulationMonitor::setTimeOut(time_out);
745 ✗ }
746 ✗ void Ida::stop()
747 {
748 ✗ SimulationMonitor::stop();
749 ✗ }
750 13348 void Ida::writeIDAOutput(const double &time, const double &h, const int &stp)
751 {
752 #ifdef RUNTIME_PROFILING
753 MEASURETIME_REGION_DEFINE(idaWriteOutputHandler, "IDAWriteOutput");
754 if(MeasureTime::getInstance() != NULL)
755 {
756 MEASURETIME_START(measuredFunctionStartValues, idaWriteOutputHandler, "IDAWriteOutput");
757 }
758 #endif
759
760
1/2
✓ Branch 0 taken 13348 times.
✗ Branch 1 not taken.
 13348 if (stp > 0)
761 {
762
1/2
✓ Branch 1 taken 13348 times.
✗ Branch 2 not taken.
 13348 if (_idasettings->getDenseOutput())
763 {
764 13348 _bWritten = false;
765 double *oldValues = NULL;
766
767 //We have to find all output-points within the last solver step
768
2/2
✓ Branch 2 taken 24999 times.
✓ Branch 3 taken 13348 times.
 38347 while (_tLastWrite + dynamic_cast<ISolverSettings*>(_idasettings)->getGlobalSettings()->gethOutput() <= time)
769 {
770
2/2
✓ Branch 0 taken 11734 times.
✓ Branch 1 taken 13265 times.
 24999 if (!_bWritten)
771 {
772 //Rescue the calculated derivatives
773
1/2
✓ Branch 1 taken 11734 times.
✗ Branch 2 not taken.
 11734 oldValues = new double[_continuous_system->getDimRHS()];
774 11734 _continuous_system->getRHS(oldValues);
775 }
776 24999 _bWritten = true;
777 24999 _tLastWrite = _tLastWrite + dynamic_cast<ISolverSettings*>(_idasettings)->getGlobalSettings()->gethOutput();
778 //Get the state vars at the output-point (interpolated)
779 24999 _idid = IDAGetDky(_idaMem, _tLastWrite, 0, _CV_yWrite);
780 24999 _time_system->setTime(_tLastWrite);
781 24999 _continuous_system->setContinuousStates(NV_DATA_S(_CV_yWrite));
782
1/2
✓ Branch 0 taken 24999 times.
✗ Branch 1 not taken.
 24999 if(_dimAE>0)
783 {
784 24999 _mixed_system->setAlgebraicDAEVars(NV_DATA_S(_CV_y)+_dimStates);
785 24999 _idid = IDAGetDky(_idaMem, _tLastWrite, 1, _CV_ypWrite);
786 24999 _continuous_system->setStateDerivatives(NV_DATA_S(_CV_ypWrite));
787 24999 _continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
788 }
789 else
790 {
791 ✗ _continuous_system->evaluateAll(IContinuous::CONTINUOUS);
792 }
793 #ifdef RUNTIME_PROFILING
794 if(MeasureTime::getInstance() != NULL)
795 {
796 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[2], idaWriteOutputHandler);
797 }
798 #endif
799 24999 SolverDefaultImplementation::writeToFile(stp, _tLastWrite, h);
800 #ifdef RUNTIME_PROFILING
801 MEASURETIME_REGION_DEFINE(idaWriteOutputHandler, "IDAWriteOutput");
802 if(MeasureTime::getInstance() != NULL)
803 {
804 (*measureTimeFunctionsArray)[2]->_sumMeasuredValues->_numCalcs--;
805 MEASURETIME_START(measuredFunctionStartValues, idaWriteOutputHandler, "IDAWriteOutput");
806 }
807 #endif
808 } //end if time -_tLastWritten
809
2/2
✓ Branch 0 taken 11734 times.
✓ Branch 1 taken 1614 times.
 13348 if (_bWritten)
810 {
811 11734 _time_system->setTime(time);
812 11734 _continuous_system->setContinuousStates(_y);
813 11734 _continuous_system->setStateDerivatives(oldValues);
814
1/2
✓ Branch 0 taken 11734 times.
✗ Branch 1 not taken.
 11734 if(_dimAE>0)
815 {
816 11734 _mixed_system->setAlgebraicDAEVars(_y+_dimStates);
817 }
818
1/2
✓ Branch 0 taken 11734 times.
✗ Branch 1 not taken.
 11734 delete[] oldValues;
819
820 }
821
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1614 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1614 else if (time == _tEnd && _tLastWrite != time)
822 {
823 ✗ _idid = IDAGetDky(_idaMem, time, 0, _CV_y);
824 ✗ _idid = IDAGetDky(_idaMem, time, 1, _CV_yp);
825 ✗ _time_system->setTime(time);
826 ✗ _continuous_system->setContinuousStates(NV_DATA_S(_CV_y));
827 ✗ if(_dimAE>0)
828 {
829 ✗ _mixed_system->setAlgebraicDAEVars(NV_DATA_S(_CV_y)+_dimStates);
830 ✗ _continuous_system->setStateDerivatives(NV_DATA_S(_CV_yp));
831 ✗ _continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
832 }
833 else
834 {
835 ✗ _continuous_system->evaluateAll(IContinuous::CONTINUOUS);
836 }
837 #ifdef RUNTIME_PROFILING
838 if(MeasureTime::getInstance() != NULL)
839 {
840 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[2], idaWriteOutputHandler);
841 }
842 #endif
843 ✗ SolverDefaultImplementation::writeToFile(stp, _tEnd, h);
844 }
845 }
846 else
847 {
848 #ifdef RUNTIME_PROFILING
849 if(MeasureTime::getInstance() != NULL)
850 {
851 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[2], idaWriteOutputHandler);
852 }
853 #endif
854 ✗ SolverDefaultImplementation::writeToFile(stp, time, h);
855 }
856 }
857 13348 }
858
859 13350 bool Ida::stateSelection()
860 {
861 13350 return SolverDefaultImplementation::stateSelection();
862 }
863 167637 int Ida::calcFunction(const double& time, const double* y, double *yp,double* res)
864 {
865 #ifdef RUNTIME_PROFILING
866 MEASURETIME_REGION_DEFINE(idaCalcFunctionHandler, "IDACalcFunction");
867 if(MeasureTime::getInstance() != NULL)
868 {
869 MEASURETIME_START(measuredFunctionStartValues, idaCalcFunctionHandler, "IDACalcFunction");
870 }
871 #endif
872
873 int returnValue = 0;
874 try
875 {
876
1/2
✓ Branch 0 taken 167637 times.
✗ Branch 1 not taken.
 167637 if(_dimAE>0)
877 {
878
1/1
✓ Branch 1 taken 167637 times.
 167637 _time_system->setTime(time);
879
1/1
✓ Branch 1 taken 167637 times.
 167637 _continuous_system->setContinuousStates(y);
880
1/1
✓ Branch 1 taken 167637 times.
 167637 _continuous_system->setStateDerivatives(yp);
881
1/1
✓ Branch 1 taken 167637 times.
 167637 _mixed_system->setAlgebraicDAEVars(y+_dimStates);
882
1/1
✓ Branch 1 taken 167637 times.
 167637 _continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
883
1/1
✓ Branch 1 taken 167637 times.
 167637 _mixed_system->getResidual(res);
884
885 }
886 else
887 {
888 ✗ _time_system->setTime(time);
889 ✗ _continuous_system->setContinuousStates(y);
890 ✗ _continuous_system->evaluateODE(IContinuous::CONTINUOUS);
891 ✗ _continuous_system->getRHS(res);
892 ✗ for(size_t i(0); i<_dimSys; ++i)
893 ✗ res[i]-=yp[i];
894
895 }
896 } //workaround until exception can be catch from c- libraries
897 ✗ catch (std::exception& ex)
898 {
899 ✗ std::string error = ex.what();
900 cerr << "IDA integration error: " << error;
901 returnValue = -1;
902 ✗ }
903
904 #ifdef RUNTIME_PROFILING
905 if(MeasureTime::getInstance() != NULL)
906 {
907 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[0], idaCalcFunctionHandler);
908 }
909 #endif
910
911 167637 return returnValue;
912 }
913
914 167637 int Ida::rhsFunctionCB(double t, N_Vector y, N_Vector ydot, N_Vector resval, void *user_data)
915 {
916
917 167637 int status = ((Ida*) user_data)->calcFunction(t, NV_DATA_S(y), NV_DATA_S(ydot),NV_DATA_S(resval));
918
919 167637 return status;
920 }
921
922 13353 void Ida::giveZeroVal(const double &t, const double *y,const double *yp,double *zeroValue)
923 {
924 #ifdef RUNTIME_PROFILING
925 MEASURETIME_REGION_DEFINE(idaEvalZeroHandler, "evaluateZeroFuncs");
926 if(MeasureTime::getInstance() != NULL)
927 {
928 MEASURETIME_START(measuredFunctionStartValues, idaEvalZeroHandler, "evaluateZeroFuncs");
929 }
930 #endif
931
932 13353 _time_system->setTime(t);
933 13353 _continuous_system->setContinuousStates(y);
934
1/2
✓ Branch 0 taken 13353 times.
✗ Branch 1 not taken.
 13353 if(_dimAE>0)
935 {
936 13353 _mixed_system->setAlgebraicDAEVars(y+_dimStates);
937 13353 _continuous_system->setStateDerivatives(yp);
938 }
939 // System aktualisieren
940 13353 _continuous_system->evaluateZeroFuncs(IContinuous::DISCRETE);
941
942 13353 _event_system->getZeroFunc(zeroValue);
943
944 #ifdef RUNTIME_PROFILING
945 if(MeasureTime::getInstance() != NULL)
946 {
947 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[3], idaEvalZeroHandler);
948 }
949 #endif
950 13353 }
951
952 13353 int Ida::zeroFunctionCB(double t, N_Vector y, N_Vector yp, double *zeroval, void *user_data)
953 {
954 13353 ((Ida*) user_data)->giveZeroVal(t, NV_DATA_S(y),NV_DATA_S(yp), zeroval);
955
956 13353 return (0);
957 }
958
959 ✗ int Ida::jacobianFunctionCB(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)
960 {
961 ✗ return ((Ida*) user_data)->calcJacobian(t,N, tmp1, tmp2, tmp3, NV_DATA_S(y), fy, Jac);
962
963 }
964
965
966 ✗ int Ida::calcJacobian(double t, long int N, N_Vector fHelp, N_Vector errorWeight, N_Vector jthCol, double* y, N_Vector fy, SUNMatrix Jac)
967 {
968 try
969 {
970 int l,g;
971 double fnorm, minInc, *f_data, *fHelp_data, *errorWeight_data, h, srur, delta_inv;
972
973 ✗ f_data = NV_DATA_S(fy);
974 ✗ errorWeight_data = NV_DATA_S(errorWeight);
975 ✗ fHelp_data = NV_DATA_S(fHelp);
976
977
978 //Get relevant info
979 ✗ _idid = IDAGetErrWeights(_idaMem, errorWeight);
980 ✗ if (_idid < 0)
981 {
982 ✗ _idid = -5;
983 ✗ throw std::invalid_argument("IDA::calcJacobian()");
984 }
985 ✗ _idid = IDAGetCurrentStep(_idaMem, &h);
986 ✗ if (_idid < 0)
987 {
988 ✗ _idid = -5;
989 ✗ throw std::invalid_argument("IDA::calcJacobian()");
990 }
991
992 srur = sqrt(UROUND);
993
994 ✗ fnorm = N_VWrmsNorm(fy, errorWeight);
995 ✗ minInc = (fnorm != 0.0) ?
996 ✗ (1000.0 * abs(h) * UROUND * N * fnorm) : 1.0;
997
998 ✗ for(int j=0;j<N;j++)
999 {
1000 ✗ _delta[j] = max(srur*abs(y[j]), minInc/errorWeight_data[j]);
1001 }
1002 ✗ for(int j=0;j<N;j++)
1003 {
1004 ✗ _deltaInv[j] = 1/_delta[j];
1005 }
1006
1007 // Calculation of the jacobian
1008
1009 ✗ if (_jacobianANonzeros != 0)
1010 {
1011 ✗ for(int color=1; color <= _maxColors; color++)
1012 {
1013 ✗ for(int k=0; k < _dimSys; k++)
1014 {
1015 ✗ if((_colorOfColumn[k] ) == color)
1016 {
1017 ✗ _ysave[k] = y[k];
1018 ✗ y[k]+= _delta[k];
1019 }
1020 }
1021
1022 ✗ calcFunction(t, y, fHelp_data,fHelp_data);
1023
1024 ✗ for (int k = 0; k < _dimSys; k++)
1025 {
1026 ✗ if((_colorOfColumn[k]) == color)
1027 {
1028 ✗ y[k] = _ysave[k];
1029
1030 ✗ int startOfColumn = k * _dimSys;
1031 ✗ for (int j = _jacobianALeadindex[k]; j < _jacobianALeadindex[k+1];j++)
1032 {
1033 ✗ l = _jacobianAIndex[j];
1034 ✗ g = l + startOfColumn;
1035 ✗ SM_DATA_D(Jac)[g] = (fHelp_data[l] - f_data[l]) * _deltaInv[k];
1036 }
1037 }
1038 }
1039 }
1040 }
1041
1042 /*
1043 //Calculation of J without colouring
1044 for (j = 0; j < N; j++)
1045 {
1046
1047
1048 //N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol);
1049
1050 _ysave[j] = y[j];
1051
1052 y[j] += _delta[j];
1053
1054 calcFunction(t, y, fHelp_data);
1055
1056 y[j] = _ysave[j];
1057
1058 delta_inv = 1.0/_delta[j];
1059 N_VLinearSum(delta_inv, fHelp, -delta_inv, fy, jthCol);
1060
1061 for(int i=0; i<_dimSys; ++i)
1062 {
1063 SM_DATA_D(Jac)[i+j*_dimSys] = NV_Ith_S(jthCol,i);
1064 }
1065
1066 //DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol);
1067 }
1068 */
1069
1070 } //workaround until exception can be catch from c- libraries
1071 ✗ catch (std::exception& ex)
1072 {
1073 ✗ std::string error = ex.what();
1074 cerr << "IDA integration error: " << error;
1075 return 1;
1076 ✗ }
1077
1078
1079 ✗ return 0;
1080 }
1081
1082
1083
1084 ✗ int Ida::reportErrorMessage(ostream& messageStream)
1085 {
1086 ✗ if (_solverStatus == ISolver::SOLVERERROR)
1087 {
1088 ✗ if (_idid == -1)
1089 messageStream << "Invalid system dimension." << std::endl;
1090 ✗ if (_idid == -2)
1091 messageStream << "Method not implemented." << std::endl;
1092 ✗ if (_idid == -3)
1093 messageStream << "No valid system/settings available." << std::endl;
1094 ✗ if (_idid == -11)
1095 messageStream << "Step size too small." << std::endl;
1096 }
1097
1098 ✗ else if (_solverStatus == ISolver::USER_STOP)
1099 {
1100 ✗ messageStream << "Simulation terminated by user at t: " << _tCurrent << std::endl;
1101 }
1102
1103 ✗ return _idid;
1104 }
1105
1106 1 void Ida::writeSimulationInfo()
1107 {
1108 long int nst, nfe, nsetups, nni, ncfn, netf;
1109 long int nfQe, netfQ;
1110 long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;
1111 long int nfQSe, netfQS;
1112
1113 int qlast, qcur;
1114 sunrealtype h0u, hlast, hcur, tcur;
1115
1116 int flag;
1117
1118 1 flag = IDAGetIntegratorStats(_idaMem, &nst, &nfe, &nsetups, &netf, &qlast, &qcur, &h0u, &hlast, &hcur, &tcur);
1119
1120 1 flag = IDAGetNonlinSolvStats(_idaMem, &nni, &ncfn);
1121
1122 ✗ LOGGER_WRITE("IDA: number steps = " + to_string(nst), LC_SOLVER, LL_INFO);
1123 ✗ LOGGER_WRITE("IDA: function evaluations 'f' = " + to_string(nfe), LC_SOLVER, LL_INFO);
1124 ✗ LOGGER_WRITE("IDA: error test failures 'netf' = " + to_string(netfS), LC_SOLVER, LL_INFO);
1125 ✗ LOGGER_WRITE("IDA: linear solver setups 'nsetups' = " + to_string(nsetups), LC_SOLVER, LL_INFO);
1126 ✗ LOGGER_WRITE("IDA: nonlinear iterations 'nni' = " + to_string(nni), LC_SOLVER, LL_INFO);
1127 ✗ LOGGER_WRITE("IDA: convergence failures 'ncfn' = " + to_string(ncfn), LC_SOLVER, LL_INFO);
1128 1 }
1129
1130 13350 int Ida::check_flag(void *flagvalue, const char *funcname, int opt)
1131 {
1132 int *errflag;
1133
1134 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
1135
1136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13350 times.
 13350 if (opt == 0 && flagvalue == NULL)
1137 {
1138 ✗ fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", funcname);
1139 ✗ return (1);
1140 }
1141
1142 /* Check if flag < 0 */
1143
1144
2/2
✓ Branch 0 taken 13348 times.
✓ Branch 1 taken 2 times.
 13350 else if (opt == 1)
1145 {
1146 errflag = (int *) flagvalue;
1147
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13348 times.
 13348 if (*errflag < 0)
1148 {
1149 ✗ fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", funcname, *errflag);
1150 ✗ return (1);
1151 }
1152 }
1153
1154 /* Check if function returned NULL pointer - no memory allocated */
1155
1156
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 else if (opt == 2 && flagvalue == NULL)
1157 {
1158 ✗ fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", funcname);
1159 ✗ return (1);
1160 }
1161
1162 return (0);
1163 }
1164
1165
1166 ✗ void Ida::errOutputIDA(int line, const char *func, const char *file, const char *msg,
1167 SUNErrCode err_code, void *err_user_data, SUNContext sunctx)
1168 {
1169 ✗ cout << "#### IDA error message #####";
1170 cout << " -> error code " << err_code
1171 ✗ << " in function " << func << " at " << file << ":" << line;
1172 /* Package level codes (IDA_* and friends) are not SUNErrCodes, so SUNGetErrMsg()
1173 only makes sense when SUNDIALS did not supply a message. */
1174 ✗ cout << " Message: " << (msg ? msg : SUNGetErrMsg(err_code));
1175 ✗ }
1176