Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 81.3% 279 / 0 / 343
Functions: 87.5% 14 / 0 / 16
Branches: 39.3% 216 / 0 / 549
OMCompiler/SimulationRuntime/cpp/Solver/DASSL/DASSL.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 solverDASSL
29 *
30 * @{
31 */
32
33
34
35 #include <Core/ModelicaDefine.h>
36 #include <Core/Modelica.h>
37 #include <Solver/DASSL/DASSL.h>
38
39 #include <Core/Utils/extension/logger.hpp>
40 #include <Core/Math/Functions.h>
41 #include <Core/Utils/extension/logger.hpp>
42
43 // Cdaskr declaration
44 extern "C" int _daskr_ddaskr_(
45 int (*res) (double *t, double *y, double *yprime, double* cj, double *delta, int *ires, double *rpar, int* ipar),
46 int *neq,
47 double *t,
48 double *y,
49 double *yprime,
50 double *tout,
51 int *info,
52 double *rtol,
53 double *atol,
54 int *idid,
55 double *rwork,
56 int *lrw,
57 int *iwork,
58 int *liw,
59 double *rpar,
60 int *ipar,
61 int (*jac) (double *t, double *y, double *yprime, double *delta, double *pd, double *cj, double *h, double *wt, double *rpar, int* ipar),
62 int (*psol) (int *neq, double *t, double *y, double *yprime, double *savr, double *wk, double *cj, double *wght, double *wp, int *iwp, double *b, double eplin, int* ier, double *rpar, int* ipar),
63 int (*rt) (int *neq, double *t, double *y, double *yp, int *nrt, double *rval, double *rpar, int* ipar),
64 int *nrt,
65 int *jroot
66 );
67
68 39 DASSL::DASSL(IMixedSystem* system, ISolverSettings* settings)
69 : SolverDefaultImplementation(system, settings)
70 39 , _info(NULL)
71 39 , _iwork(NULL)
72 39 , _iworkAcc(NULL)
73 39 , _liw(0)
74 39 , _rwork(NULL)
75 39 , _lrw(0)
76 39 , _y(NULL)
77 39 , _yPrime(NULL)
78 39 , _rtol(NULL)
79 39 , _atol(NULL)
80 39 , _idid(0)
81 39 , _zeroFound(false)
82 39 , _zeroSign(NULL)
83 39 , _tLastEvent(0.0)
84 39 , _event_n(0)
85 39 , _properties(NULL)
86 39 , _continuous_system(NULL)
87 39 , _event_system(NULL)
88 39 , _mixed_system(NULL)
89 39 , _time_system(NULL)
90 39 , _yJac(NULL)
91 39 , _dyJac(NULL)
92 39 , _fJac(NULL)
93 39 , _maxColors(0)
94 {
95 #ifdef RUNTIME_PROFILING
96 if (MeasureTime::getInstance() != NULL)
97 {
98 measureTimeFunctionsArray = new std::vector<MeasureTimeData*>(7, NULL); //0 calcFunction //1 solve ... //6 solver statistics
99 (*measureTimeFunctionsArray)[0] = new MeasureTimeData("calcFunction");
100 (*measureTimeFunctionsArray)[1] = new MeasureTimeData("solve");
101 (*measureTimeFunctionsArray)[2] = new MeasureTimeData("writeOutput");
102 (*measureTimeFunctionsArray)[3] = new MeasureTimeData("evaluateZeroFuncs");
103 (*measureTimeFunctionsArray)[4] = new MeasureTimeData("initialize");
104 (*measureTimeFunctionsArray)[5] = new MeasureTimeData("stepCompleted");
105 (*measureTimeFunctionsArray)[6] = new MeasureTimeData("solverStatistics");
106
107 MeasureTime::addResultContentBlock(system->getModelName(), "dassl", measureTimeFunctionsArray);
108 measuredFunctionStartValues = MeasureTime::getZeroValues();
109 measuredFunctionEndValues = MeasureTime::getZeroValues();
110 solveFunctionStartValues = MeasureTime::getZeroValues();
111 solveFunctionEndValues = MeasureTime::getZeroValues();
112 solverValues = new MeasureTimeValuesSolver();
113
114 delete (*measureTimeFunctionsArray)[6]->_sumMeasuredValues;
115 (*measureTimeFunctionsArray)[6]->_sumMeasuredValues = solverValues;
116 }
117 else
118 {
119 measureTimeFunctionsArray = new std::vector<MeasureTimeData*>();
120 measuredFunctionStartValues = NULL;
121 measuredFunctionEndValues = NULL;
122 solveFunctionStartValues = NULL;
123 solveFunctionEndValues = NULL;
124 solverValues = NULL;
125 }
126 #endif
127 39 }
128
129 78 DASSL::~DASSL()
130 {
131
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_info)
132 39 delete[] _info;
133
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_iwork)
134 39 delete[] _iwork;
135
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_iworkAcc)
136 39 delete[] _iworkAcc;
137
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_rwork)
138 39 delete[] _rwork;
139
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_y)
140 39 delete[] _y;
141
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_yPrime)
142 39 delete[] _yPrime;
143
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_zeroSign)
144 39 delete[] _zeroSign;
145
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_rtol)
146 39 delete[] _rtol;
147
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_atol)
148 39 delete[] _atol;
149
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_yJac)
150 39 delete[] _yJac;
151
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_dyJac)
152 39 delete[] _dyJac;
153
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (_fJac)
154 39 delete[] _fJac;
155
156 #ifdef RUNTIME_PROFILING
157 if (measuredFunctionStartValues)
158 delete measuredFunctionStartValues;
159 if (measuredFunctionEndValues)
160 delete measuredFunctionEndValues;
161 if (solveFunctionStartValues)
162 delete solveFunctionStartValues;
163 if (solveFunctionEndValues)
164 delete solveFunctionEndValues;
165 /* solverValues is owned by (*measureTimeFunctionsArray)[6] and freed by ~MeasureTimeData() */
166 #endif
167 78 }
168
169
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 void DASSL::initialize()
170 {
171 ✗ LOGGER_WRITE_BEGIN("DASSL: initialize", LC_SOLVER, LL_DEBUG);
172
173
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _properties = dynamic_cast<ISystemProperties*>(_system);
174
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _continuous_system = dynamic_cast<IContinuous*>(_system);
175
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _event_system = dynamic_cast<IEvent*>(_system);
176 39 _mixed_system = dynamic_cast<IMixedSystem*>(_system);
177
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _time_system = dynamic_cast<ITime*>(_system);
178 39 IGlobalSettings* global_settings = dynamic_cast<ISolverSettings*>(_settings)->getGlobalSettings();
179
180 39 _tLastEvent = 0.0;
181 39 _event_n = 0;
182 39 SolverDefaultImplementation::initialize();
183 39 _dimStates = _continuous_system->getDimContinuousStates();
184 39 _dimAE = _continuous_system->getDimAE();
185 39 _dimSys = _dimStates + _dimAE;
186 39 _dimZeroFunc = _event_system->getDimZeroFunc();
187
188
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 15 times.
 39 if (_dimSys == 0)
189 24 _dimSys = 1; // introduce dummy state
190
191 // Allocate memory
192
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_info)
193 ✗ delete[] _info;
194
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_iwork)
195 ✗ delete[] _iwork;
196
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_iworkAcc)
197 ✗ delete[] _iworkAcc;
198
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_rwork)
199 ✗ delete[] _rwork;
200
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_y)
201 ✗ delete[] _y;
202
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_yPrime)
203 ✗ delete[] _yPrime;
204
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_zeroSign)
205 ✗ delete[] _zeroSign;
206
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_rtol)
207 ✗ delete[] _rtol;
208
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_atol)
209 ✗ delete[] _atol;
210
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_yJac)
211 ✗ delete[] _yJac;
212
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_dyJac)
213 ✗ delete[] _dyJac;
214
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (_fJac)
215 ✗ delete[] _fJac;
216
217 39 _info = new int[20];
218 39 _liw = 40 + _dimSys;
219 39 _lrw = 60 + 9 * _dimSys + _dimSys*_dimSys + 3 * _dimZeroFunc;
220
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _iwork = new int[_liw];
221 39 _iworkAcc = new int[40];
222
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _rwork = new double[_lrw];
223
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _y = new double[_dimSys];
224
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _yPrime = new double[_dimSys];
225
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _zeroSign = new int[_dimZeroFunc];
226
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _rtol = new double[_dimSys];
227
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _atol = new double[_dimSys];
228
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _yJac = new double[_dimSys];
229
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _dyJac = new double[_dimSys];
230
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 _fJac = new double[_dimSys];
231
232 39 memset(_info, 0, 20 * sizeof(int));
233 39 memset(_iwork, 0, _liw * sizeof(int));
234 39 memset(_iworkAcc, 0, 40 * sizeof(int));
235 39 memset(_rwork, 0, _lrw * sizeof(double));
236 39 memset(_y, 0, _dimSys * sizeof(double));
237 39 memset(_yPrime, 0, _dimSys * sizeof(double));
238
239 //
240 // Setup DASSL
241 //
242
243 // Set initial values
244 39 _continuous_system->getContinuousStates(_y);
245
246 // Set tolerances
247 39 _info[1] = 1;
248
2/2
✓ Branch 0 taken 38 times.
✓ Branch 1 taken 1 time.
 39 if (_dimAE == 0) {
249 38 _atol[0] = _rtol[0] = 1.0; // in case of dummy state
250 38 _continuous_system->getNominalStates(_atol);
251
2/2
✓ Branch 0 taken 157 times.
✓ Branch 1 taken 38 times.
 195 for (int i = 0; i < _dimStates; i++) {
252
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 155 times.
 157 _atol[i] = max(_atol[i] * _settings->getATol(), 1e-10);
253
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 157 times.
 157 _rtol[i] = max(_settings->getRTol(), 1e-10);
254 }
255 }
256 else {
257
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
 3 for (int i = 0; i < _dimSys; i++) {
258 2 _atol[i] = _settings->getATol();
259 2 _rtol[i] = _settings->getRTol();
260 }
261 }
262 ✗ LOGGER_WRITE_VECTOR("atol", _atol, _dimSys, LC_SOLVER, LL_DEBUG);
263 ✗ LOGGER_WRITE_VECTOR("rtol", _rtol, _dimSys, LC_SOLVER, LL_DEBUG);
264
265 // Return after every step
266 39 _info[2] = 1;
267
268 // Use supplied Jacobian function
269
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 38 times.
 39 _info[4] = _dimAE == 0? 1: 0;
