Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 43.5% 203 / 0 / 467
Functions: 69.2% 9 / 0 / 13
Branches: 18.7% 46 / 0 / 246
OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.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 /*! \file solver_main.c
29 */
30
31 #include "omc_config.h"
32 #include "simulation/simulation_runtime.h"
33 #include "simulation/results/simulation_result.h"
34 #include "solver_main.h"
35 #include "openmodelica_func.h"
36 #include "initialization/initialization.h"
37 #include "nonlinearSystem.h"
38 #include "newtonIteration.h"
39 #include "cvode_solver.h"
40 #include "dassl.h"
41 #include "ida_solver.h"
42 #include "delay.h"
43 #include "events.h"
44 #include "util/varinfo.h"
45 #include "util/omc_strdup.h"
46 #include "model_help.h"
47 #include "simulation/solver/epsilon.h"
48 #include "simulation/solver/external_input.h"
49 #include "synchronous.h"
50 #include "linearSystem.h"
51 #include "sym_solver_ssc.h"
52 #include "gbode_main.h"
53 #include "gbode_util.h"
54 #if !defined(OMC_MINIMAL_RUNTIME)
55 #include "simulation/solver/embedded_server.h"
56 #include "simulation/solver/real_time_sync.h"
57 #endif
58 #include "simulation/simulation_input_xml.h"
59
60 #include "optimization/OptimizerInterface.h"
61
62 /*
63 * #include "dopri45.h"
64 */
65 #include "util/rtclock.h"
66 #include "util/omc_error.h"
67 #include "simulation/options.h"
68 #include <math.h>
69 #include <string.h>
70 #include <errno.h>
71 #include <float.h>
72
73 double** work_states;
74
75 typedef struct RK4_DATA
76 {
77 double** work_states;
78 int work_states_ndims;
79 const double *b;
80 const double *c;
81 double h;
82 }RK4_DATA;
83
84
85 static int euler_ex_step(DATA* data, SOLVER_INFO* solverInfo);
86 static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
87 static int sym_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
88
89 #ifdef OMC_HAVE_IPOPT
90 static int ipopt_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
91 #endif
92
93 static void writeOutputVars(char* names, DATA* data);
94
95 1 int solver_main_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
96 {
97 int retVal;
98
99
1/10
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 1 switch(solverInfo->solverMethod)
100 {
101 ✗ case S_EULER:
102 ✗ retVal = euler_ex_step(data, solverInfo);
103 ✗ if(omc_flag[FLAG_SOLVER_STEPS])
104 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
105 return retVal;
106 ✗ case S_RUNGEKUTTA:
107 ✗ retVal = rungekutta_step(data, threadData, solverInfo);
108 ✗ if(omc_flag[FLAG_SOLVER_STEPS])
109 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
110 return retVal;
111
112 #if !defined(OMC_MINIMAL_RUNTIME)
113 1 case S_DASSL:
114 1 retVal = dassl_step(data, threadData, solverInfo);
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(omc_flag[FLAG_SOLVER_STEPS])
116 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
117 return retVal;
118 #endif
119
120 #ifdef OMC_HAVE_IPOPT
121 ✗ case S_OPTIMIZATION:
122 ✗ if ((int)(data->modelData->nStates + data->modelData->nInputVars) > 0){
123 retVal = ipopt_step(data, threadData, solverInfo);
124 } else {
125 ✗ solverInfo->solverMethod = S_EULER;
126 ✗ retVal = euler_ex_step(data, solverInfo);
127 }
128 ✗ if(omc_flag[FLAG_SOLVER_STEPS])
129 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
130 return retVal;
131 #endif
132 #ifdef WITH_SUNDIALS
133 ✗ case S_IDA:
134 ✗ retVal = ida_solver_step(data, threadData, solverInfo);
135 ✗ if(omc_flag[FLAG_SOLVER_STEPS])
136 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
137 return retVal;
138 ✗ case S_CVODE:
139 ✗ retVal = cvode_solver_step(data, threadData, solverInfo);
140 ✗ if(omc_flag[FLAG_SOLVER_STEPS])
141 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
142 return retVal;
143 #endif
144 ✗ case S_SYM_SOLVER:
145 ✗ retVal = sym_solver_step(data, threadData, solverInfo);
146 ✗ if(omc_flag[FLAG_SOLVER_STEPS])
147 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
148 return retVal;
149 ✗ case S_SYM_SOLVER_SSC:
150 ✗ retVal = sym_solver_ssc_step(data, threadData, solverInfo);
151 ✗ if(omc_flag[FLAG_SOLVER_STEPS])
152 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
153 return retVal;
154 ✗ case S_GBODE:
155 ✗ retVal = gbode_main(data, threadData, solverInfo);
156 ✗ if(omc_flag[FLAG_SOLVER_STEPS])
157 ✗ data->simulationInfo->solverSteps = solverInfo->solverStats.nStepsTaken + solverInfo->solverStatsTmp.nStepsTaken;
158 return retVal;
159 ✗ default:
160 ✗ throwStreamPrint(threadData, "Unhandled case in solver_main_step.");
161 }
162 return 1;
163 }
164
165 /*! \fn initializeSolverData(DATA* data, SOLVER_INFO* solverInfo)
166 *
167 * \param [ref] [data]
168 * \param [ref] [solverInfo]
169 *
170 * This function initializes solverInfo.
171 */
172 1 int initializeSolverData(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
173 {
174 int retValue = 0;
175 int i;
176
177 1 SIMULATION_INFO *simInfo = data->simulationInfo;
178
179 /* initial solverInfo */
180 1 solverInfo->currentTime = simInfo->startTime;
181 1 solverInfo->currentStepSize = simInfo->stepSize;
182 1 solverInfo->laststep = 0;
183 1 solverInfo->solverRootFinding = 0;
184 1 solverInfo->solverNoEquidistantGrid = omc_flag[FLAG_NOEQUIDISTANT_GRID];
185 1 solverInfo->lastdesiredStep = solverInfo->currentTime + solverInfo->currentStepSize;
186 1 solverInfo->eventLst = allocList(eventListAlloc, eventListFree, eventListCopy);
187 1 solverInfo->didEventStep = 0;
188 1 solverInfo->stateEvents = 0;
189 1 solverInfo->sampleEvents = 0;
190 1 resetSolverStats(&solverInfo->solverStats);
191 1 resetSolverStats(&solverInfo->solverStatsTmp);
192
193
1/9
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1 switch (solverInfo->solverMethod)
194 {
195 case S_SYM_SOLVER:
196 case S_EULER:
197 case S_QSS: break;
198 ✗ case S_SYM_SOLVER_SSC:
199 {
200 ✗ allocateSymSolverSsc(solverInfo, data->modelData->nStates);
201 ✗ break;
202 }
203 ✗ case S_GBODE:
204 {
205 ✗ if (gbode_allocateData(data, threadData, solverInfo) != 0) {
206 ✗ throwStreamPrint(threadData, "Failed to allocate memory for generic multigrid solver.");
207 }
208 break;
209 }
210 ✗ case S_RUNGEKUTTA:
211 {
212 /* Allocate RK work arrays */
213
214 static const int rungekutta_s = 4;
215 static const double rungekutta_b[4] = { 1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0 };
216 static const double rungekutta_c[4] = { 0.0, 0.5, 0.5, 1.0 };
217
218 static const int heun_s = 2;
219 static const double heun_b[2] = { 1.0 / 2.0, 1.0 / 2.0 };
220 static const double heun_c[2] = { 0.0, 1.0 };
221
222 ✗ RK4_DATA* rungeData = (RK4_DATA*) malloc(sizeof(RK4_DATA));
223
224 ✗ rungeData->work_states_ndims = rungekutta_s;
225 ✗ rungeData->b = rungekutta_b;
226 ✗ rungeData->c = rungekutta_c;
227
228 ✗ rungeData->work_states = (double**) malloc((rungeData->work_states_ndims + 1) * sizeof(double*));
229 ✗ for (i = 0; i < rungeData->work_states_ndims + 1; i++) {
230 ✗ rungeData->work_states[i] = (double*) calloc(data->modelData->nStates, sizeof(double));
231 }
232 ✗ solverInfo->solverData = rungeData;
233 ✗ break;
234 }
235 #if !defined(OMC_MINIMAL_RUNTIME)
236 1 case S_DASSL:
237 {
238 /* Initial DASSL solver */
239 1 DASSL_DATA* dasslData = (DASSL_DATA*) malloc(sizeof(DASSL_DATA));
240 1 retValue = dassl_initial(data, threadData, solverInfo, dasslData);
241 1 solverInfo->solverData = dasslData;
242 1 break;
243 }
244 #endif
245 #ifdef OMC_HAVE_IPOPT
246 ✗ case S_OPTIMIZATION:
247 {
248 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Initializing optimizer");
249 /* solverInfo->solverData = malloc(sizeof(OptData)); */
250 ✗ break;
251 }
252 #endif
253 #ifdef WITH_SUNDIALS
254 ✗ case S_IDA:
255 {
256 IDA_SOLVER* idaData = NULL;
257 /* Allocate ida working data */
258 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Initializing IDA DAE Solver");
259 ✗ idaData = (IDA_SOLVER*) malloc(sizeof(IDA_SOLVER));
260 ✗ retValue = ida_solver_initial(data, threadData, solverInfo, idaData);
261 ✗ solverInfo->solverData = idaData;
262 ✗ break;
263 }
264 ✗ case S_CVODE:
265 {
266 CVODE_SOLVER* cvodeData = NULL;
267 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Initializing CVODE ODE Solver");
268 ✗ cvodeData = (CVODE_SOLVER*) calloc(1, sizeof(CVODE_SOLVER));
269 ✗ assertStreamPrint(threadData, cvodeData != NULL, "Out of memory");
270 ✗ retValue = cvode_solver_initial(data, threadData, solverInfo, cvodeData, 0 /* not FMI */);
271 ✗ solverInfo->solverData = cvodeData;
272 ✗ break;
273 }
274 #endif
275 ✗ default:
276 ✗ errorStreamPrint(OMC_LOG_SOLVER, 0, "Solver %s disabled on this configuration", SOLVER_METHOD_NAME[solverInfo->solverMethod]);
277 ✗ return 1;
278 }
279
280 return retValue;
281 }
282
283 /*! \fn updateSolverNominals(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
284 *
285 * \param [ref] [data]
286 * \param [ref] [threadData]
287 * \param [ref] [solverInfo]
288 *
289 * Re-read the states' nominal (and, for gbode, min and max) attributes. A
290 * nominal that is a parameter expression is only computed by
291 * updateBoundVariableAttributes, inside initializeModel, which runs after
292 * initializeSolverData because DAE mode needs the solver during initialization;
293 * until then the solver holds the modelDescription default of 1.0.
294 */
295 1 int updateSolverNominals(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
296 {
297
1/5
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1 switch (solverInfo->solverMethod)
298 {
299 ✗ case S_GBODE:
300 ✗ gbode_setVarAttributes(data, solverInfo->solverData);
301 ✗ break;
302 #if !defined(OMC_MINIMAL_RUNTIME)
303 1 case S_DASSL:
304 1 dassl_setNominals(data, solverInfo->solverData);
305 1 break;
306 #endif
307 #ifdef WITH_SUNDIALS
308 ✗ case S_IDA:
309 ✗ return ida_solver_setNominals(data, threadData, solverInfo->solverData);
310 ✗ case S_CVODE:
311 ✗ return cvode_solver_setNominals(data, threadData, solverInfo->solverData);
312 #endif
313 default:
314 break;
315 }
316
317 return 0;
318 }
319
320 /*! \fn freeSolver(DATA* data, SOLVER_INFO* solverInfo)
321 *
322 * \param [ref] [data]
323 * \param [ref] [solverInfo]
324 *
325 * This function frees solverInfo.
326 */
327 1 int freeSolverData(DATA* data, SOLVER_INFO* solverInfo)
328 {
329 int retValue = 0;
330 int i;
331
332 1 freeList(solverInfo->eventLst);
333 /* deintialize solver related workspace */
334
1/8
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1 switch (solverInfo->solverMethod)
335 {
336 case S_EULER:
337 case S_SYM_SOLVER:
338 case S_QSS: break;
339 ✗ case S_SYM_SOLVER_SSC:
340 ✗ freeSymSolverSsc(solverInfo);
341 ✗ break;
342 case S_RUNGEKUTTA:
343 /* free RK work arrays */
344 ✗ for(i = 0; i < ((RK4_DATA*)(solverInfo->solverData))->work_states_ndims + 1; i++) {
345 ✗ free(((RK4_DATA*)(solverInfo->solverData))->work_states[i]);
346 }
347 ✗ free(((RK4_DATA*)(solverInfo->solverData))->work_states);
348 ✗ free((RK4_DATA*)solverInfo->solverData);
349 ✗ break;
350 ✗ case S_GBODE:
351 ✗ gbode_freeData(data, solverInfo->solverData);
352 ✗ solverInfo->solverData = NULL;
353 ✗ break;
354 #if !defined(OMC_MINIMAL_RUNTIME)
355 1 case S_DASSL:
356 /* De-Initial DASSL solver */
357 1 dassl_deinitial(data, solverInfo->solverData);
358 1 break;
359 #endif
360 #ifdef OMC_HAVE_IPOPT
361 case S_OPTIMIZATION:
362 /* free work arrays */
363 /*destroyIpopt(solverInfo);*/
364 break;
365 #endif
366 #ifdef WITH_SUNDIALS
367 ✗ case S_IDA:
368 /* free work arrays */
369 ✗ ida_solver_deinitial(solverInfo->solverData);
370 ✗ break;
371 ✗ case S_CVODE:
372 /* free work arrays */
373 ✗ cvode_solver_deinitial(solverInfo->solverData);
374 ✗ break;
375 #endif
376 ✗ default:
377 ✗ throwStreamPrint(NULL, "Unknown solver %u encountered. Possibly leaking memory!", solverInfo->solverMethod);
378 }
379
380 1 return retValue;
381 }
382
383
384 /*! \fn initializeModel(DATA* data, const char* init_initMethod,
385 * const char* init_file, double init_time)
386 *
387 * \param [ref] [data]
388 * \param [in] [pInitMethod] user defined initialization method
389 * \param [in] [pInitFile] extra argument for initialization-method "file"
390 * \param [in] [initTime] extra argument for initialization-method "file"
391 *
392 * This function starts the initialization process of the model.
393 */
394 1 int initializeModel(DATA* data, threadData_t *threadData, const char* init_initMethod,
395 const char* init_file, double init_time)
396 {
397 1 int retValue = 0;
398 int usedLocal = 0;
399
400 1 SIMULATION_INFO *simInfo = data->simulationInfo;
401
402
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(measure_time_flag)
403 {
404 1 rt_accumulate(SIM_TIMER_PREINIT);
405 1 rt_tick(SIM_TIMER_INIT);
406 }
407
408 1 copyStartValuestoInitValues(data);
409
410 1 data->localData[0]->timeValue = simInfo->startTime;
411
412 /* read input vars */
413 1 data->callback->input_function_init(data, threadData);
414 1 externalInputUpdate(data);
415 1 data->callback->input_function_updateStartValues(data, threadData);
416 1 data->callback->input_function(data, threadData);
417
418 1 threadData->currentErrorStage = ERROR_SIMULATION;
419 /* try */
420 {
421 int success = 0;
422 #if !defined(OMC_EMCC)
423
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 OMC_TRY_INTERNAL(simulationJumpBuffer)
424 #endif
425 /* A raised error is reported by the catch below. */
426
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1 if (initialization(data, threadData, init_initMethod, init_file, init_time) && !OMC_ERROR_RAISED())
427 {
428 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Error in initialization. Storing results and exiting.\nUse -lv=LOG_INIT -w for more information.");
429 ✗ simInfo->stopTime = simInfo->startTime;
430 retValue = -1;
431 }
432
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
 1 if (!retValue && !OMC_ERROR_RAISED())
433 {
434
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (data->simulationInfo->homotopySteps == 0) {
435 1 infoStreamPrint(OMC_LOG_SUCCESS, 0, "The initialization finished successfully without homotopy method.");
436 }
437 else {
438 ✗ usedLocal = data->callback->homotopyMethod == LOCAL_EQUIDISTANT_HOMOTOPY || data->callback->homotopyMethod == LOCAL_ADAPTIVE_HOMOTOPY ;
439 ✗ infoStreamPrint(OMC_LOG_SUCCESS, 0, "The initialization finished successfully with %d %shomotopy steps.", data->simulationInfo->homotopySteps, usedLocal? "local ":"");
440 }
441 }
442
443
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { success = 1; }
444 #if !defined(OMC_EMCC)
445 1 OMC_CATCH_INTERNAL(simulationJumpBuffer)
446 #endif
447
448
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (!success)
449 {
450 ✗ retValue = -1;
451 ✗ infoStreamPrint(OMC_LOG_ASSERT, 0, "simulation terminated by an assertion at initialization");
452 }
453 }
454
455 /* adrpo: write the parameter data in the file once again after bound parameters and initialization! */
456 1 sim_result.writeParameterData(&sim_result,data,threadData);
457 1 infoStreamPrint(OMC_LOG_SOLVER, 0, "Wrote parameters to the file after initialization (for output formats that support this)");
458
459 /* Initialization complete */
460
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (measure_time_flag) {
461 1 rt_accumulate(SIM_TIMER_INIT);
462 }
463
464 1 return retValue;
465 }
466
467
468 /*! \fn finishSimulation(DATA* data, SOLVER_INFO* solverInfo)
469 *
470 * \param [ref] [data]
471 * \param [ref] [solverInfo]
472 *
473 * This function performs the last step
474 * and outputs some statistics, this this simulation terminal step.
475 */
476 1 int finishSimulation(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, const char* outputVariablesAtEnd)
477 {
478 int retValue = 0;
479 int ui;
480 double t, total100;
481
482 1 SIMULATION_INFO *simInfo = data->simulationInfo;
483
484 /* Last step with terminal()=true */
485
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
 1 if (solverInfo->currentTime >= simInfo->stopTime && solverInfo->solverMethod != S_OPTIMIZATION) {
486 1 infoStreamPrint(OMC_LOG_EVENTS_V, 0, "terminal event at stop time %g", solverInfo->currentTime);
487 1 data->simulationInfo->terminal = 1;
488 1 updateDiscreteSystem(data, threadData);
489
490 /* prevent emit if noeventemit flag is used */
491
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (!(omc_flag[FLAG_NOEVENTEMIT])) {
492 1 sim_result.emit(&sim_result, data, threadData);
493 }
494
495 1 data->simulationInfo->terminal = 0;
496 }
497
498
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (0 != strcmp("ia", data->simulationInfo->outputFormat)) {
499
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (simInfo->simulationSuccess) {
500 ✗ communicateStatus("Finished", 1, solverInfo->currentTime, solverInfo->currentStepSize);
501 } else {
502 1 communicateStatus("Simulation aborted", (solverInfo->currentTime-simInfo->startTime)/(simInfo->stopTime-simInfo->startTime), solverInfo->currentTime, 0.0);
503 }
504 }
505
506 /* we have output variables in the command line -output a,b,c */
507
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(outputVariablesAtEnd)
508 {
509 ✗ writeOutputVars(omc_strdup(outputVariablesAtEnd), data);
510 }
511
512
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(OMC_ACTIVE_STREAM(OMC_LOG_STATS))
513 {
514 1 rt_accumulate(SIM_TIMER_TOTAL);
515
516 1 infoStreamPrint(OMC_LOG_STATS, 1, "### STATISTICS ###");
517
518 1 total100 = rt_accumulated(SIM_TIMER_TOTAL)/100.0;
519
520 1 infoStreamPrint(OMC_LOG_STATS, 1, "timer");
521 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs reading init.xml", rt_accumulated(SIM_TIMER_INIT_XML));
522 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs reading info.xml", rt_accumulated(SIM_TIMER_INFO_XML));
523 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [%5.1f%%] pre-initialization", rt_accumulated(SIM_TIMER_PREINIT), rt_accumulated(SIM_TIMER_PREINIT)/total100);
524 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [%5.1f%%] initialization", rt_accumulated(SIM_TIMER_INIT), rt_accumulated(SIM_TIMER_INIT)/total100);
525 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [%5.1f%%] steps", rt_accumulated(SIM_TIMER_STEP), rt_accumulated(SIM_TIMER_STEP)/total100);
526 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [%5.1f%%] solver (excl. callbacks)", rt_accumulated(SIM_TIMER_SOLVER), rt_accumulated(SIM_TIMER_SOLVER)/total100);
527 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [%5.1f%%] creating output-file", rt_accumulated(SIM_TIMER_OUTPUT), rt_accumulated(SIM_TIMER_OUTPUT)/total100);
528 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [%5.1f%%] event-handling", rt_accumulated(SIM_TIMER_EVENT), rt_accumulated(SIM_TIMER_EVENT)/total100);
529 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [%5.1f%%] overhead", rt_accumulated(SIM_TIMER_OVERHEAD), rt_accumulated(SIM_TIMER_OVERHEAD)/total100);
530
531 1 t = rt_accumulated(SIM_TIMER_TOTAL)-rt_accumulated(SIM_TIMER_OVERHEAD)-rt_accumulated(SIM_TIMER_EVENT)-rt_accumulated(SIM_TIMER_OUTPUT)-rt_accumulated(SIM_TIMER_STEP)-rt_accumulated(SIM_TIMER_INIT)-rt_accumulated(SIM_TIMER_PREINIT)-rt_accumulated(SIM_TIMER_SOLVER);
532
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 2 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [%5.1f%%] %s", t, t/total100, S_OPTIMIZATION == solverInfo->solverMethod ? "optimization" : "simulation");
533
534 1 infoStreamPrint(OMC_LOG_STATS, 0, "%12gs [100.0%%] total", rt_accumulated(SIM_TIMER_TOTAL));
535 1 messageClose(OMC_LOG_STATS);
536
537 1 infoStreamPrint(OMC_LOG_STATS, 1, "events");
538 1 infoStreamPrint(OMC_LOG_STATS, 0, "%5ld state events", solverInfo->stateEvents);
539 1 infoStreamPrint(OMC_LOG_STATS, 0, "%5ld time events", solverInfo->sampleEvents);
540 1 messageClose(OMC_LOG_STATS);
541
542
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(S_OPTIMIZATION == solverInfo->solverMethod || /* skip solver statistics for optimization */
543 S_QSS == solverInfo->solverMethod) /* skip also for qss, since not available*/
544 {
545 }
546 else
547 {
548 /* save stats before print */
549 1 addSolverStats(&(solverInfo->solverStats), &(solverInfo->solverStatsTmp));
550
551 1 infoStreamPrint(OMC_LOG_STATS, 1, "solver: %s", SOLVER_METHOD_NAME[solverInfo->solverMethod]);
552 1 infoStreamPrint(OMC_LOG_STATS, 0, "%5d steps taken", solverInfo->solverStats.nStepsTaken);
553 1 infoStreamPrint(OMC_LOG_STATS, 0, "%5d calls of functionODE", solverInfo->solverStats.nCallsODE);
554 1 infoStreamPrint(OMC_LOG_STATS, 0, "%5d evaluations of jacobian", solverInfo->solverStats.nCallsJacobian);
555 1 infoStreamPrint(OMC_LOG_STATS, 0, "%5d error test failures", solverInfo->solverStats.nErrorTestFailures);
556 1 infoStreamPrint(OMC_LOG_STATS, 0, "%5d convergence test failures", solverInfo->solverStats.nConvergenceTestFailures);
557 1 infoStreamPrint(OMC_LOG_STATS, 0, "%gs time of jacobian evaluation", rt_accumulated(SIM_TIMER_JACOBIAN));
558
559 1 messageClose(OMC_LOG_STATS);
560 }
561
562 1 infoStreamPrint(OMC_LOG_STATS_V, 1, "function calls");
563
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (compiledInDAEMode)
564 {
565 ✗ infoStreamPrint(OMC_LOG_STATS_V, 1, "%5u calls of functionDAE", rt_ncall(SIM_TIMER_DAE));
566 ✗ infoStreamPrint(OMC_LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_DAE), rt_accumulated(SIM_TIMER_DAE)/total100);
567 ✗ messageClose(OMC_LOG_STATS_V);
568 }
569
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (data->simulationInfo->callStatistics.functionODE) {
570 1 infoStreamPrint(OMC_LOG_STATS_V, 1, "%5ld calls of functionODE", data->simulationInfo->callStatistics.functionODE);
571 1 infoStreamPrint(OMC_LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_FUNCTION_ODE), rt_accumulated(SIM_TIMER_FUNCTION_ODE)/total100);
572 1 messageClose(OMC_LOG_STATS_V);
573 }
574
575
1/2
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
 1 if (rt_ncall(SIM_TIMER_RESIDUALS)) {
576 1 infoStreamPrint(OMC_LOG_STATS_V, 1, "%5d calls of functionODE_residual", rt_ncall(SIM_TIMER_RESIDUALS));
577 1 infoStreamPrint(OMC_LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_RESIDUALS), rt_accumulated(SIM_TIMER_RESIDUALS)/total100);
578 1 messageClose(OMC_LOG_STATS_V);
579 }
580
581
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (data->simulationInfo->callStatistics.functionAlgebraics) {
582 1 infoStreamPrint(OMC_LOG_STATS_V, 1, "%5ld calls of functionAlgebraics", data->simulationInfo->callStatistics.functionAlgebraics);
583 1 infoStreamPrint(OMC_LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_ALGEBRAICS), rt_accumulated(SIM_TIMER_ALGEBRAICS)/total100);
584 1 messageClose(OMC_LOG_STATS_V);
585 }
586
587 1 infoStreamPrint(OMC_LOG_STATS_V, 1, "%5d evaluations of jacobian", rt_ncall(SIM_TIMER_JACOBIAN));
588 1 infoStreamPrint(OMC_LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_JACOBIAN), rt_accumulated(SIM_TIMER_JACOBIAN)/total100);
589 1 messageClose(OMC_LOG_STATS_V);
590
591 1 infoStreamPrint(OMC_LOG_STATS_V, 0, "%5ld calls of updateDiscreteSystem", data->simulationInfo->callStatistics.updateDiscreteSystem);
592 1 infoStreamPrint(OMC_LOG_STATS_V, 0, "%5ld calls of functionZeroCrossingsEquations", data->simulationInfo->callStatistics.functionZeroCrossingsEquations);
593
594 1 infoStreamPrint(OMC_LOG_STATS_V, 1, "%5ld calls of functionZeroCrossings", data->simulationInfo->callStatistics.functionZeroCrossings);
595 1 infoStreamPrint(OMC_LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_ZC), rt_accumulated(SIM_TIMER_ZC)/total100);
596 1 messageClose(OMC_LOG_STATS_V);
597
598 1 messageClose(OMC_LOG_STATS_V); // closes section "function calls"
599
600 1 infoStreamPrint(OMC_LOG_STATS_V, 1, "linear systems");
601
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(ui=0; ui<data->modelData->nLinearSystems; ui++)
602 ✗ printLinearSystemSolvingStatistics(data, ui, OMC_LOG_STATS_V);
603 1 messageClose(OMC_LOG_STATS_V);
604
605 1 infoStreamPrint(OMC_LOG_STATS_V, 1, "non-linear systems");
606
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(ui=0; ui<data->modelData->nNonLinearSystems; ui++)
607 ✗ printNonLinearSystemSolvingStatistics(&data->simulationInfo->nonlinearSystemData[ui], OMC_LOG_STATS_V);
608 1 messageClose(OMC_LOG_STATS_V);
609
610 1 messageClose(OMC_LOG_STATS); // closes section "### STATISTICS ###"
611 1 rt_tick(SIM_TIMER_TOTAL);
612 }
613
614 1 return retValue;
615 }
616
617 /*! \fn solver_main
618 *
619 * \param [ref] [data]
620 * \param [in] [pInitMethod] user defined initialization method
621 * \param [in] [pOptiMethod] user defined optimization method
622 * \param [in] [pInitFile] extra argument for initialization-method "file"
623 * \param [in] [initTime] extra argument for initialization-method "file"
624 * \param [in] [solverID] selects the ode solver
625 * \param [in] [outputVariablesAtEnd] ???
626 *
627 * This is the main function of the solver, it performs the simulation.
628 */
629 1 int solver_main(DATA* data, threadData_t *threadData, const char* init_initMethod, const char* init_file,
630 double init_time, int solverID, const char* outputVariablesAtEnd, const char *argv_0)
631 {
632 int i;
633 /* Read after a longjmp into the catch below. */
634 1 volatile int retVal = 1, initSolverInfo = 0;
635 unsigned int ui;
636 SOLVER_INFO solverInfo;
637 1 SIMULATION_INFO *simInfo = data->simulationInfo;
638 void *dllHandle=NULL;
639
640 1 solverInfo.solverMethod = solverID;
641
642 /* do some solver specific checks */
643 switch(solverInfo.solverMethod)
644 {
645 #ifndef OMC_HAVE_IPOPT
646 case S_OPTIMIZATION:
647 warningStreamPrint(OMC_LOG_STDOUT, 0, "Ipopt is needed but not available.");
648 return 1;
649 #endif
650 default:
651 break;
652 }
653
654 /* first initialize the model then allocate SolverData memory
655 * due to be able to use the initialized values for the integrator
656 */
657 1 simInfo->useStopTime = 1;
658
659 /* Use minStepSize if stepSize is getting too small, but
660 * allow stepSize to be zero if startTime == stopTime.
661 */
662
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1 if ((simInfo->stepSize < simInfo->minStepSize) && (simInfo->stopTime > simInfo->startTime)){
663 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The step-size %g is too small. Adjust the step-size to %g.", simInfo->stepSize, simInfo->minStepSize);
664 ✗ simInfo->stepSize = simInfo->minStepSize;
665 ✗ simInfo->numSteps = round((simInfo->stopTime - simInfo->startTime)/simInfo->stepSize);
666 }
667 /* Check step size is not larger then stopTime-startTime, up to 6 decimals
668 * Ignored when linearizing model */
669
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
 1 if (!data->modelData->create_linearmodel && simInfo->stepSize > (simInfo->stopTime - simInfo->startTime + 1e-7)) {
670 ✗ warningStreamPrint(OMC_LOG_STDOUT, 1, "Integrator step size greater than length of experiment");
671 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "start time: %f, stop time: %f, integrator step size: %f",simInfo->startTime, simInfo->stopTime, simInfo->stepSize);
672 ✗ messageCloseWarning(OMC_LOG_STDOUT);
673 }
674 #if !defined(OMC_EMCC)
675
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 OMC_TRY_INTERNAL(simulationJumpBuffer)
676 #endif
677
678 /* initialize external input structure */
679 1 externalInputallocate(data);
680 /* set tolerance for ZeroCrossings */
681 1 setZCtol(fmin(data->simulationInfo->stepSize, data->simulationInfo->tolerance));
682 1 omc_alloc_interface.collect_a_little();
683
684 /* initialize solver data */
685 /* For the DAEmode we need to initialize solverData before the initialization,
686 * since the solver is used to obtain consistent values also via updateDiscreteSystem
687 */
688 1 retVal = initializeSolverData(data, threadData, &solverInfo);
689 1 initSolverInfo = 1;
690
691 /* initialize all parts of the model */
692
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (0 == retVal){
693 1 retVal = initializeModel(data, threadData, init_initMethod, init_file, init_time);
694 1 omc_alloc_interface.collect_a_little();
695 }
696
697 /* the nominal values are only final now */
698
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (0 == retVal){
699 1 retVal = updateSolverNominals(data, threadData, &solverInfo);
700 }
701
702 #if !defined(OMC_MINIMAL_RUNTIME)
703 1 dllHandle = embedded_server_load_functions(omc_flagValue[FLAG_EMBEDDED_SERVER]);
704 1 omc_real_time_sync_init(threadData, data);
705 int port = 4841;
706 /* If an embedded server is specified */
707
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (dllHandle != NULL) {
708 ✗ if (omc_flag[FLAG_EMBEDDED_SERVER_PORT]) {
709 ✗ port = atoi(omc_flagValue[FLAG_EMBEDDED_SERVER_PORT]);
710 /* In case of a bad conversion, don't spawn a server on port 0...*/
711 ✗ if (port == 0) {
712 port = 4841;
713 }
714 }
715 }
716 1 data->embeddedServerState = embedded_server_init(data, data->localData[0]->timeValue, solverInfo.currentStepSize, argv_0, omc_real_time_sync_update, port);
717 /* If an embedded server is specified */
718
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (dllHandle != NULL) {
719 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "The embedded server is initialized.");
720 }
721 1 wait_for_step(data->embeddedServerState);
722 #endif
723
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(0 == retVal) {
724 1 retVal = -1;
725 /* if the model has no time changing variables skip the main loop*/
726
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(data->modelData->nVariablesReal == 0 &&
727 ✗ data->modelData->nVariablesInteger == 0 &&
728 ✗ data->modelData->nVariablesBoolean == 0 &&
729 ✗ data->modelData->nVariablesString == 0 ) {
730 /* prevent emit if noeventemit flag is used */
731 ✗ if (!(omc_flag[FLAG_NOEVENTEMIT])) {
732 ✗ sim_result.emit(&sim_result, data, threadData);
733 }
734
735 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "The model has no time changing variables, no integration will be performed.");
736 ✗ solverInfo.currentTime = simInfo->stopTime;
737 ✗ data->localData[0]->timeValue = simInfo->stopTime;
738 ✗ overwriteOldSimulationData(data);
739 ✗ retVal = finishSimulation(data, threadData, &solverInfo, outputVariablesAtEnd);
740
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 } else if(S_QSS == solverInfo.solverMethod) {
741 /* starts the simulation main loop - special solvers */
742 ✗ sim_result.emit(&sim_result,data,threadData);
743
744 /* overwrite the whole ring-buffer with initialized values */
745 ✗ overwriteOldSimulationData(data);
746
747 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Start numerical integration (startTime: %g, stopTime: %g)", simInfo->startTime, simInfo->stopTime);
748 ✗ retVal = data->callback->performQSSSimulation(data, threadData, &solverInfo);
749 ✗ omc_alloc_interface.collect_a_little();
750
751 /* terminate the simulation */
752 ✗ finishSimulation(data, threadData, &solverInfo, outputVariablesAtEnd);
753 ✗ omc_alloc_interface.collect_a_little();
754 } else {
755 /* starts the simulation main loop - standard solver interface */
756
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(omc_flag[FLAG_SOLVER_STEPS])
757 ✗ data->simulationInfo->solverSteps = 0;
758
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(solverInfo.solverMethod != S_OPTIMIZATION) {
759 1 sim_result.emit(&sim_result,data,threadData);
760 }
761
762 /* overwrite the whole ring-buffer with initialized values */
763 1 overwriteOldSimulationData(data);
764
765 /* store all values for non-dassl event search */
766 1 storeOldValues(data);
767
768 1 infoStreamPrint(OMC_LOG_SOLVER, 0, "Start numerical solver from %g to %g", simInfo->startTime, simInfo->stopTime);
769 1 retVal = data->callback->performSimulation(data, threadData, &solverInfo);
770 1 omc_alloc_interface.collect_a_little();
771 /* terminate the simulation */
772 //if (solverInfo.solverMethod == S_SYM_SOLVER_SSC) data->callback->symbolicInlineSystems(data, threadData, 0, 2);
773 1 finishSimulation(data, threadData, &solverInfo, outputVariablesAtEnd);
774 1 omc_alloc_interface.collect_a_little();
775 }
776 }
777
778
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (data->real_time_sync.enabled) {
779 ✗ int tMaxLate=0;
780 ✗ const char *unit = prettyPrintNanoSec(data->real_time_sync.maxLate, &tMaxLate);
781 ✗ infoStreamPrint(OMC_LOG_RT, 0, "Maximum real-time latency was (positive=missed dealine, negative is slack): %d %s", tMaxLate, unit);
782 }
783 #if !defined(OMC_MINIMAL_RUNTIME)
784 1 embedded_server_deinit(data->embeddedServerState);
785 1 embedded_server_unload_functions(dllHandle);
786 #endif
787
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
788
789 #if !defined(OMC_EMCC)
790 1 OMC_CATCH_INTERNAL(simulationJumpBuffer)
791 #endif
792
793 /* free external input data */
794 1 externalInputFree(data);
795
796 /* free SolverInfo memory */
797
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (initSolverInfo)
798 {
799 1 freeSolverData(data, &solverInfo);
800 }
801
802
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (!retVal)
803 1 infoStreamPrint(OMC_LOG_SUCCESS, 0, "The simulation finished successfully.");
804
805 1 return retVal;
806 }
807
808 /*************************************** EULER_EXP *********************************/
809 ✗ static int euler_ex_step(DATA* data, SOLVER_INFO* solverInfo)
810 {
811 int i;
812 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
813 ✗ SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
814 ✗ modelica_real* stateDer = sDataOld->realVars + data->modelData->nStates;
815
816 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
817
818 ✗ solverInfo->currentTime = sDataOld->timeValue + solverInfo->currentStepSize;
819
820 ✗ for(i = 0; i < data->modelData->nStates; i++)
821 {
822 ✗ sData->realVars[i] = sDataOld->realVars[i] + stateDer[i] * solverInfo->currentStepSize;
823 }
824 ✗ sData->timeValue = solverInfo->currentTime;
825
826 /* save stats */
827 /* steps */
828 ✗ solverInfo->solverStatsTmp.nStepsTaken += 1;
829 /* function ODE evaluation is done directly after this function */
830 ✗ solverInfo->solverStatsTmp.nCallsODE += 1;
831
832 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
833
834 ✗ return 0;
835 }
836
837 /*************************************** SYM_SOLVER *********************************/
838 ✗ static int sym_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo){
839 int retVal = 0, i;
840
841 ✗ modelica_integer nStates = data->modelData->nStates;
842 ✗ SIMULATION_DATA *sData = data->localData[0];
843 ✗ SIMULATION_DATA *sDataOld = data->localData[1];
844 ✗ modelica_real* stateDer = sDataOld->realVars + data->modelData->nStates;
845
846 ✗ if (solverInfo->currentStepSize >= DASSL_STEP_EPS){
847 /* time */
848 ✗ solverInfo->currentTime = sDataOld->timeValue + solverInfo->currentStepSize;
849 ✗ sData->timeValue = solverInfo->currentTime;
850 /* update dt */
851 ✗ data->simulationInfo->inlineData->dt = solverInfo->currentStepSize;
852 /* copy old states to workspace */
853 ✗ memcpy(data->simulationInfo->inlineData->algOldVars, sDataOld->realVars, nStates * sizeof(double));
854 ✗ memcpy(sData->realVars, sDataOld->realVars, nStates * sizeof(double));
855
856 /* read input vars */
857 ✗ externalInputUpdate(data);
858 ✗ data->callback->input_function(data, threadData);
859 ✗ retVal = data->callback->symbolicInlineSystems(data, threadData);
860
861 ✗ if(retVal != 0){
862 return -1;
863 }
864
865
866 /* update der(x) */
867 ✗ for(i=0; i<nStates; ++i)
868 {
869 ✗ stateDer[i] = (sData->realVars[i]-data->simulationInfo->inlineData->algOldVars[i])/solverInfo->currentStepSize;
870 }
871
872 /* save stats */
873 /* steps */
874 ✗ solverInfo->solverStatsTmp.nStepsTaken += 1;
875 /* function ODE evaluation is done directly after this */
876 ✗ solverInfo->solverStatsTmp.nCallsODE += 1;
877 }
878 else
879 /* in case desired step size is too small */
880 {
881 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Desired step to small try next one");
882 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Interpolate linear");
883
884 /* explicit euler step*/
885 ✗ for(i = 0; i < nStates; i++)
886 {
887 ✗ sData->realVars[i] = sDataOld->realVars[i] + stateDer[i] * solverInfo->currentStepSize;
888 }
889 ✗ sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
890 ✗ solverInfo->currentTime = sData->timeValue;
891 }
892 return retVal;
893 }
894
895 /*************************************** RK4 ***********************************/
896 ✗ static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
897 {
898 ✗ RK4_DATA *rk = ((RK4_DATA*)(solverInfo->solverData));
899 ✗ double** k = rk->work_states;
900 double sum;
901 int i,j;
902 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
903 ✗ SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
904 ✗ modelica_real* stateDer = sData->realVars + data->modelData->nStates;
905 ✗ modelica_real* stateDerOld = sDataOld->realVars + data->modelData->nStates;
906
907 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
908
909 ✗ solverInfo->currentTime = sDataOld->timeValue + solverInfo->currentStepSize;
910
911 /* We calculate k[0] before returning from this function.
912 * We only want to calculate f() 4 times per call */
913 ✗ memcpy(k[0], stateDerOld, data->modelData->nStates*sizeof(modelica_real));
914
915 ✗ for (j = 1; j < rk->work_states_ndims; j++)
916 {
917 ✗ for(i = 0; i < data->modelData->nStates; i++)
918 {
919 ✗ sData->realVars[i] = sDataOld->realVars[i] + solverInfo->currentStepSize * rk->c[j] * k[j - 1][i];
920 }
921 ✗ sData->timeValue = sDataOld->timeValue + rk->c[j] * solverInfo->currentStepSize;
922 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
923 /* read input vars */
924 ✗ externalInputUpdate(data);
925 ✗ data->callback->input_function(data, threadData);
926 /* eval ode equations */
927 ✗ data->callback->functionODE(data, threadData);
928 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
929 ✗ memcpy(k[j], stateDer, data->modelData->nStates*sizeof(modelica_real));
930
931 }
932
933 ✗ for(i = 0; i < data->modelData->nStates; i++)
934 {
935 sum = 0;
936 ✗ for(j = 0; j < rk->work_states_ndims; j++)
937 {
938 ✗ sum = sum + rk->b[j] * k[j][i];
939 }
940 ✗ sData->realVars[i] = sDataOld->realVars[i] + solverInfo->currentStepSize * sum;
941 }
942 ✗ sData->timeValue = solverInfo->currentTime;
943
944 /* save stats */
945 /* steps */
946 ✗ solverInfo->solverStatsTmp.nStepsTaken += 1;
947 /* function ODE evaluation is done directly after this */
948 ✗ solverInfo->solverStatsTmp.nCallsODE += rk->work_states_ndims+1;
949 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
950
951 ✗ return 0;
952 }
953
954 /*************************************** Run Ipopt for optimization ***********************************/
955 #if defined(OMC_HAVE_IPOPT)
956 static int ipopt_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
957 {
958 int cJ, res;
959
960 ✗ cJ = threadData->currentErrorStage;
961 ✗ threadData->currentErrorStage = ERROR_OPTIMIZE;
962 ✗ res = runOptimizer(data, threadData, solverInfo);
963 ✗ threadData->currentErrorStage = cJ;
964 return res;
965 }
966 #endif
967
968 ✗ static void writeOutputVars(char* names, DATA* data)
969 {
970 int i = 0;
971 char *p = NULL;
972 /* fix names to contain | instead of , for splitting */
973 ✗ parseVariableStr(names);
974 ✗ p = strtok(names, "!");
975
976 ✗ fprintf(stdout, "time=%.20g", data->localData[0]->timeValue);
977
978 ✗ while(p)
979 {
980 ✗ for(i = 0; i < data->modelData->nVariablesReal; i++)
981 ✗ if(!strcmp(p, data->modelData->realVarsData[i].info.name))
982 ✗ fprintf(stdout, ",%s=%.20g", p, (data->localData[0])->realVars[i]);
983 ✗ for(i = 0; i < data->modelData->nVariablesInteger; i++)
984 ✗ if(!strcmp(p, data->modelData->integerVarsData[i].info.name))
985 ✗ fprintf(stdout, ",%s=" OMC_INT_FORMAT, p, (data->localData[0])->integerVars[i]);
986 ✗ for(i = 0; i < data->modelData->nVariablesBoolean; i++)
987 ✗ if(!strcmp(p, data->modelData->booleanVarsData[i].info.name))
988 ✗ fprintf(stdout, ",%s=%i", p, (data->localData[0])->booleanVars[i]);
989 ✗ for(i = 0; i < data->modelData->nVariablesString; i++)
990 ✗ if(!strcmp(p, data->modelData->stringVarsData[i].info.name))
991 ✗ fprintf(stdout, ",%s=\"%s\"", p, omc_string_data((data->localData[0])->stringVars[i]));
992
993 ✗ for(i = 0; i < data->modelData->nAliasReal; i++)
994 ✗ if(!strcmp(p, data->modelData->realAlias[i].info.name))
995 {
996 ✗ if(data->modelData->realAlias[i].negate)
997 ✗ fprintf(stdout, ",%s=%.20g", p, -(data->localData[0])->realVars[data->modelData->realAlias[i].nameID]);
998 else
999 ✗ fprintf(stdout, ",%s=%.20g", p, (data->localData[0])->realVars[data->modelData->realAlias[i].nameID]);
1000 }
1001 ✗ for(i = 0; i < data->modelData->nAliasInteger; i++)
1002 ✗ if(!strcmp(p, data->modelData->integerAlias[i].info.name))
1003 {
1004 ✗ if(data->modelData->integerAlias[i].negate)
1005 ✗ fprintf(stdout, ",%s=" OMC_INT_FORMAT, p, -(data->localData[0])->integerVars[data->modelData->integerAlias[i].nameID]);
1006 else
1007 ✗ fprintf(stdout, ",%s=" OMC_INT_FORMAT, p, (data->localData[0])->integerVars[data->modelData->integerAlias[i].nameID]);
1008 }
1009 ✗ for(i = 0; i < data->modelData->nAliasBoolean; i++)
1010 ✗ if(!strcmp(p, data->modelData->booleanAlias[i].info.name))
1011 {
1012 ✗ if(data->modelData->booleanAlias[i].negate)
1013 ✗ fprintf(stdout, ",%s=%i", p, -(data->localData[0])->booleanVars[data->modelData->booleanAlias[i].nameID]);
1014 else
1015 ✗ fprintf(stdout, ",%s=%i", p, (data->localData[0])->booleanVars[data->modelData->booleanAlias[i].nameID]);
1016 }
1017 ✗ for(i = 0; i < data->modelData->nAliasString; i++)
1018 ✗ if(!strcmp(p, data->modelData->stringAlias[i].info.name))
1019 ✗ fprintf(stdout, ",%s=\"%s\"", p, omc_string_data((data->localData[0])->stringVars[data->modelData->stringAlias[i].nameID]));
1020
1021 /* parameters */
1022 ✗ for(i = 0; i < data->modelData->nParametersReal; i++)
1023 ✗ if(!strcmp(p, data->modelData->realParameterData[i].info.name))
1024 ✗ fprintf(stdout, ",%s=%.20g", p, data->simulationInfo->realParameter[i]);
1025
1026 ✗ for(i = 0; i < data->modelData->nParametersInteger; i++)
1027 ✗ if(!strcmp(p, data->modelData->integerParameterData[i].info.name))
1028 ✗ fprintf(stdout, ",%s=" OMC_INT_FORMAT, p, data->simulationInfo->integerParameter[i]);
1029
1030 ✗ for(i = 0; i < data->modelData->nParametersBoolean; i++)
1031 ✗ if(!strcmp(p, data->modelData->booleanParameterData[i].info.name))
1032 ✗ fprintf(stdout, ",%s=%i", p, data->simulationInfo->booleanParameter[i]);
1033
1034 ✗ for(i = 0; i < data->modelData->nParametersString; i++)
1035 ✗ if(!strcmp(p, data->modelData->stringParameterData[i].info.name))
1036 ✗ fprintf(stdout, ",%s=\"%s\"", p, omc_string_data(data->simulationInfo->stringParameter[i]));
1037
1038 /* move to next */
1039 ✗ p = strtok(NULL, "!");
1040 }
1041 ✗ fprintf(stdout, "\n"); fflush(stdout);
1042 ✗ }
1043
1044 /**
1045 * @brief Set all solver stats to zero.
1046 *
1047 * @param stats Pointer to solver stats.
1048 */
1049 2 void resetSolverStats(SOLVERSTATS* stats) {
1050 2 stats->nStepsTaken = 0;
1051 2 stats->nCallsODE = 0;
1052 2 stats->nCallsJacobian = 0;
1053 2 stats->nErrorTestFailures = 0;
1054 2 stats->nConvergenceTestFailures = 0;
1055 2 }
1056
1057 /**
1058 * @brief Add two solver statistics.
1059 *
1060 * destStats += addStats
1061 *
1062 * @param destStats Pointer to solver stats to add stats to.
1063 * On return has result of addition.
1064 * @param addStats Pointer to solver stats to add.
1065 */
1066 1 void addSolverStats(SOLVERSTATS* destStats, SOLVERSTATS* addStats) {
1067 1 destStats->nStepsTaken += addStats->nStepsTaken;
1068 1 destStats->nCallsODE += addStats->nCallsODE;
1069 1 destStats->nCallsJacobian += addStats->nCallsJacobian;
1070 1 destStats->nErrorTestFailures += addStats->nErrorTestFailures;
1071 1 destStats->nConvergenceTestFailures += addStats->nConvergenceTestFailures;
1072 1 }
1073