Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/solver_main.c
Warning:line 1037, column 30
The expression is an uninitialized value. The computed value will also be garbage

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2014, Open Source Modelica Consortium (OSMC),
5 * c/o Linköpings universitet, Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
11 * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
14 * ACCORDING TO RECIPIENTS CHOICE.
15 *
16 * The OpenModelica software and the OSMC (Open Source Modelica Consortium)
17 * Public License (OSMC-PL) are obtained from OSMC, either from the above
18 * address, from the URLs: http://www.openmodelica.org or
19 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica
20 * distribution. GNU version 3 is obtained from:
21 * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from:
22 * http://www.opensource.org/licenses/BSD-3-Clause.
23 *
24 * This program is distributed WITHOUT ANY WARRANTY; without even the implied
25 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS
26 * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE
27 * CONDITIONS OF OSMC-PL.
28 *
29 */
30
31 /*! \file solver_main.c
32 */
33
34#include "omc_config.h"
35#include "simulation/simulation_runtime.h"
36#include "simulation/results/simulation_result.h"
37#include "solver_main.h"
38#include "openmodelica_func.h"
39#include "initialization/initialization.h"
40#include "nonlinearSystem.h"
41#include "newtonIteration.h"
42#include "dassl.h"
43#include "ida_solver.h"
44#include "delay.h"
45#include "events.h"
46#include "util/parallel_helper.h"
47#include "util/varinfo.h"
48#include "radau.h"
49#include "model_help.h"
50#include "meta/meta_modelica.h"
51#include "simulation/solver/epsilon.h"
52#include "simulation/solver/external_input.h"
53#include "linearSystem.h"
54#include "sym_solver_ssc.h"
55#include "irksco.h"
56#if !defined(OMC_MINIMAL_RUNTIME)
57#include "simulation/solver/embedded_server.h"
58#include "simulation/solver/real_time_sync.h"
59#endif
60#include "simulation/simulation_input_xml.h"
61
62#include "optimization/OptimizerInterface.h"
63
64/*
65 * #include "dopri45.h"
66 */
67#include "util/rtclock.h"
68#include "util/omc_error.h"
69#include "simulation/options.h"
70#include <math.h>
71#include <string.h>
72#include <errno(*__errno_location ()).h>
73#include <float.h>
74
75double** work_states;
76
77typedef struct RK4_DATA
78{
79 double** work_states;
80 int work_states_ndims;
81 const double *b;
82 const double *c;
83 double h;
84}RK4_DATA;
85
86
87static int euler_ex_step(DATA* data, SOLVER_INFO* solverInfo);
88static int rungekutta_step_ssc(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
89static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
90static int sym_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
91
92static int radau_lobatto_step(DATA* data, SOLVER_INFO* solverInfo);
93
94#ifdef WITH_IPOPT
95static int ipopt_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
96#endif
97
98static void writeOutputVars(char* names, DATA* data);
99
100int solver_main_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
101{
102 TRACE_PUSH
103 int retVal;
104
105 switch(solverInfo->solverMethod)
1
Control jumps to 'case S_SYM_SOLVER:' at line 165
106 {
107 case S_EULER:
108 retVal = euler_ex_step(data, solverInfo);
109 if(omc_flag[FLAG_SOLVER_STEPS])
110 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
111 TRACE_POP
112 return retVal;
113 case S_HEUN:
114 case S_RUNGEKUTTA:
115 retVal = rungekutta_step(data, threadData, solverInfo);
116 if(omc_flag[FLAG_SOLVER_STEPS])
117 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
118 TRACE_POP
119 return retVal;
120
121#if !defined(OMC_MINIMAL_RUNTIME)
122 case S_DASSL:
123 retVal = dassl_step(data, threadData, solverInfo);
124 if(omc_flag[FLAG_SOLVER_STEPS])
125 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
126 TRACE_POP
127 return retVal;
128#endif
129
130#ifdef WITH_IPOPT
131 case S_OPTIMIZATION:
132 if ((int)(data->modelData->nStates + data->modelData->nInputVars) > 0){
133 retVal = ipopt_step(data, threadData, solverInfo);
134 } else {
135 solverInfo->solverMethod = S_EULER;
136 retVal = euler_ex_step(data, solverInfo);
137 }
138 if(omc_flag[FLAG_SOLVER_STEPS])
139 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
140 TRACE_POP
141 return retVal;
142#endif
143#ifdef WITH_SUNDIALS
144 case S_IMPEULER:
145 case S_TRAPEZOID:
146 case S_IMPRUNGEKUTTA:
147 retVal = radau_lobatto_step(data, solverInfo);
148 if(omc_flag[FLAG_SOLVER_STEPS])
149 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
150 TRACE_POP
151 return retVal;
152 case S_IDA:
153 retVal = ida_solver_step(data, threadData, solverInfo);
154 if(omc_flag[FLAG_SOLVER_STEPS])
155 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
156 TRACE_POP
157 return retVal;
158#endif
159 case S_ERKSSC:
160 retVal = rungekutta_step_ssc(data, threadData, solverInfo);
161 if(omc_flag[FLAG_SOLVER_STEPS])
162 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
163 TRACE_POP
164 return retVal;
165 case S_SYM_SOLVER:
166 retVal = sym_solver_step(data, threadData, solverInfo);
2
Calling 'sym_solver_step'
167 if(omc_flag[FLAG_SOLVER_STEPS])
168 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
169 return retVal;
170 case S_SYM_SOLVER_SSC:
171 retVal = sym_solver_ssc_step(data, threadData, solverInfo);
172 if(omc_flag[FLAG_SOLVER_STEPS])
173 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
174 return retVal;
175 case S_IRKSCO:
176 {
177 retVal = irksco_midpoint_rule(data, threadData, solverInfo);
178 if(omc_flag[FLAG_SOLVER_STEPS])
179 data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
180 return retVal;
181 }
182
183 }
184
185 TRACE_POP
186 return 1;
187}
188
189/*! \fn initializeSolverData(DATA* data, SOLVER_INFO* solverInfo)
190 *
191 * \param [ref] [data]
192 * \param [ref] [solverInfo]
193 *
194 * This function initializes solverInfo.
195 */
196int initializeSolverData(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
197{
198 TRACE_PUSH
199 int retValue = 0;
200 int i;
201
202 SIMULATION_INFO *simInfo = data->simulationInfo;
203
204 /* initial solverInfo */
205 solverInfo->currentTime = simInfo->startTime;
206 solverInfo->currentStepSize = simInfo->stepSize;
207 solverInfo->laststep = 0;
208 solverInfo->solverRootFinding = 0;
209 solverInfo->solverNoEquidistantGrid = 0;
210 solverInfo->lastdesiredStep = solverInfo->currentTime + solverInfo->currentStepSize;
211 solverInfo->eventLst = allocList(sizeof(long));
212 solverInfo->didEventStep = 0;
213 solverInfo->stateEvents = 0;
214 solverInfo->sampleEvents = 0;
215 solverInfo->solverStats = (unsigned int*) calloc(numStatistics, sizeof(unsigned int));
216 solverInfo->solverStatsTmp = (unsigned int*) calloc(numStatistics, sizeof(unsigned int));
217
218 /* if FLAG_NOEQUIDISTANT_GRID is set, choose integrator step method */
219 if (omc_flag[FLAG_NOEQUIDISTANT_GRID])
220 {
221 solverInfo->integratorSteps = 1; /* TRUE */
222 }
223 else
224 {
225 solverInfo->integratorSteps = 0;
226 }
227
228 /* set tolerance for ZeroCrossings */
229 /* TODO: Check this! */
230 /* setZCtol(fmin(simInfo->stepSize, simInfo->tolerance)); */
231
232 switch (solverInfo->solverMethod)
233 {
234 case S_SYM_SOLVER:
235 case S_EULER: break;
236 case S_SYM_SOLVER_SSC:
237 {
238 allocateSymSolverSsc(solverInfo, data->modelData->nStates);
239 break;
240 }
241 case S_IRKSCO:
242 {
243 allocateIrksco(solverInfo, data->modelData->nStates, data->modelData->nZeroCrossings);
244 break;
245 }
246 case S_ERKSSC:
247 case S_RUNGEKUTTA:
248 case S_HEUN:
249 {
250 /* Allocate RK work arrays */
251
252 static const int rungekutta_s = 4;
253 static const double rungekutta_b[4] = { 1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0 };
254 static const double rungekutta_c[4] = { 0.0, 0.5, 0.5, 1.0 };
255
256 static const int heun_s = 2;
257 static const double heun_b[2] = { 1.0 / 2.0, 1.0 / 2.0 };
258 static const double heun_c[2] = { 0.0, 1.0 };
259
260 RK4_DATA* rungeData = (RK4_DATA*) malloc(sizeof(RK4_DATA));
261
262 if (solverInfo->solverMethod==S_RUNGEKUTTA) {
263 rungeData->work_states_ndims = rungekutta_s;
264 rungeData->b = rungekutta_b;
265 rungeData->c = rungekutta_c;
266 } else if (solverInfo->solverMethod==S_HEUN) {
267 rungeData->work_states_ndims = heun_s;
268 rungeData->b = heun_b;
269 rungeData->c = heun_c;
270 } else if (solverInfo->solverMethod==S_ERKSSC) {
271 rungeData->h = (omc_flag[FLAG_INITIAL_STEP_SIZE]) ? atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]) : solverInfo->currentStepSize;
272 rungeData->work_states_ndims = 5;
273 } else {
274 throwStreamPrint(threadData, "Unknown RK solver");
275 }
276
277 rungeData->work_states = (double**) malloc((rungeData->work_states_ndims + 1) * sizeof(double*));
278 for (i = 0; i < rungeData->work_states_ndims + 1; i++) {
279 rungeData->work_states[i] = (double*) calloc(data->modelData->nStates, sizeof(double));
280 }
281 solverInfo->solverData = rungeData;
282 break;
283 }
284 case S_QSS: break;
285#if !defined(OMC_MINIMAL_RUNTIME)
286 case S_DASSL:
287 {
288 /* Initial DASSL solver */
289 DASSL_DATA* dasslData = (DASSL_DATA*) malloc(sizeof(DASSL_DATA));
290 retValue = dassl_initial(data, threadData, solverInfo, dasslData);
291 solverInfo->solverData = dasslData;
292 break;
293 }
294#endif
295#ifdef WITH_IPOPT
296 case S_OPTIMIZATION:
297 {
298 infoStreamPrint(LOG_SOLVER, 0, "Initializing optimizer");
299 /* solverInfo->solverData = malloc(sizeof(OptData)); */
300 break;
301 }
302#endif
303#ifdef WITH_SUNDIALS
304 case S_IMPEULER:
305 case S_TRAPEZOID:
306 case S_IMPRUNGEKUTTA:
307 {
308 int usedImpRKOrder = DEFAULT_IMPRK_ORDER5;
309 if (solverInfo->solverMethod == S_IMPEULER)
310 usedImpRKOrder = 1;
311 if (solverInfo->solverMethod == S_TRAPEZOID)
312 usedImpRKOrder = 2;
313
314 /* Check the order if set */
315 if (omc_flag[FLAG_IMPRK_ORDER])
316 {
317 usedImpRKOrder = atoi(omc_flagValue[FLAG_IMPRK_ORDER]);
318 if (usedImpRKOrder>6 || usedImpRKOrder<1)
319 {
320 warningStreamPrint(LOG_STDOUT, 0, "Selected order %d is out of range[1-6]. Use default order %d", usedImpRKOrder, DEFAULT_IMPRK_ORDER5);
321 usedImpRKOrder = DEFAULT_IMPRK_ORDER5;
322 }
323 }
324
325 /* Allocate implicit Runge-Kutta methods */
326 infoStreamPrint(LOG_SOLVER, 0, "Initializing Runge-Kutta method with order %d", usedImpRKOrder);
327 solverInfo->solverData = calloc(1, sizeof(KINODE));
328 allocateKinOde(data, threadData, solverInfo, usedImpRKOrder);
329 break;
330 }
331 case S_IDA:
332 {
333 IDA_SOLVER* idaData = NULL((void*)0);
334 /* Allocate ida working data */
335 infoStreamPrint(LOG_SOLVER, 0, "Initializing IDA DAE Solver");
336 idaData = (IDA_SOLVER*) malloc(sizeof(IDA_SOLVER));
337 retValue = ida_solver_initial(data, threadData, solverInfo, idaData);
338 solverInfo->solverData = idaData;
339 break;
340 }
341#endif
342 default:
343 errorStreamPrint(LOG_SOLVER, 0, "Solver %s disabled on this configuration", SOLVER_METHOD_NAME[solverInfo->solverMethod]);
344 TRACE_POP
345 return 1;
346 }
347
348 TRACE_POP
349 return retValue;
350}
351
352/*! \fn freeSolver(DATA* data, SOLVER_INFO* solverInfo)
353 *
354 * \param [ref] [data]
355 * \param [ref] [solverInfo]
356 *
357 * This function frees solverInfo.
358 */
359int freeSolverData(DATA* data, SOLVER_INFO* solverInfo)
360{
361 int retValue = 0;
362 int i;
363
364 freeList(solverInfo->eventLst);
365 /* free solver statistics */
366 free(solverInfo->solverStats);
367 free(solverInfo->solverStatsTmp);
368 /* deintialize solver related workspace */
369 if (solverInfo->solverMethod == S_SYM_SOLVER_SSC)
370 {
371 freeSymSolverSsc(solverInfo);
372 }
373 else if(solverInfo->solverMethod == S_RUNGEKUTTA || solverInfo->solverMethod == S_HEUN
374 || solverInfo->solverMethod == S_ERKSSC)
375 {
376 /* free RK work arrays */
377 for(i = 0; i < ((RK4_DATA*)(solverInfo->solverData))->work_states_ndims + 1; i++)
378 free(((RK4_DATA*)(solverInfo->solverData))->work_states[i]);
379 free(((RK4_DATA*)(solverInfo->solverData))->work_states);
380 free((RK4_DATA*)solverInfo->solverData);
381 }
382 else if (solverInfo->solverMethod == S_IRKSCO)
383 {
384 freeIrksco(solverInfo);
385 }
386#if !defined(OMC_MINIMAL_RUNTIME)
387 else if(solverInfo->solverMethod == S_DASSL)
388 {
389 /* De-Initial DASSL solver */
390 dassl_deinitial(solverInfo->solverData);
391 }
392#endif
393#ifdef WITH_IPOPT
394 else if(solverInfo->solverMethod == S_OPTIMIZATION)
395 {
396 /* free work arrays */
397 /*destroyIpopt(solverInfo);*/
398 }
399#endif
400#ifdef WITH_SUNDIALS
401 else if(solverInfo->solverMethod == S_IMPEULER ||
402 solverInfo->solverMethod == S_TRAPEZOID ||
403 solverInfo->solverMethod == S_IMPRUNGEKUTTA)
404 {
405 /* free work arrays */
406 freeKinOde(data, solverInfo);
407 }
408 else if(solverInfo->solverMethod == S_IDA)
409 {
410 /* free work arrays */
411 ida_solver_deinitial(solverInfo->solverData);
412 }
413#endif
414 {
415 /* free other solver memory */
416 }
417
418 return retValue;
419}
420
421
422/*! \fn initializeModel(DATA* data, const char* init_initMethod,
423 * const char* init_file, double init_time)
424 *
425 * \param [ref] [data]
426 * \param [in] [pInitMethod] user defined initialization method
427 * \param [in] [pInitFile] extra argument for initialization-method "file"
428 * \param [in] [initTime] extra argument for initialization-method "file"
429 *
430 * This function starts the initialization process of the model .
431 */
432int initializeModel(DATA* data, threadData_t *threadData, const char* init_initMethod,
433 const char* init_file, double init_time)
434{
435 TRACE_PUSH
436 int retValue = 0;
437
438 SIMULATION_INFO *simInfo = data->simulationInfo;
439
440 if(measure_time_flag)
441 {
442 rt_accumulate(SIM_TIMER_PREINIT6);
443 rt_tick(SIM_TIMER_INIT1);
444 }
445
446 copyStartValuestoInitValues(data);
447
448 /* read input vars */
449 data->callback->input_function_init(data, threadData);
450 externalInputUpdate(data);
451 data->callback->input_function_updateStartValues(data, threadData);
452 data->callback->input_function(data, threadData);
453
454 data->localData[0]->timeValue = simInfo->startTime;
455
456 threadData->currentErrorStage = ERROR_SIMULATION;
457 /* try */
458 {
459 int success = 0;
460#if !defined(OMC_EMCC)
461 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
462#endif
463 if(initialization(data, threadData, init_initMethod, init_file, init_time))
464 {
465 warningStreamPrint(LOG_STDOUT, 0, "Error in initialization. Storing results and exiting.\nUse -lv=LOG_INIT -w for more information.");
466 simInfo->stopTime = simInfo->startTime;
467 retValue = -1;
468 }
469 if (!retValue)
470 {
471 if (data->simulationInfo->homotopySteps == 0)
472 infoStreamPrint(LOG_SUCCESS, 0, "The initialization finished successfully without homotopy method.");
473 else
474 infoStreamPrint(LOG_SUCCESS, 0, "The initialization finished successfully with %d homotopy steps.", data->simulationInfo->homotopySteps);
475 }
476
477 success = 1;
478#if !defined(OMC_EMCC)
479 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
480#endif
481
482 if (!success)
483 {
484 retValue = -1;
485 infoStreamPrint(LOG_ASSERT, 0, "simulation terminated by an assertion at initialization");
486 }
487 }
488
489 /* adrpo: write the parameter data in the file once again after bound parameters and initialization! */
490 sim_result.writeParameterData(&sim_result,data,threadData);
491 infoStreamPrint(LOG_SOLVER, 0, "Wrote parameters to the file after initialization (for output formats that support this)");
492
493 /* Initialization complete */
494 if (measure_time_flag) {
495 rt_accumulate(SIM_TIMER_INIT1);
496 }
497
498 TRACE_POP
499 return retValue;
500}
501
502
503/*! \fn finishSimulation(DATA* data, SOLVER_INFO* solverInfo)
504 *
505 * \param [ref] [data]
506 * \param [ref] [solverInfo]
507 *
508 * This function performs the last step
509 * and outputs some statistics, this this simulation terminal step.
510 */
511int finishSimulation(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, const char* outputVariablesAtEnd)
512{
513 TRACE_PUSH
514
515 int retValue = 0;
516 int ui;
517 double t, total100;
518
519 SIMULATION_INFO *simInfo = data->simulationInfo;
520
521 /* Last step with terminal()=true */
522 if(solverInfo->currentTime >= simInfo->stopTime && solverInfo->solverMethod != S_OPTIMIZATION) {
523 infoStreamPrint(LOG_EVENTS_V, 0, "terminal event at stop time %g", solverInfo->currentTime);
524 data->simulationInfo->terminal = 1;
525 updateDiscreteSystem(data, threadData);
526
527 /* prevent emit if noeventemit flag is used */
528 if (!(omc_flag[FLAG_NOEVENTEMIT])) {
529 sim_result.emit(&sim_result, data, threadData);
530 }
531
532 data->simulationInfo->terminal = 0;
533 }
534
535 if (0 != strcmp("ia", data->simulationInfo->outputFormat)) {
536 communicateStatus("Finished", 1, solverInfo->currentTime, solverInfo->currentStepSize);
537 }
538
539 /* we have output variables in the command line -output a,b,c */
540 if(outputVariablesAtEnd)
541 {
542 writeOutputVars(strdup(outputVariablesAtEnd), data);
543 }
544
545 if(ACTIVE_STREAM(LOG_STATS)(useStream[LOG_STATS]))
546 {
547 rt_accumulate(SIM_TIMER_TOTAL0);
548
549 infoStreamPrint(LOG_STATS, 1, "### STATISTICS ###");
550
551 total100 = rt_accumulated(SIM_TIMER_TOTAL0)/100.0;
552
553 infoStreamPrint(LOG_STATS, 1, "timer");
554 infoStreamPrint(LOG_STATS, 0, "%12gs reading init.xml", rt_accumulated(SIM_TIMER_INIT_XML13));
555 infoStreamPrint(LOG_STATS, 0, "%12gs reading info.xml", rt_accumulated(SIM_TIMER_INFO_XML14));
556 infoStreamPrint(LOG_STATS, 0, "%12gs [%5.1f%%] pre-initialization", rt_accumulated(SIM_TIMER_PREINIT6), rt_accumulated(SIM_TIMER_PREINIT6)/total100);
557 infoStreamPrint(LOG_STATS, 0, "%12gs [%5.1f%%] initialization", rt_accumulated(SIM_TIMER_INIT1), rt_accumulated(SIM_TIMER_INIT1)/total100);
558 infoStreamPrint(LOG_STATS, 0, "%12gs [%5.1f%%] steps", rt_accumulated(SIM_TIMER_STEP2), rt_accumulated(SIM_TIMER_STEP2)/total100);
559 infoStreamPrint(LOG_STATS, 0, "%12gs [%5.1f%%] solver (excl. callbacks)", rt_accumulated(SIM_TIMER_SOLVER12), rt_accumulated(SIM_TIMER_SOLVER12)/total100);
560 infoStreamPrint(LOG_STATS, 0, "%12gs [%5.1f%%] creating output-file", rt_accumulated(SIM_TIMER_OUTPUT3), rt_accumulated(SIM_TIMER_OUTPUT3)/total100);
561 infoStreamPrint(LOG_STATS, 0, "%12gs [%5.1f%%] event-handling", rt_accumulated(SIM_TIMER_EVENT4), rt_accumulated(SIM_TIMER_EVENT4)/total100);
562 infoStreamPrint(LOG_STATS, 0, "%12gs [%5.1f%%] overhead", rt_accumulated(SIM_TIMER_OVERHEAD7), rt_accumulated(SIM_TIMER_OVERHEAD7)/total100);
563
564 t = rt_accumulated(SIM_TIMER_TOTAL0)-rt_accumulated(SIM_TIMER_OVERHEAD7)-rt_accumulated(SIM_TIMER_EVENT4)-rt_accumulated(SIM_TIMER_OUTPUT3)-rt_accumulated(SIM_TIMER_STEP2)-rt_accumulated(SIM_TIMER_INIT1)-rt_accumulated(SIM_TIMER_PREINIT6)-rt_accumulated(SIM_TIMER_SOLVER12);
565 infoStreamPrint(LOG_STATS, 0, "%12gs [%5.1f%%] %s", t, t/total100, S_OPTIMIZATION == solverInfo->solverMethod ? "optimization" : "simulation");
566
567 infoStreamPrint(LOG_STATS, 0, "%12gs [100.0%%] total", rt_accumulated(SIM_TIMER_TOTAL0));
568 messageClose(LOG_STATS);
569
570 infoStreamPrint(LOG_STATS, 1, "events");
571 infoStreamPrint(LOG_STATS, 0, "%5ld state events", solverInfo->stateEvents);
572 infoStreamPrint(LOG_STATS, 0, "%5ld time events", solverInfo->sampleEvents);
573 messageClose(LOG_STATS);
574
575 if(S_OPTIMIZATION == solverInfo->solverMethod || /* skip solver statistics for optimization */
576 S_QSS == solverInfo->solverMethod) /* skip also for qss, since not available*/
577 {
578 }
579 else
580 {
581 /* save stats before print */
582 for(ui=0; ui<numStatistics; ui++)
583 solverInfo->solverStats[ui] += solverInfo->solverStatsTmp[ui];
584
585 infoStreamPrint(LOG_STATS, 1, "solver: %s", SOLVER_METHOD_NAME[solverInfo->solverMethod]);
586 infoStreamPrint(LOG_STATS, 0, "%5d steps taken", solverInfo->solverStats[0]);
587
588 infoStreamPrint(LOG_STATS, 0, "%5d calls of functionODE", solverInfo->solverStats[1]);
589
590 infoStreamPrint(LOG_STATS, 0, "%5d evaluations of jacobian", solverInfo->solverStats[2]);
591
592 infoStreamPrint(LOG_STATS, 0, "%5d error test failures", solverInfo->solverStats[3]);
593 infoStreamPrint(LOG_STATS, 0, "%5d convergence test failures", solverInfo->solverStats[4]);
594 infoStreamPrint(LOG_STATS, 0, "%gs time of jacobian evaluation", rt_accumulated(SIM_TIMER_JACOBIAN5));
595#ifdef USE_PARJAC
596 infoStreamPrint(LOG_STATS, 0, "%i OpenMP-threads used for jacobian evaluation", omc_get_max_threads());
597 int chunk_size;
598 omp_sched_t kind;
599 omp_get_schedule(&kind, &chunk_size);
600 infoStreamPrint(LOG_STATS, 0, "Schedule: %i Chunk Size: %i", kind, chunk_size);
601#endif
602
603 messageClose(LOG_STATS);
604 }
605
606 infoStreamPrint(LOG_STATS_V, 1, "function calls");
607 if (compiledInDAEMode)
608 {
609 infoStreamPrint(LOG_STATS_V, 1, "%5ld calls of functionDAE", rt_ncall(SIM_TIMER_DAE15));
610 infoStreamPrint(LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_DAE15), rt_accumulated(SIM_TIMER_DAE15)/total100);
611 messageClose(LOG_STATS_V);
612 }
613 if (data->simulationInfo->callStatistics.functionODE) {
614 infoStreamPrint(LOG_STATS_V, 1, "%5ld calls of functionODE", data->simulationInfo->callStatistics.functionODE);
615 infoStreamPrint(LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_FUNCTION_ODE8), rt_accumulated(SIM_TIMER_FUNCTION_ODE8)/total100);
616 messageClose(LOG_STATS_V);
617 }
618
619 if (rt_ncall(SIM_TIMER_RESIDUALS9)) {
620 infoStreamPrint(LOG_STATS_V, 1, "%5d calls of functionODE_residual", rt_ncall(SIM_TIMER_RESIDUALS9));
621 infoStreamPrint(LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_RESIDUALS9), rt_accumulated(SIM_TIMER_RESIDUALS9)/total100);
622 messageClose(LOG_STATS_V);
623 }
624
625 if (data->simulationInfo->callStatistics.functionAlgebraics) {
626 infoStreamPrint(LOG_STATS_V, 1, "%5d calls of functionAlgebraics", data->simulationInfo->callStatistics.functionAlgebraics);
627 infoStreamPrint(LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_ALGEBRAICS10), rt_accumulated(SIM_TIMER_ALGEBRAICS10)/total100);
628 messageClose(LOG_STATS_V);
629 }
630
631 infoStreamPrint(LOG_STATS_V, 1, "%5d evaluations of jacobian", rt_ncall(SIM_TIMER_JACOBIAN5));
632 infoStreamPrint(LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_JACOBIAN5), rt_accumulated(SIM_TIMER_JACOBIAN5)/total100);
633 messageClose(LOG_STATS_V);
634
635 infoStreamPrint(LOG_STATS_V, 0, "%5ld calls of updateDiscreteSystem", data->simulationInfo->callStatistics.updateDiscreteSystem);
636 infoStreamPrint(LOG_STATS_V, 0, "%5ld calls of functionZeroCrossingsEquations", data->simulationInfo->callStatistics.functionZeroCrossingsEquations);
637
638 infoStreamPrint(LOG_STATS_V, 1, "%5ld calls of functionZeroCrossings", data->simulationInfo->callStatistics.functionZeroCrossings);
639 infoStreamPrint(LOG_STATS_V, 0, "%12gs [%5.1f%%]", rt_accumulated(SIM_TIMER_ZC11), rt_accumulated(SIM_TIMER_ZC11)/total100);
640 messageClose(LOG_STATS_V);
641
642 messageClose(LOG_STATS_V);
643
644 infoStreamPrint(LOG_STATS_V, 1, "linear systems");
645 for(ui=0; ui<data->modelData->nLinearSystems; ui++)
646 printLinearSystemSolvingStatistics(data, ui, LOG_STATS_V);
647 messageClose(LOG_STATS_V);
648
649 infoStreamPrint(LOG_STATS_V, 1, "non-linear systems");
650 for(ui=0; ui<data->modelData->nNonLinearSystems; ui++)
651 printNonLinearSystemSolvingStatistics(data, ui, LOG_STATS_V);
652 messageClose(LOG_STATS_V);
653
654 messageClose(LOG_STATS);
655 rt_tick(SIM_TIMER_TOTAL0);
656 }
657
658 TRACE_POP
659 return retValue;
660}
661
662/*! \fn solver_main
663 *
664 * \param [ref] [data]
665 * \param [in] [pInitMethod] user defined initialization method
666 * \param [in] [pOptiMethod] user defined optimization method
667 * \param [in] [pInitFile] extra argument for initialization-method "file"
668 * \param [in] [initTime] extra argument for initialization-method "file"
669 * \param [in] [solverID] selects the ode solver
670 * \param [in] [outputVariablesAtEnd] ???
671 *
672 * This is the main function of the solver, it performs the simulation.
673 */
674int solver_main(DATA* data, threadData_t *threadData, const char* init_initMethod, const char* init_file,
675 double init_time, int solverID, const char* outputVariablesAtEnd, const char *argv_0)
676{
677 TRACE_PUSH
678
679 int i, retVal = 1, initSolverInfo = 0;
680 unsigned int ui;
681 SOLVER_INFO solverInfo;
682 SIMULATION_INFO *simInfo = data->simulationInfo;
683 void *dllHandle=NULL((void*)0);
684
685 solverInfo.solverMethod = solverID;
686
687 /* do some solver specific checks */
688 switch(solverInfo.solverMethod)
689 {
690#ifndef WITH_SUNDIALS
691 case S_IMPEULER:
692 case S_TRAPEZOID:
693 case S_IMPRUNGEKUTTA:
694 warningStreamPrint(LOG_STDOUT, 0, "Sundial/kinsol is needed but not available. Please choose other solver.");
695 TRACE_POP
696 return 1;
697#endif
698
699#ifndef WITH_IPOPT
700 case S_OPTIMIZATION:
701 warningStreamPrint(LOG_STDOUT, 0, "Ipopt is needed but not available.");
702 TRACE_POP
703 return 1;
704#endif
705
706 }
707
708 /* first initialize the model then allocate SolverData memory
709 * due to be able to use the initialized values for the integrator
710 */
711 simInfo->useStopTime = 1;
712
713 /* Use minStepSize if stepSize is getting too small, but
714 * allow stepSize to be zero if startTime == stopTime.
715 */
716 if ((simInfo->stepSize < simInfo->minStepSize) && (simInfo->stopTime > simInfo->startTime)){
717 warningStreamPrint(LOG_STDOUT, 0, "The step-size %g is too small. Adjust the step-size to %g.", simInfo->stepSize, simInfo->minStepSize);
718 simInfo->stepSize = simInfo->minStepSize;
719 simInfo->numSteps = round((simInfo->stopTime - simInfo->startTime)/simInfo->stepSize);
720 }
721#if !defined(OMC_EMCC)
722 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
723#endif
724
725 /* initialize external input structure */
726 externalInputallocate(data);
727 /* set tolerance for ZeroCrossings */
728 setZCtol(fmin(data->simulationInfo->stepSize, data->simulationInfo->tolerance));
729 omc_alloc_interface.collect_a_little();
730
731 /* initialize solver data */
732 /* For the DAEmode we need to initialize solverData before the initialization,
733 * since the solver is used to obtain consistent values also via updateDiscreteSystem
734 */
735 retVal = initializeSolverData(data, threadData, &solverInfo);
736 initSolverInfo = 1;
737
738 /* initialize all parts of the model */
739 if (0 == retVal){
740 retVal = initializeModel(data, threadData, init_initMethod, init_file, init_time);
741 omc_alloc_interface.collect_a_little();
742 }
743
744#if !defined(OMC_MINIMAL_RUNTIME)
745 dllHandle = embedded_server_load_functions(omc_flagValue[FLAG_EMBEDDED_SERVER]);
746 omc_real_time_sync_init(threadData, data);
747 int port = 4841;
748 /* If an embedded server is specified */
749 if (dllHandle != NULL((void*)0)) {
750 if (omc_flag[FLAG_EMBEDDED_SERVER_PORT]) {
751 port = atoi(omc_flagValue[FLAG_EMBEDDED_SERVER_PORT]);
752 /* In case of a bad conversion, don't spawn a server on port 0...*/
753 if (port == 0) {
754 port = 4841;
755 }
756 }
757 }
758 data->embeddedServerState = embedded_server_init(data, data->localData[0]->timeValue, solverInfo.currentStepSize, argv_0, omc_real_time_sync_update, port);
759 /* If an embedded server is specified */
760 if (dllHandle != NULL((void*)0)) {
761 infoStreamPrint(LOG_STDOUT, 0, "The embedded server is initialized.");
762 }
763 wait_for_step(data->embeddedServerState);
764#endif
765 if(0 == retVal) {
766 retVal = -1;
767 /* if the model has no time changing variables skip the main loop*/
768 if(data->modelData->nVariablesReal == 0 &&
769 data->modelData->nVariablesInteger == 0 &&
770 data->modelData->nVariablesBoolean == 0 &&
771 data->modelData->nVariablesString == 0 ) {
772 /* prevent emit if noeventemit flag is used */
773 if (!(omc_flag[FLAG_NOEVENTEMIT])) {
774 sim_result.emit(&sim_result, data, threadData);
775 }
776
777 infoStreamPrint(LOG_SOLVER, 0, "The model has no time changing variables, no integration will be performed.");
778 solverInfo.currentTime = simInfo->stopTime;
779 data->localData[0]->timeValue = simInfo->stopTime;
780 overwriteOldSimulationData(data);
781 retVal = finishSimulation(data, threadData, &solverInfo, outputVariablesAtEnd);
782 } else if(S_QSS == solverInfo.solverMethod) {
783 /* starts the simulation main loop - special solvers */
784 sim_result.emit(&sim_result,data,threadData);
785
786 /* overwrite the whole ring-buffer with initialized values */
787 overwriteOldSimulationData(data);
788
789 infoStreamPrint(LOG_SOLVER, 0, "Start numerical integration (startTime: %g, stopTime: %g)", simInfo->startTime, simInfo->stopTime);
790 retVal = data->callback->performQSSSimulation(data, threadData, &solverInfo);
791 omc_alloc_interface.collect_a_little();
792
793 /* terminate the simulation */
794 finishSimulation(data, threadData, &solverInfo, outputVariablesAtEnd);
795 omc_alloc_interface.collect_a_little();
796 } else {
797 /* starts the simulation main loop - standard solver interface */
798 if(omc_flag[FLAG_SOLVER_STEPS])
799 data->simulationInfo->solverSteps = 0;
800 if(solverInfo.solverMethod != S_OPTIMIZATION) {
801 sim_result.emit(&sim_result,data,threadData);
802 }
803
804 /* overwrite the whole ring-buffer with initialized values */
805 overwriteOldSimulationData(data);
806
807 /* store all values for non-dassl event search */
808 storeOldValues(data);
809
810 infoStreamPrint(LOG_SOLVER, 0, "Start numerical solver from %g to %g", simInfo->startTime, simInfo->stopTime);
811 retVal = data->callback->performSimulation(data, threadData, &solverInfo);
812 omc_alloc_interface.collect_a_little();
813 /* terminate the simulation */
814 //if (solverInfo.solverMethod == S_SYM_SOLVER_SSC) data->callback->symbolicInlineSystems(data, threadData, 0, 2);
815 finishSimulation(data, threadData, &solverInfo, outputVariablesAtEnd);
816 omc_alloc_interface.collect_a_little();
817 }
818 }
819
820 if (data->real_time_sync.enabled) {
821 int tMaxLate=0;
822 const char *unit = prettyPrintNanoSec(data->real_time_sync.maxLate, &tMaxLate);
823 infoStreamPrint(LOG_RT, 0, "Maximum real-time latency was (positive=missed dealine, negative is slack): %d %s", tMaxLate, unit);
824 }
825#if !defined(OMC_MINIMAL_RUNTIME)
826 embedded_server_deinit(data->embeddedServerState);
827 embedded_server_unload_functions(dllHandle);
828#endif
829
830#if !defined(OMC_EMCC)
831 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
832#endif
833
834 /* free external input data */
835 externalInputFree(data);
836
837 /* free SolverInfo memory */
838 if (initSolverInfo)
839 {
840 freeSolverData(data, &solverInfo);
841 }
842
843 if (!retVal)
844 infoStreamPrint(LOG_SUCCESS, 0, "The simulation finished successfully.");
845
846 TRACE_POP
847 return retVal;
848}
849
850/*************************************** EULER_EXP *********************************/
851static int euler_ex_step(DATA* data, SOLVER_INFO* solverInfo)
852{
853 int i;
854 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
855 SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
856 modelica_real* stateDer = sDataOld->realVars + data->modelData->nStates;
857
858 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
859
860 solverInfo->currentTime = sDataOld->timeValue + solverInfo->currentStepSize;
861
862 for(i = 0; i < data->modelData->nStates; i++)
863 {
864 sData->realVars[i] = sDataOld->realVars[i] + stateDer[i] * solverInfo->currentStepSize;
865 }
866 sData->timeValue = solverInfo->currentTime;
867
868 /* save stats */
869 /* steps */
870 solverInfo->solverStatsTmp[0] += 1;
871 /* function ODE evaluation is done directly after this function */
872 solverInfo->solverStatsTmp[1] += 1;
873
874 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
875
876 return 0;
877}
878
879/*************************************** EXP_RK_SSC *********************************/
880static int rungekutta_step_ssc(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
881{
882 // see.
883 // Solving Stiff Systems of ODEs by Explicit Methods with Conformed Stability
884 // Domains
885 // 2016 9th EUROSIM Congress on Modelling and Simulation
886 // A. E. Novikov...
887
888 RK4_DATA *rk = ((RK4_DATA*)(solverInfo->solverData));
889 const double Atol = data->simulationInfo->tolerance;
890 const double Rtol = data->simulationInfo->tolerance;
891 const int nx = data->modelData->nStates;
892 modelica_real h = rk->h;
893 double** k = rk->work_states;
894 int j, i;
895 double sum;
896 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
897 SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
898 modelica_real* stateDer = sData->realVars + nx;
899 modelica_real* stateDerOld = sDataOld->realVars + nx;
900 double t = sDataOld->timeValue;
901 const double targetTime = t + solverInfo->currentStepSize;
902 const short isMaxStepSizeSet = (short) omc_flagValue[FLAG_MAX_STEP_SIZE];
903 const double maxStepSize = isMaxStepSizeSet ? atof(omc_flagValue[FLAG_MAX_STEP_SIZE]) : -1;
904#if defined(_MSC_VER)
905 /* handle stupid compilers */
906 double *x0 = (double*)malloc(nx*sizeof(double));
907#else
908 double x0[nx];
909#endif
910
911 if(t + h > targetTime)
912 h = solverInfo->currentStepSize;
913
914 memcpy(k[0], stateDerOld, nx*sizeof(double));
915 memcpy(x0, sDataOld->realVars, nx*sizeof(double));
916 while(t < targetTime && h > 0){
917 //printf("\nfrom %g to %g by %g\n", t, targetTime, h);
918 for(j = 0; j < 5; ++j){
919
920 if(j > 0)
921 memcpy(k[j], stateDer, nx*sizeof(double));
922
923 switch(j){
924 case 0:
925 //yn + k1/3
926 for(i = 0; i < nx; ++i)
927 sData->realVars[i] = x0[i] + h * k[0][i]/3.0;
928 sData->timeValue = t + h/3.0;
929 break;
930
931 case 1:
932 //yn + 1/6(k1 + k2)
933 for(i = 0; i < nx; ++i)
934 sData->realVars[i] = x0[i] + h/6.0 * (k[0][i] + k[1][i]);
935 sData->timeValue = t + h/3.0;
936 break;
937
938 case 2:
939 //yn + 1/8*(k1 + 3*k3)
940 for(i = 0; i < nx; ++i)
941 sData->realVars[i] = x0[i] + h/8.0 * (k[0][i] + 3*k[2][i]);
942 sData->timeValue = t + h/2.0;
943 break;
944
945 case 3:
946 //yn + 1/2*(k1 - 3*k3 + 4*k4)
947 for(i = 0; i < nx; ++i)
948 sData->realVars[i] = x0[i] + h/2.0 * (k[0][i] - 3*k[2][i] + 4*k[3][i]);
949 sData->timeValue = t + h;
950 break;
951
952 case 4:
953 //yn + 1/6*(k1 + 4*k3 + k5)
954 for(i = 0; i < nx; ++i){
955 sData->realVars[i] = x0[i] + h/6.0 * (k[0][i] + 4*k[3][i] + k[4][i]);
956 }
957 sData->timeValue = t + h;
958 break;
959
960 }
961 //f(yn + ...)
962 /* read input vars */
963 externalInputUpdate(data);
964 data->callback->input_function(data, threadData);
965 /* eval ode equations */
966 data->callback->functionODE(data, threadData);
967 }
968 t += h;
969 sData->timeValue = t;
970 solverInfo->currentTime = t;
971
972 /* save stats */
973 /* steps */
974 solverInfo->solverStatsTmp[0] += 1;
975 /* function ODE evaluation is done directly after this */
976 solverInfo->solverStatsTmp[1] += 4;
977
978
979 //stepsize
980 for(i = 0, sum = 0.0; i < nx; ++i)
981 sum = fmax(fabs(k[0][i] + 4*k[3][i] -(4.5*k[2][i] + k[4][i]))/(fabs(k[4][i]) + fabs(k[2][i]) + fabs(k[3][i])+ + fabs(k[0][i]) + Atol), sum);
982 sum *= 2.0/30.0;
983
984
985 h = fmin(0.9*fmax(pow(sum,1/4.0)/(Atol ), 1e-12)*h + 1e-12, (targetTime - h));
986 if(isMaxStepSizeSet && h > maxStepSize) h = maxStepSize;
987 if (h > 0) rk->h = h;
988
989 if(t < targetTime){
990 memcpy(x0, sData->realVars, nx*sizeof(double));
991 memcpy(k[0], k[4], nx*sizeof(double));
992 sim_result.emit(&sim_result, data, threadData);
993 }
994 }
995
996 //assert(sData->timeValue == targetTime);
997 //assert(solverInfo->currentTime == targetTime);
998#if defined(_MSC_VER)
999 /* handle stupid compilers */
1000 free(x0);
1001#endif
1002
1003 return 0;
1004}
1005
1006
1007/*************************************** SYM_SOLVER *********************************/
1008static int sym_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo){
1009 int retVal,i,j;
3
'j' declared without an initial value
1010
1011 modelica_integer nStates = data->modelData->nStates;
1012 SIMULATION_DATA *sData = data->localData[0];
1013 SIMULATION_DATA *sDataOld = data->localData[1];
1014 modelica_real* stateDer = sDataOld->realVars + data->modelData->nStates;
1015
1016 if (solverInfo->currentStepSize >= DASSL_STEP_EPS){
4
Taking true branch
1017 /* time */
1018 solverInfo->currentTime = sDataOld->timeValue + solverInfo->currentStepSize;
1019 sData->timeValue = solverInfo->currentTime;
1020 /* update dt */
1021 data->simulationInfo->inlineData->dt = solverInfo->currentStepSize;
1022 /* copy old states to workspace */
1023 memcpy(data->simulationInfo->inlineData->algOldVars, sDataOld->realVars, nStates * sizeof(double));
1024 memcpy(sData->realVars, sDataOld->realVars, nStates * sizeof(double));
1025
1026 /* read input vars */
1027 externalInputUpdate(data);
1028 data->callback->input_function(data, threadData);
1029 retVal = data->callback->symbolicInlineSystems(data, threadData);
1030
1031 if(retVal != 0){
5
Assuming 'retVal' is equal to 0
6
Taking false branch
1032 return -1;
1033 }
1034
1035
1036 /* update der(x) */
1037 for(i=0; i<nStates; ++i, ++j)
7
Assuming 'i' is < 'nStates'
8
Loop condition is true. Entering loop body
9
The expression is an uninitialized value. The computed value will also be garbage
1038 {
1039 stateDer[i] = (sData->realVars[i]-data->simulationInfo->inlineData->algOldVars[i])/solverInfo->currentStepSize;
1040 }
1041
1042 /* save stats */
1043 /* steps */
1044 solverInfo->solverStatsTmp[0] += 1;
1045 /* function ODE evaluation is done directly after this */
1046 solverInfo->solverStatsTmp[1] += 1;
1047 }
1048 else
1049 /* in case desired step size is too small */
1050 {
1051 infoStreamPrint(LOG_SOLVER, 0, "Desired step to small try next one");
1052 infoStreamPrint(LOG_SOLVER, 0, "Interpolate linear");
1053
1054 /* explicit euler step*/
1055 for(i = 0; i < nStates; i++)
1056 {
1057 sData->realVars[i] = sDataOld->realVars[i] + stateDer[i] * solverInfo->currentStepSize;
1058 }
1059 sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
1060 solverInfo->currentTime = sData->timeValue;
1061 }
1062 return retVal;
1063}
1064
1065/*************************************** RK4 ***********************************/
1066static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
1067{
1068 RK4_DATA *rk = ((RK4_DATA*)(solverInfo->solverData));
1069 double** k = rk->work_states;
1070 double sum;
1071 int i,j;
1072 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
1073 SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
1074 modelica_real* stateDer = sData->realVars + data->modelData->nStates;
1075 modelica_real* stateDerOld = sDataOld->realVars + data->modelData->nStates;
1076
1077 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1078
1079 solverInfo->currentTime = sDataOld->timeValue + solverInfo->currentStepSize;
1080
1081 /* We calculate k[0] before returning from this function.
1082 * We only want to calculate f() 4 times per call */
1083 memcpy(k[0], stateDerOld, data->modelData->nStates*sizeof(modelica_real));
1084
1085 for (j = 1; j < rk->work_states_ndims; j++)
1086 {
1087 for(i = 0; i < data->modelData->nStates; i++)
1088 {
1089 sData->realVars[i] = sDataOld->realVars[i] + solverInfo->currentStepSize * rk->c[j] * k[j - 1][i];
1090 }
1091 sData->timeValue = sDataOld->timeValue + rk->c[j] * solverInfo->currentStepSize;
1092 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1093 /* read input vars */
1094 externalInputUpdate(data);
1095 data->callback->input_function(data, threadData);
1096 /* eval ode equations */
1097 data->callback->functionODE(data, threadData);
1098 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1099 memcpy(k[j], stateDer, data->modelData->nStates*sizeof(modelica_real));
1100
1101 }
1102
1103 for(i = 0; i < data->modelData->nStates; i++)
1104 {
1105 sum = 0;
1106 for(j = 0; j < rk->work_states_ndims; j++)
1107 {
1108 sum = sum + rk->b[j] * k[j][i];
1109 }
1110 sData->realVars[i] = sDataOld->realVars[i] + solverInfo->currentStepSize * sum;
1111 }
1112 sData->timeValue = solverInfo->currentTime;
1113
1114 /* save stats */
1115 /* steps */
1116 solverInfo->solverStatsTmp[0] += 1;
1117 /* function ODE evaluation is done directly after this */
1118 solverInfo->solverStatsTmp[1] += rk->work_states_ndims+1;
1119 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1120
1121 return 0;
1122}
1123
1124/*************************************** Run Ipopt for optimization ***********************************/
1125#if defined(WITH_IPOPT)
1126static int ipopt_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
1127{
1128 int cJ, res;
1129
1130 cJ = threadData->currentErrorStage;
1131 threadData->currentErrorStage = ERROR_OPTIMIZE;
1132 res = runOptimizer(data, threadData, solverInfo);
1133 threadData->currentErrorStage = cJ;
1134 return res;
1135}
1136#endif
1137
1138#if defined(WITH_SUNDIALS) && !defined(OMC_MINIMAL_RUNTIME)
1139/*************************************** Radau/Lobatto ***********************************/
1140int radau_lobatto_step(DATA* data, SOLVER_INFO* solverInfo)
1141{
1142 if(kinsolOde(solverInfo) == 0)
1143 {
1144 solverInfo->currentTime += solverInfo->currentStepSize;
1145 return 0;
1146 }
1147 return -1;
1148}
1149#endif
1150
1151
1152static void writeOutputVars(char* names, DATA* data)
1153{
1154 int i = 0;
1155 char *p = NULL((void*)0);
1156 /* fix names to contain | instead of , for splitting */
1157 parseVariableStr(names);
1158 p = strtok(names, "!");
1159
1160 fprintf(stdoutstdout, "time=%.20g", data->localData[0]->timeValue);
1161
1162 while(p)
1163 {
1164 for(i = 0; i < data->modelData->nVariablesReal; i++)
1165 if(!strcmp(p, data->modelData->realVarsData[i].info.name))
1166 fprintf(stdoutstdout, ",%s=%.20g", p, (data->localData[0])->realVars[i]);
1167 for(i = 0; i < data->modelData->nVariablesInteger; i++)
1168 if(!strcmp(p, data->modelData->integerVarsData[i].info.name))
1169 fprintf(stdoutstdout, ",%s=%li", p, (data->localData[0])->integerVars[i]);
1170 for(i = 0; i < data->modelData->nVariablesBoolean; i++)
1171 if(!strcmp(p, data->modelData->booleanVarsData[i].info.name))
1172 fprintf(stdoutstdout, ",%s=%i", p, (data->localData[0])->booleanVars[i]);
1173 for(i = 0; i < data->modelData->nVariablesString; i++)
1174 if(!strcmp(p, data->modelData->stringVarsData[i].info.name))
1175 fprintf(stdoutstdout, ",%s=\"%s\"", p, MMC_STRINGDATA((data->localData[0])->stringVars[i])(((struct mmc_string*)((void*)((char*)((data->localData[0]
)->stringVars[i]) - 3)))->data)
);
1176
1177 for(i = 0; i < data->modelData->nAliasReal; i++)
1178 if(!strcmp(p, data->modelData->realAlias[i].info.name))
1179 {
1180 if(data->modelData->realAlias[i].negate)
1181 fprintf(stdoutstdout, ",%s=%.20g", p, -(data->localData[0])->realVars[data->modelData->realAlias[i].nameID]);
1182 else
1183 fprintf(stdoutstdout, ",%s=%.20g", p, (data->localData[0])->realVars[data->modelData->realAlias[i].nameID]);
1184 }
1185 for(i = 0; i < data->modelData->nAliasInteger; i++)
1186 if(!strcmp(p, data->modelData->integerAlias[i].info.name))
1187 {
1188 if(data->modelData->integerAlias[i].negate)
1189 fprintf(stdoutstdout, ",%s=%li", p, -(data->localData[0])->integerVars[data->modelData->integerAlias[i].nameID]);
1190 else
1191 fprintf(stdoutstdout, ",%s=%li", p, (data->localData[0])->integerVars[data->modelData->integerAlias[i].nameID]);
1192 }
1193 for(i = 0; i < data->modelData->nAliasBoolean; i++)
1194 if(!strcmp(p, data->modelData->booleanAlias[i].info.name))
1195 {
1196 if(data->modelData->booleanAlias[i].negate)
1197 fprintf(stdoutstdout, ",%s=%i", p, -(data->localData[0])->booleanVars[data->modelData->booleanAlias[i].nameID]);
1198 else
1199 fprintf(stdoutstdout, ",%s=%i", p, (data->localData[0])->booleanVars[data->modelData->booleanAlias[i].nameID]);
1200 }
1201 for(i = 0; i < data->modelData->nAliasString; i++)
1202 if(!strcmp(p, data->modelData->stringAlias[i].info.name))
1203 fprintf(stdoutstdout, ",%s=\"%s\"", p, MMC_STRINGDATA((data->localData[0])->stringVars[data->modelData->stringAlias[i].nameID])(((struct mmc_string*)((void*)((char*)((data->localData[0]
)->stringVars[data->modelData->stringAlias[i].nameID
]) - 3)))->data)
);
1204
1205 /* parameters */
1206 for(i = 0; i < data->modelData->nParametersReal; i++)
1207 if(!strcmp(p, data->modelData->realParameterData[i].info.name))
1208 fprintf(stdoutstdout, ",%s=%.20g", p, data->simulationInfo->realParameter[i]);
1209
1210 for(i = 0; i < data->modelData->nParametersInteger; i++)
1211 if(!strcmp(p, data->modelData->integerParameterData[i].info.name))
1212 fprintf(stdoutstdout, ",%s=%li", p, data->simulationInfo->integerParameter[i]);
1213
1214 for(i = 0; i < data->modelData->nParametersBoolean; i++)
1215 if(!strcmp(p, data->modelData->booleanParameterData[i].info.name))
1216 fprintf(stdoutstdout, ",%s=%i", p, data->simulationInfo->booleanParameter[i]);
1217
1218 for(i = 0; i < data->modelData->nParametersString; i++)
1219 if(!strcmp(p, data->modelData->stringParameterData[i].info.name))
1220 fprintf(stdoutstdout, ",%s=\"%s\"", p, MMC_STRINGDATA(data->simulationInfo->stringParameter[i])(((struct mmc_string*)((void*)((char*)(data->simulationInfo
->stringParameter[i]) - 3)))->data)
);
1221
1222 /* move to next */
1223 p = strtok(NULL((void*)0), "!");
1224 }
1225 fprintf(stdoutstdout, "\n"); fflush(stdoutstdout);
1226}