Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 122
Functions: 0.0% 0 / 0 / 8
Branches: 0.0% 0 / 0 / 96
OMCompiler/SimulationRuntime/c/simulation/solver/delay.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 #if !defined(OMC_NDELAY_EXPRESSIONS) || OMC_NDELAY_EXPRESSIONS>0
29
30 /*! \file delay.c
31 */
32
33 #include "delay.h"
34 #include "epsilon.h"
35 #include "../../util/omc_error.h"
36 #include "../../util/ringbuffer.h"
37 #include "../../openmodelica.h"
38
39 #include <stdio.h>
40 #include <stdlib.h>
41
42 /* Private prototypes */
43 void printDelayBuffer(void* data, int stream, void* elemPointer);
44
45
46 /**
47 * @brief Find row with greatest time that is smaller than or equal to 'time'.
48 *
49 * @param[in] time Time value to search for.
50 * @param[in] delayStruct Ringbuffer with stored delay values.
51 * Looks like a matrix with columns of type TIME_AND_VALUE.
52 * @return int Row with maximum time value smaller equal to time.
53 */
54 ✗ static int findTime(double time, RINGBUFFER *delayStruct)
55 {
56 ✗ int end = ringBufferLength(delayStruct);
57 int pos = 0;
58 double curTime;
59 TIME_AND_VALUE* bufferElem;
60
61 /* Check if ring buffer is valid */
62 ✗ assertStreamPrint(NULL, ringBufferLength(delayStruct) > 0, "delay: In function findTime\nEmpty ring buffer.");
63 ✗ bufferElem = getRingData(delayStruct, pos);
64 ✗ curTime = bufferElem->t;
65
66 /* If searched time is smaller then first element return first position */
67 ✗ if (time < curTime) {
68 return pos;
69 }
70
71 /* Search for time starting at begin of ring buffer */
72 ✗ while (pos < end-1) {
73 ✗ pos++;
74 ✗ bufferElem = getRingData(delayStruct, pos);
75 ✗ curTime = bufferElem->t;
76
77 ✗ if (curTime > time) {
78 pos--;
79 // Found time in previous step
80 break;
81 }
82 }
83 ✗ assertStreamPrint(NULL, pos < end, "delay: In function findTime\nCould not find time");
84
85 return pos;
86 }
87
88
89 /**
90 * @brief Look for events between `oldTime` and `newTime`
91 *
92 * @param[in] time Time value to search for.
93 * @param[in] delayStruct Ringbuffer with stored delay values.
94 * Looks like a matrix with columns of type TIME_AND_VALUE.
95 * @return modelica_boolean Boolean indicating if an event was found.
96 */
97 ✗ static modelica_boolean searchEvent(double time, RINGBUFFER *delayStruct)
98 {
99 ✗ int end = ringBufferLength(delayStruct);
100 int pos = 0;
101 double curTime, prevTime;
102 TIME_AND_VALUE* bufferElem;
103 modelica_boolean foundEvent = FALSE;
104
105 ✗ bufferElem = getRingData(delayStruct, pos);
106 ✗ curTime = bufferElem->t;
107
108 /* If searched time is smaller then first element we have no event */
109 ✗ if (time < curTime) {
110 return FALSE;
111 }
112
113 /* Search for time starting at begin of ring buffer */
114 ✗ while (pos < end-1) {
115 ✗ pos++;
116 ✗ bufferElem = getRingData(delayStruct, pos);
117 prevTime = curTime;
118 ✗ curTime = bufferElem->t;
119
120 /* Check for an event */
121 ✗ if (fabs(prevTime-curTime) < 1e-12) {
122 foundEvent = TRUE;
123 break;
124 }
125 ✗ if (curTime > time) {
126 break;
127 }
128 }
129
130 ✗ if (foundEvent) {
131 ✗ printRingBuffer(delayStruct, OMC_LOG_DEBUG, printDelayBuffer);
132 }
133
134 return foundEvent;
135 }
136
137
138 /**
139 * @brief Store expression value in delay.
140 *
141 * @param data Storing all simulation/model data.
142 * @param threadData Used for error handling.
143 * @param exprNumber Index of delay.
144 * @param exprValue Value to store in delay ringbuffer.
145 * @param delayTime Time to delay expValue.
146 * @param delayMax Maximum allowed delay time, defaults to delayTime.
147 */
148 ✗ void storeDelayedExpression(DATA* data, threadData_t *threadData, int exprNumber, double exprValue, double delayTime, double delayMax)
149 {
150 ✗ RINGBUFFER* delayStruct = data->simulationInfo->delayStructure[exprNumber];
151 int row;
152 ✗ int length = ringBufferLength(delayStruct);
153 ✗ double time = data->localData[0]->timeValue;
154 TIME_AND_VALUE tpl;
155 TIME_AND_VALUE* lastElem;
156
157 ✗ assertStreamPrint(threadData, exprNumber < data->modelData->nDelayExpressions, "storeDelayedExpression: invalid expression number %d", exprNumber);
158 ✗ assertStreamPrint(threadData, 0 <= exprNumber, "storeDelayedExpression: invalid expression number %d", exprNumber);
159
160 /* Check if time is greater equal then last stored time in delay structure */
161 /* ph: Is this needed because of event search? */
162 ✗ if (length > 0) {
163 ✗ lastElem = getRingData(delayStruct, length-1);
164 ✗ while (time < lastElem->t && length > 0) {
165 ✗ removeLastRingData(delayStruct,1);
166 ✗ length = ringBufferLength(delayStruct);
167 ✗ if (length > 0) {
168 ✗ lastElem = getRingData(delayStruct, length-1);
169 }
170 }
171 }
172
173 /* Check if (time,value) pair is already saved
174 * This should happen after an event was found and the event iteration finished */
175 ✗ if (length > 0) {
176 ✗ if (fabs(lastElem->t-time) < 1e-10 && fabs(lastElem->value-exprValue) < 1e-10) {
177 /* Dequeue no longer needed values from ring buffer */
178 ✗ row = findTime(time-delayTime+1e-10, delayStruct);
179 ✗ if (row > 0) {
180 ✗ dequeueNFirstRingDatas(delayStruct, row);
181 }
182 ✗ return;
183 }
184 }
185
186 /* Append expression value to delay ring buffer */
187 ✗ tpl.t = time;
188 ✗ tpl.value = exprValue;
189 ✗ appendRingData(delayStruct, &tpl);
190
191 /* Dequeue no longer needed values from ring buffer */
192 ✗ row = findTime(time-delayTime+DBL_EPSILON, delayStruct);
193 ✗ if (row > 0 && !searchEvent(time-delayTime+DBL_EPSILON, delayStruct)) {
194 ✗ dequeueNFirstRingDatas(delayStruct, row);
195 }
196
197 /* Debug print */
198 ✗ infoStreamPrint(OMC_LOG_DELAY, 0, "storeDelayed[%d] (%g,%g) position=%d", exprNumber, time, exprValue, ringBufferLength(delayStruct));
199 ✗ printRingBuffer(delayStruct, OMC_LOG_DELAY, printDelayBuffer);
200 }
201
202
203 /**
204 * @brief Evaluate delay expression.
205 *
206 * @param data Pointer to data.
207 * @param threadData Pointer to thread data.
208 * @param exprNumber Index of delay expression.
209 * @param exprValue Value of delay expression.
210 * @param delayTime Amount of time exprValue should be delayed.
211 * @param delayMax Maximum time to delay exprValue.
212 * @return double Return delayed value.
213 */
214 ✗ double delayImpl(DATA* data, threadData_t *threadData, int exprNumber, double exprValue, double delayTime, double delayMax)
215 {
216 ✗ RINGBUFFER* delayStruct = data->simulationInfo->delayStructure[exprNumber];
217 double timeStamp;
218 double time0, time1, value0, value1;
219 double timedif;
220 double dt1;
221 double retVal;
222 int i;
223 ✗ int length = ringBufferLength(delayStruct);
224 ✗ double time = data->localData[0]->timeValue;
225
226 ✗ infoStreamPrint(OMC_LOG_DELAY, 0, "delayImpl: exprNumber = %d, exprValue = %g, time = %g, delayTime = %g", exprNumber, exprValue, time, delayTime);
227
228 /* Check for errors */
229 ✗ assertStreamPrint(threadData, 0 <= exprNumber, "invalid exprNumber = %d", exprNumber);
230 ✗ assertStreamPrint(threadData, exprNumber < data->modelData->nDelayExpressions, "invalid exprNumber = %d", exprNumber);
231 ✗ assertStreamPrint(threadData, delayTime >= 0, "Negative delay requested: delayTime = %g", delayTime);
232 ✗ assertStreamPrint(threadData, delayTime >= DASSL_STEP_EPS, "delayImpl: delayTime is zero or too small.\n" \
233 "OpenModelica doesn't support delay operator with zero delay time.");
234 ✗ assertStreamPrint(threadData, delayTime <= delayMax, "Too large delay requested: delayTime = %g, delayMax = %g", delayTime, delayMax);
235
236 /* Return expression value before simulation start */
237 ✗ if (time <= data->simulationInfo->startTime) {
238 return exprValue;
239 }
240
241 /* Empty delay buffer at initialization phase */
242 ✗ if (length == 0) {
243 ✗ infoStreamPrint(OMC_LOG_EVENTS, 0, "delayImpl: Missing initial value, using argument value %g instead.", exprValue);
244 ✗ return exprValue;
245 }
246
247 /* Return oldest element in ring buffer */
248 ✗ if (time <= data->simulationInfo->startTime + delayTime) {
249 ✗ return ((TIME_AND_VALUE*)getRingData(delayStruct, 0))->value;
250 }
251 /* return expr(time-delayTime) */
252 else {
253 ✗ timeStamp = time - delayTime;
254
255 /* find the row for the lower limit */
256 ✗ if (timeStamp > ((TIME_AND_VALUE*)getRingData(delayStruct, length - 1))->t) {
257 /* delay between the last accepted time step and the current time */
258 ✗ time0 = ((TIME_AND_VALUE*)getRingData(delayStruct, length - 1))->t;
259 ✗ value0 = ((TIME_AND_VALUE*)getRingData(delayStruct, length - 1))->value;
260 time1 = time;
261 value1 = exprValue;
262 } else {
263 ✗ i = findTime(timeStamp, delayStruct);
264 ✗ assertStreamPrint(threadData, i < length, "%d = i < length = %d", i, length);
265 ✗ time0 = ((TIME_AND_VALUE*)getRingData(delayStruct, i))->t;
266 ✗ value0 = ((TIME_AND_VALUE*)getRingData(delayStruct, i))->value;
267
268 /* was it the last value? */
269 ✗ if (i+1 == length) {
270 return value0;
271 }
272 ✗ time1 = ((TIME_AND_VALUE*)getRingData(delayStruct, i+1))->t;
273 ✗ value1 = ((TIME_AND_VALUE*)getRingData(delayStruct, i+1))->value;
274 }
275 /* Return left value */
276 ✗ if (time0 == timeStamp) {
277 return value0;
278 }
279 /* Return right value */
280 ✗ else if (time1 == timeStamp) {
281 return value1;
282 }
283 /* linear interpolation */
284 /* FIXME instead of linear, do the same interpolation order as the integrator */
285 else {
286 ✗ timedif = time1 - time0;
287 ✗ dt1 = timeStamp - time0;
288 /* Exact when value0 == value1, unlike (value0*dt0 + value1*dt1)/timedif */
289 ✗ retVal = value0 + (value1 - value0) * (dt1 / timedif);
290 ✗ return retVal;
291 }
292 }
293 }
294
295
296 /**
297 * @brief Returns value of zero crossing at current simulation time.
298 *
299 * @param data Storing all simulation/model data.
300 * @param threadData Used for error handling.
301 * @param exprNumber Index of delay.
302 * @param relationIndex Index of relation used for zero crossing.
303 * @param delayTime Time to delay expValue.
304 * @return double Value of zeroCrossing at current simulation time.
305 */
306 ✗ double delayZeroCrossing(DATA* data, threadData_t *threadData, unsigned int exprNumber, unsigned int relationIndex, double delayTime)
307 {
308 ✗ RINGBUFFER* delayStruct = data->simulationInfo->delayStructure[exprNumber];
309 ✗ double zeroCrossingValue = data->simulationInfo->zeroCrossingsPre[relationIndex];
310 ✗ double time = data->localData[0]->timeValue;
311
312 ✗ if (ringBufferLength(delayStruct) == 0) {
313 return zeroCrossingValue;
314 }
315
316 /* Flip sign of ZC if an event was found */
317 ✗ if (searchEvent(time - delayTime, delayStruct)) {
318 ✗ return -zeroCrossingValue;
319 } else {
320 return zeroCrossingValue;
321 }
322 }
323
324
325 /**
326 * @brief Print transported quantity data to stream.
327 *
328 * Prints tuple (time, value).
329 *
330 * @param data Void pointer to bufferElemData.
331 * Will be casted to TIME_AND_VALUE*.
332 * @param stream Stream of OMC_LOG_STREAM type.
333 * @param elemPointer Address of element storing this data.
334 */
335 ✗ void printDelayBuffer(void* data, int stream, void* elemPointer)
336 {
337 TIME_AND_VALUE* bufferElemData = (TIME_AND_VALUE*) data;
338 ✗ infoStreamPrint(stream, 0, "%p: (%e,%e)", elemPointer, bufferElemData->t, bufferElemData->value);
339 ✗ }
340
341 #endif
342
343 /**
344 * @brief The delay buffers as flat words, for an FMU state.
345 *
346 * Per delay expression its length n, then n (time, value) pairs.
347 *
348 * @param data Runtime data struct.
349 * @param out Receives the words, or NULL to only count them.
350 * @return Number of words.
351 */
352 ✗ size_t delayStateWords(DATA* data, double* out)
353 {
354 size_t k = 0;
355 long i;
356 int j, n;
357
358 ✗ if (!data->simulationInfo->delayStructure) {
359 return 0;
360 }
361 ✗ for (i = 0; i < data->modelData->nDelayExpressions; i++) {
362 ✗ RINGBUFFER* rb = data->simulationInfo->delayStructure[i];
363 ✗ n = ringBufferLength(rb);
364 ✗ if (out) out[k] = n;
365 ✗ k++;
366 ✗ for (j = 0; j < n; j++) {
367 ✗ TIME_AND_VALUE* e = (TIME_AND_VALUE*) getRingData(rb, j);
368 ✗ if (out) {
369 ✗ out[k] = e->t;
370 ✗ out[k+1] = e->value;
371 }
372 ✗ k += 2;
373 }
374 }
375 return k;
376 }
377
378 /**
379 * @brief Inverse of delayStateWords.
380 *
381 * @param data Runtime data struct.
382 * @param w Words delayStateWords wrote.
383 * @param len Number of words available.
384 * @return Number of words read, or -1 if they are not what delayStateWords wrote.
385 */
386 ✗ long setDelayStateWords(DATA* data, const double* w, size_t len)
387 {
388 size_t k = 0;
389 long i;
390 int j, n;
391
392 ✗ if (!data->simulationInfo->delayStructure) {
393 return 0;
394 }
395 ✗ for (i = 0; i < data->modelData->nDelayExpressions; i++) {
396 ✗ RINGBUFFER* rb = data->simulationInfo->delayStructure[i];
397 ✗ if (k >= len) return -1;
398 ✗ n = (int) w[k++];
399 ✗ if (n < 0 || k + 2 * (size_t) n > len) return -1;
400 ✗ removeLastRingData(rb, ringBufferLength(rb));
401 ✗ for (j = 0; j < n; j++) {
402 TIME_AND_VALUE e;
403 ✗ e.t = w[k];
404 ✗ e.value = w[k+1];
405 ✗ k += 2;
406 ✗ appendRingData(rb, &e);
407 }
408 }
409 ✗ return (long) k;
410 }
411