270 39 _maxColors = _system->getAMaxColors();
271
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 39 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 39 if (_system->isAnalyticJacobianGenerated() && _continuous_system->getDimContinuousStates() > 0)
272 {
273 ✗ LOGGER_WRITE("Jacobian size " + to_string(_dimSys) + ", generated symbolically", LC_SOLVER, LL_DEBUG);
274 }
275
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 26 times.
 39 else if (_maxColors > 0)
276 {
277 ✗ LOGGER_WRITE("Jacobian size " + to_string(_dimSys) + " with " + to_string(_maxColors) + " colors", LC_SOLVER, LL_DEBUG);
278 }
279 else
280 {
281 ✗ LOGGER_WRITE("Jacobian size " + to_string(_dimSys) + ", dense numerical", LC_SOLVER, LL_DEBUG);
282 }
283
284 // Max step size
285 39 _info[6] = 1;
286 39 _rwork[1] = _settings->getGlobalSettings()->gethOutput();
287
288 // Initial step size
289 39 _info[7] = 1;
290 39 _rwork[2] = _settings->gethInit();
291
292 // Adapt tolerances for zero crossings and end time to output interval
293
4/4
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 31 times.
✓ Branch 4 taken 37 times.
✓ Branch 5 taken 2 times.
 47 double tol = min(1e-6, max(1e-6 * _settings->getGlobalSettings()->gethOutput(), DBL_EPSILON));
