Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 1.6% 7 / 0 / 440
Functions: 13.3% 2 / 1 / 16
Branches: 0.6% 2 / 0 / 340
OMCompiler/SimulationRuntime/c/simulation/solver/spatialDistribution.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 spatialDistribution.c
31 */
32
33 #include "spatialDistribution.h"
34 #include "../../util/omc_error.h"
35 #include "../../util/ringbuffer.h"
36 #include "../../openmodelica.h"
37 #include "epsilon.h"
38
39 #include <math.h>
40 #include <stdio.h>
41 #include <stdlib.h>
42
43
44 /**
45 * @brief Describing value z(x,t).
46 *
47 * See Modelica specification 3.7.2.2 spatialDistribution for details
48 * on transported quantity z(x,t).
49 * https://specification.modelica.org/maint/3.4/Ch3.html#spatialdistribution
50 */
51 typedef struct TRANSPORTED_QUANTITY_DATA {
52 double position; /* position x */
53 double value; /* transported quantity at position x */
54 } TRANSPORTED_QUANTITY_DATA;
55
56 /**
57 * @brief Saving an event at given position.
58 *
59 * The zero crossing function will return 0 on this event,
60 * zeroCrossValue until next event position
61 * and -1*zeroCrossValue before this event.
62 */
63 typedef struct TRANSPORTED_EVENT_DATA {
64 double position; /* position x */
65 double zeroCrossValue; /* Value of zero crossing at position x
66 * Either +1 or -1 */
67 } TRANSPORTED_EVENT_DATA;
68
69
70 /* Private function prototypes */
71 double interpolateTransportedQuantity(threadData_t *threadData, const TRANSPORTED_QUANTITY_DATA* leftData, const TRANSPORTED_QUANTITY_DATA* rightData, const double interpolationPos);
72 double extrapolateTransportedQuantity(threadData_t *threadData, const TRANSPORTED_QUANTITY_DATA* leftData, const TRANSPORTED_QUANTITY_DATA* rightData, const double extrapolationPos);
73 void addNewNodeSpatialDistribution(threadData_t *threadData, SPATIAL_DISTRIBUTION_DATA* spatialDistribution, int isPositiveVelocity, double position, double value, int isEvent);
74 int findOppositeEndSpatialDistribution(threadData_t *threadData, SPATIAL_DISTRIBUTION_DATA* spatialDistribution, double in0, double in1, double posX, int isPositiveVelocity, double* eventPreValue, double* outValue);
75 int pruneSpatialDistribution(threadData_t *threadData, SPATIAL_DISTRIBUTION_DATA* spatialDistribution, int isPositiveVelocity);
76
77
78 /* SPATIAL_EPS and SPATIAL_ZERO_DELTA_X are absolute, positions and values are
79 * not: scale them, clamped at 1 so nothing gets tighter than unscaled. */
80 static const double SPATIAL_EPS_ULPS = 8.0;
81
82 static double spatialScale(double a, double b) {
83 ✗ double scaleA = fabs(a);
84 ✗ double scaleB = fabs(b);
85 ✗ double scale = (scaleA > scaleB) ? scaleA : scaleB;
86 ✗ return (scale > 1.0) ? scale : 1.0;
87 }
88
89 static double spatialPosEps(double posA, double posB) {
90 ✗ return SPATIAL_EPS_ULPS * SPATIAL_EPS * spatialScale(posA, posB);
91 }
92
93 static double spatialValEps(double valA, double valB) {
94 ✗ return SPATIAL_EPS_ULPS * SPATIAL_EPS * spatialScale(valA, valB);
95 }
96
97 /* Never below the resolution of the position coordinate. */
98 static double spatialZeroDeltaX(double posA, double posB) {
99 double posEps = spatialPosEps(posA, posB);
100 ✗ return (SPATIAL_ZERO_DELTA_X > posEps) ? SPATIAL_ZERO_DELTA_X : posEps;
101 }
102
103 // ############################################################################
104 //
105 // Section for allocating/ deallocating spatial distribution data
106 //
107 // ############################################################################
108
109
110 /**
111 * @brief Allocates memory for spatial distribution structs.
112 *
113 * Returns pointer to array with allocated spatial distribution structs.
114 * To free memroy call freeSpatialDistribution.
115 *
116 * @param nSpatialDistributions Number of spacial distributions to be allocated.
117 * @return SPATIAL_DISTRIBUTION_DATA* Array with allocated spatial distributions.
118 */
119 1 SPATIAL_DISTRIBUTION_DATA* allocSpatialDistribution(unsigned int nSpatialDistributions) {
120 /* Debug info */
121 1 infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "Allocating memory for %i spatial distribution(s).", nSpatialDistributions);
122
123 /* Variables */
124 int i;
125 SPATIAL_DISTRIBUTION_DATA* spatialDistributionData;
126
127
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (nSpatialDistributions==0) {
128 return NULL;
129 }
130
131 ✗ spatialDistributionData = (SPATIAL_DISTRIBUTION_DATA*) calloc(nSpatialDistributions, sizeof(SPATIAL_DISTRIBUTION_DATA));
132
133 ✗ for(i=0; i<nSpatialDistributions; i++) {
134 ✗ spatialDistributionData[i].index = i;
135 ✗ spatialDistributionData[i].isInitialized = 0 /* false */;
136 ✗ spatialDistributionData[i].startPosXSet = 0 /* false */;
137 ✗ spatialDistributionData[i].startPosX = 0.0 /* false */;
138 ✗ spatialDistributionData[i].oldPosX = 0.0;
139 ✗ spatialDistributionData[i].transportedQuantity = allocDoubleEndedList(sizeof(TRANSPORTED_QUANTITY_DATA)); /* empty double ended list */
140 ✗ spatialDistributionData[i].storedEvents = allocDoubleEndedList(sizeof(TRANSPORTED_EVENT_DATA)); /* empty double ended list */
141 ✗ spatialDistributionData[i].lastStoredEventValue = 0;
142 ✗ spatialDistributionData[i].nWarningsRemovedEvents = 0;
143 ✗ spatialDistributionData[i].nWarningsOutputEvents = 0;
144 }
145
146 return spatialDistributionData;
147 }
148
149
150 /**
151 * @brief Frees array of spatial distributions.
152 *
153 * @param spatialDistributionData Array with spatial distribution of length nSpatialDistributions.
154 * @param nSpatialDistributions Length of spatialDistributionData.
155 */
156 1 void freeSpatialDistribution(SPATIAL_DISTRIBUTION_DATA* spatialDistributionData, unsigned int nSpatialDistributions) {
157 /* Debug info */
158 1 infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "Freeing %i spatial distribution(s).", nSpatialDistributions);
159
160 /* Variables */
161 int i;
162
163
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(i=0; i<nSpatialDistributions; i++) {
164 ✗ freeDoubleEndedList(spatialDistributionData[i].transportedQuantity);
165 ✗ freeDoubleEndedList(spatialDistributionData[i].storedEvents);
166 }
167 1 }
168
169
170 /**
171 * @brief Initializes transportedQuantity of single spacial distribution.
172 *
173 * Spatial distribution array data->simulationInfo->spatialDistributionData has
174 * to be allocated before using allocSpatialDistribution.
175 *
176 * @param data Data
177 * @param threadData threadDate for error handling
178 * @param index Index of spatial distribution, has to match position data->simulationInfo->spatialDistributionData[index].
179 * @param initialPoints Array with initial points.
180 * Is ordered from 0.0 = initialPoints[0] < initialPoints[i] < initialPoints[length] = 1.0
181 * @param initialValues Array with initial values at initial points.
182 * @param length Length of arrays initialPoints and initialValues.
183 */
184 ✗ void initSpatialDistribution(DATA* data, threadData_t* threadData, unsigned int index, real_array* initialPoints, real_array* initialValues, unsigned int length) {
185 /* Debug info */
186 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 1, "Initializing spatial distributions (index=%i)", index);
187
188 /* Variables */
189 int i;
190 SPATIAL_DISTRIBUTION_DATA* spatialDistributionData;
191 DOUBLE_ENDED_LIST* transportedQuantityList;
192 TRANSPORTED_QUANTITY_DATA tmpData;
193 TRANSPORTED_EVENT_DATA eventData;
194 int numSamePos = 0;
195 double lastZeroCrossValue = -1;
196 ✗ modelica_real* initPnts = (modelica_real *) initialPoints->data;
197 ✗ modelica_real* initVals = (modelica_real *) initialValues->data;
198
199 /* Error checking */
200 ✗ if (fabs(initPnts[0]) > SPATIAL_EPS ) {
201 ✗ errorStreamPrint(OMC_LOG_STDOUT, 1, "Initialization of spatial distribution with index %i failed.", index);
202 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "initialPoints[0] = %e is not zero.", initPnts[0]);
203 ✗ messageClose(OMC_LOG_STDOUT);
204 ✗ omc_throw_function(threadData);
205 }
206 ✗ else if (fabs(initPnts[length-1] - 1.0) > SPATIAL_EPS) {
207 ✗ errorStreamPrint(OMC_LOG_STDOUT, 1, "Initialization of spatial distribution with index %i failed.", index);
208 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "initialPoints[end] = %e is not one.", initPnts[length-1]);
209 ✗ messageClose(OMC_LOG_STDOUT);
210 ✗ omc_throw_function(threadData);
211 }
212 ✗ for (i=0; i<length-1; i++) {
213 ✗ if (initPnts[i] > initPnts[i+1]) {
214 ✗ errorStreamPrint(OMC_LOG_STDOUT, 1, "Initialization of spatial distribution with index %i failed.", index);
215 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "initialPoints[%i] > initialPoints[%i]", i, i+1);
216 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "%f > %f", initPnts[i], initPnts[i+1]);
217 ✗ messageClose(OMC_LOG_STDOUT);
218 ✗ omc_throw_function(threadData);
219 }
220 }
221 ✗ spatialDistributionData = &(data->simulationInfo->spatialDistributionData[index]);
222 ✗ assertStreamPrint(threadData, 1 != spatialDistributionData->isInitialized, "SpatialDistribution was allready allocated!");
223
224 /* Initialize quantity list */
225 ✗ transportedQuantityList = spatialDistributionData->transportedQuantity;
226 ✗ for (i=0; i<length-1; i++) {
227 ✗ tmpData.position = initPnts[i];
228 ✗ tmpData.value = initVals[i];
229 ✗ pushBackDoubleEndedList(transportedQuantityList, (const void*) &tmpData);
230 ✗ if (initPnts[i] == initPnts[i+1]) {
231 ✗ numSamePos += 1;
232 ✗ if (numSamePos > 1) {
233 ✗ errorStreamPrint(OMC_LOG_STDOUT, 1, "Initialization of spatial distribution with index %i failed.", index);
234 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "initialPoints[%i] = initialPoints[%i] = initialPoints[%i]", i-1, i, i+1);
235 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Only events with one pre-value and one value are allowed.");
236 ✗ messageClose(OMC_LOG_STDOUT);
237 ✗ omc_throw_function(threadData);
238 }
239 ✗ eventData.position = initPnts[i];
240 ✗ lastZeroCrossValue = lastZeroCrossValue*(-1);
241 ✗ eventData.zeroCrossValue = lastZeroCrossValue;
242 ✗ pushBackDoubleEndedList(spatialDistributionData->storedEvents, (const void*) &eventData);
243 } else {
244 numSamePos = 0;
245 }
246 }
247 ✗ tmpData.position = initPnts[length-1];
248 ✗ tmpData.value = initVals[length-1];
249 ✗ pushBackDoubleEndedList(transportedQuantityList, (const void*) &tmpData);
250
251 ✗ spatialDistributionData->isInitialized = 1 /* true */;
252
253 /* Debug info */
254 ✗ doubleEndedListPrint(transportedQuantityList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
255 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "List of events");
256 ✗ doubleEndedListPrint(spatialDistributionData->storedEvents, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
257 ✗ messageClose(OMC_LOG_SPATIALDISTR);
258 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "Finished initializing spatial distribution (index=%i)", index);
259 ✗ }
260
261
262 // ############################################################################
263 //
264 // Section for evaluating spatialDistribution operator
265 //
266 // ############################################################################
267
268
269 /**
270 * @brief Shift posX so that the operator internally starts at x = 0.
271 *
272 * The spatialDistribution operator only depends on the change of the spatial
273 * coordinate x (the transport distance), not on its absolute value, and its
274 * initial profile (initialPoints/initialValues) is stored assuming x(t0) = 0.
275 * The value of x at the very first call is captured once and subtracted from
276 * every subsequent posX, so a model where x has a nonzero start value behaves
277 * exactly like one starting at x = 0 (and no longer triggers a spurious
278 * "x got reinitialized during an event" error at the initial event).
279 *
280 * @param spatialDistribution Spatial distribution to shift for.
281 * @param posX Value of position x.
282 * @return double posX relative to its value at the first call.
283 */
284 static double shiftToStartPosX(SPATIAL_DISTRIBUTION_DATA* spatialDistribution, double posX) {
285 ✗ if (!spatialDistribution->startPosXSet) {
286 ✗ spatialDistribution->startPosX = posX;
287 ✗ spatialDistribution->startPosXSet = 1 /* true */;
288 }
289 ✗ return posX - spatialDistribution->startPosX;
290 }
291
292 ✗ static void warnStepSizeTooBig(DATA* data, unsigned long* nDisplayed, const char* what, unsigned int index, int nEvents) {
293 ✗ unsigned long maxWarnDisplays = data->simulationInfo->maxWarnDisplays;
294
295 ✗ if (++*nDisplayed > maxWarnDisplays || !OMC_ACTIVE_WARNING_STREAM(OMC_LOG_STDOUT)) {
296 return;
297 }
298 ✗ warningStreamPrint(OMC_LOG_STDOUT, 1, "%s more then one event from spatialDistribution. Step size to big!", what);
299 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "time: %f, spatialDistribution index: %i, number of events: %i", data->localData[0]->timeValue, index, nEvents);
300 ✗ messageCloseWarning(OMC_LOG_STDOUT);
301 ✗ if (*nDisplayed == maxWarnDisplays) {
302 ✗ warningStreamPrintLimitReached(OMC_LOG_STDOUT, 0, maxWarnDisplays);
303 }
304 }
305
306
307 /**
308 * @brief Store spatial distribution data for an accepted step.
309 *
310 * @param data Data
311 * @param threadData Thread data for error handling
312 * @param index Index of spatial distribution.
313 * @param in0 First input to spatial distribution.
314 * @param in1 Second input to spatial distribution
315 * @param posX Value of position x.
316 * @param isPositiveVelocity Boolean describing if velocity v is positive (>=0).
317 * Velocity v is `v:=der(x)`.
318 */
319 ✗ void storeSpatialDistribution(DATA* data, threadData_t *threadData, unsigned int index, double in0, double in1, double posX, int isPositiveVelocity) {
320 /* Variables */
321 SPATIAL_DISTRIBUTION_DATA* spatialDistribution;
322 DOUBLE_ENDED_LIST* transportedQuantityList;
323 DOUBLE_ENDED_LIST* storedEventsList;
324 int walkedOverEvents = 0;
325 double deltaX, realDirection;
326
327 /* Access spatialDistribution */
328 ✗ spatialDistribution = &(data->simulationInfo->spatialDistributionData[index]);
329 ✗ transportedQuantityList = spatialDistribution->transportedQuantity;
330 ✗ storedEventsList = spatialDistribution->storedEvents;
331
332 /* Shift x so the operator starts at x = 0 (only the change of x matters) */
333 posX = shiftToStartPosX(spatialDistribution, posX);
334
335 /* Debug log */
336 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 1, "Calling storeSpatialDistribution (index=%i, time=%e)", index, data->localData[0]->timeValue);
337 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "spatialDistribution(%f, %f, %f, %s)", in0, in1, posX, isPositiveVelocity?"true":"false");
338 ✗ doubleEndedListPrint(transportedQuantityList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
339 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "List of events");
340 ✗ doubleEndedListPrint(storedEventsList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
341
342 ✗ if (data->simulationInfo->discreteCall) {
343 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Discrete call of storeSpatialDistribution");
344 ✗ omc_throw_function(threadData);
345 }
346
347 /* Get deltaX */
348 ✗ deltaX = spatialDistribution->oldPosX - posX;
349 ✗ if (deltaX > 0) {
350 realDirection = 1 /* positive */;
351 ✗ } else if (deltaX < 0) {
352 realDirection = -1 /* negative */;
353 deltaX = -deltaX;
354 } else {
355 realDirection = 0 /* standing still */;
356 }
357
358 /* deltaX = oldPosX - posX, so realDirection > 0 means x decreased */
359 ✗ if (realDirection > 0) {
360 isPositiveVelocity = 0;
361 ✗ } else if (realDirection < 0) {
362 isPositiveVelocity = 1;
363 }
364 /* Event nodes go onto the edge position: prune computes edge positions as
365 * edge +/- 1, so -posX can be an ulp on the wrong side of the edge. */
366 ✗ if (isPositiveVelocity) {
367 ✗ TRANSPORTED_QUANTITY_DATA* front = (TRANSPORTED_QUANTITY_DATA*) firstDataDoubleEndedList(transportedQuantityList);
368 ✗ if (fabs(-posX - front->position) < spatialPosEps(-posX, front->position)) {
369 ✗ if (fabs(front->value - in0) > spatialValEps(front->value, in0)) {
370 ✗ addNewNodeSpatialDistribution(threadData, spatialDistribution, isPositiveVelocity, front->position, in0, 1 /* true */);
371 }
372 } else {
373 ✗ addNewNodeSpatialDistribution(threadData, spatialDistribution, isPositiveVelocity, -posX, in0, 0 /* false */);
374 }
375 } else {
376 ✗ TRANSPORTED_QUANTITY_DATA* last = (TRANSPORTED_QUANTITY_DATA*) lastDataDoubleEndedList(transportedQuantityList);
377 ✗ if (fabs(-posX+1 - last->position) < spatialPosEps(-posX+1, last->position)) {
378 ✗ if (fabs(last->value - in1) > spatialValEps(last->value, in1)) {
379 ✗ addNewNodeSpatialDistribution(threadData, spatialDistribution, isPositiveVelocity, last->position, in1, 1 /* true */);
380 }
381 } else {
382 ✗ addNewNodeSpatialDistribution(threadData, spatialDistribution, isPositiveVelocity, -posX+1, in1, 0 /* false */);
383 }
384 }
385
386 /* Remove nodes that droppen of spatial distribution */
387 ✗ walkedOverEvents = pruneSpatialDistribution(threadData, spatialDistribution, isPositiveVelocity);
388 ✗ if (walkedOverEvents > 1) {
389 ✗ warnStepSizeTooBig(data, &spatialDistribution->nWarningsRemovedEvents, "Removed", index, walkedOverEvents);
390 }
391
392 /* Update oldPosX */
393 ✗ spatialDistribution->oldPosX = posX;
394 ✗ messageClose(OMC_LOG_SPATIALDISTR);
395 ✗ return;
396 }
397
398
399 /**
400 * @brief Evaluate spatialDistribution operator.
401 *
402 * (out0, out1) = spatialDistribution (in0, in1, posX, isPositiveVelocity)
403 * If an event was outputted integrator needs to iterate.
404 * Doesn't store in0 or in1 because this function doesn't know if the step will be accepted.
405 *
406 * @param data Data
407 * @param threadData Thread data for error handling
408 * @param index Index of spatial distribution.
409 * @param in0 First input to spatial distribution.
410 * @param in1 Second input to spatial distribution
411 * @param posX Value of position x.
412 * @param isPositiveVelocity Boolean describing if velocity v is positive (>=0).
413 * Velocity v is `v:=der(x)`.
414 * @param out1 Second output of spatial distribution.
415 * @return double out0, first output of spatial distribution.
416 */
417 ✗ double spatialDistribution(DATA* data, threadData_t *threadData, unsigned int index, double in0, double in1, double posX, int isPositiveVelocity, double* out1) {
418 /* Variables */
419 SPATIAL_DISTRIBUTION_DATA* spatialDistribution;
420 DOUBLE_ENDED_LIST* transportedQuantityList;
421 DOUBLE_ENDED_LIST_NODE* firstNode;
422 DOUBLE_ENDED_LIST_NODE* lastNode;
423 TRANSPORTED_QUANTITY_DATA* firstNodeData;
424 TRANSPORTED_QUANTITY_DATA* secondNodeData;
425 TRANSPORTED_QUANTITY_DATA* lastNodeData;
426 TRANSPORTED_QUANTITY_DATA* forelastNodeData;
427 int walkedOverEvents;
428 int realDirection;
429 int jumped = 0;
430 double deltaX;
431 ✗ double eventPreValue = NAN; /* only written if an event is walked over */
432 double outValue;
433 double out0; /* First output variable */
434 double out1Val; /* Second output variable, only written to *out1 if out1 != NULL */
435
436 /* Access spatialDistribution */
437 ✗ spatialDistribution = &(data->simulationInfo->spatialDistributionData[index]);
438 ✗ transportedQuantityList = spatialDistribution->transportedQuantity;
439
440 /* Shift x so the operator starts at x = 0 (only the change of x matters) */
441 posX = shiftToStartPosX(spatialDistribution, posX);
442
443 /* Debug log */
444 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 1, "Calling spatialDistribution (index=%i, time=%e)", index, data->localData[0]->timeValue);
445 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "(out0,out1) = spatialDistribution(in0=%f, in1=%f, x=%f, isPositiveVelocity=%s)", in0, in1, posX, isPositiveVelocity?"true":"false");
446 ✗ doubleEndedListPrint(transportedQuantityList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
447
448 /* Get deltaX */
449 ✗ deltaX = spatialDistribution->oldPosX - posX;
450 ✗ if (deltaX > 0) {
451 realDirection = 1 /* positive */;
452 ✗ } else if (deltaX < 0) {
453 realDirection = -1 /* negative */;
454 ✗ deltaX = -deltaX;
455 } else {
456 realDirection = 0 /* standing still */;
457 }
458
459 ✗ if (deltaX > spatialZeroDeltaX(spatialDistribution->oldPosX, posX) &&
460 ✗ ((isPositiveVelocity && realDirection > 0) || (!isPositiveVelocity && realDirection < 0))) {
461 ✗ isPositiveVelocity = !isPositiveVelocity;
462 jumped = 1 /* true */;
463 }
464
465 /* Check if x was reinitialized */
466 ✗ if (deltaX > spatialZeroDeltaX(spatialDistribution->oldPosX, posX) && data->simulationInfo->discreteCall) {
467 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "x got reinitialized during an event at time %f. OpenModelica can't handle that.", data->localData[0]->timeValue);
468 ✗ omc_throw_function(threadData);
469 }
470
471 /* Special case: Zero progress */
472 ✗ if (deltaX < spatialPosEps(spatialDistribution->oldPosX, posX)) {
473 ✗ firstNodeData = (TRANSPORTED_QUANTITY_DATA*) firstDataDoubleEndedList(transportedQuantityList);
474 ✗ lastNodeData = (TRANSPORTED_QUANTITY_DATA*) lastDataDoubleEndedList(transportedQuantityList);
475 ✗ out0 = firstNodeData->value;
476 ✗ out1Val = lastNodeData->value;
477 ✗ if (out1 != NULL) {
478 ✗ *out1 = out1Val;
479 }
480 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "(out0,out1) = (%f, %f)", out0, out1Val);
481 ✗ messageClose(OMC_LOG_SPATIALDISTR);
482 ✗ return out0;
483 }
484
485 /* Get value of ou0/out1 by walkling over list */
486 ✗ walkedOverEvents = findOppositeEndSpatialDistribution(threadData, spatialDistribution, in0, in1, posX, isPositiveVelocity, &eventPreValue, &outValue);
487
488 /* Handle events that would come out of spatialDistribution */
489 ✗ if (walkedOverEvents > 1) {
490 ✗ warnStepSizeTooBig(data, &spatialDistribution->nWarningsOutputEvents, "Need to output", index, walkedOverEvents);
491 }
492 ✗ if (walkedOverEvents>0 && !data->simulationInfo->discreteCall && !isnan(eventPreValue)) {
493 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "Found event in spatial distribution at time %f", data->localData[0]->timeValue);
494 ✗ outValue = eventPreValue;
495 }
496
497 /* Extrapolate return values to break up quasi-loop with inputs */
498 ✗ firstNodeData = (TRANSPORTED_QUANTITY_DATA*) firstDataDoubleEndedList(transportedQuantityList);
499 ✗ secondNodeData = dataDoubleEndedList(getNextNodeDoubleEndedList(getFirstNodeDoubleEndedList(transportedQuantityList)));
500 ✗ lastNodeData = (TRANSPORTED_QUANTITY_DATA*) lastDataDoubleEndedList(transportedQuantityList);
501 ✗ forelastNodeData = dataDoubleEndedList(getPreviousNodeDoubleEndedList(getLastNodeDoubleEndedList(transportedQuantityList)));
502 /* jumped only suppresses the extrapolation: in0/in1 must not be used, the
503 * velocity sign that selects between them is what is in doubt here. */
504 ✗ if (isPositiveVelocity) {
505 ✗ if (!jumped && deltaX > spatialPosEps(spatialDistribution->oldPosX, posX) &&
506 ✗ fabs(firstNodeData->position-secondNodeData->position) > spatialPosEps(firstNodeData->position, secondNodeData->position)) {
507 ✗ out0 = extrapolateTransportedQuantity(threadData, firstNodeData, secondNodeData, -posX);
508 } else {
509 ✗ out0 = firstNodeData->value;
510 }
511 ✗ out1Val = outValue;
512 } else {
513 ✗ out0 = outValue;
514 ✗ if (!jumped && deltaX > spatialPosEps(spatialDistribution->oldPosX, posX) &&
515 ✗ fabs(forelastNodeData->position-lastNodeData->position) > spatialPosEps(forelastNodeData->position, lastNodeData->position)) {
516 ✗ out1Val = extrapolateTransportedQuantity(threadData, forelastNodeData, lastNodeData, -posX+1);
517 } else {
518 ✗ out1Val = lastNodeData->value;
519 }
520 }
521
522 ✗ if (out1 != NULL) {
523 ✗ *out1 = out1Val;
524 }
525
526 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "(out0,out1) = (%f, %f)", out0, out1Val);
527 ✗ messageClose(OMC_LOG_SPATIALDISTR);
528 ✗ return out0;
529 }
530
531
532 // ############################################################################
533 //
534 // Section for evaluating spatialDistribution zero-crossing function
535 //
536 // ############################################################################
537
538
539 /**
540 * @brief Returns value of zero crossing at position x.
541 *
542 * zeroCross(x):= -1 if there are no events or before the first event.
543 * Otherwise zeroCross(x):=(-1)*zeroCross(x_E), where x_E is the position of the nearest event with bigger position.
544 * If there is no event with bigger position zeroCross(x):=zeroCross(x_E) where x_E is the event with the biggest position.
545 *
546 * @param data Data
547 * @param threadData threadDate for error handling
548 * @param index Index of spatial distribution, has to match position data->simulationInfo->spatialDistributionData[index].
549 * @param posX Value of position x.
550 * @param isPositiveVelocity Unused
551 * @return double Value of zeroCrossing at position posX.
552 *
553 * Event positions are compared with the absolute SPATIAL_EPS: a wider tolerance
554 * flips the value before the event is reached, leaving no sign change to find.
555 */
556 ✗ double spatialDistributionZeroCrossing(DATA* data, threadData_t *threadData, unsigned int index, unsigned int relationIndex, double posX, int isPositiveVelocity) {
557 /* Variables */
558 SPATIAL_DISTRIBUTION_DATA* spatialDistribution;
559 DOUBLE_ENDED_LIST* storedEventsList;
560 DOUBLE_ENDED_LIST_NODE* currentNode;
561 TRANSPORTED_EVENT_DATA* currentNodeData;
562 double zeroCrossingValue = -1;
563 double prevPosition, prevValue;
564
565 /* Access spatialDistribution */
566 ✗ spatialDistribution = &(data->simulationInfo->spatialDistributionData[index]);
567 ✗ storedEventsList = spatialDistribution->storedEvents;
568
569 /* Shift x so the operator starts at x = 0 (only the change of x matters).
570 * Do NOT capture the start position here: the zero-crossing function is
571 * evaluated unconditionally by the solver, also while the operator is frozen
572 * inside an inactive if-branch. Capturing the start position here would mark
573 * the operator as started too early and make the guarded storeSpatialDistribution/
574 * spatialDistribution calls see a spurious jump in x (#16099). */
575 ✗ if (spatialDistribution->startPosXSet) {
576 ✗ posX = posX - spatialDistribution->startPosX;
577 } else {
578 /* Not started: its zero point will be the x of its first call, so report the
579 * value it will have then (x = 0 in the operator's own coordinate). Using the
580 * x the inactive branch is not reading would let activating the branch flip
581 * the crossing, reporting an event for a discontinuity that has not moved. */
582 posX = 0.0;
583 }
584
585 ✗ if (doubleEndedListLen(storedEventsList) == 0) {
586 ✗ zeroCrossingValue = data->simulationInfo->zeroCrossingsPre[relationIndex];
587 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "spatialDistributionZeroCrossing(%e) = %e (no stored events, returning previous value)", posX, zeroCrossingValue);
588 ✗ return zeroCrossingValue;
589 }
590
591 ✗ if (isPositiveVelocity) {
592 ✗ currentNode = getLastNodeDoubleEndedList(storedEventsList);
593 ✗ currentNodeData = dataDoubleEndedList(currentNode);
594 // -posX+1 is behind last event
595 ✗ if (currentNodeData->position < -posX+1 ) {
596 ✗ zeroCrossingValue = -currentNodeData->zeroCrossValue;
597 } else {
598 ✗ while (currentNode != NULL) {
599 // Am I on an event?
600 ✗ if (fabs(currentNodeData->position+posX-1) <= SPATIAL_EPS) {
601 ✗ zeroCrossingValue = -currentNodeData->zeroCrossValue;
602 ✗ break;
603 }
604
605 prevPosition = currentNodeData->position;
606 ✗ prevValue = currentNodeData->zeroCrossValue;
607 ✗ currentNode = getPreviousNodeDoubleEndedList(currentNode);
608 // Did I walk over the first element in the list?
609 ✗ if (currentNode==NULL) {
610 zeroCrossingValue = prevValue; /* prevValue value of first list element */
611 break;
612 }
613 ✗ currentNodeData = dataDoubleEndedList(currentNode);
614
615 // Are we between two events?
616 ✗ if (currentNodeData->position < -posX+1 && -posX+1 < prevPosition) {
617 zeroCrossingValue = prevValue;
618 break;
619 }
620 }
621 }
622 } else {
623 ✗ currentNode = getFirstNodeDoubleEndedList(storedEventsList);
624 ✗ currentNodeData = dataDoubleEndedList(currentNode);
625 // -posX is before first event
626 ✗ if (currentNodeData->position > -posX ) {
627 ✗ zeroCrossingValue = currentNodeData->zeroCrossValue;
628 } else {
629 ✗ while (currentNode != NULL) {
630 // Am I on an event?
631 ✗ if (fabs(currentNodeData->position+posX) <= SPATIAL_EPS) {
632 ✗ zeroCrossingValue = -currentNodeData->zeroCrossValue;
633 ✗ break;
634 }
635
636 prevPosition = currentNodeData->position;
637 ✗ prevValue = currentNodeData->zeroCrossValue;
638 ✗ currentNode = getNextNodeDoubleEndedList(currentNode);
639 // Did I walk over the first element in the list?
640 ✗ if (currentNode==NULL) {
641 ✗ zeroCrossingValue = -prevValue; /* prevValue value of first list element */
642 ✗ break;
643 }
644 ✗ currentNodeData = dataDoubleEndedList(currentNode);
645
646 // Are we between two events?
647 ✗ if (currentNodeData->position > -posX && -posX > prevPosition) {
648 ✗ zeroCrossingValue = -prevValue;
649 ✗ break;
650 }
651 }
652 }
653 }
654
655
656 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "List of events for spatialDistributionZeroCrossing(%e) = %e", posX, zeroCrossingValue);
657 ✗ doubleEndedListPrint(storedEventsList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
658
659 ✗ return zeroCrossingValue;
660 }
661
662
663 // ############################################################################
664 //
665 // Section for "small" helper functions
666 //
667 // ############################################################################
668
669
670 /**
671 * @brief Linear interpolation between left and right (position,value) pair at given position.
672 *
673 * leftData->position < interpolationPos < rightData->position must hold.
674 *
675 * @param leftData Left (position,value) pair
676 * @param rightData Right (position,value) pair
677 * @param interpolationPos Position where to interpolate.
678 * @return double Interpolated value
679 */
680 ✗ double interpolateTransportedQuantity(threadData_t *threadData, const TRANSPORTED_QUANTITY_DATA* leftData, const TRANSPORTED_QUANTITY_DATA* rightData, const double interpolationPos) {
681 double leftPosition, rightPosition;
682 double leftValue, rightValue;
683 double distPos;
684 double interpolatedValue;
685
686 ✗ leftPosition = leftData->position;
687 ✗ leftValue = leftData->value;
688 ✗ rightPosition = rightData->position;
689 ✗ rightValue = rightData->value;
690 ✗ distPos = rightPosition - leftPosition;
691
692 ✗ assertStreamPrint(threadData, distPos > 0, "interpolateTransportedQuantity: wrong order or same position!");
693
694 ✗ interpolatedValue = leftValue * ((rightPosition-interpolationPos)/distPos)
695 ✗ + rightValue * ((interpolationPos-leftPosition)/distPos);
696
697 ✗ return interpolatedValue;
698 }
699
700
701 /**
702 * @brief Linear extrapolation at given position.
703 *
704 * @param leftData Left (position,value) pair
705 * @param rightData Right (position,value) pair
706 * @param extrapolationPos Position where to interpolate.
707 * @return double Extrapolated value.
708 */
709 ✗ double extrapolateTransportedQuantity(threadData_t *threadData, const TRANSPORTED_QUANTITY_DATA* leftData, const TRANSPORTED_QUANTITY_DATA* rightData, const double extrapolationPos) {
710 double leftPosition, rightPosition;
711 double leftValue, rightValue;
712 double distPos;
713 double extrapolatedValue;
714
715 ✗ leftPosition = leftData->position;
716 ✗ leftValue = leftData->value;
717 ✗ rightPosition = rightData->position;
718 ✗ rightValue = rightData->value;
719 ✗ distPos = rightPosition - leftPosition;
720
721 ✗ assertStreamPrint(threadData, distPos > 0, "interpolateTransportedQuantity: wrong order or same position!");
722
723 ✗ extrapolatedValue = leftValue + (rightValue-leftValue)/(distPos) * (extrapolationPos - leftPosition);
724 ✗ return extrapolatedValue;
725 }
726
727
728 /**
729 * @brief Adding new pair (position, value) to front or back of spatial distribution.
730 *
731 * For positive velocity add at frond, else at back.
732 * If this node is an event node add an event to stored events list as well.
733 *
734 * @param transportedQuantityList Double ended list representing spatial distribution.
735 * @param front Boolean value if node should be added at the front (true) or the end (false).
736 * @param position Position of new node.
737 * @param value Value of new node.
738 * @param isEvent Boolean value if new node is an event node.
739 */
740 ✗ void addNewNodeSpatialDistribution(threadData_t *threadData, SPATIAL_DISTRIBUTION_DATA* spatialDistribution, int front, double position, double value, int isEvent) {
741 /* Variables */
742 ✗ DOUBLE_ENDED_LIST* transportedQuantityList = spatialDistribution->transportedQuantity;
743 ✗ DOUBLE_ENDED_LIST* storedEventsList = spatialDistribution->storedEvents;
744 TRANSPORTED_QUANTITY_DATA newNodeData;
745 TRANSPORTED_EVENT_DATA newEventNodeData;
746
747 /* New node */
748 ✗ newNodeData.position = position;
749 ✗ newNodeData.value = value;
750 ✗ newEventNodeData.position = position;
751
752 /* Add node to transported quantity list */
753 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "Adding (%e,%e) at %s.", newNodeData.position, newNodeData.value, front?"front":"back");
754 ✗ if (front) {
755 // Make sure new first node is smaller then previous first node
756 ✗ TRANSPORTED_QUANTITY_DATA* oldFront = (TRANSPORTED_QUANTITY_DATA*) firstDataDoubleEndedList(transportedQuantityList);
757 ✗ assertStreamPrint(threadData, position<=oldFront->position, "New front position is not smaller then previous first node.");
758 ✗ pushFrontDoubleEndedList(transportedQuantityList, (const void*) &newNodeData);
759 } else {
760 // Make sure new first node is smaller then previous first node
761 ✗ TRANSPORTED_QUANTITY_DATA* oldEnd = (TRANSPORTED_QUANTITY_DATA*) lastDataDoubleEndedList(transportedQuantityList);
762 ✗ assertStreamPrint(threadData, position>=oldEnd->position, "New end position is not bigger then previous last node.");
763 ✗ pushBackDoubleEndedList(transportedQuantityList, (const void*) &newNodeData);
764 }
765
766 /* Add event to stored event list */
767 ✗ if (isEvent == 1) {
768 ✗ if (front) {
769 ✗ if (doubleEndedListLen(storedEventsList) == 0) {
770 ✗ if (spatialDistribution->lastStoredEventValue==0) {
771 ✗ newEventNodeData.zeroCrossValue = 1;
772 } else {
773 ✗ newEventNodeData.zeroCrossValue = -spatialDistribution->lastStoredEventValue;
774 }
775 } else {
776 // Make sure new first node is smaller then previous first node
777 ✗ TRANSPORTED_EVENT_DATA* oldEventFront = (TRANSPORTED_EVENT_DATA*) firstDataDoubleEndedList(storedEventsList);
778 ✗ assertStreamPrint(threadData, position<=oldEventFront->position, "New front position is not smaller then previous first event node.");
779 ✗ newEventNodeData.zeroCrossValue = oldEventFront->zeroCrossValue*(-1);
780 }
781 ✗ pushFrontDoubleEndedList(storedEventsList, (const void*) &newEventNodeData);
782 } else {
783 ✗ if (doubleEndedListLen(storedEventsList) == 0) {
784 ✗ newEventNodeData.zeroCrossValue = 1;
785 } else {
786 // Make sure new first node is smaller then previous first node
787 ✗ TRANSPORTED_EVENT_DATA* oldEventEnd = (TRANSPORTED_EVENT_DATA*) lastDataDoubleEndedList(storedEventsList);
788 ✗ assertStreamPrint(threadData, position>=oldEventEnd->position, "New end position is not bigger then previous last event node.");
789 ✗ newEventNodeData.zeroCrossValue = oldEventEnd->zeroCrossValue*(-1);
790 }
791 ✗ pushBackDoubleEndedList(storedEventsList, (const void*) &newEventNodeData);
792 }
793 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "Adding event (%e,%e) at %s.", newEventNodeData.position, newEventNodeData.zeroCrossValue, front?"front":"back");
794 }
795
796 /* Debug prints */
797 ✗ doubleEndedListPrint(transportedQuantityList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
798 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "List of events");
799 ✗ doubleEndedListPrint(storedEventsList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
800 ✗ }
801
802
803 /**
804 * @brief Gets value from opposite end of list.
805 *
806 * @param transportedQuantityList Double ended list containing spatial distribution.
807 * @param isPositiveVelocity Boolean describing if velocity v is positive (>=0).
808 * Velocity v is `v:=der(x)`.
809 * @param eventPreValue On output containing value of first/last node before event.
810 * This value is only written when function returned 1 or greater.
811 * @return int Return number of events that were encountered.
812 */
813 ✗ int findOppositeEndSpatialDistribution(threadData_t *threadData, SPATIAL_DISTRIBUTION_DATA* spatialDistribution, double in0, double in1, double posX, int isPositiveVelocity, double* eventPreValue, double* outValue) {
814 /* Variables */
815 ✗ DOUBLE_ENDED_LIST* transportedQuantityList = spatialDistribution->transportedQuantity;
816 ✗ DOUBLE_ENDED_LIST* storedEventsList = spatialDistribution->storedEvents;
817 DOUBLE_ENDED_LIST_NODE* currentNode;
818 DOUBLE_ENDED_LIST_NODE* firstNode;
819 DOUBLE_ENDED_LIST_NODE* lastNode;
820 DOUBLE_ENDED_LIST_NODE* prevVisitedNode;
821 TRANSPORTED_QUANTITY_DATA* currentNodeData;
822 TRANSPORTED_QUANTITY_DATA* prevVisitedNodeData;
823 TRANSPORTED_QUANTITY_DATA* firstNodeData;
824 TRANSPORTED_QUANTITY_DATA* lastNodeData;
825 TRANSPORTED_QUANTITY_DATA tempData;
826 double edgeNodePosition;
827 double readPosition;
828 double currentDistance;
829 int walkedOverEvents = 0;
830
831 /* Step 0
832 * Check if we are still in spatialDistribution intervall or if deltaX > 1
833 */
834 ✗ firstNode = getFirstNodeDoubleEndedList(transportedQuantityList);
835 ✗ firstNodeData = firstDataDoubleEndedList(transportedQuantityList);
836 ✗ lastNode = getLastNodeDoubleEndedList(transportedQuantityList);
837 ✗ lastNodeData = lastDataDoubleEndedList(transportedQuantityList);
838 ✗ if (isPositiveVelocity) {
839 ✗ if (-posX+1 < firstNodeData->position) {
840 // We need to interpolate (-posX,in0) <-> (-posX+1,out1) <-> (firstNodeData->position, firstNodeData->value)
841 // ^
842 // |
843 ✗ tempData.position = -posX;
844 ✗ tempData.value = in0;
845 ✗ *outValue = interpolateTransportedQuantity(threadData, &tempData, firstNodeData, -posX + 1);
846 ✗ *eventPreValue = *outValue;
847 ✗ return doubleEndedListLen(storedEventsList);
848 }
849 } else {
850 ✗ if (-posX > lastNodeData->position) {
851 // We need to interpolate (lastNodeData->position,lastNodeData->value) <-> (-posX, out0) <-> (-posX+1,in1)
852 // ^
853 // |
854 ✗ tempData.position = -posX+1;
855 ✗ tempData.value = in1;
856 ✗ *outValue = interpolateTransportedQuantity(threadData, lastNodeData, &tempData, -posX);
857 ✗ *eventPreValue = *outValue;
858 ✗ return doubleEndedListLen(storedEventsList);
859 }
860 }
861
862 /* Step 1
863 * Walk from the opposite end to the read position: z(1,t) sits at -posX+1 and
864 * z(0,t) at -posX, for the x of this call. Clamped to the stored profile in
865 * case x moved backwards since the last stored step.
866 */
867 ✗ if (isPositiveVelocity) {
868 ✗ edgeNodePosition = firstNodeData->position;
869 ✗ readPosition = -posX+1;
870 ✗ if (readPosition > lastNodeData->position) {
871 readPosition = lastNodeData->position;
872 }
873 currentNode = lastNode;
874 } else {
875 ✗ edgeNodePosition = lastNodeData->position;
876 ✗ readPosition = -posX;
877 ✗ if (readPosition < firstNodeData->position) {
878 readPosition = firstNodeData->position;
879 }
880 currentNode = firstNode;
881 }
882 ✗ currentNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(currentNode);
883
884 ✗ currentDistance = fabs(currentNodeData->position - edgeNodePosition);
885 ✗ if (currentDistance + spatialPosEps(currentNodeData->position, edgeNodePosition) < 1) {
886 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Error for spatialDistribution in function findOppositeEndSpatialDistribution.\nThis case should not be possible. Please open a bug report about it.");
887 ✗ omc_throw_function(threadData);
888 return walkedOverEvents;
889 }
890
891 /* Move to neighbor */
892 prevVisitedNode = currentNode;
893 ✗ prevVisitedNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(prevVisitedNode);
894
895 ✗ while (currentNode != NULL) {
896 ✗ if (isPositiveVelocity) {
897 ✗ currentNode = getPreviousNodeDoubleEndedList(currentNode);
898 } else {
899 ✗ currentNode = getNextNodeDoubleEndedList(currentNode);
900 }
901 ✗ if(currentNode == NULL) {
902 break;
903 }
904 ✗ currentNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(currentNode);
905
906 /* Check for event:
907 * Current node position equal to previous visited node position
908 */
909 ✗ if (fabs(prevVisitedNodeData->position - currentNodeData->position) < spatialPosEps(prevVisitedNodeData->position, currentNodeData->position)) {
910 ✗ *eventPreValue = prevVisitedNodeData->value;
911 ✗ walkedOverEvents += 1;
912 }
913
914 /* Check if the read position is passed */
915 ✗ currentDistance = fabs(currentNodeData->position - readPosition);
916 ✗ if (currentDistance > spatialPosEps(currentNodeData->position, readPosition) &&
917 (isPositiveVelocity ? currentNodeData->position < readPosition
918 : currentNodeData->position > readPosition)) {
919 break;
920 } else {
921 prevVisitedNode = currentNode;
922 ✗ prevVisitedNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(prevVisitedNode);
923 }
924 }
925
926 /* Step 2
927 * Interpolate at the read position.
928 */
929 ✗ if (currentNode == NULL) {
930 /* Read position is at or beyond the far end of the stored profile */
931 ✗ if (isPositiveVelocity) {
932 ✗ *outValue = firstNodeData->value;
933 } else {
934 ✗ *outValue = lastNodeData->value;
935 }
936 } else {
937 ✗ if (isPositiveVelocity) {
938 ✗ *outValue = interpolateTransportedQuantity(threadData, currentNodeData, prevVisitedNodeData, readPosition);
939 } else {
940 ✗ *outValue = interpolateTransportedQuantity(threadData, prevVisitedNodeData, currentNodeData, readPosition);
941 }
942 }
943
944 return walkedOverEvents;
945 }
946
947
948
949 /**
950 * @brief Remove nodes until distance between first and last element is 1.
951 *
952 * @param transportedQuantityList Double ended list containing spatial distribution.
953 * @param isPositiveVelocity Boolean describing if velocity v is positive (>=0).
954 * Velocity v is `v:=der(x)`.
955 * @param eventPreValue On output containing value of first/last node before event.
956 * This value is only written when function returned 1 or greater.
957 * @return int Return number of events that were encountered.
958 */
959 ✗ int pruneSpatialDistribution(threadData_t *threadData, SPATIAL_DISTRIBUTION_DATA* spatialDistribution, int isPositiveVelocity) {
960 /* Variables */
961 ✗ DOUBLE_ENDED_LIST* transportedQuantityList = spatialDistribution->transportedQuantity;
962 ✗ DOUBLE_ENDED_LIST* storedEventsList = spatialDistribution->storedEvents;
963 DOUBLE_ENDED_LIST_NODE* edgeNode;
964 DOUBLE_ENDED_LIST_NODE* currentNode;
965 DOUBLE_ENDED_LIST_NODE* prevVisitedNode;
966 TRANSPORTED_QUANTITY_DATA* edgeNodeData;
967 TRANSPORTED_QUANTITY_DATA* currentNodeData;
968 TRANSPORTED_QUANTITY_DATA* prevVisitedNodeData;
969 TRANSPORTED_EVENT_DATA* eventData;
970 int walkedOverEvents = 0;
971 int i;
972 double currentDistance;
973
974 /* Step 1
975 * Walk over list, starting from opposite side of edgeNode,
976 * until distance between currentNode and edgeNode < 1.
977 */
978 ✗ if (isPositiveVelocity) {
979 ✗ edgeNode = getFirstNodeDoubleEndedList(transportedQuantityList);
980 ✗ currentNode = getLastNodeDoubleEndedList(transportedQuantityList);
981 } else {
982 ✗ edgeNode = getLastNodeDoubleEndedList(transportedQuantityList);
983 ✗ currentNode = getFirstNodeDoubleEndedList(transportedQuantityList);
984 }
985 ✗ edgeNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(edgeNode);
986 ✗ currentNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(currentNode);
987
988 ✗ currentDistance = fabs(currentNodeData->position - edgeNodeData->position);
989 ✗ if (currentDistance + spatialPosEps(currentNodeData->position, edgeNodeData->position) < 1) {
990 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Error for spatialDistribution in function pruneSpatialDistribution.\nThis case should not be possible. Please open a bug reoprt about it.");
991 ✗ omc_throw_function(threadData);
992 }
993
994 /* Move to neighbor */
995 prevVisitedNode = currentNode;
996 ✗ prevVisitedNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(prevVisitedNode);
997
998 ✗ while (currentNode != edgeNode) {
999 ✗ if (isPositiveVelocity) {
1000 ✗ currentNode = getPreviousNodeDoubleEndedList(currentNode);
1001 } else {
1002 ✗ currentNode = getNextNodeDoubleEndedList(currentNode);
1003 }
1004 ✗ currentNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(currentNode);
1005
1006 /* Check for event:
1007 * Current node position equal to previous visited node position
1008 */
1009 ✗ if (fabs(prevVisitedNodeData->position - currentNodeData->position) < spatialPosEps(prevVisitedNodeData->position, currentNodeData->position)) {
1010 ✗ walkedOverEvents += 1;
1011 }
1012
1013 /* Check if distance between currentNode and edgeNode is < 1 */
1014 ✗ currentDistance = fabs(currentNodeData->position - edgeNodeData->position);
1015 ✗ if (currentDistance + spatialPosEps(currentNodeData->position, edgeNodeData->position) < 1) {
1016 break;
1017 } else {
1018 prevVisitedNode = currentNode;
1019 ✗ prevVisitedNodeData = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(prevVisitedNode);
1020 }
1021 }
1022
1023 /* Step 2
1024 * Interpolate at edgeNode->position +/- 1.
1025 */
1026 ✗ if (currentDistance + spatialPosEps(currentNodeData->position, edgeNodeData->position) < 1) {
1027 ✗ if (isPositiveVelocity) {
1028 ✗ prevVisitedNodeData->value = interpolateTransportedQuantity(threadData, currentNodeData, prevVisitedNodeData, edgeNodeData->position + 1);
1029 ✗ prevVisitedNodeData->position = edgeNodeData->position + 1;
1030 } else {
1031 ✗ prevVisitedNodeData->value = interpolateTransportedQuantity(threadData, prevVisitedNodeData, currentNodeData, edgeNodeData->position - 1);
1032 ✗ prevVisitedNodeData->position = edgeNodeData->position - 1;
1033 }
1034 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "Interpolate at %s", isPositiveVelocity?"end":"front");
1035 }
1036
1037 /* Step 3
1038 * Remove all nodes that have a distance to edge > 1.
1039 */
1040 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "Removing nodes %s node %p", isPositiveVelocity?"after":"before", (void*)prevVisitedNode);
1041 ✗ if (isPositiveVelocity) {
1042 ✗ clearAfterNodeDoubleEndedList(transportedQuantityList, prevVisitedNode);
1043 } else {
1044 ✗ clearBeforeNodeDoubleEndedList(transportedQuantityList, prevVisitedNode);
1045 }
1046
1047 /* Step 4
1048 * Remove all events that are outside spatial distribution [leftEdge-SPATIAL_ZERO_DELTA_X, rightEdge+SPATIAL_ZERO_DELTA_X]
1049 */
1050 ✗ if (doubleEndedListLen(storedEventsList) > 0) {
1051 ✗ if (isPositiveVelocity) {
1052 ✗ eventData = lastDataDoubleEndedList(storedEventsList);
1053 ✗ while (edgeNodeData->position+1 + spatialZeroDeltaX(edgeNodeData->position, eventData->position) < eventData->position) {
1054 ✗ spatialDistribution->lastStoredEventValue = eventData->zeroCrossValue;
1055 ✗ removeLastDoubleEndedList(storedEventsList);
1056 ✗ if (doubleEndedListLen(storedEventsList) == 0) {
1057 break;
1058 } else {
1059 ✗ eventData = lastDataDoubleEndedList(storedEventsList);
1060 }
1061 }
1062 } else {
1063 ✗ eventData = firstDataDoubleEndedList(storedEventsList);
1064 ✗ while (edgeNodeData->position-1 - spatialZeroDeltaX(edgeNodeData->position, eventData->position) > eventData->position) {
1065 ✗ spatialDistribution->lastStoredEventValue = eventData->zeroCrossValue;
1066 ✗ removeFirstDoubleEndedList(storedEventsList);
1067 ✗ if (doubleEndedListLen(storedEventsList) == 0) {
1068 break;
1069 } else {
1070 ✗ eventData = firstDataDoubleEndedList(storedEventsList);
1071 }
1072 }
1073 }
1074 }
1075
1076 /* Debug prints */
1077 ✗ doubleEndedListPrint(transportedQuantityList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
1078 ✗ infoStreamPrint(OMC_LOG_SPATIALDISTR, 0, "List of events");
1079 ✗ doubleEndedListPrint(storedEventsList, OMC_LOG_SPATIALDISTR, &printTransportedQuantity);
1080
1081 ✗ return walkedOverEvents;
1082 }
1083
1084
1085 /**
1086 * @brief Print transported quantity data to stream.
1087 *
1088 * Prints tuple (position, value).
1089 *
1090 * @param data Void pointer to transportedQuantityData.
1091 * Will be casted to TRANSPORTED_QUANTITY_DATA*.
1092 * @param stream Stream of OMC_LOG_STREAM type.
1093 * @param nodePointer Address of node storing this data.
1094 */
1095 ✗ void printTransportedQuantity(void* data, int stream, void* nodePointer) {
1096 TRANSPORTED_QUANTITY_DATA* transportedQuantityData = (TRANSPORTED_QUANTITY_DATA*) data;
1097 ✗ infoStreamPrint(stream, 0, "%p: (%e,%e)", nodePointer, transportedQuantityData->position, transportedQuantityData->value);
1098 ✗ }
1099
1100
1101 //#endif
1102
1103 /**
1104 * @brief The spatialDistribution operators as flat words, for an FMU state.
1105 *
1106 * Per operator its scalar fields, then the transported quantity and the stored
1107 * events, each led by its length.
1108 *
1109 * @param data Runtime data struct.
1110 * @param out Receives the words, or NULL to only count them.
1111 * @return Number of words.
1112 */
1113 ✗ size_t spatialDistributionStateWords(DATA* data, double* out)
1114 {
1115 size_t k = 0;
1116 unsigned int i;
1117 DOUBLE_ENDED_LIST_NODE* node;
1118
1119 ✗ if (!data->simulationInfo->spatialDistributionData) {
1120 return 0;
1121 }
1122 ✗ for (i = 0; i < data->modelData->nSpatialDistributions; i++) {
1123 ✗ SPATIAL_DISTRIBUTION_DATA* s = &data->simulationInfo->spatialDistributionData[i];
1124 ✗ if (out) {
1125 ✗ out[k] = s->isInitialized;
1126 ✗ out[k+1] = s->oldPosX;
1127 ✗ out[k+2] = s->startPosXSet;
1128 ✗ out[k+3] = s->startPosX;
1129 ✗ out[k+4] = s->lastStoredEventValue;
1130 ✗ out[k+5] = (double) s->nWarningsRemovedEvents;
1131 ✗ out[k+6] = (double) s->nWarningsOutputEvents;
1132 ✗ out[k+7] = doubleEndedListLen(s->transportedQuantity);
1133 }
1134 ✗ k += 8;
1135 ✗ for (node = getFirstNodeDoubleEndedList(s->transportedQuantity); node; node = getNextNodeDoubleEndedList(node)) {
1136 ✗ TRANSPORTED_QUANTITY_DATA* q = (TRANSPORTED_QUANTITY_DATA*) dataDoubleEndedList(node);
1137 ✗ if (out) {
1138 ✗ out[k] = q->position;
1139 ✗ out[k+1] = q->value;
1140 }
1141 ✗ k += 2;
1142 }
1143 ✗ if (out) out[k] = doubleEndedListLen(s->storedEvents);
1144 ✗ k++;
1145 ✗ for (node = getFirstNodeDoubleEndedList(s->storedEvents); node; node = getNextNodeDoubleEndedList(node)) {
1146 ✗ TRANSPORTED_EVENT_DATA* e = (TRANSPORTED_EVENT_DATA*) dataDoubleEndedList(node);
1147 ✗ if (out) {
1148 ✗ out[k] = e->position;
1149 ✗ out[k+1] = e->zeroCrossValue;
1150 }
1151 ✗ k += 2;
1152 }
1153 }
1154 return k;
1155 }
1156
1157 /**
1158 * @brief Inverse of spatialDistributionStateWords.
1159 *
1160 * @param data Runtime data struct.
1161 * @param w Words spatialDistributionStateWords wrote.
1162 * @param len Number of words available.
1163 * @return Number of words read, or -1 if they are not what it wrote.
1164 */
1165 ✗ long setSpatialDistributionStateWords(DATA* data, const double* w, size_t len)
1166 {
1167 size_t k = 0;
1168 unsigned int i;
1169 int j, n;
1170
1171 ✗ if (!data->simulationInfo->spatialDistributionData) {
1172 return 0;
1173 }
1174 ✗ for (i = 0; i < data->modelData->nSpatialDistributions; i++) {
1175 ✗ SPATIAL_DISTRIBUTION_DATA* s = &data->simulationInfo->spatialDistributionData[i];
1176 ✗ if (k + 8 > len) return -1;
1177 ✗ s->isInitialized = w[k] != 0;
1178 ✗ s->oldPosX = w[k+1];
1179 ✗ s->startPosXSet = w[k+2] != 0;
1180 ✗ s->startPosX = w[k+3];
1181 ✗ s->lastStoredEventValue = (int) w[k+4];
1182 ✗ s->nWarningsRemovedEvents = (unsigned long) w[k+5];
1183 ✗ s->nWarningsOutputEvents = (unsigned long) w[k+6];
1184 ✗ n = (int) w[k+7];
1185 k += 8;
1186 ✗ if (n < 0 || k + 2 * (size_t) n + 1 > len) return -1;
1187 ✗ clearDoubleEndedList(s->transportedQuantity);
1188 ✗ for (j = 0; j < n; j++) {
1189 TRANSPORTED_QUANTITY_DATA q;
1190 ✗ q.position = w[k];
1191 ✗ q.value = w[k+1];
1192 ✗ k += 2;
1193 ✗ pushBackDoubleEndedList(s->transportedQuantity, &q);
1194 }
1195 ✗ n = (int) w[k++];
1196 ✗ if (n < 0 || k + 2 * (size_t) n > len) return -1;
1197 ✗ clearDoubleEndedList(s->storedEvents);
1198 ✗ for (j = 0; j < n; j++) {
1199 TRANSPORTED_EVENT_DATA e;
1200 ✗ e.position = w[k];
1201 ✗ e.zeroCrossValue = w[k+1];
1202 ✗ k += 2;
1203 ✗ pushBackDoubleEndedList(s->storedEvents, &e);
1204 }
1205 }
1206 ✗ return (long) k;
1207 }
1208