Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 453
Functions: 0.0% 0 / 0 / 13
Branches: 0.0% 0 / 0 / 240
OMCompiler/SimulationRuntime/c/simulation/solver/cvode_solver.c
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 /* Standard C headers */
29 #include <float.h>
30 #include <math.h>
31 #include <string.h>
32 #include <stdio.h>
33 #include <stdlib.h>
34
35 #include "cvode_solver.h"
36
37 /* OMC headers */
38 #include "../../util/context.h"
39 #include "../options.h"
40 #include "../solver/external_input.h"
41 #include "../arrayIndex.h"
42 #include "model_help.h"
43 #include "omc_math.h"
44
45 #include "../../util/omc_error.h"
46 #include "../../gc/omc_gc.h"
47
48 #include "dassl.h"
49 #include "epsilon.h"
50 #ifndef OMC_FMI_RUNTIME
51 #include "../jacobian_util.h"
52 #include "sundials_util.h"
53 #endif
54
55
56 #ifdef WITH_SUNDIALS
57
58 #define CVODE_LMM_MAX 2
59 const char *CVODE_LMM_NAME[CVODE_LMM_MAX + 1] = {
60 "undefined",
61 "CV_ADAMS", /* 1 */
62 "CV_BDF" /* 2 */
63 };
64
65 const char *CVODE_LMM_DESC[CVODE_LMM_MAX + 1] = {
66 "undefined",
67 "Adams-Moulton linear multistep method. Use together with CV_ITER_FIXED_POINT for nonstiff problems.",
68 "BDF linear multistep method. Use together with CV_ITER_NEWTON for stiff problems. Default option."};
69
70 #define CVODE_ITER_MAX 2
71 const char *CVODE_ITER_NAME[CVODE_ITER_MAX + 1] = {
72 "undefined",
73 "CV_ITER_FIXED_POINT", /* 1 */
74 "CV_ITER_NEWTON" /* 2 */
75 };
76
77 const char *CVODE_ITER_DESC[CVODE_ITER_MAX + 1] = {
78 "undefined",
79 "Nonlinear system solution through fixed-point iterations",
80 "Nonlinear system solution through Newton iterations"
81 };
82
83 /* Internal function prototypes */
84 int cvodeRightHandSideODEFunction(sunrealtype time, N_Vector y, N_Vector ydot, void *userData);
85 void cvodeGetConfig(CVODE_CONFIG *config, threadData_t *threadData, sunbooleantype isFMI);
86
87 /**
88 * @brief Computes the ODE right-hand side for a given value of the independent variable t and state vector y
89 *
90 * @param time is the current value of the independent variable
91 * @param y is the current value of the dependent variable vector, y(t).
92 * @param ydot is the output vector f(t, y).
93 * @param userData user data containing CVODE_SOLVER
94 * @return int
95 */
96 ✗ int cvodeRightHandSideODEFunction(sunrealtype time, N_Vector y, N_Vector ydot, void *userData)
97 {
98 /* Variables */
99 CVODE_SOLVER *cvodeData;
100 DATA *data;
101 threadData_t *threadData;
102 long int i;
103 ✗ int success = 0, retVal = 0;
104 int saveJumpState;
105
106 /* Access userData */
107 cvodeData = (CVODE_SOLVER *)userData;
108 ✗ data = cvodeData->simData->data;
109 ✗ threadData = cvodeData->simData->threadData;
110
111 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "### eval cvodeRightHandSideODEFunction ###");
112
113 /* TODO: Add scaling of y and ydot */
114
115 ✗ if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
116 {
117 ✗ setContext(data, time, CONTEXT_ODE);
118 }
119 /* Set time */
120 ✗ data->localData[0]->timeValue = time;
121
122 ✗ saveJumpState = threadData->currentErrorStage;
123 ✗ threadData->currentErrorStage = ERROR_INTEGRATOR;
124
125 /* try */
126 #if !defined(OMC_EMCC)
127 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
128 #endif
129
130 /*
131 fix issue https://github.com/OpenModelica/OpenModelica/issues/13582
132 Update y*/
133 ✗ for (i = 0; i < cvodeData->N; i++)
134 {
135 ✗ data->localData[0]->realVars[i] = NV_Ith_S(y, i);
136 }
137
138 /* Debug print for states (input) */
139 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER_V))
140 {
141 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "y at time=%f", time);
142 ✗ for (i = 0; i < cvodeData->N; i++)
143 {
144 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, "y[%ld] = %e", i, NV_Ith_S(y, i));
145 }
146 ✗ messageClose(OMC_LOG_SOLVER_V);
147 }
148
149 /* Read input vars (exclude from timer) */
150 ✗ if (measure_time_flag)
151 ✗ rt_accumulate(SIM_TIMER_SOLVER);
152 #ifndef OMC_FMI_RUNTIME
153 ✗ externalInputUpdate(data);
154 ✗ data->callback->input_function(data, threadData);
155 #endif
156 ✗ if (measure_time_flag)
157 ✗ rt_tick(SIM_TIMER_SOLVER);
158
159 /* eval function ODE (exclude from timer) */
160 ✗ if (measure_time_flag)
161 ✗ rt_accumulate(SIM_TIMER_SOLVER);
162 ✗ data->callback->functionODE(data, threadData);
163 ✗ if (measure_time_flag)
164 ✗ rt_tick(SIM_TIMER_SOLVER);
165
166 /* Update ydot */
167 ✗ for (i = 0; i < cvodeData->N; i++)
168 {
169 ✗ NV_Ith_S(ydot, i) = data->localData[0]->realVars[cvodeData->N + i];
170 }
171
172 /* Debug print for derived states (output) */
173 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER_V))
174 {
175 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "ydot at time=%f", time);
176 ✗ for (i = 0; i < cvodeData->N; i++)
177 {
178 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, "ydot[%ld] = %e", i, NV_Ith_S(ydot, i));
179 }
180 ✗ messageClose(OMC_LOG_SOLVER_V);
181 }
182
183 /* TODO: Scale result */
184
185 /* catch */
186 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { success = 1; }
187 #if !defined(OMC_EMCC)
188 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
189 #endif
190
191 ✗ for (i = 0; success && i < cvodeData->N; i++)
192 {
193 ✗ success = isfinite(NV_Ith_S(ydot, i));
194 }
195
196 ✗ if (!success)
197 {
198 ✗ retVal = 1; /* Recoverable error, reduce step size and retry */
199 #ifndef OMC_FMI_RUNTIME
200 /* At the start point fall back to the derivatives DASSL and IDA start from:
201 * the model can be singular at exactly that point. */
202 ✗ if (cvodeData->fStart != NULL && time == cvodeData->startTime
203 ✗ && memcmp(N_VGetArrayPointer(y), cvodeData->yStart, cvodeData->N * sizeof(double)) == 0)
204 {
205 ✗ memcpy(N_VGetArrayPointer(ydot), cvodeData->fStart, cvodeData->N * sizeof(double));
206 retVal = 0;
207 }
208 #endif
209 }
210
211 ✗ threadData->currentErrorStage = saveJumpState;
212
213 ✗ if (data->simulationInfo->currentContext == CONTEXT_ODE)
214 {
215 ✗ unsetContext(data);
216 }
217 ✗ messageClose(OMC_LOG_SOLVER_V);
218 ✗ if (measure_time_flag)
219 ✗ rt_accumulate(SIM_TIMER_SOLVER);
220
221 ✗ return retVal;
222 }
223
224
225 #ifndef OMC_FMI_RUNTIME
226 /**
227 * @brief Colored numerical Jacobian J = df/dy in the sparse matrix Jac.
228 *
229 * @param t Independent variable (time).
230 * @param y Dependent variable vector, restored on return.
231 * @param fy Current value of f(t,y).
232 * @param Jac Output Jacobian.
233 * @param cvodeData CVODE solver data.
234 * @return int 0 on success, 1 if a perturbed f could not be evaluated.
235 */
236 ✗ static int jacColoredNumericalSparse(double t, N_Vector y, N_Vector fy, SUNMatrix Jac, CVODE_SOLVER *cvodeData)
237 {
238 ✗ DATA *data = cvodeData->simData->data;
239 ✗ const SPARSE_PATTERN *sp = getJacobianCscPattern(getSymbolicOdeJacobian(data));
240 ✗ double *states = N_VGetArrayPointer(y);
241 ✗ double *f = N_VGetArrayPointer(fy);
242 ✗ double *fProbe = N_VGetArrayPointer(cvodeData->fProbe);
243 ✗ double *abstol = N_VGetArrayPointer(cvodeData->absoluteTolerance);
244 ✗ double *ysave = cvodeData->ysave;
245 ✗ double *delta_hh = cvodeData->delta_hh;
246 ✗ double rtol = data->simulationInfo->tolerance;
247 double h, hf;
248 long int i, ii;
249 unsigned int nth;
250 int retVal = 0;
251
252 ✗ CVodeGetCurrentStep(cvodeData->cvode_mem, &h);
253 ✗ setContext(data, t, CONTEXT_JACOBIAN);
254
255 ✗ for (i = 0; i < sp->maxColors && retVal == 0; i++)
256 {
257 ✗ for (ii = 0; ii < cvodeData->N; ii++)
258 {
259 ✗ if (sp->colorCols[ii] - 1 == i)
260 {
261 ✗ hf = h * f[ii];
262 /* abstol is nominal*rtol */
263 ✗ delta_hh[ii] = numericalJacobianStep(states[ii], hf, rtol * fabs(states[ii]) + abstol[ii],
264 ✗ cvodeData->jacNominalFactor * abstol[ii] / rtol);
265 ✗ delta_hh[ii] = (hf >= 0 ? delta_hh[ii] : -delta_hh[ii]);
266 ✗ delta_hh[ii] = (states[ii] + delta_hh[ii]) - states[ii];
267 ✗ ysave[ii] = states[ii];
268 ✗ states[ii] += delta_hh[ii];
269 ✗ delta_hh[ii] = 1. / delta_hh[ii];
270 }
271 }
272
273 ✗ retVal = cvodeRightHandSideODEFunction(t, y, cvodeData->fProbe, cvodeData);
274 ✗ increaseJacContext(data);
275
276 ✗ for (ii = 0; ii < cvodeData->N; ii++)
277 {
278 ✗ if (sp->colorCols[ii] - 1 == i)
279 {
280 ✗ for (nth = sp->leadindex[ii]; retVal == 0 && nth < sp->leadindex[ii + 1]; nth++)
281 {
282 ✗ setJacElementSundialsSparse(sp->index[nth], ii, nth, (fProbe[sp->index[nth]] - f[sp->index[nth]]) * delta_hh[ii], Jac, cvodeData->N);
283 }
284 ✗ states[ii] = ysave[ii];
285 }
286 }
287 }
288 ✗ setSundialsSparseColPtrs(sp, Jac);
289
290 ✗ unsetContext(data);
291 ✗ return retVal == 0 ? 0 : 1;
292 }
293
294 /**
295 * @brief Colored symbolical Jacobian J = df/dy in the sparse matrix Jac.
296 *
297 * The model already holds the point: CVODE evaluated f(t,y) last.
298 *
299 * @param t Independent variable (time).
300 * @param Jac Output Jacobian.
301 * @param cvodeData CVODE solver data.
302 * @return int 0 on success, 1 if the model raised an error.
303 */
304 ✗ static int jacColoredSymbolicalSparse(double t, SUNMatrix Jac, CVODE_SOLVER *cvodeData)
305 {
306 ✗ DATA *data = cvodeData->simData->data;
307 ✗ threadData_t *threadData = cvodeData->simData->threadData;
308 ✗ JACOBIAN *jac = getSymbolicOdeJacobian(data);
309 ✗ int saveJumpState, success = 0;
310
311 ✗ SUNMatZero(Jac);
312 ✗ setContext(data, t, CONTEXT_SYM_JACOBIAN);
313 ✗ saveJumpState = threadData->currentErrorStage;
314 ✗ threadData->currentErrorStage = ERROR_INTEGRATOR;
315
316 #if !defined(OMC_EMCC)
317 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
318 #endif
319 ✗ setSundialsSparsePattern(jac, Jac);
320 ✗ evalJacobian(data, threadData, jac, NULL, SM_DATA_S(Jac), FALSE);
321 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { success = 1; }
322 #if !defined(OMC_EMCC)
323 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
324 #endif
325
326 ✗ threadData->currentErrorStage = saveJumpState;
327 ✗ unsetContext(data);
328 ✗ return success ? 0 : 1;
329 }
330
331 /**
332 * @brief CVLsJacFn: J = df/dy for the sparse linear solver.
333 *
334 * @param t Independent variable (time).
335 * @param y Dependent variable vector.
336 * @param fy Current value of f(t,y).
337 * @param Jac Output Jacobian.
338 * @param user_data CVODE solver data.
339 * @param tmp1 Unused work space.
340 * @param tmp2 "
341 * @param tmp3 "
342 * @return int 0 on success, positive value for a recoverable error.
343 */
344 ✗ static int callSparseJacobian(double t, N_Vector y, N_Vector fy,
345 SUNMatrix Jac, void *user_data,
346 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
347 {
348 CVODE_SOLVER *cvodeData = (CVODE_SOLVER *)user_data;
349 ✗ JACOBIAN_METHOD method = cvodeData->config.jacobianMethod;
350 int retVal;
351
352 ✗ if (measure_time_flag)
353 ✗ rt_accumulate(SIM_TIMER_SOLVER);
354 ✗ rt_tick(SIM_TIMER_JACOBIAN);
355
356 ✗ if (method == COLOREDSYMJAC || method == COLOREDSYMJACADJ || method == BICOLOREDSYMJAC)
357 {
358 ✗ retVal = jacColoredSymbolicalSparse(t, Jac, cvodeData);
359 }
360 else
361 {
362 ✗ retVal = jacColoredNumericalSparse(t, y, fy, Jac, cvodeData);
363 }
364
365 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC))
366 {
367 ✗ sundialsPrintSparseMatrix(Jac, "CVODE-Solver: Matrix A", OMC_LOG_JAC);
368 }
369
370 ✗ rt_accumulate(SIM_TIMER_JACOBIAN);
371 ✗ if (measure_time_flag)
372 ✗ rt_tick(SIM_TIMER_SOLVER);
373
374 ✗ return retVal;
375 }
376 #endif /* OMC_FMI_RUNTIME */
377
378 /**
379 * @brief Root function for CVODE
380 *
381 * @param time Current time.
382 * @param y State vector.
383 * @param gout Zero crossing array.
384 * @param userData User data.
385 * @return int Will return 0 on success.
386 */
387 ✗ int rootsFunctionCVODE(double time, N_Vector y, double *gout, void *userData)
388 {
389 CVODE_SOLVER *cvodeData = (CVODE_SOLVER *)userData;
390 ✗ DATA *data = (DATA *)(((CVODE_USERDATA *)cvodeData->simData)->data);
391 ✗ threadData_t *threadData = (threadData_t *)(((CVODE_USERDATA *)((CVODE_SOLVER *)userData)->simData)->threadData);
392
393 int saveJumpState;
394
395 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "### eval rootsFunctionCVODE ###");
396
397 ✗ if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
398 {
399 ✗ setContext(data, time, CONTEXT_EVENTS);
400 }
401
402 /* TODO: re-scale cvodeData->y to evaluate the equations */
403
404 ✗ saveJumpState = threadData->currentErrorStage;
405 ✗ threadData->currentErrorStage = ERROR_EVENTSEARCH;
406
407 ✗ data->localData[0]->timeValue = time;
408
409 /* Read input vars (exclude from timer) */
410 ✗ if (measure_time_flag)
411 ✗ rt_accumulate(SIM_TIMER_SOLVER);
412 #ifndef OMC_FMI_RUNTIME
413 ✗ externalInputUpdate(data);
414 ✗ data->callback->input_function(data, threadData);
415 #endif
416 /* eval needed equations (exclude from timer) */
417 ✗ data->callback->function_ZeroCrossingsEquations(data, threadData);
418 ✗ data->callback->function_ZeroCrossings(data, threadData, gout);
419 ✗ if (measure_time_flag)
420 ✗ rt_tick(SIM_TIMER_SOLVER);
421
422 ✗ threadData->currentErrorStage = saveJumpState;
423
424 /* TODO: scale data again */
425
426 ✗ if (data->simulationInfo->currentContext == CONTEXT_EVENTS)
427 {
428 ✗ unsetContext(data);
429 }
430 ✗ messageClose(OMC_LOG_SOLVER_V);
431 ✗ if (measure_time_flag)
432 ✗ rt_tick(SIM_TIMER_SOLVER);
433
434 ✗ return 0;
435 }
436
437 /**
438 * @brief Get settings for CVODE from user flags.
439 *
440 * If the user didn't provide any flags following settings will be chosen:
441 * config->lmm = CV_BDF
442 * config->iter = CV_ITER_NEWTON
443 *
444 * @param cvodeData CVODE solver data struckt
445 * @param threadData Thread data for error handling
446 */
447 ✗ void cvodeGetConfig(CVODE_CONFIG *config, threadData_t *threadData, sunbooleantype isFMI)
448 {
449 /* Variables */
450 int i;
451
452 /* ### Options for CVodeCreate ### */
453
454 /* Set linear multistep method */
455 ✗ if (omc_flag[FLAG_CVODE_LMM])
456 {
457 ✗ if (strcmp((const char *)omc_flagValue[FLAG_CVODE_LMM], CVODE_LMM_NAME[CV_ADAMS]) == 0)
458 {
459 ✗ config->lmm = CV_ADAMS;
460 }
461 ✗ else if (strcmp((const char *)omc_flagValue[FLAG_CVODE_LMM], CVODE_LMM_NAME[CV_BDF]) == 0)
462 {
463 ✗ config->lmm = CV_BDF;
464 }
465 else
466 {
467 ✗ if (OMC_ACTIVE_WARNING_STREAM(OMC_LOG_SOLVER))
468 {
469 ✗ warningStreamPrint(OMC_LOG_SOLVER, 1, "Unrecognized linear multistep method %s for CVODE, current options are:", (const char *)omc_flagValue[FLAG_CVODE_LMM]);
470 ✗ for (i = 1; i <= CVODE_LMM_MAX; ++i)
471 {
472 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "%s [%s]", CVODE_LMM_NAME[i], CVODE_LMM_DESC[i]);
473 }
474 ✗ messageClose(OMC_LOG_SOLVER);
475 }
476 ✗ throwStreamPrint(threadData, "Unrecognized linear multistep method %s for CVODE.", (const char *)omc_flagValue[FLAG_CVODE_LMM]);
477 }
478 }
479 else /* No user provided flag */
480 {
481 ✗ config->lmm = CV_BDF;
482 }
483
484 /* Set nonlinear solver iteration type */
485 ✗ if (omc_flag[FLAG_CVODE_ITER])
486 {
487 ✗ if (strcmp((const char *)omc_flagValue[FLAG_CVODE_ITER], CVODE_ITER_NAME[CV_ITER_FIXED_POINT]) == 0)
488 {
489 ✗ config->iter = CV_ITER_FIXED_POINT;
490 }
491 ✗ else if (strcmp((const char *)omc_flagValue[FLAG_CVODE_ITER], CVODE_ITER_NAME[CV_ITER_NEWTON]) == 0)
492 {
493 ✗ config->iter = CV_ITER_NEWTON;
494 }
495 else
496 {
497 ✗ if (OMC_ACTIVE_WARNING_STREAM(OMC_LOG_SOLVER))
498 {
499 ✗ warningStreamPrint(OMC_LOG_SOLVER, 1, "Unrecognized type of nonlinear solver iteration %s for CVODE, current options are:", (const char *)omc_flagValue[FLAG_CVODE_ITER]);
500 ✗ for (i = 1; i <= CVODE_ITER_MAX; ++i)
501 {
502 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "%s [%s]", CVODE_ITER_NAME[i], CVODE_ITER_DESC[i]);
503 }
504 ✗ messageClose(OMC_LOG_SOLVER);
505 }
506 ✗ throwStreamPrint(threadData, "Unrecognized type of nonlinear solver iteration %s for CVODE.", (const char *)omc_flagValue[FLAG_CVODE_ITER]);
507 }
508 }
509 else /* No user provided flag */
510 {
511 ✗ if (config->lmm == CV_ADAMS)
512 {
513 ✗ config->iter = CV_ITER_FIXED_POINT;
514 }
515 else
516 {
517 ✗ config->iter = CV_ITER_NEWTON;
518 }
519 }
520
521 /* Check for compability of lmn and iter */
522 ✗ if ((config->lmm == CV_ADAMS && config->iter != CV_ITER_FIXED_POINT) ||
523 ✗ (config->lmm == CV_BDF && config->iter != CV_ITER_NEWTON))
524 {
525 ✗ if (OMC_ACTIVE_WARNING_STREAM(OMC_LOG_SOLVER))
526 {
527 ✗ warningStreamPrint(OMC_LOG_SOLVER, 1, "Combination of %s and %s not recommended.", CVODE_LMM_NAME[config->lmm], CVODE_ITER_NAME[config->iter]);
528 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "Use simflags %s and %s to set.", FLAG_NAME[FLAG_CVODE_LMM], FLAG_NAME[FLAG_CVODE_ITER]);
529 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "Use (CV_BDF, CV_ITER_NEWTON) for stiff problems (Default) or");
530 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "Use (CV_ADAMS, CV_ITER_FIXED_POINT) for nonstiff problems.");
531 ✗ messageClose(OMC_LOG_SOLVER);
532 }
533 }
534 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE linear multistep method %s", CVODE_LMM_NAME[config->lmm]);
535 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE maximum integration order %s", CVODE_ITER_NAME[config->iter]);
536
537 /* if FLAG_NOEQUIDISTANT_GRID is set, choose ida step method */
538 ✗ if (omc_flag[FLAG_NOEQUIDISTANT_GRID])
539 {
540 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "Ignoring user supplied flag \"%s\", using equidistant time grid.", omc_flagValue[FLAG_NOEQUIDISTANT_GRID]);
541 }
542 ✗ config->internalSteps = FALSE; // TODO: Setting not used yet
543 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE use equidistant time grid %s", config->internalSteps ? "NO" : "YES");
544
545 /* Set jacobian method, see cvode_solver_initial */
546 #ifdef OMC_FMI_RUNTIME
547 if (omc_flag[FLAG_JACOBIAN])
548 {
549 warningStreamPrint(OMC_LOG_SOLVER, 0, "Ignoring user supplied flag \"%s\", using internal dense Jacobian of CVODE.", omc_flagValue[FLAG_JACOBIAN]);
550 }
551 #endif
552 ✗ config->jacobianMethod = INTERNALNUMJAC;
553
554 /* Maximum absolute step size */
555 /* TODO: Check flags FLAG_NOEQUIDISTANT_OUT_FREQ, FLAG_NOEQUIDISTANT_OUT_TIME */
556 ✗ config->maxStepSize = 0.0; /* default value a.k.a. no maximum step size */
557
558 /* Initial step size */
559 ✗ if (omc_flag[FLAG_INITIAL_STEP_SIZE])
560 {
561 ✗ config->initStepSize = atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]);
562 ✗ assertStreamPrint(threadData, config->initStepSize >= DASSL_STEP_EPS, "Selected initial step size %e is too small.", config->initStepSize);
563 }
564 else
565 {
566 ✗ config->initStepSize = 0.0; /* use default */
567 }
568
569 /* Maximum integration order */
570 ✗ if (omc_flag[FLAG_MAX_ORDER])
571 {
572 ✗ config->maxOrderLinearMultistep = atoi(omc_flagValue[FLAG_MAX_ORDER]);
573 }
574 ✗ else if (config->lmm == CV_ADAMS)
575 {
576 ✗ config->maxOrderLinearMultistep = 12 /* From ADAMS_Q_MAX */;
577 }
578 ✗ else if (config->lmm == CV_BDF)
579 {
580 ✗ config->maxOrderLinearMultistep = 5 /* From BDF_Q_MAX */;
581 }
582 else
583 {
584 ✗ throwStreamPrint(threadData, "Unrecognized linear multistep method. Can't set maximum order.");
585 }
586 /* Maximum number of nonlinear convergence failures */
587 /* TODO: Add a user flag */
588 ✗ config->maxConvFailPerStep = 10;
589
590 /* Use BDF stability limit detection */
591 /* TODO: Add a user flag */
592 ✗ if (config->lmm == CV_BDF)
593 {
594 ✗ config->BDFStabDetect = TRUE;
595 }
596 else
597 {
598 ✗ config->BDFStabDetect = FALSE;
599 }
600
601 ✗ if(omc_flag[FLAG_NO_ROOTFINDING] || isFMI)
602 {
603 ✗ config->solverRootFinding = FALSE;
604 }
605 else
606 {
607 ✗ config->solverRootFinding = TRUE;
608 }
609 ✗ }
610
611 /**
612 * @brief Read the states' nominal values into the absolute tolerances.
613 *
614 * Re-read by updateSolverNominals once initialization has computed the nominals
615 * that are parameter expressions.
616 *
617 * @param data Runtime data struct
618 * @param threadData Thread data for error handling
619 * @param cvodeData CVODE solver data struct with absoluteTolerance allocated.
620 * @return int Return 0 on success.
621 */
622 ✗ int cvode_solver_setNominals(DATA *data, threadData_t *threadData, CVODE_SOLVER *cvodeData)
623 {
624 int flag;
625 long int i;
626 ✗ double *abstol = N_VGetArrayPointer_Serial(cvodeData->absoluteTolerance);
627
628 ✗ for (i = 0; i < cvodeData->N; ++i)
629 {
630 ✗ const modelica_real nominal = getNominalFromScalarIdx(data->simulationInfo, data->modelData, VAR_KIND_STATE, i);
631 ✗ abstol[i] = fmax(fabs(nominal), 1e-32) * data->simulationInfo->tolerance;
632 }
633 ✗ flag = CVodeSVtolerances(cvodeData->cvode_mem, data->simulationInfo->tolerance, cvodeData->absoluteTolerance);
634 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSVtolerances");
635
636 ✗ return 0;
637 }
638
639 /**
640 * @brief Allocate memory, initialize and set configurations for CVODE solver
641 *
642 * @param data Runtime data struct
643 * @param threadData Thread data for error handling
644 * @param solverInfo Information about main solver. Unused at the moment.
645 * @param cvodeData CVODE solver data struct.
646 * @return int Return 0 on success.
647 */
648 ✗ int cvode_solver_initial(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, CVODE_SOLVER *cvodeData, int isFMI)
649 {
650 /* Variables */
651 int flag;
652 int i;
653 double *abstol_tmp;
654 #ifndef OMC_FMI_RUNTIME
655 const SPARSE_PATTERN *cscPattern;
656 #else
657 JACOBIAN *jacobian;
658 #endif
659
660 /* Log cvode_initial */
661 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, "### Start initialize of CVODE solver ###");
662
663 /* Set simData */
664 ✗ cvodeData->simData = (CVODE_USERDATA *)malloc(sizeof(CVODE_USERDATA));
665 ✗ cvodeData->simData->data = data;
666 ✗ cvodeData->simData->threadData = threadData;
667
668 ✗ cvodeData->isInitialized = FALSE;
669
670 /* Get CVODE settings from user flags */
671 ✗ cvodeGetConfig(&(cvodeData->config), threadData, isFMI);
672
673 /* Create the SUNDIALS context every other SUNDIALS object is created with */
674 ✗ flag = SUNContext_Create(SUN_COMM_NULL, &cvodeData->sunctx);
675 ✗ assertStreamPrint(threadData, flag == SUN_SUCCESS, "SUNDIALS_ERROR: SUNContext_Create failed.");
676 ✗ sundialsSilenceLogger(cvodeData->sunctx);
677
678 /* Set error handler */
679 ✗ flag = SUNContext_PushErrHandler(cvodeData->sunctx, sundialsErrorHandlerFunction, cvodeData);
680 ✗ assertStreamPrint(threadData, flag == SUN_SUCCESS, "SUNDIALS_ERROR: SUNContext_PushErrHandler failed.");
681
682 /* Initialize states */
683 ✗ cvodeData->N = (long int)data->modelData->nStates;
684 ✗ cvodeData->y = N_VMake_Serial(cvodeData->N, (sunrealtype *)data->localData[0]->realVars, cvodeData->sunctx);
685 ✗ assertStreamPrint(threadData, NULL != cvodeData->y, "SUNDIALS_ERROR: N_VMake_Serial failed - returned NULL pointer.");
686
687 /* Allocate CVODE memory block */
688 ✗ cvodeData->cvode_mem = CVodeCreate(cvodeData->config.lmm, cvodeData->sunctx);
689 ✗ assertStreamPrint(threadData, NULL != cvodeData->cvode_mem, "CVODE_ERROR: CVodeCreate failed - returned NULL pointer.");
690
691 ✗ if (measure_time_flag)
692 {
693 ✗ rt_tick(SIM_TIMER_SOLVER); /* Maybe use SIM_TIMER_OVERHEAD instead? */
694 }
695
696 /* Provide problem and solution specifications, allocate internal memory and initializes CVODE */
697 ✗ flag = CVodeInit(cvodeData->cvode_mem,
698 cvodeRightHandSideODEFunction,
699 ✗ data->simulationInfo->startTime,
700 cvodeData->y);
701 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeInit");
702
703 /* Set CVODE relative and absolute error tolerances */
704 ✗ abstol_tmp = (double *)calloc(cvodeData->N, sizeof(double)); /* Is freed with `free(NV_DATA_S(cvodeData->absoluteTolerance));` */
705 ✗ assertStreamPrint(threadData, abstol_tmp != NULL, "Out of memory.");
706 ✗ cvodeData->absoluteTolerance = N_VMake_Serial(cvodeData->N, abstol_tmp, cvodeData->sunctx);
707 ✗ assertStreamPrint(threadData, NULL != cvodeData->absoluteTolerance, "SUNDIALS_ERROR: N_VMake_Serial failed - returned NULL pointer.");
708 ✗ cvode_solver_setNominals(data, threadData, cvodeData);
709 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE Using relative error tolerance %e", data->simulationInfo->tolerance);
710
711 /* Provide cvodeData as user data */
712 ✗ flag = CVodeSetUserData(cvodeData->cvode_mem, cvodeData);
713 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetUserData");
714
715 /* Set linear solver used by CVODE: KLU over the ODE Jacobian's sparsity pattern,
716 * dense with CVODE's internal difference quotient without one. */
717 ✗ cvodeData->y_linSol = N_VNew_Serial(cvodeData->N, cvodeData->sunctx);
718 #ifndef OMC_FMI_RUNTIME
719 ✗ cvodeData->config.jacobianMethod = getRequestedJacobianMethod(threadData);
720 ✗ cscPattern = getJacobianCscPattern(initSymbolicOdeJacobian(data, threadData, &cvodeData->config.jacobianMethod, FALSE));
721 ✗ if (cvodeData->config.jacobianMethod == SYMJAC)
722 {
723 ✗ cvodeData->config.jacobianMethod = COLOREDSYMJAC;
724 }
725 ✗ else if (cvodeData->config.jacobianMethod == NUMJAC)
726 {
727 ✗ cvodeData->config.jacobianMethod = COLOREDNUMJAC;
728 }
729 ✗ if (cscPattern == NULL)
730 {
731 ✗ cvodeData->config.jacobianMethod = INTERNALNUMJAC;
732 }
733 #else
734 jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
735 data->callback->initialAnalyticJacobianA(data, threadData, jacobian);
736 #endif
737
738 ✗ if (cvodeData->config.jacobianMethod == INTERNALNUMJAC)
739 {
740 ✗ cvodeData->J = SUNDenseMatrix(cvodeData->N, cvodeData->N, cvodeData->sunctx);
741 ✗ cvodeData->linSol = SUNLinSol_Dense(cvodeData->y_linSol, cvodeData->J, cvodeData->sunctx);
742 ✗ assertStreamPrint(threadData, NULL != cvodeData->linSol, "##CVODE## SUNLinSol_Dense failed.");
743 ✗ flag = CVodeSetLinearSolver(cvodeData->cvode_mem, cvodeData->linSol, cvodeData->J);
744 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CVLS_FLAG, "CVodeSetLinearSolver");
745 ✗ flag = CVodeSetJacFn(cvodeData->cvode_mem, NULL);
746 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CVLS_FLAG, "CVodeSetJacFn");
747 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE Using dense internal linear solver SUNLinSol_Dense.");
748 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE Use internal dense numeric jacobian method.");
749 }
750 #ifndef OMC_FMI_RUNTIME
751 else
752 {
753 /* Room for the diagonal CVODE's I - gamma*J adds */
754 ✗ cvodeData->J = SUNSparseMatrix(cvodeData->N, cvodeData->N, cscPattern->nnz + cvodeData->N, SUN_CSC_MAT, cvodeData->sunctx);
755 ✗ cvodeData->linSol = SUNLinSol_KLU(cvodeData->y_linSol, cvodeData->J, cvodeData->sunctx);
756 ✗ assertStreamPrint(threadData, NULL != cvodeData->linSol, "##CVODE## SUNLinSol_KLU failed.");
757 ✗ flag = CVodeSetLinearSolver(cvodeData->cvode_mem, cvodeData->linSol, cvodeData->J);
758 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CVLS_FLAG, "CVodeSetLinearSolver");
759 ✗ flag = CVodeSetJacFn(cvodeData->cvode_mem, callSparseJacobian);
760 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CVLS_FLAG, "CVodeSetJacFn");
761 ✗ cvodeData->fProbe = N_VNew_Serial(cvodeData->N, cvodeData->sunctx);
762 ✗ cvodeData->ysave = (double *)malloc(cvodeData->N * sizeof(double));
763 ✗ cvodeData->delta_hh = (double *)malloc(cvodeData->N * sizeof(double));
764 ✗ assertStreamPrint(threadData, cvodeData->ysave != NULL && cvodeData->delta_hh != NULL, "Out of memory.");
765 ✗ cvodeData->jacNominalFactor = omc_flag[FLAG_JACOBIAN_NOMINAL_FACTOR]
766 ✗ ? atof(omc_flagValue[FLAG_JACOBIAN_NOMINAL_FACTOR]) : 1.0;
767 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE Using sparse linear solver SUNLinSol_KLU.");
768 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE Use sparse Jacobian method %s", JACOBIAN_METHOD_NAME[cvodeData->config.jacobianMethod]);
769 }
770 #endif
771
772 /* Set optional non-linear solver module */
773 ✗ switch (cvodeData->config.iter)
774 {
775 ✗ case CV_ITER_FIXED_POINT:
776 ✗ cvodeData->y_nonLinSol = N_VNew_Serial(cvodeData->N, cvodeData->sunctx);
777 ✗ cvodeData->nonLinSol = SUNNonlinSol_FixedPoint(cvodeData->y_nonLinSol, cvodeData->N /* Num acceleration vectors for Anderson's method, m <= dimension*/, cvodeData->sunctx);
778 ✗ assertStreamPrint(threadData, NULL != cvodeData->nonLinSol, "##CVODE## SUNNonlinSol_FixedPoint failed.");
779 ✗ flag = CVodeSetNonlinearSolver(cvodeData->cvode_mem, cvodeData->nonLinSol);
780 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetNonlinearSolver");
781 ✗ break;
782 ✗ case CV_ITER_NEWTON:
783 /* Default option, no allocation needed */
784 ✗ cvodeData->y_nonLinSol = NULL;
785 ✗ cvodeData->nonLinSol = NULL;
786 ✗ break;
787 ✗ case CV_ITER_MAX:
788 ✗ throwStreamPrint(threadData, "##CVODE## Non-linear solver method not set.");
789 ✗ default:
790 ✗ throwStreamPrint(threadData, "##CVODE## Unknown non-linear solver method %s.", CVODE_ITER_NAME[cvodeData->config.iter]);
791 }
792
793 /* Set root finding function */
794 ✗ if (cvodeData->config.solverRootFinding)
795 {
796 ✗ solverInfo->solverRootFinding = 1;
797 ✗ flag = CVodeRootInit(cvodeData->cvode_mem, data->modelData->nZeroCrossings, rootsFunctionCVODE);
798 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeRootInit");
799 }
800 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE uses internal root finding method %s", solverInfo->solverRootFinding ? "YES" : "NO");
801
802 /* ### Set optional settings ### */
803 /* Maximum absolute step size */
804 ✗ flag = CVodeSetMaxStep(cvodeData->cvode_mem, cvodeData->config.maxStepSize);
805 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetMaxStep");
806 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE maximum absolut step size %g", cvodeData->config.maxStepSize);
807
808 /* Initial step size */
809 ✗ flag = CVodeSetInitStep(cvodeData->cvode_mem, cvodeData->config.initStepSize);
810 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetInitStep");
811 ✗ if (cvodeData->config.initStepSize == 0)
812 {
813 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE initial step size is set automatically");
814 }
815 else
816 {
817 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE initial step size %g", cvodeData->config.initStepSize);
818 }
819
820 /* Maximum integration order */
821 ✗ flag = CVodeSetMaxOrd(cvodeData->cvode_mem, cvodeData->config.maxOrderLinearMultistep);
822 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetMaxOrd");
823 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE maximum integration order %d", cvodeData->config.maxOrderLinearMultistep);
824
825 /* Maximum number of nonlinear convergence failures */
826 ✗ flag = CVodeSetMaxConvFails(cvodeData->cvode_mem, cvodeData->config.maxConvFailPerStep);
827 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetMaxConvFails");
828 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE maximum number of nonlinear convergence failures permitted during one step %d", cvodeData->config.maxConvFailPerStep);
829
830 /* BDF stability limit detection */
831 ✗ flag = CVodeSetStabLimDet(cvodeData->cvode_mem, cvodeData->config.BDFStabDetect);
832 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetStabLimDet");
833 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "CVODE BDF stability limit detection algorithm %s", cvodeData->config.BDFStabDetect ? "ON" : "OFF");
834
835 /* TODO: Add stuff in cvodeGetConfig for this */
836 ✗ flag = CVodeSetMaxNonlinIters(cvodeData->cvode_mem, 5); /* Maximum number of iterations */
837 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetMaxNonlinIters");
838 ✗ flag = CVodeSetMaxErrTestFails(cvodeData->cvode_mem, 100); /* Maximum number of error test failures */
839 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetMaxErrTestFails");
840 ✗ flag = CVodeSetMaxNumSteps(cvodeData->cvode_mem, 1000); /* Maximum number of steps */
841 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetMaxNumSteps");
842
843 /* Log cvode_initial */
844 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, "### Finished initialize of CVODE solver successfully ###");
845
846 ✗ if (measure_time_flag)
847 {
848 ✗ rt_clear(SIM_TIMER_SOLVER); /* Initialization should not add to this timer... */
849 }
850
851 ✗ return 0;
852 }
853
854 /**
855 * @brief Reinitialize CVODE solver
856 * Provide required problem specifications and reinitialize CVODE.
857 * If scaling is used y will be scaled accordingly.
858 *
859 * @param data Runtime data struct.
860 * @param threadData Thread data for error handling.
861 * @param solverInfo Information about main solver. Unused at the moment.
862 * @param cvodeData CVODE solver data struckt.
863 * @return int Return 0 on success.
864 */
865 ✗ int cvode_solver_reinit(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, CVODE_SOLVER *cvodeData)
866 {
867 /* Variables */
868 int flag, i;
869
870 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Re-initialized CVODE Solver");
871
872 /* Calculate matrix for residual scaling */
873 /* TODO: Add scaling */
874
875 ✗ flag = CVodeReInit(cvodeData->cvode_mem,
876 solverInfo->currentTime,
877 cvodeData->y);
878 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeReInit");
879
880 /* Calculate matrix for residual scaling */
881 /* TODO: Add rescaling */
882
883 ✗ return 0;
884 }
885
886 /**
887 * @brief Deinitialize CVODE data
888 *
889 * @param cvodeData
890 * @return int Return 0 on success.
891 */
892 ✗ int cvode_solver_deinitial(CVODE_SOLVER *cvodeData)
893 {
894 /* Free work arrays */
895 ✗ N_VDestroy_Serial(cvodeData->y);
896 ✗ free(NV_DATA_S(cvodeData->absoluteTolerance));
897 ✗ N_VDestroy_Serial(cvodeData->absoluteTolerance);
898
899 /* Free linear solver data */
900 ✗ N_VDestroy_Serial(cvodeData->y_linSol);
901 ✗ SUNMatDestroy(cvodeData->J);
902 ✗ SUNLinSolFree(cvodeData->linSol);
903 #ifndef OMC_FMI_RUNTIME
904 ✗ if (cvodeData->fProbe)
905 {
906 ✗ N_VDestroy_Serial(cvodeData->fProbe);
907 }
908 ✗ free(cvodeData->ysave);
909 ✗ free(cvodeData->delta_hh);
910 ✗ free(cvodeData->yStart);
911 ✗ free(cvodeData->fStart);
912 ✗ freeSymbolicOdeJacobian(cvodeData->simData->data);
913 #endif
914
915 /* Free non-linear solver data */
916 ✗ N_VDestroy_Serial(cvodeData->y_nonLinSol);
917 ✗ SUNNonlinSolFree(cvodeData->nonLinSol);
918
919 /* Free CVODE internal data */
920 ✗ CVodeFree(&cvodeData->cvode_mem);
921
922 ✗ SUNContext_Free(&cvodeData->sunctx);
923 ✗ free(cvodeData->simData);
924
925 #ifdef OMC_FMI_RUNTIME
926 cvodeData->freeSolverMemory(cvodeData);
927 #else
928 ✗ free(cvodeData);
929 #endif
930
931 /* Log cvode_solver_deinitial */
932 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "### Finished deinitialization of CVODE solver successfully ###");
933 ✗ return 0;
934 }
935
936 /**
937 * @brief Save solver statistics.
938 *
939 * If flag OMC_LOG_SOLVER_V is provided even more statistics will be collected.
940 *
941 * @param cvode_mem Pointer to CVODE memory block.
942 * @param solverStats Pointer to solverStats of solverInfo.
943 * @param threadData Thread data for error handling.
944 */
945 ✗ void cvode_save_statistics(void *cvode_mem, SOLVERSTATS *solverStats, threadData_t *threadData)
946 {
947 /* Variables */
948 long int tmp1, tmp2;
949 double dtmp;
950 int flag;
951
952 /* Get number of internal steps taken by CVODE */
953 ✗ tmp1 = 0;
954 ✗ flag = CVodeGetNumSteps(cvode_mem, &tmp1);
955 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeGetNumSteps");
956 ✗ solverStats->nStepsTaken = tmp1;
957
958 /* Get number of right hand side evaluations */
959 /* TODO: Is it okay to count number of rhs evaluations instead of residual evaluations? */
960 ✗ tmp1 = 0;
961 ✗ flag = CVodeGetNumRhsEvals(cvode_mem, &tmp1);
962 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeGetNumRhsEvals");
963 ✗ solverStats->nCallsODE = tmp1;
964
965 /* Get number of Jacobian evaluations */
966 ✗ tmp1 = 0;
967 ✗ flag = CVodeGetNumJacEvals(cvode_mem, &tmp1);
968 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CVLS_FLAG, "CVodeGetNumJacEvals");
969 ✗ solverStats->nCallsJacobian = tmp1;
970
971 /* Get number of local error test failures */
972 ✗ tmp1 = 0;
973 ✗ flag = CVodeGetNumErrTestFails(cvode_mem, &tmp1);
974 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeGetNumErrTestFails");
975 ✗ solverStats->nErrorTestFailures = tmp1;
976
977 /* Get number of nonlinear convergence failures */
978 ✗ tmp1 = 0;
979 ✗ flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &tmp1);
980 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeGetNumNonlinSolvConvFails");
981 ✗ solverStats->nConvergenceTestFailures = tmp1;
982
983 /* Get even more statistics */
984 ✗ if (omc_useStream[OMC_LOG_SOLVER_V])
985 {
986 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "### CVODEStats ###");
987 /* Nonlinear stats */
988 ✗ tmp1 = tmp2 = 0;
989 ✗ flag = CVodeGetNonlinSolvStats(cvode_mem, &tmp1, &tmp2);
990 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, " ## Cumulative number of nonlinear iterations performed: %ld", tmp1);
991 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, " ## Cumulative number of nonlinear convergence failures that have occurred: %ld", tmp2);
992
993 /* Others stats */
994 ✗ flag = CVodeGetTolScaleFactor(cvode_mem, &dtmp);
995 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, " ## Suggested scaling factor for user tolerances: %g", dtmp);
996
997 ✗ flag = CVodeGetNumLinSolvSetups(cvode_mem, &tmp1);
998 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, " ## Number of calls made to the linear solver setup function: %ld", tmp1);
999
1000 ✗ messageClose(OMC_LOG_SOLVER_V);
1001 }
1002 ✗ }
1003
1004 /**
1005 * @brief DASSL's and IDA's first step, min(0.001*tdist, 0.5/||der||) in the
1006 * weighted RMS norm.
1007 *
1008 * CVODE's own estimate differences f over the step, which after an event
1009 * straddles the discontinuity and comes out tiny: an ideal diode on the edge of
1010 * conducting then switches back within it, restart after restart.
1011 *
1012 * @param data Runtime data struct, holding the post-event states and derivatives.
1013 * @param cvodeData CVODE solver data struct.
1014 * @param tdist Distance to the next output point.
1015 * @return double Initial step size.
1016 */
1017 ✗ static double cvodeRestartStep(DATA *data, CVODE_SOLVER *cvodeData, double tdist)
1018 {
1019 ✗ const double *states = data->localData[0]->realVars;
1020 ✗ const double *ders = states + cvodeData->N;
1021 ✗ const double *abstol = N_VGetArrayPointer(cvodeData->absoluteTolerance);
1022 ✗ const double rtol = data->simulationInfo->tolerance;
1023 double sum = 0.0, w, norm, h;
1024 long int i;
1025
1026 ✗ for (i = 0; i < cvodeData->N; i++)
1027 {
1028 ✗ w = ders[i] / (rtol * fabs(states[i]) + abstol[i]);
1029 ✗ sum += w * w;
1030 }
1031 ✗ norm = sqrt(sum / fmax(cvodeData->N, 1));
1032 ✗ h = 0.001 * fabs(tdist);
1033 ✗ return norm * h > 0.5 ? 0.5 / norm : h;
1034 }
1035
1036 /**
1037 * @brief Main CVODE function to make a step.
1038 *
1039 * Integrates on current time interval.
1040 *
1041 * @param data Runtime data struct
1042 * @param threadData Thread data for error handling
1043 * @param cvodeData CVODE solver data struct.
1044 * @return int Returns 0 on success and return flag from CVode else.
1045 */
1046 ✗ int cvode_solver_step(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
1047 {
1048 /* Variabes */
1049 int saveJumpState;
1050 int flag;
1051 ✗ int retVal = 0;
1052 ✗ int finished = FALSE;
1053 double tout = 0;
1054
1055 CVODE_SOLVER *cvodeData;
1056 SIMULATION_DATA *simulationData;
1057 SIMULATION_INFO *simulationInfo;
1058
1059 /* Measure time */
1060 ✗ if (measure_time_flag)
1061 ✗ rt_tick(SIM_TIMER_SOLVER);
1062
1063 /* Access data */
1064 ✗ cvodeData = (CVODE_SOLVER *)solverInfo->solverData;
1065 ✗ simulationData = data->localData[0];
1066 ✗ simulationInfo = data->simulationInfo;
1067
1068 /* Set work array */
1069 ✗ N_VSetArrayPointer(data->localData[0]->realVars, cvodeData->y);
1070
1071 /* Reinitialize after event or at first call to cvode_solver_step() */
1072 ✗ if (solverInfo->didEventStep || !cvodeData->isInitialized)
1073 {
1074 #ifndef OMC_FMI_RUNTIME
1075 ✗ if (!cvodeData->isInitialized)
1076 {
1077 ✗ cvodeData->startTime = solverInfo->currentTime;
1078 ✗ cvodeData->yStart = (double *)malloc(cvodeData->N * sizeof(double));
1079 ✗ cvodeData->fStart = (double *)malloc(cvodeData->N * sizeof(double));
1080 ✗ assertStreamPrint(threadData, cvodeData->yStart != NULL && cvodeData->fStart != NULL, "Out of memory.");
1081 ✗ memcpy(cvodeData->yStart, simulationData->realVars, cvodeData->N * sizeof(double));
1082 ✗ memcpy(cvodeData->fStart, simulationData->realVars + cvodeData->N, cvodeData->N * sizeof(double));
1083 }
1084 #endif
1085 ✗ cvode_solver_reinit(data, threadData, solverInfo, cvodeData);
1086 ✗ cvodeData->isInitialized = TRUE;
1087 }
1088
1089 ✗ saveJumpState = threadData->currentErrorStage;
1090 ✗ threadData->currentErrorStage = ERROR_INTEGRATOR;
1091
1092 /* Try */
1093 #if !defined(OMC_EMCC)
1094 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1095 #endif
1096
1097 /* Check current step size */
1098 ✗ if (solverInfo->currentStepSize < DASSL_STEP_EPS)
1099 {
1100 ✗ throwStreamPrint(threadData, "##CVODE## Desired step to small!");
1101 infoStreamPrint(OMC_LOG_SOLVER, 0, "Interpolate constant");
1102
1103 /* Constant extrapolation */
1104 /* TODO: Interpolate linear solution */
1105 simulationData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
1106 if (measure_time_flag)
1107 rt_accumulate(SIM_TIMER_SOLVER);
1108 data->callback->functionODE(data, threadData);
1109 solverInfo->currentTime = simulationData->timeValue;
1110
1111 return 0;
1112 }
1113
1114 /* CVODE may step past tout and interpolates back to it, but never past the
1115 * next time event */
1116 ✗ tout = solverInfo->currentTime + solverInfo->currentStepSize;
1117 ✗ if (simulationInfo->nextSampleEvent < DBL_MAX)
1118 {
1119 ✗ CVodeSetStopTime(cvodeData->cvode_mem, fmax(simulationInfo->nextSampleEvent, tout));
1120 }
1121 else
1122 {
1123 ✗ CVodeClearStopTime(cvodeData->cvode_mem);
1124 }
1125
1126 ✗ if (solverInfo->didEventStep && !omc_flag[FLAG_INITIAL_STEP_SIZE])
1127 {
1128 ✗ flag = CVodeSetInitStep(cvodeData->cvode_mem, cvodeRestartStep(data, cvodeData, tout - solverInfo->currentTime));
1129 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_CV_FLAG, "CVodeSetInitStep");
1130 }
1131 /* Integrator loop */
1132 do
1133 {
1134 ✗ infoStreamPrint(OMC_LOG_SOLVER, 1, "##CVODE## new step from %.15g to %.15g", solverInfo->currentTime, tout);
1135
1136 /* Read input vars (exclude from timer) */
1137 ✗ if (measure_time_flag)
1138 ✗ rt_accumulate(SIM_TIMER_SOLVER);
1139 #ifndef OMC_FMI_RUNTIME
1140 ✗ externalInputUpdate(data);
1141 ✗ data->callback->input_function(data, threadData);
1142 #endif
1143 ✗ if (measure_time_flag)
1144 ✗ rt_tick(SIM_TIMER_SOLVER);
1145
1146 /* TODO: Add scaling */
1147
1148 /* Call CVODE integrator */
1149 ✗ flag = CVode(cvodeData->cvode_mem,
1150 tout,
1151 cvodeData->y,
1152 ✗ &(solverInfo->currentTime),
1153 CV_NORMAL);
1154
1155 /* Error handling */
1156 ✗ if ((flag == CV_SUCCESS || flag == CV_TSTOP_RETURN) && solverInfo->currentTime >= tout)
1157 {
1158 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##CVODE## step done to time = %.15g", solverInfo->currentTime);
1159 finished = TRUE;
1160 }
1161 ✗ else if (flag == CV_ROOT_RETURN)
1162 {
1163 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##CVODE## root found at time = %.15g", solverInfo->currentTime);
1164 finished = TRUE;
1165 }
1166 ✗ else if (flag == CV_TOO_MUCH_WORK)
1167 {
1168 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "##CVODE## has done too much work with small steps at time = %.15g", solverInfo->currentTime);
1169 }
1170 else
1171 {
1172 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "##CVODE## %d error occurred at time = %.15g", flag, solverInfo->currentTime);
1173 finished = TRUE;
1174 retVal = flag;
1175 }
1176
1177 /* Closing new step message */
1178 ✗ messageClose(OMC_LOG_SOLVER); // TODO make sure this is called even if something in between fails
1179
1180 /* Set time to current time */
1181 ✗ simulationData->timeValue = solverInfo->currentTime;
1182 ✗ } while (!finished && !OMC_ERROR_RAISED());
1183
1184 /* Catch */
1185 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
1186 #if !defined(OMC_EMCC)
1187 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1188 #endif
1189 ✗ threadData->currentErrorStage = saveJumpState;
1190
1191 /* If a state event occured no sample event needs to be activated */
1192 ✗ if (simulationInfo->sampleActivated && solverInfo->currentTime < simulationInfo->nextSampleEvent)
1193 {
1194 ✗ simulationInfo->sampleActivated = 0 /* false */;
1195 }
1196
1197 /* Save statistics */
1198 ✗ cvode_save_statistics(cvodeData->cvode_mem, &solverInfo->solverStatsTmp, threadData);
1199
1200 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##CVODE## Finished Integrator step.");
1201 /* Measure time */
1202 ✗ if (measure_time_flag)
1203 ✗ rt_accumulate(SIM_TIMER_SOLVER);
1204
1205 ✗ return retVal;
1206 }
1207
1208 #ifdef OMC_FMI_RUNTIME
1209
1210 /**
1211 * @brief Integration step with CVODE for fmi2DoStep
1212 *
1213 * @param comp Pointer to FMU component.
1214 * @param tNext Next desired time step for integrator to end.
1215 * @param states States vector.
1216 * @return int Returns 0 on success and -1 else.
1217 */
1218 int cvode_solver_fmi_step(ModelInstance *comp, double tNext, double* states)
1219 {
1220 DATA* data = comp->fmuData;
1221 threadData_t* threadData = comp->threadData;
1222 SOLVER_INFO* solverInfo = comp->solverInfo;
1223 /* Variables */
1224 int flag;
1225 int retVal = 0;
1226
1227 CVODE_SOLVER *cvodeData;
1228
1229 cvodeData = (CVODE_SOLVER*) solverInfo->solverData;
1230 solverInfo->currentTime = data->localData[0]->timeValue;
1231
1232 N_VSetArrayPointer(states, cvodeData->y);
1233 if (solverInfo->didEventStep || !cvodeData->isInitialized) // TODO Save if we have had an event
1234 {
1235 cvode_solver_reinit(data, threadData, solverInfo, cvodeData);
1236 cvodeData->isInitialized = TRUE;
1237 }
1238 flag = CVodeSetStopTime(cvodeData->cvode_mem, tNext);
1239 if (flag < 0) {
1240 FILTERED_LOG(comp, fmi2Fatal, LOG_STATUSFATAL, "fmi2DoStep: ##CVODE## CVodeSetStopTime failed with flag %i.", flag)
1241 return -1;
1242 }
1243 flag = CVode(cvodeData->cvode_mem,
1244 tNext,
1245 cvodeData->y,
1246 &(solverInfo->currentTime),
1247 CV_NORMAL);
1248 /* Error handling */
1249 if ((flag == CV_SUCCESS || flag == CV_TSTOP_RETURN) && solverInfo->currentTime >= tNext)
1250 {
1251 FILTERED_LOG(comp, fmi2OK, LOG_ALL, "fmi2DoStep:##CVODE## step done to time = %.15g.", comp->solverInfo->currentTime)
1252 }
1253 else
1254 {
1255 FILTERED_LOG(comp, fmi2Fatal, LOG_STATUSFATAL, "fmi2DoStep: ##CVODE## %d error occurred at time = %.15g.", flag, solverInfo->currentTime)
1256 return -1;
1257 }
1258
1259 return 0;
1260 }
1261
1262 #endif /* OMC_FMI_RUNTIME */
1263
1264 #else /* WITH_SUNDIALS */
1265
1266 int cvode_solver_initial(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, CVODE_SOLVER *cvodeData, int isFMI)
1267 {
1268 #ifdef OMC_FMI_RUNTIME
1269 printf("##CVODE## SUNDIALS not available in FMU. See OpenModelica command line flag \"--fmiFlags\" from \"omc --help\" on how to enable CVODE in FMUs.\n");
1270 return -1;
1271 #else
1272 throwStreamPrint(threadData, "##CVODE## SUNDIALS not available. Reconfigure omc with SUNDIALS.\n");
1273 #endif
1274 }
1275
1276 int cvode_solver_deinitial(CVODE_SOLVER *cvodeData)
1277 {
1278 #ifdef OMC_FMI_RUNTIME
1279 printf("##CVODE## SUNDIALS not available in FMU. See OpenModelica command line flag \"--fmiFlags\" from \"omc --help\" on how to enable CVODE in FMUs.\n");
1280 return -1;
1281 #else
1282 throwStreamPrint(NULL, "##CVODE## SUNDIALS not available. Reconfigure omc with SUNDIALS.\n");
1283 #endif
1284 }
1285
1286 int cvode_solver_step(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
1287 {
1288 #ifdef OMC_FMI_RUNTIME
1289 printf("##CVODE## SUNDIALS not available in FMU. See OpenModelica command line flag \"--fmiFlags\" from \"omc --help\" on how to enable CVODE in FMUs.\n");
1290 return -1;
1291 #else
1292 throwStreamPrint(threadData, "##CVODE## SUNDIALS not available. Reconfigure omc with SUNDIALS.\n");
1293 #endif
1294 }
1295
1296 #ifdef OMC_FMI_RUNTIME
1297 int cvode_solver_fmi_step(ModelInstance *comp, double tNext, double* states)
1298 {
1299 printf("##CVODE## SUNDIALS not available in FMU. See OpenModelica command line flag \"--fmiFlags\" from \"omc --help\" on how to enable CVODE in FMUs.\n");
1300 return -1;
1301 }
1302 #endif
1303
1304 #endif /* #ifdef WITH_SUNDIALS */
1305