294 39 _event_system->setZeroTol(tol);
295 39 _settings->setEndTimeTol(tol);
296
297 ✗ LOGGER_WRITE_END(LC_SOLVER, LL_DEBUG);
298 39 }
299
300 3799 void DASSL::solve(const SOLVERCALL action)
301 {
302 3799 bool writeEventOutput = (_settings->getGlobalSettings()->getOutputPointType() == OPT_ALL);
303 3799 bool writeOutput = !(_settings->getGlobalSettings()->getOutputPointType() == OPT_NONE);
304
305 #ifdef RUNTIME_PROFILING
306 MEASURETIME_REGION_DEFINE(dasslSolveFunctionHandler, "solve");
307 if (MeasureTime::getInstance() != NULL)
308 MEASURETIME_START(solveFunctionStartValues, dasslSolveFunctionHandler, "solve");
309 #endif
310
311
2/4
✓ Branch 0 taken 3799 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3799 times.
 3799 if (!_settings || !_system)
312 ✗ throw ModelicaSimulationError(SOLVER, "DASSL::solve missing system or settings");
313
314 // prepare solver and system
315
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 3760 times.
 3799 if ((action & RECORDCALL) && (action & FIRST_CALL))
316 {
317 #ifdef RUNTIME_PROFILING
318 MEASURETIME_REGION_DEFINE(dasslInitializeHandler, "DASSLInitialize");
319 if (MeasureTime::getInstance() != NULL)
320 MEASURETIME_START(measuredFunctionStartValues, dasslInitializeHandler, "DASSLInitialize");
321 #endif
322
323 39 initialize();
324
325 #ifdef RUNTIME_PROFILING
326 if (MeasureTime::getInstance() != NULL)
327 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[4], dasslInitializeHandler);
328 #endif
329
330
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (writeOutput)
331 39 writeToFile(_accStps, _tCurrent, _h);
332
333 39 return;
334 }
335
336
2/2
✓ Branch 0 taken 1908 times.
✓ Branch 1 taken 1852 times.
 3760 if ((action & RECORDCALL) && !(action & FIRST_CALL))
337 {
338 1908 writeToFile(_accStps, _tCurrent, _h);
339 1908 return;
340 }
341
342 // recored new state after time event
343
2/2
✓ Branch 0 taken 1821 times.
✓ Branch 1 taken 31 times.
 1852 if (action & RECALL)
344 {
345 1821 _firstStep = true;
346
1/2
✓ Branch 0 taken 1821 times.
✗ Branch 1 not taken.
 1821 if (writeOutput || writeEventOutput)
347 1821 writeToFile(_accStps, _tCurrent, _h);
348 1821 _continuous_system->getContinuousStates(_y);
349 }
350
351 // solver shall continue
352 1852 _solverStatus = ISolver::CONTINUE;
353
354
3/4
✓ Branch 0 taken 1852 times.
✓ Branch 1 taken 1852 times.
✓ Branch 2 taken 1852 times.
✗ Branch 3 not taken.
 3704 while ((_solverStatus & ISolver::CONTINUE) && !_interrupt)
