Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 8.2% 13 / 0 / 158
Functions: 25.0% 3 / 0 / 12
Branches: 5.9% 6 / 0 / 102
OMCompiler/SimulationRuntime/c/simulation/solver/events.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 #include "events.h"
29 #include "../../util/omc_error.h"
30 #include "../options.h"
31 #include "../../simulation_data.h"
32 #include "../../openmodelica.h" /* for modelica types */
33 #include "../../openmodelica_func.h" /* for modelica functions */
34 #include "solver_main.h"
35 #include "model_help.h"
36 #include "external_input.h"
37 #include "epsilon.h"
38
39 #include <math.h>
40 #include <stdio.h>
41 #include <stdlib.h>
42 #include <string.h>
43
44 int maxBisectionIterations = 0;
45 void bisection(DATA* data, threadData_t *threadData, double*, double*, double*, double*, LIST*, LIST*);
46 void saveZeroCrossingsAfterEvent(DATA *data, threadData_t *threadData);
47
48 /*! \fn checkForSampleEvent
49 *
50 * \param [ref] [data]
51 * \param [ref] [solverInfo]
52 * \return indicates if a time event is occurred or not.
53 *
54 * Function check if a sample expression should be activated
55 * before next step and sets then the next step size to the
56 * time event.
57 *
58 */
59 1 void checkForSampleEvent(DATA *data, SOLVER_INFO* solverInfo)
60 {
61 1 double nextTimeStep = solverInfo->currentTime + solverInfo->currentStepSize;
62 /* A short run would otherwise lose output points to the pulled-in step. */
63 1 double eps = fmin(SAMPLE_EPS, 1e-9 * fabs(data->simulationInfo->stopTime - data->simulationInfo->startTime));
64
65
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1 if ((data->simulationInfo->nextSampleEvent <= nextTimeStep + eps) && (data->simulationInfo->nextSampleEvent >= solverInfo->currentTime))
66 {
67 const double t = solverInfo->currentTime, te = data->simulationInfo->nextSampleEvent;
68 ✗ double h = te - t;
69 int i;
70 /* A relation `time >= te` only switches if the step lands on te exactly. */
71 ✗ for (i = 0; i < 4 && t + h < te; i++) h = nextafter(h, INFINITY);
72 ✗ for (i = 0; i < 4 && t + h > te; i++) h = nextafter(h, 0.0);
73 ✗ solverInfo->currentStepSize = h;
74 ✗ data->simulationInfo->sampleActivated = 1;
75 ✗ infoStreamPrint(OMC_LOG_EVENTS_V, 0, "Adjust step-size to %.15g at time %.15g to get next sample event at %.15g", solverInfo->currentStepSize, solverInfo->currentTime, data->simulationInfo->nextSampleEvent );
76 }
77 1 }
78
79 /*! \fn checkForStateEvent
80 *
81 * \param [ref] [data]
82 * \param [ref] [eventList]
83 *
84 * This function checks for events in interval=[oldTime, timeValue]
85 * If a zero crossing function cause a sign change, root finding
86 * process will start
87 */
88 1 int checkForStateEvent(DATA* data, LIST *eventList)
89 {
90 long i=0;
91
92
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(i=0; i<data->modelData->nZeroCrossings; i++)
93 {
94 int *eq_indexes;
95
96 // Check if sign of zero crossing changed
97 ✗ if(sign(data->simulationInfo->zeroCrossings[i]) != sign(data->simulationInfo->zeroCrossingsPre[i]))
98 {
99 ✗ listPushFront(eventList, &(data->simulationInfo->zeroCrossingIndex[i]));
100 }
101 }
102
103
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if(listLen(eventList) > 0)
104 {
105 ✗ return 1;
106 }
107
108 return 0;
109 }
110
111 /*! \fn checkEvents
112 *
113 * This function check if a time event or a state event should
114 * processed. If sample and state event have the same event-time
115 * then time events are prioritize, since they handle also
116 * state event. It returns 1 if state event is before time event
117 * then it de-activate the time events.
118 *
119 * \param [ref] [data]
120 * \param [ref] [threadData]
121 * \param [ref] [eventLst]
122 * \param [in] [useRootFinding]
123 * \param [out] [eventTime]
124 * \return 0: no event; 1: time event; 2: state event
125 */
126 1 int checkEvents(DATA* data, threadData_t *threadData, LIST* eventLst, modelica_boolean useRootFinding, double *eventTime)
127 {
128 1 int found = checkForStateEvent(data, eventLst);
129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(found && useRootFinding)
130 {
131 ✗ *eventTime = findRoot(data, threadData, eventLst, data->simulationInfo->timeValueOld, data->simulationInfo->realVarsOld, data->localData[0]->timeValue, data->localData[0]->realVars);
132 }
133
134
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if(data->simulationInfo->sampleActivated == 1)
135 {
136 return 1;
137 }
138
139
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if(listLen(eventLst) > 0)
140 {
141 ✗ return 2;
142 }
143
144 return 0;
145 }
146
147 /* Chattering: numEvents state events in a row within less than the step size
148 * and intervalFraction*(stopTime-startTime). numEventLimit is the largest
149 * numEvents. */
150 static const struct { int numEvents; double intervalFraction; } chatteringLimits[] = {
151 {1000, 1e-6},
152 {100, 1e-9}
153 };
154
155 ✗ static double chatteringTimeLimit(SIMULATION_INFO *simulationInfo, double intervalFraction)
156 {
157 ✗ double interval = simulationInfo->stopTime - simulationInfo->startTime;
158 ✗ if (interval > 0 && isfinite(interval)) {
159 ✗ return fmin(simulationInfo->stepSize, intervalFraction * interval);
160 }
161 ✗ return simulationInfo->stepSize;
162 }
163
164 /*! \fn handleEvents
165 *
166 * \param [ref] [data]
167 * \param [ref] [eventList]
168 * \param [in] [eventTime]
169 *
170 * This handles all zero crossing events from event list at event time
171 */
172 ✗ void handleEvents(DATA* data, threadData_t *threadData, LIST* eventLst, double *eventTime, SOLVER_INFO* solverInfo)
173 {
174 ✗ double time = data->localData[0]->timeValue;
175 long i;
176 LIST_NODE* it;
177
178 /* time event */
179 ✗ if(data->simulationInfo->sampleActivated)
180 {
181 ✗ storePreValues(data);
182
183 /* activate time event */
184 ✗ for(i=0; i<data->modelData->nSamples; ++i) {
185 ✗ if(data->simulationInfo->nextSampleTimes[i] <= time + SAMPLE_EPS)
186 {
187 ✗ data->simulationInfo->samples[i] = 1;
188 ✗ infoStreamPrint(OMC_LOG_EVENTS, 0, "[%ld] sample(%g, %g)", data->modelData->samplesInfo[i].index, data->modelData->samplesInfo[i].start, data->modelData->samplesInfo[i].interval);
189 }
190 }
191
192 ✗ solverInfo->sampleEvents++;
193 }
194 /* state event */
195 ✗ if(listLen(eventLst)>0)
196 {
197 ✗ CHATTERING_INFO *chattering = &data->simulationInfo->chatteringInfo;
198 ✗ data->localData[0]->timeValue = *eventTime;
199 /* time = data->localData[0]->timeValue; */
200
201 ✗ if (omc_useStream[OMC_LOG_EVENTS])
202 {
203 ✗ for(it = listFirstNode(eventLst); it; it = listNextNode(it))
204 {
205 ✗ long ix = *((long*) listNodeData(it));
206 int *eq_indexes;
207 ✗ const char *exp_str = data->callback->zeroCrossingDescription(ix,&eq_indexes);
208 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_EVENTS, omc_dummyFileInfo, 0, eq_indexes, "[%ld] %s", ix+1, exp_str);
209 }
210 }
211
212 ✗ solverInfo->stateEvents++;
213 ✗ if (chattering->stateEventsInARow < chattering->numEventLimit) {
214 ✗ chattering->stateEventsInARow++;
215 }
216 ✗ chattering->lastTimes[chattering->currentIndex]=time;
217
218 ✗ for (i = 0; !chattering->messageEmitted && i < (long) (sizeof(chatteringLimits)/sizeof(chatteringLimits[0])); i++)
219 {
220 ✗ int numEvents = chatteringLimits[i].numEvents;
221 double t0, limit;
222 ✗ if (chattering->stateEventsInARow < numEvents) {
223 ✗ continue;
224 }
225 ✗ t0 = chattering->lastTimes[(chattering->currentIndex + chattering->numEventLimit - (numEvents-1)) % chattering->numEventLimit];
226 ✗ limit = chatteringTimeLimit(data->simulationInfo, chatteringLimits[i].intervalFraction);
227 ✗ if (time - t0 < limit)
228 {
229 ✗ long ix = *((long*) listNodeData(listFirstNode(eventLst)));
230 int *eq_indexes;
231 ✗ const char *exp_str = data->callback->zeroCrossingDescription(ix,&eq_indexes);
232 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_STDOUT, omc_dummyFileInfo, 0, eq_indexes, "Chattering detected around time %.12g..%.12g (%d state events in a row with a total time delta less than %.12g, the smaller of the step size and %g times the simulation interval). This can be a performance bottleneck. Use -lv LOG_EVENTS for more information. The zero-crossing was: %s", t0, time, numEvents, limit, chatteringLimits[i].intervalFraction, exp_str);
233 ✗ chattering->messageEmitted = 1;
234 ✗ if (omc_flag[FLAG_ABORT_SLOW])
235 {
236 ✗ throwStreamPrintWithEquationIndexes(threadData, omc_dummyFileInfo, eq_indexes, "Aborting simulation due to chattering being detected and the simulation flags requesting we do not continue further.");
237 }
238 }
239 }
240 ✗ chattering->currentIndex = (chattering->currentIndex+1) % chattering->numEventLimit;
241
242 ✗ listClear(eventLst);
243 } else {
244 ✗ data->simulationInfo->chatteringInfo.stateEventsInARow = 0;
245 }
246
247 /* update the whole system */
248 ✗ updateDiscreteSystem(data, threadData);
249 ✗ if (OMC_ERROR_RAISED()) {
250 return;
251 }
252 ✗ saveZeroCrossingsAfterEvent(data, threadData);
253 /*sim_result_emit(data);*/
254
255 ✗ updateNextSampleEvent(data, threadData);
256 ✗ data->simulationInfo->sampleActivated = 0;
257 }
258
259 /*! \fn findRoot
260 *
261 * \param [ref] [data]
262 * \param [ref] [threadData]
263 * \param [ref] [eventList]
264 * \param [in] [time_left]
265 * \param [in] [values_left]
266 * \param [in] [time_right]
267 * \param [in] [values_right]
268 * \return: first event of interval [time_left, time_right]
269 */
270 ✗ double findRoot(DATA* data, threadData_t* threadData, LIST* eventList, double time_left, double* values_left, double time_right, double* values_right)
271 {
272 LIST_NODE* it;
273 fortran_integer i=0;
274 ✗ LIST *tmpEventList = allocList(eventListAlloc, eventListFree, eventListCopy);
275
276 /* static work arrays */
277 ✗ double *states_left = data->simulationInfo->states_left;
278 ✗ double *states_right = data->simulationInfo->states_right;
279
280
281 /* write states to work arrays */
282 ✗ memcpy(states_left, values_left, data->modelData->nStates * sizeof(double));
283 ✗ memcpy(states_right, values_right, data->modelData->nStates * sizeof(double));
284
285 ✗ for(it=listFirstNode(eventList); it; it=listNextNode(it))
286 {
287 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "search for current event. Events in list: %ld", *((long*)listNodeData(it)));
288 }
289
290 /* Search for event time and event_id with bisection method */
291 ✗ bisection(data, threadData, &time_left, &time_right, states_left, states_right, tmpEventList, eventList);
292
293 /* what happens here? */
294 ✗ if(listLen(tmpEventList) == 0)
295 {
296 ✗ double value = fabs(data->simulationInfo->zeroCrossings[*((long*) listFirstData(eventList))]);
297 ✗ for(it = listFirstNode(eventList); it; it = listNextNode(it))
298 {
299 ✗ double fvalue = fabs(data->simulationInfo->zeroCrossings[*((long*) listNodeData(it))]);
300 ✗ if(value > fvalue)
301 {
302 value = fvalue;
303 }
304 }
305 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "Minimum value: %e", value);
306 ✗ for(it = listFirstNode(eventList); it; it = listNextNode(it))
307 {
308 ✗ if(value == fabs(data->simulationInfo->zeroCrossings[*((long*) listNodeData(it))]))
309 {
310 ✗ listPushBack(tmpEventList, listNodeData(it));
311 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "added tmp event : %ld", *((long*) listNodeData(it)));
312 }
313 }
314 }
315
316 ✗ listClear(eventList);
317
318 ✗ while(listLen(tmpEventList) > 0)
319 {
320 ✗ long event_id = *((long*)listFirstData(tmpEventList));
321 ✗ listPushFrontNodeNoCopy(eventList, listPopFrontNode(tmpEventList));
322 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "Event id: %ld", event_id);
323 }
324
325 ✗ data->localData[0]->timeValue = time_left;
326 ✗ memcpy(data->localData[0]->realVars, states_left, data->modelData->nStates * sizeof(double));
327
328 /* determined continuous system */
329 ✗ data->callback->updateContinuousSystem(data, threadData);
330 ✗ updateRelationsPre(data);
331 /*sim_result_emit(data);*/
332
333 ✗ data->localData[0]->timeValue = time_right;
334 ✗ memcpy(data->localData[0]->realVars, states_right, data->modelData->nStates * sizeof(double));
335
336 ✗ freeList(tmpEventList);
337
338 ✗ return time_right;
339 }
340
341 /*! \fn bisection
342 *
343 * \param [ref] [data]
344 * \param [ref] [a]
345 * \param [ref] [b]
346 * \param [ref] [states_a]
347 * \param [ref] [states_b]
348 * \param [ref] [eventListTmp]
349 * \param [in] [eventList]
350 *
351 * Method to find root in interval [oldTime, timeValue]
352 */
353 ✗ void bisection(DATA* data, threadData_t *threadData, double* a, double* b, double* states_a, double* states_b, LIST *tmpEventList, LIST *eventList)
354 {
355 ✗ double TTOL = MINIMAL_STEP_SIZE + MINIMAL_STEP_SIZE*fabs(*b-*a); /* absTol + relTol*abs(b-a) */
356 double c;
357 long i=0;
358 /* n >= log(2)/log(2) + log(|b-a|/TOL)/log(2)*/
359 ✗ unsigned int n = maxBisectionIterations > 0 ? maxBisectionIterations : 1 + ceil(log(fabs(*b - *a)/TTOL)/log(2));
360
361 ✗ memcpy(data->simulationInfo->zeroCrossingsBackup, data->simulationInfo->zeroCrossings, data->modelData->nZeroCrossings * sizeof(modelica_real));
362
363 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "bisection method starts in interval [%e, %e]", *a, *b);
364 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "TTOL is set to %e and maximum number of intersections %d.", TTOL, n);
365
366 ✗ while(fabs(*b - *a) > MINIMAL_STEP_SIZE && n-- > 0)
367 {
368 ✗ c = 0.5 * (*a + *b);
369 ✗ data->localData[0]->timeValue = c;
370
371 /*calculates states at time c */
372 ✗ for(i=0; i < data->modelData->nStates; i++)
373 {
374 ✗ data->localData[0]->realVars[i] = 0.5*(states_a[i] + states_b[i]);
375 }
376
377 /*calculates Values dependents on new states*/
378 /* read input vars */
379 ✗ externalInputUpdate(data);
380 ✗ data->callback->input_function(data, threadData);
381 /* eval needed equations*/
382 ✗ data->callback->function_ZeroCrossingsEquations(data, threadData);
383
384 ✗ data->callback->function_ZeroCrossings(data, threadData, data->simulationInfo->zeroCrossings);
385
386 ✗ if(checkZeroCrossings(data, tmpEventList, eventList)) /* If Zerocrossing in left Section */
387 {
388 ✗ memcpy(states_b, data->localData[0]->realVars, data->modelData->nStates * sizeof(modelica_real));
389 ✗ *b = c;
390 ✗ memcpy(data->simulationInfo->zeroCrossingsBackup, data->simulationInfo->zeroCrossings, data->modelData->nZeroCrossings * sizeof(modelica_real));
391 }
392 else /*else Zerocrossing in right Section */
393 {
394 ✗ memcpy(states_a, data->localData[0]->realVars, data->modelData->nStates * sizeof(modelica_real));
395 ✗ *a = c;
396 ✗ memcpy(data->simulationInfo->zeroCrossingsPre, data->simulationInfo->zeroCrossings, data->modelData->nZeroCrossings * sizeof(modelica_real));
397 ✗ memcpy(data->simulationInfo->zeroCrossings, data->simulationInfo->zeroCrossingsBackup, data->modelData->nZeroCrossings * sizeof(modelica_real));
398 }
399 }
400 ✗ }
401
402 /*! \fn checkZeroCrossings
403 *
404 * Function checks for an event list on events
405 *
406 * \param [ref] [data]
407 * \param [ref] [eventListTmp]
408 * \param [in] [eventList]
409 * \return boolean value
410 */
411 ✗ int checkZeroCrossings(DATA *data, LIST *tmpEventList, LIST *eventList)
412 {
413 LIST_NODE *it;
414
415 ✗ listClear(tmpEventList);
416 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "bisection checks for condition changes");
417
418 ✗ for(it=listFirstNode(eventList); it; it=listNextNode(it))
419 {
420 /* found event in left section */
421 ✗ if((data->simulationInfo->zeroCrossings[*((long*) listNodeData(it))] == -1 &&
422 ✗ data->simulationInfo->zeroCrossingsPre[*((long*) listNodeData(it))] == 1) ||
423 ✗ (data->simulationInfo->zeroCrossings[*((long*) listNodeData(it))] == 1 &&
424 ✗ data->simulationInfo->zeroCrossingsPre[*((long*) listNodeData(it))] == -1))
425 {
426 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "%ld changed from %s to current %s",
427 ✗ *((long*) listNodeData(it)),
428 ✗ (data->simulationInfo->zeroCrossingsPre[*((long*) listNodeData(it))] > 0) ? "TRUE" : "FALSE",
429 ✗ (data->simulationInfo->zeroCrossings[*((long*) listNodeData(it))] > 0) ? "TRUE" : "FALSE");
430 ✗ listPushFront(tmpEventList, listNodeData(it));
431 }
432 }
433
434 ✗ if(listLen(tmpEventList) > 0)
435 {
436 ✗ return 1; /* event in left section */
437 }
438
439 return 0; /* event in right section */
440 }
441
442 /*! \fn saveZeroCrossingsAfterEvent
443 *
444 * Function saves all zero-crossing values as pre(zero-crossing)
445 *
446 * \param [ref] [data]
447 */
448 ✗ void saveZeroCrossingsAfterEvent(DATA *data, threadData_t *threadData)
449 {
450 long i=0;
451
452 ✗ infoStreamPrint(OMC_LOG_ZEROCROSSINGS, 0, "save all zerocrossings after an event at time=%g", data->localData[0]->timeValue); /* ??? */
453
454 ✗ data->callback->function_ZeroCrossings(data, threadData, data->simulationInfo->zeroCrossings);
455 ✗ for(i=0; i<data->modelData->nZeroCrossings; i++) {
456 ✗ data->simulationInfo->zeroCrossingsPre[i] = data->simulationInfo->zeroCrossings[i];
457 }
458 ✗ }
459
460 /**
461 * @brief Allocate memory for eventList elements.
462 *
463 * @param data Unused.
464 * @return void* Allocated memory for LIST_NODE data.
465 */
466 ✗ void* eventListAlloc(const void* data) {
467 ✗ void* newElem = malloc(sizeof(long));
468 ✗ assertStreamPrint(NULL, newElem != NULL, "eventListAlloc: Out of memory");
469 ✗ return newElem;
470 }
471
472 /**
473 * @brief Free memory allocated with eventListAlloc.
474 *
475 * @param data Void pointer, representing index for new list element.
476 */
477 ✗ void eventListFree(void* data) {
478 ✗ free(data);
479 ✗ }
480
481 /**
482 * @brief Copy data of eventList elements.
483 *
484 * @param dest Void pointer of destination data, representing long index.
485 * @param src Void pointer of source data, representing long index.
486 */
487 ✗ void eventListCopy(void* dest, const void* src) {
488 long* dest_event = (long*) dest;
489 ✗ *dest_event = *((long*) src);
490 ✗ }
491