355 {
356 // call solver
357 1852 DASSLCore();
358 }
359
360 // not successful and not interruped by user
361
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1852 times.
 1852 if (_solverStatus == ISolver::SOLVERERROR)
362 {
363 ✗ throw ModelicaSimulationError(SOLVER, "DASSL: solve failed with idid = " + to_string(_idid));
364 }
365
366 1852 _firstCall = false;
367
368 #ifdef RUNTIME_PROFILING
369 if (MeasureTime::getInstance() != NULL)
370 {
371 MEASURETIME_END(solveFunctionStartValues, solveFunctionEndValues, (*measureTimeFunctionsArray)[1], dasslSolveFunctionHandler);
372
373 // DASKR reports its statistics in the integer work array, see writeSimulationInfo().
374 // _iworkAcc holds the counts accumulated over previous restarts, _iwork those of the current one.
375 unsigned long long nst = _iworkAcc[10] + _iwork[10]; // steps taken
376 unsigned long long nre = _iworkAcc[11] + _iwork[11]; // residual evaluations
377 unsigned long long netf = _iworkAcc[13] + _iwork[13]; // error test failures
378
379 MeasureTimeValuesSolver solverVals = MeasureTimeValuesSolver(nre, netf);
380 (*measureTimeFunctionsArray)[6]->_sumMeasuredValues->_numCalcs += nst;
381 (*measureTimeFunctionsArray)[6]->_sumMeasuredValues->add(&solverVals);
382 }
383 #endif
384 }
385
386 2848 bool DASSL::isInterrupted()
387 {
388
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2848 times.
 2848 if (_interrupt)
389 {
390 ✗ _solverStatus = DONE;
391 ✗ return true;
392 }
393 else
394 {
395 return false;
396 }
397 }
398
399
1/2
✓ Branch 0 taken 1852 times.
✗ Branch 1 not taken.
 1852 void DASSL::DASSLCore()
400 {
401 ✗ LOGGER_WRITE_BEGIN("DASSL: solve at t = " + to_string(_tCurrent), LC_SOLVER, LL_DEBUG);
402
403 1852 bool writeEventOutput = (_settings->getGlobalSettings()->getOutputPointType() == OPT_ALL);
404 1852 bool writeOutput = !(_settings->getGlobalSettings()->getOutputPointType() == OPT_NONE);
405
406 1852 _info[0] = 0; // (re-)start dassl
407
408
3/4
✓ Branch 0 taken 310220 times.
✓ Branch 1 taken 1852 times.
✓ Branch 2 taken 310220 times.
✗ Branch 3 not taken.
 312072 while ((_solverStatus & ISolver::CONTINUE) && !_interrupt)
409 {
410
2/2
✓ Branch 0 taken 3276 times.
✓ Branch 1 taken 306944 times.
 310220 if (_info[0] == 0)
411 {
412 // accumulate previous solver stats upon restart
413
2/2
✓ Branch 0 taken 85176 times.
✓ Branch 1 taken 3276 times.
 88452 for (int i = 10; i <= 35; i++)
414 85176 _iworkAcc[i] += _iwork[i];
415 }
416
417
2/2
✓ Branch 1 taken 310159 times.
✓ Branch 2 taken 61 times.
 310220 if (_tEnd - _tCurrent > _settings->getEndTimeTol())
418 310159 _daskr_ddaskr_(_res, &_dimSys, &_tCurrent, _y, _yPrime, &_tEnd, _info,
419 _rtol, _atol, &_idid, _rwork, &_lrw, _iwork, &_liw, /*rpar*/NULL,
420 (int *)this, _jac, /*psol*/NULL, _rt, &_dimZeroFunc, _zeroSign);
421 else
422 61 _idid = 3; // daskr would return -33 as end time too close
423
424
2/2
✓ Branch 0 taken 3276 times.
✓ Branch 1 taken 306944 times.
 310220 if (_idid != 1)
425 ✗ LOGGER_WRITE("proceed to t = " + to_string(_tCurrent) + ", idid = " + to_string(_idid), LC_SOLVER, LL_DEBUG);
426
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 310220 times.
 310220 if (_idid < 0)
427 {
428 ✗ _rejStps ++;
429 ✗ _solverStatus = ISolver::SOLVERERROR;
430 ✗ break;
431 }
432
2/2
✓ Branch 0 taken 1852 times.
✓ Branch 1 taken 308368 times.
 310220 else if (1 < _idid && _idid < 4)
433 1852 _solverStatus = DONE;
434
435 try
436 {
437 // complete step for system and check for terminate
438
2/3
✓ Branch 1 taken 310220 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 310220 times.
 310220 if (_continuous_system->stepCompleted(_tCurrent))
439 ✗ _solverStatus = DONE;
440
441
1/2
✓ Branch 0 taken 310220 times.
✗ Branch 1 not taken.
 310220 if (writeOutput)
442 {
443
2/2
✓ Branch 0 taken 1852 times.
✓ Branch 1 taken 308368 times.
 310220 if (_idid == 3)
444
1/1
✓ Branch 1 taken 1852 times.
 1852 _time_system->setTime(_tEnd); // interpolated time point
445
1/1
✓ Branch 1 taken 310220 times.
 310220 _continuous_system->setContinuousStates(_y);
446
2/2
✓ Branch 0 taken 570 times.
✓ Branch 1 taken 309650 times.
 310220 if (_dimAE > 0) {
447
1/1
✓ Branch 1 taken 570 times.
 570 _mixed_system->setAlgebraicDAEVars(_y + _dimStates);
448
1/1
✓ Branch 1 taken 570 times.
 570 _continuous_system->setStateDerivatives(_yPrime);
449 }
450
1/1
✓ Branch 1 taken 310220 times.
 310220 _continuous_system->evaluateAll(IContinuous::CONTINUOUS);
451
1/1
✓ Branch 1 taken 310220 times.
 310220 writeToFile(_accStps, _tCurrent, _h);
452 }
453
454 #ifdef RUNTIME_PROFILING
455 MEASURETIME_REGION_DEFINE(dasslStepCompletedHandler, "DASSLStepCompleted");
456 if (MeasureTime::getInstance() != NULL)
457 MEASURETIME_START(measuredFunctionStartValues, dasslStepCompletedHandler, "DASSLStepCompleted");
458 if (MeasureTime::getInstance() != NULL)
459 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[5], dasslStepCompletedHandler);
460 #endif
461
462 // Perform state selection
463
1/1
✓ Branch 1 taken 310220 times.
 310220 bool state_selection = stateSelection();
464
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 310220 times.
 310220 if (state_selection) {
465 ✗ _continuous_system->getContinuousStates(_y);
466 }
467 310220 _zeroFound = false;
468
469 // Check for found root
470
4/5
✓ Branch 0 taken 1424 times.
✓ Branch 1 taken 308796 times.
✓ Branch 3 taken 1424 times.
✓ Branch 5 taken 1424 times.
✗ Branch 6 not taken.
 310220 if (_idid == 5 && !isInterrupted())
471 {
472
1/2
✓ Branch 0 taken 1424 times.
✗ Branch 1 not taken.
 1424 _zeros ++;
473 ✗ LOGGER_WRITE_VECTOR("jroot", _zeroSign, _dimZeroFunc, LC_SOLVER, LL_DEBUG);
474 // DASSL sets _tCurrent to the time where the first event occurred
475 1424 double _abs = fabs(_tLastEvent - _tCurrent);
476 1424 _zeroFound = true;
477
478
3/4
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 1402 times.
✓ Branch 2 taken 22 times.
✗ Branch 3 not taken.
 1424 if (_abs < 1e-3 && _event_n == 0)
479 {
480 22 _tLastEvent = _tCurrent;
481 22 _event_n++;
482 }
483
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1402 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1402 else if ((_abs < 1e-3) && (_event_n >= 1 && _event_n < 500))
484 {
485 ✗ _event_n++;
486 }
487
1/2
✓ Branch 0 taken 1402 times.
✗ Branch 1 not taken.
 1402 else if (_abs >= 1e-3)
488 {
489 //restart event counter
490 1402 _tLastEvent = _tCurrent;
491 1402 _event_n = 0;
492 }
493 else
494 {
495 ✗ _solverStatus = ISolver::SOLVERERROR;
496 ✗ break;
497 }
498
499 // DASSL has interpolated the states at time _tCurrent
500
1/1
✓ Branch 1 taken 1424 times.
 1424 _time_system->setTime(_tCurrent);
501
502 // To get steep steps in the result file, two value points (P1 and P2) are added
503 //
504 // Y | (P2) X...........
505 // | :
506 // | :
507 // |........X (P1)
508 // |---------------------------------->
509 // | ^ t
510 // _tCurrent
511
512 // Write the values of (P1) if not done via writeOutput above
513
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1424 times.
 1424 if (writeEventOutput && !writeOutput)
514 {
515 try
516 {
517 ✗ _continuous_system->evaluateAll(IContinuous::CONTINUOUS);
518 }
519 ✗ catch (std::exception& ex)
520 {
521 // if a zero crossing was dected before the event iteration was called and evalutateAll throws an error
522 // for this time step the event iteration evaluates the system with corrected values.
523 ✗ }
524
525 ✗ writeToFile(_accStps, _tCurrent, _h);
526 }
527
528
2/2
✓ Branch 0 taken 29127 times.
✓ Branch 1 taken 1424 times.
 30551 for (int i = 0; i < _dimZeroFunc; i++)
529 29127 _events[i] = (_zeroSign[i] != 0);
530
531
3/3
✓ Branch 1 taken 1424 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 1421 times.
 1424 if (_mixed_system->handleSystemEvents(_events))
532 {
533 // State variables were reinitialized, thus we have to give these values to dassl
534
1/1
✓ Branch 1 taken 3 times.
 3 _continuous_system->getContinuousStates(_y);
535 }
536 }
537
538
5/7
✓ Branch 0 taken 308796 times.
✓ Branch 1 taken 1424 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 308796 times.
✓ Branch 5 taken 1424 times.
✓ Branch 7 taken 1424 times.
✗ Branch 8 not taken.
 310220 if ((_zeroFound || state_selection) && !isInterrupted())
539 {
540 // Write the values of (P2)
541
1/2
✓ Branch 0 taken 1424 times.
✗ Branch 1 not taken.
 1424 if (writeEventOutput)
542 {
543 // If we want to write the event-results, we should evaluate the whole system again
544
1/1
✓ Branch 1 taken 1424 times.
 1424 _continuous_system->evaluateAll(IContinuous::CONTINUOUS);
545
1/1
✓ Branch 1 taken 1424 times.
 1424 writeToFile(_accStps, _tCurrent, _h);
546 }
547
548 1424 _info[0] = 0; // restart dassl
549
550 // Check for event at end time
551
2/3
✓ Branch 1 taken 1424 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1424 times.
 1424 if (_tEnd - _tCurrent <= _settings->getEndTimeTol())
552 ✗ _solverStatus = DONE;
553
2/3
✓ Branch 1 taken 1424 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1424 times.
 1424 if (_continuous_system->stepCompleted(_tCurrent))
554 ✗ _solverStatus = DONE;
555 }
556
557 310220 _accStps ++;
558 310220 _tLastSuccess = _tCurrent;
559 }
560 ✗ catch (const std::exception& ex)
561 {
562 ✗ LOGGER_WRITE("DASSL: failed step at t = " + to_string(_tCurrent) + ": " + ex.what(), LC_SOLVER, LL_ERROR);
563 ✗ _solverStatus = ISolver::SOLVERERROR;
564 break;
565 ✗ }
566 }
567
568 ✗ LOGGER_WRITE_END(LC_SOLVER, LL_DEBUG);
569 1852 }
570
571 ✗ void DASSL::setTimeOut(unsigned int time_out)
572 {
573 ✗ SimulationMonitor::setTimeOut(time_out);
574 ✗ }
575
576 ✗ void DASSL::stop()
577 {
578 ✗ SimulationMonitor::stop();
579 ✗ }
580
581 310298 bool DASSL::stateSelection()
582 {
583 310298 return SolverDefaultImplementation::stateSelection();
584 }
585
586 484480 int DASSL::_res(double *t, double *y, double *yp,
587 double *cj, double *delta, int *ires, double *rpar, int *ipar)
588 {
589 484480 int success = ((DASSL *)ipar)->calcFunction(*t, y, yp, delta);
590
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 484480 times.
 484480 if (!success)
591 ✗ *ires = -1;
592 484480 return 0;
593 }
594
595 1763372 int DASSL::calcFunction(const double& time, const double* y, const double *yp, double* f)
596 {
597 int success = 0;
598
599 #ifdef RUNTIME_PROFILING
600 MEASURETIME_REGION_DEFINE(dasslCalcFunctionHandler, "DASSLCalcFunction");
601 if (MeasureTime::getInstance() != NULL)
602 {
603 MEASURETIME_START(measuredFunctionStartValues, dasslCalcFunctionHandler, "DASSLCalcFunction");
604 }
605 #endif
606
607 try
608 {
609 1763372 f[0] = 0.0; // in case of dummy state
610
1/1
✓ Branch 1 taken 1763372 times.
 1763372 _time_system->setTime(time);
611
1/1
✓ Branch 1 taken 1763372 times.
 1763372 _continuous_system->setContinuousStates(y);
612
2/2
✓ Branch 0 taken 1762660 times.
✓ Branch 1 taken 712 times.
 1763372 if (_dimAE == 0) {
613
1/1
✓ Branch 1 taken 1762660 times.
 1762660 _continuous_system->evaluateODE(IContinuous::CONTINUOUS);
614
1/1
✓ Branch 1 taken 1762660 times.
 1762660 _continuous_system->getRHS(f);
615
2/2
✓ Branch 0 taken 28705418 times.
✓ Branch 1 taken 1762660 times.
 30468078 for (int i = 0; i < _dimStates; i++)
616 28705418 f[i] -= yp[i];
617 }
618 else {
619
1/1
✓ Branch 1 taken 712 times.
 712 _mixed_system->setAlgebraicDAEVars(y + _dimStates);
620
1/1
✓ Branch 1 taken 712 times.
 712 _continuous_system->setStateDerivatives(yp);
621
1/1
✓ Branch 1 taken 712 times.
 712 _continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
622
1/1
✓ Branch 1 taken 712 times.
 712 _mixed_system->getResidual(f);
623 }
624 success = 1;
625 }
626 ✗ catch (std::exception & ex)
627 {
628 ✗ LOGGER_WRITE("DASSL: failed evaluation of residual at t = " + to_string(_tCurrent) + ": " + ex.what(), LC_SOLVER, LL_DEBUG);
629 ✗ }
630
631 #ifdef RUNTIME_PROFILING
632 if (MeasureTime::getInstance() != NULL)
633 {
634 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[0], dasslCalcFunctionHandler);
635 }
636 #endif
637
638 1763372 return success;
639 }
640
641 302299 int DASSL::_rt(int *neq, double *t, double *y, double *yp,
642 int *nrt, double *rval, double *rpar, int *ipar)
643 {
644 302299 int success = ((DASSL *)ipar)->calcRoots(*t, y, rval);
645
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 302299 times.
 302299 if (!success)
646 ✗ memset(rval, 0, *nrt * sizeof(double));
647 302299 return 0;
648 }
649
650 302299 int DASSL::calcRoots(double t, const double *y, double *zeroValue)
651 {
652 int success = 0;
653
654 #ifdef RUNTIME_PROFILING
655 MEASURETIME_REGION_DEFINE(dasslEvalZeroHandler, "evaluateZeroFuncs");
656 if (MeasureTime::getInstance() != NULL)
657 {
658 MEASURETIME_START(measuredFunctionStartValues, dasslEvalZeroHandler, "evaluateZeroFuncs");
659 }
660 #endif
661
662 try
663 {
664
1/1
✓ Branch 1 taken 302299 times.
 302299 _time_system->setTime(t);
665
1/1
✓ Branch 1 taken 302299 times.
 302299 _continuous_system->setContinuousStates(y);
666
1/1
✓ Branch 1 taken 302299 times.
 302299 _continuous_system->evaluateZeroFuncs(IContinuous::DISCRETE);
667
1/1
✓ Branch 1 taken 302299 times.
 302299 _event_system->getZeroFunc(zeroValue);
668 success = 1;
669 }
670 ✗ catch (std::exception & ex)
671 {
672 ✗ LOGGER_WRITE("DASSL: failed evaluation of roots at t = " + to_string(_tCurrent) + ": " + ex.what(), LC_SOLVER, LL_WARNING);
673 ✗ }
674
675 #ifdef RUNTIME_PROFILING
676 if (MeasureTime::getInstance() != NULL)
677 {
678 MEASURETIME_END(measuredFunctionStartValues, measuredFunctionEndValues, (*measureTimeFunctionsArray)[3], dasslEvalZeroHandler);
679 }
680 #endif
681
682 302299 return success;
683 }
684
685 118715 int DASSL::_jac(double *t, double *y, double *yp, double *delta,
686 double *pd, double *cj, double *h, double *wt, double *rpar, int *ipar)
687 {
688 118715 int success = ((DASSL *)ipar)->calcJacobian(*t, y, yp, delta, pd, *cj, *h, wt);
689 118715 int n = ((DASSL *)ipar)->_dimSys;
690
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 118715 times.
 118715 if (!success)
691 ✗ memset(pd, 0, n * n * sizeof(double));
692 else
693
2/2
✓ Branch 0 taken 1906639 times.
✓ Branch 1 taken 118715 times.
 2025354 for (int i = 0; i < n; i++)
694 1906639 pd[i*n + i] -= *cj;
695 118715 return 0;
696 }
697
698 118715 int DASSL::calcJacobian(double t, double *y, double *yp, double *delta,
699 double *pd, double cj, double h, double *wt)
700 {
701 int success = 0;
702
703 try
704 {
705
2/7
✓ Branch 1 taken 118715 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 118715 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 118715 if (_system->isAnalyticJacobianGenerated() && _continuous_system->getDimContinuousStates() > 0)
706 {
707 ✗ memcpy(pd, &_system->getJacobian().data()[0], _dimSys * _dimSys * sizeof(double));
708 }
709 else
710 {
711
2/2
✓ Branch 0 taken 1906639 times.
✓ Branch 1 taken 118715 times.
 2025354 for (int j = 0; j < _dimSys; j++)
712 {
713
4/4
✓ Branch 0 taken 86354 times.
✓ Branch 1 taken 1820285 times.
✓ Branch 2 taken 1313805 times.
✓ Branch 3 taken 592834 times.
 3813278 _dyJac[j] = max(1e-10, 1e-8 * max(max(abs(y[j]), abs(h * yp[j])), abs(1.0 / wt[j])));
714 1906639 _dyJac[j] = y[j] + _dyJac[j];
715 1906639 _dyJac[j] -= y[j];
716 1906639 _yJac[j] = y[j];
717 }
718
719
2/2
✓ Branch 0 taken 116699 times.
✓ Branch 1 taken 2016 times.
 118715 if (_maxColors > 0) // colored numerical
720 {
721
2/2
✓ Branch 0 taken 1276734 times.
✓ Branch 1 taken 116699 times.
 1393433 for (int color = 1; color <= _maxColors; color++)
722 {
723
3/3
✓ Branch 1 taken 1276734 times.
✓ Branch 3 taken 1904481 times.
✓ Branch 4 taken 1276734 times.
 3181215 for (int j: _system->getAColumnsOfColor(color))
724 {
725 1904481 _yJac[j] += _dyJac[j];
726 }
727
728
1/1
✓ Branch 1 taken 1276734 times.
 1276734 calcFunction(t, _yJac, yp, _fJac);
729
730
3/3
✓ Branch 1 taken 1276734 times.
✓ Branch 3 taken 1904481 times.
✓ Branch 4 taken 1276734 times.
 3181215 for (int j: _system->getAColumnsOfColor(color))
731 {
732 1904481 int startOfColumn = j * _dimSys;
733
3/3
✓ Branch 1 taken 1904481 times.
✓ Branch 3 taken 7779495 times.
✓ Branch 4 taken 1904481 times.
 9683976 for (int i: _system->getADependenciesOfColumn(j))
734 {
735 7779495 pd[startOfColumn + i] = (_fJac[i] - delta[i]) / _dyJac[j];
736 }
737 1904481 _yJac[j] = y[j];
738 }
739 }
740 }
741 else // dense numerical
742 {
743
2/2
✓ Branch 0 taken 2158 times.
✓ Branch 1 taken 2016 times.
 4174 for (int j = 0; j < _dimSys; j++)
744 {
745 2158 _yJac[j] += _dyJac[j];
746
747
1/1
✓ Branch 1 taken 2158 times.
 2158 calcFunction(t, _yJac, yp, _fJac);
748
749 2158 int startOfColumn = j * _dimSys;
750
2/2
✓ Branch 0 taken 2584 times.
✓ Branch 1 taken 2158 times.
 4742 for (int i = 0; i < _dimSys; i++)
751 {
752 2584 pd[startOfColumn + i] = (_fJac[i] - delta[i]) / _dyJac[j];
753 }
754
755 2158 _yJac[j] = y[j];
756 }
757 }
758 }
759 success = 1;
760 }
761 ✗ catch (std::exception & ex)
762 {
763 ✗ LOGGER_WRITE("DASSL: failed evaluation of Jacobian at t = " + to_string(_tCurrent) + ": " + ex.what(), LC_SOLVER, LL_WARNING);
764 ✗ }
765
766 118715 return success;
767 }
768
769 31 void DASSL::writeSimulationInfo()
770 {
771 // don't write before memory has been initialized
772
2/4
✓ Branch 0 taken 31 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 31 times.
✗ Branch 3 not taken.
 31 if (_rwork == NULL || _iwork == NULL)
773 return;
774
775
2/2
✓ Branch 0 taken 806 times.
✓ Branch 1 taken 31 times.
 837 for (int i = 10; i <= 35; i++)
776 806 _iworkAcc[i] += _iwork[i];
777
778
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: steps taken nst = " + to_string(_iworkAcc[10]), LC_SOLVER, LL_INFO);
779
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: residual evaluations nre = " + to_string(_iworkAcc[11]), LC_SOLVER, LL_INFO);
780
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: jacobian evaluations nje = " + to_string(_iworkAcc[12]), LC_SOLVER, LL_INFO);
781
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: root evaluations nrt = " + to_string(_iworkAcc[35]), LC_SOLVER, LL_INFO);
782
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: error test failures netf = " + to_string(_iworkAcc[13]), LC_SOLVER, LL_INFO);
783
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: nonlinear convergence failures ncfn = " + to_string(_iworkAcc[14]), LC_SOLVER, LL_INFO);
784 //LOGGER_WRITE("DASSL: linear convergence failures ncfl = " + to_string(_iworkAcc[15]), LC_SOLVER, LL_INFO);
785
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: nonlinear iterations nni = " + to_string(_iworkAcc[18]), LC_SOLVER, LL_INFO);
786 //LOGGER_WRITE("DASSL: linear iterations nli = " + to_string(_iworkAcc[19]), LC_SOLVER, LL_INFO);
787 //LOGGER_WRITE("DASSL: preconditioning calls nps = " + to_string(_iworkAcc[20]), LC_SOLVER, LL_INFO);
788
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: last evaluation time t = " + to_string(_rwork[3]), LC_SOLVER, LL_INFO);
789 //LOGGER_WRITE("DASSL: next step size h = " + to_string(_rwork[2]), LC_SOLVER, LL_INFO);
790
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: last step size h = " + to_string(_rwork[6]), LC_SOLVER, LL_INFO);
791 //LOGGER_WRITE("DASSL: next used order k = " + to_string(_iwork[6]), LC_SOLVER, LL_INFO);
792
3/6
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 LOGGER_WRITE("DASSL: last used order k = " + to_string(_iwork[7]), LC_SOLVER, LL_INFO);
793 }
794
795 /** @} */ // end of solverDASSL
796