Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 49.0% 253 / 0 / 516
Functions: 52.6% 10 / 0 / 19
Branches: 33.1% 81 / 0 / 245
OMCompiler/SimulationRuntime/c/simulation/solver/dassl.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 #include <float.h>
28 #include <math.h>
29 #include <string.h>
30 #include <setjmp.h>
31 #include <time.h>
32
33 #include "openmodelica.h"
34 #include "openmodelica_func.h"
35 #include "simulation_data.h"
36
37 #include "gc/omc_gc.h"
38 #include "util/context.h"
39 #include "simulation/jacobian_util.h"
40 #include "util/omc_error.h"
41
42 #include "../arrayIndex.h"
43 #include "epsilon.h"
44 #include "external_input.h"
45 #include "model_help.h"
46 #include "omc_math.h"
47 #include "simulation/options.h"
48 #include "simulation/results/simulation_result.h"
49 #include "simulation/simulation_runtime.h"
50 #include "solver_main.h"
51
52 #include "dassl.h"
53
54 #define UNUSED(x) (void)(x) /* Surpress compiler warnings for unused function input */
55 #define DASSL_FIRST_STEP_RESTARTS 3
56
57 #ifdef __cplusplus
58 extern "C" {
59 #endif
60
61 /* experimental flag for SKF TLM Master Solver Interface
62 * - it's used with -noEquidistantTimeGrid flag.
63 * - it's set to 1 if the continuous system is evaluated
64 * when dassl finished a step, otherwise it's 0.
65 */
66 int RHSFinalFlag;
67
68 /* provides a dummy Jacobian to be used with DASSL */
69 ✗ static int dummy_Jacobian(double *t, double *y, double *yprime,double *deltaD,
70 double *delta, double *cj, double *h, double *wt,
71 double *rpar, int* ipar) {
72 ✗ return 0;
73 }
74
75 /* provides a dummy zero crossing function to be used with DASSL */
76 ✗ static int dummy_zeroCrossing(int *neqm, double *t, double *y, double *yp,
77 int *ng, double *gout, double *rpar, int* ipar) {
78 ✗ return 0;
79 }
80
81 /* provides a dumm precondition function to be used with DASSL */
82 ✗ static int dummy_precondition(int *neq, double *t, double *y, double *yprime,
83 double *savr, double *pwk, double *cj,
84 double *wt, double *wp, int *iwp, double *b,
85 double eplin, int* ires, double *rpar, int* ipar){
86 ✗ return 0;
87 }
88
89 /* Function prototypes */
90 static int callJacobian(double *t, double *y, double *yprime, double *deltaD,
91 double *pd, double *cj, double *h, double *wt,
92 double *rpar, int* ipar);
93
94 int jacA_num(double *t, double *y, double *yprime, double *deltaD,
95 double *pd, double *cj, double *h, double *wt,
96 double *rpar, int* ipar);
97
98 int jacA_numColored(double *t, double *y, double *yprime,
99 double *deltaD, double *pd, double *cj, double *h,
100 double *wt, double *rpar, int* ipar);
101
102 int jacA_sym(double *t, double *y, double *yprime, double *deltaD,
103 double *pd, double *cj, double *h, double *wt,
104 double *rpar, int* ipar);
105
106 int jacA_symColored(double *t, double *y, double *yprime,
107 double *deltaD, double *pd, double *cj, double *h,
108 double *wt, double *rpar, int* ipar);
109
110 void DDASKR(
111 int (*res) (double *t, double *y, double *yprime, double* cj, double *delta, int *ires, double *rpar, int* ipar),
112 int *neq,
113 double *t,
114 double *y,
115 double *yprime,
116 double *tout,
117 int *info,
118 double *rtol,
119 double *atol,
120 int *idid,
121 double *rwork,
122 int *lrw,
123 int *iwork,
124 int *liw,
125 double *rpar,
126 int *ipar,
127 int (*jac) (double *t, double *y, double *yprime, double *deltaD, double *delta, double *cj, double *h, double *wt, double *rpar, int* ipar),
128 int (*psol) (int *neq, double *t, double *y, double *yprime, double *savr, double *pwk, double *cj, double *wt, double *wp, int *iwp, double *b, double eplin, int* ires, double *rpar, int* ipar),
129 int (*g) (int *neqm, double *t, double *y, double *yp, int *ng, double *gout, double *rpar, int* ipar),
130 int *ng,
131 int *jroot
132 );
133
134 static int continue_DASSL(int* idid, double* tolarence);
135 static int dasslStuck(DASSL_DATA* dasslData, double t);
136
137 /* function for calculating state values on residual form */
138 static int functionODE_residual(double *t, double *y, double *yd, double* cj,
139 double *delta, int *ires, double *rpar, int *ipar);
140
141 /* function for calculating zeroCrossings */
142 static int function_ZeroCrossingsDASSL(int *neqm, double *t, double *y,
143 double *yp, int *ng, double *gout,
144 double *rpar, int* ipar);
145
146
147 /*
148 * \brief Read the states' nominal values into the absolute tolerances.
149 *
150 * Re-read by updateSolverNominals once initialization has computed the nominals
151 * that are parameter expressions.
152 */
153 2 void dassl_setNominals(DATA* data, DASSL_DATA *dasslData)
154 {
155 int i;
156 char name[2048];
157 const array_index_t *ix;
158 const STATIC_REAL_DATA *var;
159
160 2 infoStreamPrint(OMC_LOG_SOLVER, 1, "The relative tolerance is %g. Following absolute tolerances are used for the states: ", data->simulationInfo->tolerance);
161
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 for(i=0; i<dasslData->N; ++i)
162 {
163 4 const modelica_real nominal = getNominalFromScalarIdx(data->simulationInfo, data->modelData, VAR_KIND_STATE, i);
164 4 dasslData->nominal[i] = fmax(fabs(nominal), 1e-32);
165 4 dasslData->rtol[i] = data->simulationInfo->tolerance;
166 4 dasslData->atol[i] = data->simulationInfo->tolerance * dasslData->nominal[i];
167
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER_V)) {
168 4 ix = &data->simulationInfo->realVarsReverseIndex[i];
169 4 var = &data->modelData->realVarsData[ix->array_idx];
170 4 printArrayElementName(name, sizeof(name), var->info.name, &var->dimension, ix->dim_idx, FALSE);
171 4 infoStreamPrint(OMC_LOG_SOLVER_V, 0, "%d. %s -> %g", i+1, name, dasslData->atol[i]);
172 }
173 }
174 2 messageClose(OMC_LOG_SOLVER);
175 2 }
176
177 /*
178 * \brief Configure DASSL solver
179 *
180 * Allocate memory for intern data of `dasslData`.
181 * Configures DASSL:
182 * - Set relative and absolute tolerance
183 * - Set maximum step size, initial step size
184 * - Set maximum integration order
185 * - Set time grid
186 * - Set method for jacobian computation
187 * - Set root finding method
188 * - Set event handling and restart option
189 *
190 */
191 1 int dassl_initial(DATA* data, threadData_t *threadData,
192 SOLVER_INFO* solverInfo, DASSL_DATA *dasslData)
193 {
194 /* work arrays for DASSL */
195 unsigned int i;
196 long N;
197 SIMULATION_DATA tmpSimData = {0};
198
199 1 dasslData->residualFunction = functionODE_residual;
200 1 N = data->modelData->nStates;
201
202 1 dasslData->N = N;
203
204 1 RHSFinalFlag = 0;
205
206 1 dasslData->liw = 40 + N;
207 1 dasslData->lrw = 60 + ((maxOrder + 4) * N) + (N * N) + (3*data->modelData->nZeroCrossings);
208 1 dasslData->rwork = (double*) calloc(dasslData->lrw, sizeof(double));
209
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 assertStreamPrint(threadData, 0 != dasslData->rwork,"out of memory");
210 1 dasslData->iwork = (int*) calloc(dasslData->liw, sizeof(int));
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 assertStreamPrint(threadData, 0 != dasslData->iwork,"out of memory");
212 1 dasslData->ng = (int) data->modelData->nZeroCrossings;
213 1 dasslData->jroot = (int*) calloc(data->modelData->nZeroCrossings, sizeof(int));
214 1 dasslData->rpar = (double**) malloc(3*sizeof(double*));
215 1 dasslData->ipar = (int*) malloc(sizeof(int));
216 1 dasslData->ipar[0] = OMC_ACTIVE_STREAM(OMC_LOG_JAC);
217 assertStreamPrint(threadData, 0 != dasslData->ipar,"out of memory");
218 1 dasslData->atol = (double*) malloc(N*sizeof(double));
219 1 dasslData->rtol = (double*) malloc(N*sizeof(double));
220 1 dasslData->nominal = (double*) malloc(N*sizeof(double));
221 1 dasslData->info = (int*) calloc(infoLength, sizeof(int));
222
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 assertStreamPrint(threadData, 0 != dasslData->info,"out of memory");
223
224 1 dasslData->idid = 0;
225 1 dasslData->tinySteps = 0;
226
227 1 dasslData->ysave = (double*) malloc(N*sizeof(double));
228 1 dasslData->delta_hh = (double*) malloc(N*sizeof(double));
229 1 dasslData->newdelta = (double*) malloc(N*sizeof(double));
230 1 dasslData->stateDer = (double*) calloc(N, sizeof(double));
231 1 dasslData->states = (double*) malloc(N*sizeof(double));
232
233 1 data->simulationInfo->currentContext = CONTEXT_ALGEBRAIC;
234
235 /* ### start configuration of dassl ### */
236 1 infoStreamPrint(OMC_LOG_SOLVER, 1, "Configuration of the dassl code:");
237
238
239
240 2 dasslData->jacNominalFactor = omc_flag[FLAG_JACOBIAN_NOMINAL_FACTOR]
241
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 ? atof(omc_flagValue[FLAG_JACOBIAN_NOMINAL_FACTOR]) : 1.0;
242
243 /* set nominal values of the states for absolute tolerances */
244 1 dasslData->info[1] = 1;
245 1 dassl_setNominals(data, dasslData);
246
247
248 /* let dassl return at every internal step */
249 1 dasslData->info[2] = 1;
250
251
252 /* define maximum step size dassl is allowed to go */
253
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (omc_flag[FLAG_MAX_STEP_SIZE])
254 {
255 ✗ double maxStepSize = atof(omc_flagValue[FLAG_MAX_STEP_SIZE]);
256
257 ✗ assertStreamPrint(threadData, maxStepSize >= DASSL_STEP_EPS, "Selected maximum step size %e is too small.", maxStepSize);
258
259 ✗ dasslData->rwork[1] = maxStepSize;
260 ✗ dasslData->info[6] = 1;
261 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "maximum step size %g", dasslData->rwork[1]);
262 }
263 else
264 {
265 1 infoStreamPrint(OMC_LOG_SOLVER, 0, "maximum step size not set");
266 }
267
268
269 /* define initial step size, which is dassl is used every time it restarts */
270
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (omc_flag[FLAG_INITIAL_STEP_SIZE])
271 {
272 ✗ double initialStepSize = atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]);
273
274 ✗ assertStreamPrint(threadData, initialStepSize >= DASSL_STEP_EPS, "Selected initial step size %e is too small.", initialStepSize);
275
276 ✗ dasslData->rwork[2] = initialStepSize;
277 ✗ dasslData->info[7] = 1;
278 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "initial step size %g", dasslData->rwork[2]);
279 }
280 else
281 {
282 1 infoStreamPrint(OMC_LOG_SOLVER, 0, "initial step size not set");
283 }
284
285
286 /* define maximum integration order of dassl */
287
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (omc_flag[FLAG_MAX_ORDER])
288 {
289 ✗ int maxOrder = atoi(omc_flagValue[FLAG_MAX_ORDER]);
290
291 ✗ assertStreamPrint(threadData, maxOrder >= 1 && maxOrder <= 5, "Selected maximum order %d is out of range (1-5).", maxOrder);
292
293 ✗ dasslData->iwork[2] = maxOrder;
294 ✗ dasslData->info[8] = 1;
295 }
296
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 infoStreamPrint(OMC_LOG_SOLVER, 0, "maximum integration order %d", dasslData->info[8]?dasslData->iwork[2]:maxOrder);
297
298
299 /* if FLAG_NOEQUIDISTANT_GRID is set, choose dassl step method */
300
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (omc_flag[FLAG_NOEQUIDISTANT_GRID])
301 {
302 ✗ dasslData->dasslSteps = 1; /* TRUE */
303 ✗ solverInfo->solverNoEquidistantGrid = TRUE;
304 }
305 else
306 {
307 1 dasslData->dasslSteps = 0; /* FALSE */
308 }
309
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 2 infoStreamPrint(OMC_LOG_SOLVER, 0, "use equidistant time grid %s", dasslData->dasslSteps?"NO":"YES");
310
311 /* check if Flags FLAG_NOEQUIDISTANT_OUT_FREQ or FLAG_NOEQUIDISTANT_OUT_TIME are set */
312
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (dasslData->dasslSteps){
313 ✗ if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ])
314 {
315 ✗ dasslData->dasslStepsFreq = atoi(omc_flagValue[FLAG_NOEQUIDISTANT_OUT_FREQ]);
316 }
317 ✗ else if (omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME])
318 {
319 ✗ dasslData->dasslStepsTime = atof(omc_flagValue[FLAG_NOEQUIDISTANT_OUT_TIME]);
320 ✗ dasslData->rwork[1] = dasslData->dasslStepsTime;
321 ✗ dasslData->info[6] = 1;
322 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "maximum step size %g", dasslData->rwork[1]);
323 } else {
324 ✗ dasslData->dasslStepsFreq = 1;
325 ✗ dasslData->dasslStepsTime = 0.0;
326 }
327
328 ✗ if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ] && omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]){
329 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The flags are \"noEquidistantOutputFrequency\" "
330 "and \"noEquidistantOutputTime\" are in opposition "
331 "to each other. The flag \"noEquidistantOutputFrequency\" superiors.");
332 }
333 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "as the output frequency control is used: %d", dasslData->dasslStepsFreq);
334 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "as the output frequency time step control is used: %f", dasslData->dasslStepsTime);
335 }
336
337 /* Choose and initialize the ODE Jacobian. The mapping from the `-jacobian` flag to
338 * the forward / adjoint / bidirectional Jacobian is shared with IDA and GBODE. */
339 1 dasslData->dasslJacobian = getRequestedJacobianMethod(threadData);
340 1 JACOBIAN* jacobian = initSymbolicOdeJacobian(data, threadData, &dasslData->dasslJacobian, FALSE);
341
342 /* default use a user sub-routine for JAC */
343 1 dasslData->info[4] = 1;
344
345 /* set up the appropriate function pointer */
346
1/7
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 1 switch (dasslData->dasslJacobian){
347 1 case COLOREDNUMJAC:
348 1 data->simulationInfo->jacobianEvals = jacobian->sparsePattern->maxColors;
349 1 dasslData->jacobianFunction = jacA_numColored;
350 1 break;
351 ✗ case BICOLOREDSYMJAC:
352 ✗ data->simulationInfo->jacobianEvals = jacobian->sparsePattern->maxColors
353 ✗ + jacobian->adjointJacobian->sparsePattern->maxColors;
354 ✗ dasslData->jacobianFunction = jacA_symColored;
355 ✗ break;
356 ✗ case COLOREDSYMJAC:
357 case COLOREDSYMJACADJ:
358 ✗ data->simulationInfo->jacobianEvals = jacobian->sparsePattern->maxColors;
359 ✗ dasslData->jacobianFunction = jacA_symColored;
360 ✗ break;
361 ✗ case SYMJAC:
362 ✗ dasslData->jacobianFunction = jacA_sym;
363 ✗ break;
364 ✗ case NUMJAC:
365 ✗ dasslData->jacobianFunction = jacA_num;
366 ✗ break;
367 ✗ case INTERNALNUMJAC:
368 ✗ dasslData->jacobianFunction = dummy_Jacobian;
369 /* no user sub-routine for JAC */
370 ✗ dasslData->info[4] = 0;
371 ✗ break;
372 ✗ default:
373 ✗ throwStreamPrint(threadData,"unrecognized jacobian calculation method %s", (const char*)omc_flagValue[FLAG_JACOBIAN]);
374 break;
375 }
376 1 infoStreamPrint(OMC_LOG_SOLVER, 0, "jacobian is calculated by %s", JACOBIAN_METHOD_DESC[dasslData->dasslJacobian]);
377
378 /* if FLAG_NO_ROOTFINDING is set, choose dassl with out internal root finding */
379
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(omc_flag[FLAG_NO_ROOTFINDING])
380 {
381 ✗ dasslData->dasslRootFinding = 0;
382 ✗ dasslData->zeroCrossingFunction = dummy_zeroCrossing;
383 ✗ dasslData->ng = 0;
384 }
385 else
386 {
387 1 solverInfo->solverRootFinding = 1;
388 1 dasslData->dasslRootFinding = 1;
389 1 dasslData->zeroCrossingFunction = function_ZeroCrossingsDASSL;
390 }
391
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 infoStreamPrint(OMC_LOG_SOLVER, 0, "dassl uses internal root finding method %s", dasslData->dasslRootFinding?"YES":"NO");
392
393
394 /* if FLAG_NO_RESTART is set, choose dassl step method */
395
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (omc_flag[FLAG_NO_RESTART])
396 {
397 ✗ dasslData->dasslAvoidEventRestart = 1; /* TRUE */
398 }
399 else
400 {
401 1 dasslData->dasslAvoidEventRestart = 0; /* FALSE */
402 }
403
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 2 infoStreamPrint(OMC_LOG_SOLVER, 0, "dassl performs an restart after an event occurs %s", dasslData->dasslAvoidEventRestart?"NO":"YES");
404
405 /* ### end configuration of dassl ### */
406
407 1 messageClose(OMC_LOG_SOLVER);
408 1 return 0;
409 }
410
411
412 /*
413 * \brief Deallocates `DASSL_DATA`
414 */
415 1 int dassl_deinitial(DATA* data, DASSL_DATA *dasslData)
416 {
417 unsigned int i;
418
419 /* free work arrays for DASSL */
420 1 free(dasslData->rwork);
421 1 free(dasslData->iwork);
422 1 free(dasslData->jroot);
423 1 free(dasslData->rpar);
424 1 free(dasslData->ipar);
425 1 free(dasslData->atol);
426 1 free(dasslData->rtol);
427 1 free(dasslData->nominal);
428 1 free(dasslData->info);
429 1 free(dasslData->ysave);
430 1 free(dasslData->delta_hh);
431 1 free(dasslData->newdelta);
432 1 free(dasslData->stateDer);
433 1 free(dasslData->states);
434
435 /* Free Jacobians */
436 1 freeSymbolicOdeJacobian(data);
437
438 1 free(dasslData);
439
440 1 return 0;
441 }
442
443 /* \fn printCurrentStatesVector(int logLevel, double* y, DATA* data, double time)
444 *
445 * \param [in] [logLevel]
446 * \param [in] [states]
447 * \param [in] [data]
448 * \param [in] [time]
449 *
450 * This function outputs states vector.
451 *
452 */
453 39 int printCurrentStatesVector(int logLevel, double* states, DATA* data, double time)
454 {
455 int i;
456 39 infoStreamPrint(logLevel, 1, "states at time=%g", time);
457
2/2
✓ Branch 0 taken 78 times.
✓ Branch 1 taken 39 times.
 117 for(i=0;i<data->modelData->nStates;++i)
458 {
459 78 infoStreamPrint(logLevel, 0, "%d. %s = %g", i+1, data->modelData->realVarsData[i].info.name, states[i]);
460 }
461 39 messageClose(logLevel);
462
463 39 return 0;
464 }
465
466 /* \fn printVector(int logLevel, double* y, DATA* data, double time)
467 *
468 * \param [in] [logLevel]
469 * \param [in] [name]
470 * \param [in] [vec]
471 * \param [in] [size]
472 * \param [in] [time]
473 *
474 * This function outputs a vector of size
475 *
476 */
477 78 int printVector(int logLevel, const char* name, double* vec, int n, double time)
478 {
479 int i;
480 78 infoStreamPrint(logLevel, 1, "%s at time=%g", name, time);
481
2/2
✓ Branch 0 taken 156 times.
✓ Branch 1 taken 78 times.
 234 for(i=0; i<n; ++i)
482 {
483 156 infoStreamPrint(logLevel, 0, "%d. %g", i+1, vec[i]);
484 }
485 78 messageClose(logLevel);
486
487 78 return 0;
488 }
489
490 ✗ int printJacobianMatrix(int logLevel, const char* name, double* matrix, DATA* data, int n, double time)
491 {
492 int row, col;
493
494 ✗ infoStreamPrint(logLevel, 1, "%s at time=%g", name, time);
495 ✗ for (col = 0; col < n; ++col)
496 {
497 ✗ const char* colName = data->modelData->realVarsData[col].info.name;
498 ✗ for (row = 0; row < n; ++row)
499 {
500 ✗ const char* rowName = data->modelData->realVarsData[row].info.name;
501 ✗ const int idx = col * n + row;
502 ✗ infoStreamPrint(logLevel, 0,
503 "J(row=%d:'%s', col=%d:'%s') = %.16g [flat=%d]",
504 ✗ row, rowName, col, colName, matrix[idx], idx);
505 }
506 }
507 ✗ messageClose(logLevel);
508
509 ✗ return 0;
510 }
511
512
513 /**********************************************************************************************
514 * DASSL with synchronous treating of when equation
515 * - without integrated ZeroCrossing method.
516 * + ZeroCrossing are handled outside DASSL.
517 * + if no event occurs outside DASSL performs a warm-start
518 **********************************************************************************************/
519 1 int dassl_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
520 {
521 1 double tout = 0;
522 double tStepStart = 0;
523 int i = 0;
524 unsigned int ui = 0;
525 1 int retVal = 0;
526 int saveJumpState;
527 static unsigned int dasslStepsOutputCounter = 1;
528 1 int return_from_small_step = 0;
529 1 int firstStepRestarts = 0;
530
531 1 DASSL_DATA *dasslData = (DASSL_DATA*) solverInfo->solverData;
532
533 1 SIMULATION_DATA *sData = data->localData[0];
534 1 SIMULATION_DATA *sDataOld = data->localData[1];
535
536 1 modelica_real* states = sData->realVars;
537 1 modelica_real* stateDer = dasslData->stateDer;
538
539
540 MODEL_DATA *mData = (MODEL_DATA*) data->modelData;
541
542
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
543
544 1 memcpy(stateDer, data->localData[1]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
545
546 1 dasslData->rpar[0] = (double*) (void*) data;
547 1 dasslData->rpar[1] = (double*) (void*) dasslData;
548 1 dasslData->rpar[2] = (double*) (void*) threadData;
549
550 1 saveJumpState = threadData->currentErrorStage;
551 1 threadData->currentErrorStage = ERROR_INTEGRATOR;
552
553 /* try */
554 #if !defined(OMC_EMCC)
555
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 OMC_TRY_INTERNAL(simulationJumpBuffer)
556 #endif
557
558
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 assertStreamPrint(threadData, 0 != dasslData->rpar, "could not passed to DDASKR");
559
560 /* If an event is triggered and processed restart dassl. */
561
3/6
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
 1 if(!dasslData->dasslAvoidEventRestart && (solverInfo->didEventStep || 0 == dasslData->idid))
562 {
563 /* obtain reset */
564 1 dasslData->info[0] = 0;
565 1 dasslData->idid = 0;
566
567 }
568
569 /* Calculate steps until TOUT is reached */
570
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (dasslData->dasslSteps)
571 {
572 /* If dasslsteps is selected, the dassl run to stopTime or next sample event */
573 ✗ if (data->simulationInfo->nextSampleEvent < data->simulationInfo->stopTime)
574 {
575 ✗ tout = data->simulationInfo->nextSampleEvent;
576 }
577 else
578 {
579 ✗ tout = data->simulationInfo->stopTime;
580 }
581 }
582 else
583 {
584 1 tout = solverInfo->currentTime + solverInfo->currentStepSize;
585 }
586
587 /* Never step past the next time event: what the model computes beyond it
588 * continues the left limit and is no part of the solution. */
589
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (data->simulationInfo->nextSampleEvent < DBL_MAX &&
590 ✗ (0 == dasslData->info[0] || data->simulationInfo->nextSampleEvent >= dasslData->rwork[3]))
591 {
592 ✗ dasslData->info[3] = 1;
593 ✗ dasslData->rwork[0] = fmax(data->simulationInfo->nextSampleEvent, tout);
594 }
595 else
596 {
597 1 dasslData->info[3] = 0;
598 }
599
600 /* Check that tout is not less than timeValue
601 * else will dassl get in trouble. If that is the case we skip the current step.
602 Also check if step size is smaller than DASSL_STEP_EPS or DASSL_STEP_EPS times simulation interval */
603
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if ((solverInfo->currentStepSize < DASSL_STEP_EPS) ||
604
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 (solverInfo->currentStepSize < DASSL_STEP_EPS*(data->simulationInfo->stopTime - data->simulationInfo->startTime)) )
605 {
606 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "Desired step size %e too small.", solverInfo->currentStepSize);
607 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "Interpolate linear");
608
609 /*euler step*/
610 ✗ for(i = 0; i < data->modelData->nStates; i++)
611 {
612 ✗ sData->realVars[i] = sDataOld->realVars[i] + stateDer[i] * solverInfo->currentStepSize;
613 }
614 ✗ sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
615 ✗ data->callback->functionODE(data, threadData);
616 ✗ solverInfo->currentTime = sData->timeValue;
617
618 ✗ return_from_small_step = 1;
619 }
620 else
621 {
622 do
623 {
624 21 infoStreamPrint(OMC_LOG_DASSL, 1, "new step at time = %.15g", solverInfo->currentTime);
625
626 /* rhs final flag is FALSE during for dassl evaluation */
627 21 RHSFinalFlag = 0;
628
629
1/2
✓ Branch 0 taken 21 times.
✗ Branch 1 not taken.
 21 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
630 /* read input vars */
631 21 externalInputUpdate(data);
632 21 data->callback->input_function(data, threadData);
633
1/2
✓ Branch 0 taken 21 times.
✗ Branch 1 not taken.
 21 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
634
635 21 tStepStart = solverInfo->currentTime;
636 21 DDASKR(dasslData->residualFunction, (int*) &dasslData->N,
637 &solverInfo->currentTime, states, stateDer, &tout,
638 dasslData->info, dasslData->rtol, dasslData->atol, &dasslData->idid,
639 dasslData->rwork, &dasslData->lrw, dasslData->iwork, &dasslData->liw,
640 21 (double*) (void*) dasslData->rpar, dasslData->ipar, callJacobian, dummy_precondition,
641 dasslData->zeroCrossingFunction, (int*) &dasslData->ng, dasslData->jroot);
642 21 dasslData->info[7] = omc_flag[FLAG_INITIAL_STEP_SIZE] ? 1 : 0;
643
644 /* A step landing on TSTOP returns there even past TOUT; called again from
645 * the old T, DDASKR interpolates back to TOUT. */
646
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 21 if (dasslData->idid == 2 && solverInfo->currentTime > tout)
647 {
648 ✗ solverInfo->currentTime = tStepStart;
649 ✗ dasslData->idid = 1;
650 ✗ messageClose(OMC_LOG_DASSL);
651 ✗ continue;
652 }
653
654 /* closing new step message */
655 21 messageClose(OMC_LOG_DASSL);
656
657 /* set ringbuffer time to current time */
658 21 sData->timeValue = solverInfo->currentTime;
659
660 /* rhs final flag is TRUE during for output evaluation */
661 21 RHSFinalFlag = 1;
662
663
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
 21 if(dasslData->idid == -1)
664 {
665 ✗ fflush(stderr);
666 ✗ fflush(stdout);
667 ✗ warningStreamPrint(OMC_LOG_DASSL, 0, "A large amount of work has been expended.(About 500 steps). Trying to continue ...");
668 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "DASSL will try again...");
669 ✗ dasslData->info[0] = 1; /* try again */
670 ✗ if (solverInfo->currentTime <= data->simulationInfo->stopTime)
671 ✗ continue;
672 }
673
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 21 else if(dasslData->idid == -7 && dasslData->iwork[10] == 0 && firstStepRestarts < DASSL_FIRST_STEP_RESTARTS)
674 {
675 /* DASKR gives up after ten corrector failures, each quartering H, so a
676 * first step needing a smaller H is never taken: restart from the H it
677 * reached (RWORK(3)). */
678 ✗ firstStepRestarts++;
679 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "The corrector could not converge on the first step. Restarting with initial step size %g.", dasslData->rwork[2]);
680 ✗ dasslData->info[0] = 0;
681 ✗ dasslData->info[7] = 1;
682 ✗ dasslData->idid = 1;
683 ✗ continue;
684 }
685
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
 21 else if(dasslData->idid < 0)
686 {
687 ✗ fflush(stderr);
688 ✗ fflush(stdout);
689 ✗ retVal = continue_DASSL(&dasslData->idid, &data->simulationInfo->tolerance);
690 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "can't continue. time = %f", sData->timeValue);
691 break;
692 }
693
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 else if(dasslStuck(dasslData, solverInfo->currentTime))
694 {
695 retVal = -1;
696 break;
697 }
698
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
 21 if(dasslData->idid == 5)
699 {
700 ✗ threadData->currentErrorStage = ERROR_EVENTSEARCH;
701 }
702
703 /* emit step, if dasslsteps is selected */
704
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
 21 if (dasslData->dasslSteps)
705 {
706 ✗ if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ]){
707 /* output every n-th time step */
708 ✗ if (dasslStepsOutputCounter >= dasslData->dasslStepsFreq){
709 ✗ dasslStepsOutputCounter = 1; /* next line set it to one */
710 ✗ break;
711 }
712 ✗ dasslStepsOutputCounter++;
713 ✗ } else if (omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]){
714 /* output when time>=k*timeValue */
715 ✗ if (solverInfo->currentTime > dasslStepsOutputCounter * dasslData->dasslStepsTime){
716 ✗ dasslStepsOutputCounter++;
717 ✗ break;
718 }
719 } else {
720 break;
721 }
722 }
723
724
3/4
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 21 } while(dasslData->idid == 1 && !OMC_ERROR_RAISED());
725
726 1 states = dasslData->states;
727 }
728
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
729
730 #if !defined(OMC_EMCC)
731 1 OMC_CATCH_INTERNAL(simulationJumpBuffer)
732 #endif
733 1 threadData->currentErrorStage = saveJumpState;
734
735
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (return_from_small_step) {
736 /* need to do this outside the try-catch macro */
737 return 0;
738 }
739
740 /* if a state event occurs than no sample event does need to be activated */
741
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1 if (data->simulationInfo->sampleActivated && solverInfo->currentTime < data->simulationInfo->nextSampleEvent)
742 {
743 ✗ data->simulationInfo->sampleActivated = 0;
744 }
745
746
747 /* save dassl stats */
748 // TODO: Who thought this is an acceptable way to log stats in iwork? Never heard of structs?
749 1 solverInfo->solverStatsTmp.nStepsTaken = dasslData->iwork[10];
750 1 solverInfo->solverStatsTmp.nCallsODE = dasslData->iwork[11];
751 1 solverInfo->solverStatsTmp.nCallsJacobian = dasslData->iwork[12];
752 1 solverInfo->solverStatsTmp.nErrorTestFailures = dasslData->iwork[13];
753 1 solverInfo->solverStatsTmp.nConvergenceTestFailures = dasslData->iwork[14];
754
755
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if(OMC_ACTIVE_STREAM(OMC_LOG_DASSL))
756 {
757 ✗ infoStreamPrint(OMC_LOG_DASSL, 1, "dassl call statistics: ");
758 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "value of idid: %d", (int)dasslData->idid);
759 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "current time value: %0.4g", solverInfo->currentTime);
760 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "current integration time value: %0.4g", dasslData->rwork[3]);
761 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "step size H to be attempted on next step: %0.4g", dasslData->rwork[2]);
762 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "step size used on last successful step: %0.4g", dasslData->rwork[6]);
763 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "the order of the method used on the last step: %d", dasslData->iwork[7]);
764 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "the order of the method to be attempted on the next step: %d", dasslData->iwork[8]);
765 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "number of steps taken so far: %d", solverInfo->solverStatsTmp.nStepsTaken);
766 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "number of calls of functionODE() : %d", solverInfo->solverStatsTmp.nCallsODE);
767 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "number of calculation of jacobian : %d", solverInfo->solverStatsTmp.nCallsJacobian);
768 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "total number of convergence test failures: %d", solverInfo->solverStatsTmp.nConvergenceTestFailures);
769 ✗ infoStreamPrint(OMC_LOG_DASSL, 0, "total number of error test failures: %d", solverInfo->solverStatsTmp.nErrorTestFailures);
770 ✗ messageClose(OMC_LOG_DASSL);
771 }
772
773 1 infoStreamPrint(OMC_LOG_DASSL, 0, "Finished DASSL step.");
774
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
775
776 1 return retVal;
777 }
778
779 #define DASSL_STUCK_STEPS 1000
780
781 /* A run of accepted steps that are each only a few hundred ulp of time long.
782 * DASKR accepts them, so without this the simulation never ends. */
783 21 static int dasslStuck(DASSL_DATA* dasslData, double t)
784 {
785 21 const double tiny = 1000 * DBL_EPSILON * fmax(fabs(t), 1.0);
786 const char *suppressed;
787
788
1/2
✓ Branch 0 taken 21 times.
✗ Branch 1 not taken.
 21 if (dasslData->rwork[6] >= tiny) {
789 21 dasslData->tinySteps = 0;
790 21 return 0;
791 }
792 ✗ if (0 == dasslData->tinySteps++) {
793 ✗ omc_clear_last_suppressed_error();
794 }
795 ✗ if (dasslData->tinySteps < DASSL_STUCK_STEPS) {
796 return 0;
797 }
798 ✗ errorStreamPrint(OMC_LOG_STDOUT, 1, "The integrator is stuck at time %.15g: its last %d steps were each shorter than %g, too short to move time forward. The model is probably singular or discontinuous here.", t, DASSL_STUCK_STEPS, tiny);
799 ✗ suppressed = omc_last_suppressed_error();
800 ✗ if (suppressed[0]) {
801 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "The last error a nonlinear solver recovered from: %s", suppressed);
802 }
803 ✗ messageClose(OMC_LOG_STDOUT);
804 ✗ return 1;
805 }
806
807 ✗ static int continue_DASSL(int* idid, double* atol)
808 {
809 int retValue = -1;
810
811 ✗ switch(*idid)
812 {
813 case 1:
814 case 2:
815 case 3:
816 /* 1-4 means success */
817 break;
818 ✗ case -1:
819 ✗ warningStreamPrint(OMC_LOG_DASSL, 0, "A large amount of work has been expended.(About 500 steps). Trying to continue ...");
820 retValue = 1; /* adrpo: try to continue */
821 ✗ break;
822 ✗ case -2:
823 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The error tolerances are too stringent");
824 retValue = -2;
825 ✗ break;
826 ✗ case -3:
827 /* wbraun: don't throw at this point let the solver handle it */
828 /* throwStreamPrint("DDASKR: THE LAST STEP TERMINATED WITH A NEGATIVE IDID value"); */
829 retValue = -3;
830 ✗ break;
831 ✗ case -6:
832 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "DDASSL had repeated error test failures on the last attempted step.");
833 retValue = -6;
834 ✗ break;
835 ✗ case -7:
836 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The corrector could not converge.");
837 retValue = -7;
838 ✗ break;
839 ✗ case -8:
840 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The matrix of partial derivatives is singular.");
841 retValue = -8;
842 ✗ break;
843 ✗ case -9:
844 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The corrector could not converge. There were repeated error test failures in this step.");
845 retValue = -9;
846 ✗ break;
847 ✗ case -10:
848 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "A Modelica assert prevents the integrator to continue. For more information use -lv LOG_SOLVER");
849 retValue = -10;
850 ✗ break;
851 ✗ case -11:
852 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "IRES equal to -2 was encountered and control is being returned to the calling program.");
853 retValue = -11;
854 ✗ break;
855 ✗ case -12:
856 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "DDASSL failed to compute the initial YPRIME.");
857 retValue = -12;
858 ✗ break;
859 ✗ case -33:
860 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The code has encountered trouble from which it cannot recover.");
861 retValue = -33;
862 ✗ break;
863 }
864
865 ✗ return retValue;
866 }
867
868 /**
869 * @brief ODE residual function.
870 *
871 * Compute difference between old and new state derivatives.
872 *
873 * @param t Independent variable (time).
874 * @param y Array with state variables, size dasslData->N.
875 * @param yd Array with state derivatives, size dasslData->N.
876 * @param cj Unused, specified by DDASKR interface but can be ignored.
877 * @param delta Output: state derivatives - yd
878 * @param ires If not successfull set to -1 on exit.
879 * @param rpar Struct storing user data.
880 * Type {DATA*, DASSL_DATA*, threadData_t*}
881 * TODO: Why is this of type double* and not void*? I guess DDASKR needs it in this specific format...
882 * @param ipar Unused, specified by DDASKR interface.
883 * @return int Return 0.
884 */
885 39 static int functionODE_residual(double *t, double *y, double *yd, double* cj,
886 double *delta, int *ires, double *rpar, int *ipar)
887 {
888 UNUSED(cj); UNUSED(ipar); /* Silence compíler warnings */
889
890 39 DATA* data = (DATA*)((double**)rpar)[0];
891 DASSL_DATA* dasslData = (DASSL_DATA*)((double**)rpar)[1];
892 39 threadData_t *threadData = (threadData_t*)((double**)rpar)[2];
893
894 double timeBackup;
895 long i;
896 int saveJumpState;
897 39 int success = 0;
898
899
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
900
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (measure_time_flag) rt_tick(SIM_TIMER_RESIDUALS);
901
902
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 14 times.
 39 if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
903 {
904 25 setContext(data, *t, CONTEXT_ODE);
905 }
906 39 printCurrentStatesVector(OMC_LOG_DASSL_STATES, y, data, *t);
907 39 printVector(OMC_LOG_DASSL_STATES, "yd", yd, data->modelData->nStates, *t);
908
909 39 timeBackup = data->localData[0]->timeValue;
910 39 data->localData[0]->timeValue = *t;
911
912 39 saveJumpState = threadData->currentErrorStage;
913 39 threadData->currentErrorStage = ERROR_INTEGRATOR;
914
915 /* try */
916 #if !defined(OMC_EMCC)
917
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 OMC_TRY_INTERNAL(simulationJumpBuffer)
918 #endif
919
920 /* read input vars */
921 39 externalInputUpdate(data);
922 /* eval input vars */
923 39 data->callback->input_function(data, threadData);
924
925 /* Compute state derivatives */
926 // TODO: Why is y not used to update states in data->localData[0]->realVars before computing state derivatives?
927 39 data->callback->functionODE(data, threadData);
928
929 /* Difference between old and currend state derivatives */
930
2/2
✓ Branch 0 taken 78 times.
✓ Branch 1 taken 39 times.
 117 for(i=0; i < data->modelData->nStates; i++)
931 {
932 78 delta[i] = data->localData[0]->realVars[data->modelData->nStates + i] - yd[i];
933 }
934 39 printVector(OMC_LOG_DASSL_STATES, "dd", delta, data->modelData->nStates, *t);
935
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { success = 1; }
936 #if !defined(OMC_EMCC)
937 39 OMC_CATCH_INTERNAL(simulationJumpBuffer)
938 #endif
939
940
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
 39 if (!success) {
941 ✗ *ires = -1;
942 }
943
944 39 threadData->currentErrorStage = saveJumpState;
945
946 39 data->localData[0]->timeValue = timeBackup;
947
948
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 14 times.
 39 if (data->simulationInfo->currentContext == CONTEXT_ODE){
949 25 unsetContext(data);
950 }
951
952
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (measure_time_flag) rt_accumulate(SIM_TIMER_RESIDUALS);
953
1/2
✓ Branch 0 taken 39 times.
✗ Branch 1 not taken.
 39 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
954
955 39 return 0;
956 }
957
958 ✗ static int function_ZeroCrossingsDASSL(int *neqm, double *t, double *y, double *yp,
959 int *ng, double *gout, double *rpar, int* ipar)
960 {
961 ✗ DATA* data = (DATA*)(void*)((double**)rpar)[0];
962 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
963 ✗ threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
964
965 double timeBackup;
966 int saveJumpState;
967
968 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
969 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_EVENT);
970
971 ✗ if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
972 {
973 ✗ setContext(data, *t, CONTEXT_EVENTS);
974 }
975
976 ✗ saveJumpState = threadData->currentErrorStage;
977 ✗ threadData->currentErrorStage = ERROR_EVENTSEARCH;
978
979 ✗ timeBackup = data->localData[0]->timeValue;
980 ✗ data->localData[0]->timeValue = *t;
981
982 /* read input vars */
983 ✗ externalInputUpdate(data);
984 ✗ data->callback->input_function(data, threadData);
985 /* eval needed equations*/
986 ✗ data->callback->function_ZeroCrossingsEquations(data, threadData);
987
988 ✗ data->callback->function_ZeroCrossings(data, threadData, gout);
989
990 ✗ threadData->currentErrorStage = saveJumpState;
991 ✗ data->localData[0]->timeValue = timeBackup;
992
993 ✗ if (data->simulationInfo->currentContext == CONTEXT_EVENTS){
994 ✗ unsetContext(data);
995 }
996
997 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_EVENT);
998 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
999
1000 ✗ return 0;
1001 }
1002
1003 /* \fn jacA_symColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
1004 double *rpar, int* ipar)
1005 *
1006 * This function calculates the Jacobian matrix symbolically, exploiting the coloring.
1007 *
1008 * It handles all three symbolic evaluation modes transparently, since evalJacobian()
1009 * dispatches on the properties of the selected Jacobian:
1010 * coloredSymbolical -> forward (column) evaluation
1011 * coloredSymbolicalAdjoint -> adjoint (row) evaluation
1012 * bicoloredSymbolical -> bidirectional (column + row) evaluation
1013 */
1014 ✗ int jacA_symColored(double *t, double *y, double *yprime, double *delta,
1015 double *matrixA, double *cj, double *h, double *wt,
1016 double *rpar, int *ipar)
1017 {
1018 ✗ DATA* data = (DATA*)(void*)((double**)rpar)[0];
1019 ✗ threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
1020 ✗ JACOBIAN* jac = getSymbolicOdeJacobian(data);
1021 ✗ evalJacobian(data, threadData, jac, NULL, matrixA, TRUE);
1022
1023 ✗ return 0;
1024 }
1025
1026 /* \fn jacA_sym(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
1027 double *rpar, int* ipar)
1028 *
1029 *
1030 * This function calculates symbolically the jacobian matrix.
1031 */
1032 ✗ int jacA_sym(double *t, double *y, double *yprime, double *delta,
1033 double *matrixA, double *cj, double *h, double *wt, double *rpar,
1034 int *ipar)
1035 {
1036
1037 ✗ DATA* data = (DATA*)(void*)((double**)rpar)[0];
1038 ✗ threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
1039
1040 ✗ const int index = data->callback->INDEX_JAC_A;
1041 ✗ JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[index]);
1042 ✗ unsigned int columns = jac->sizeCols;
1043 ✗ unsigned int rows = jac->sizeRows;
1044 unsigned int i, j;
1045
1046 /* Evaluate constant equations if available */
1047 ✗ if (jac->constantEqns != NULL) {
1048 ✗ jac->constantEqns(data, threadData, jac, NULL);
1049 }
1050
1051 ✗ for(i=0; i < columns; i++)
1052 {
1053 ✗ jac->seedVars[i] = 1.0;
1054 ✗ data->callback->functionJacA_column(data, threadData, jac, NULL);
1055
1056 ✗ for(j = 0; j < rows; j++)
1057 {
1058 ✗ matrixA[i*columns+j] = jac->resultVars[j];
1059 }
1060
1061 ✗ jac->seedVars[i] = 0.0;
1062 }
1063
1064 ✗ return 0;
1065 }
1066
1067 /**
1068 * @brief Calculate Jacobian matrix numericaly.
1069 *
1070 * Calculate Jacobian matrix using forward finite differences.
1071 *
1072 * @param t Independent variable (time).
1073 * @param y Array with state variables, size dasslData->N.
1074 * @param yprime Array with state derivatives, size dasslData->N.
1075 * @param delta Previous f(t,y) - dy, Array of size dasslData->N.
1076 * @param matrixA On output contains values of Jacobian matrix
1077 * J = (∂F)/(∂y).
1078 * Array of size dasslData->N*dasslData->N, storing matrix in row-major order.
1079 * @param cj Specified by library interface. Given to residualFunction, which ignores it.
1080 * @param h Step size of DASSL solver.
1081 * @param wt Array with error weights, size dasslData->N.
1082 * @param rpar Struct storing user data.
1083 * Type: {DATA*, DASSL_DATA*, threadData_t*}
1084 * @param ipar Specified by library interface. Given to residualFunction, which ignores it.
1085 * @return int Return 0.
1086 */
1087 ✗ int jacA_num(double *t, double *y, double *yprime, double *delta,
1088 double *matrixA, double *cj, double *h, double *wt, double *rpar,
1089 int *ipar)
1090 {
1091 ✗ DATA* data = (DATA*)(void*)((double**)rpar)[0];
1092 ✗ DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
1093 threadData_t* threadData = (threadData_t*)(void*)((double**)rpar)[2];
1094
1095 double delta_hh, delta_hhh, deltaInv;
1096 double ysave;
1097 int ires;
1098 int col, row;
1099
1100 /* set context for the start values extrapolation of non-linear algebraic loops */
1101 ✗ setContext(data, *t, CONTEXT_JACOBIAN);
1102
1103 ✗ for(col=dasslData->N-1; col >= 0; col--)
1104 {
1105 ✗ delta_hhh = *h * yprime[col];
1106 ✗ delta_hh = numericalJacobianStep(y[col], delta_hhh, fabs(1. / wt[col]),
1107 ✗ dasslData->jacNominalFactor * dasslData->nominal[col]);
1108 ✗ delta_hh = (delta_hhh >= 0 ? delta_hh : -delta_hh);
1109 ✗ delta_hh = y[col] + delta_hh - y[col]; // Due to floating-point arithmetic rounding errors can result in: delta_hh != y[i] + delta_hh - y[i]
1110 ✗ deltaInv = 1. / delta_hh;
1111 ysave = y[col];
1112 ✗ y[col] += delta_hh;
1113
1114 ✗ (*dasslData->residualFunction)(t, y, yprime, cj, dasslData->newdelta, &ires, rpar, ipar);
1115 // TODO: What if residualFunction failed (ires=-1)?
1116
1117 ✗ increaseJacContext(data);
1118
1119 ✗ for(row = dasslData->N-1; row >= 0 ; row--)
1120 {
1121 ✗ matrixA[col*dasslData->N + row] = (dasslData->newdelta[row] - delta[row]) * deltaInv;
1122 // -I*cj will be added in callJacobian()
1123 }
1124 ✗ y[col] = ysave;
1125 }
1126
1127 ✗ return 0;
1128 }
1129
1130 /**
1131 * @brief Calculate colored Jacobian matrix numericaly.
1132 *
1133 * Calculate Jacobian matrix using forward finite differences and use coloring.
1134 *
1135 * @param t Independent variable (time).
1136 * @param y Array with state variables, size dasslData->N.
1137 * @param yprime Array with state derivatives, size dasslData->N.
1138 * @param delta Previous f(t,y) - dy, Array of size dasslData->N.
1139 * @param matrixA On output contains values of Jacobian matrix
1140 * J = (∂F)/(∂y).
1141 * Array of size dasslData->N*dasslData->N, storing matrix in row-major order.
1142 * @param cj Specified by library interface. Given to residualFunction, which ignores it.
1143 * @param h Step size of DASSL solver.
1144 * @param wt Array with error weights, size dasslData->N.
1145 * @param rpar Struct storing user data.
1146 * Type: {DATA*, DASSL_DATA*, threadData_t*}
1147 * @param ipar Specified by library interface. Given to residualFunction, which ignores it.
1148 * @return int Return 0.
1149 */
1150 14 int jacA_numColored(double *t, double *y, double *yprime, double *delta,
1151 double *matrixA, double *cj, double *h, double *wt,
1152 double *rpar, int *ipar)
1153 {
1154
1155 14 DATA* data = (DATA*)(void*)((double**)rpar)[0];
1156 14 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
1157 threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
1158
1159 14 const int index = data->callback->INDEX_JAC_A;
1160 14 JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
1161
1162 double delta_hhh;
1163 int ires;
1164 14 double* delta_hh = dasslData->delta_hh;
1165 14 double* ysave = dasslData->ysave;
1166
1167 unsigned int i,j,l,k,ii;
1168
1169 /* set context for the start values extrapolation of non-linear algebraic loops */
1170 14 setContext(data, *t, CONTEXT_JACOBIAN);
1171
1172
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 14 times.
 28 for(i = 0; i < jacobian->sparsePattern->maxColors; i++)
1173 {
1174
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 14 times.
 42 for(ii=0; ii < jacobian->sizeCols; ii++)
1175 {
1176
1/2
✓ Branch 0 taken 28 times.
✗ Branch 1 not taken.
 28 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
1177 {
1178 28 delta_hhh = *h * yprime[ii];
1179 28 delta_hh[ii] = numericalJacobianStep(y[ii], delta_hhh, fabs(1./wt[ii]),
1180 28 dasslData->jacNominalFactor * dasslData->nominal[ii]);
1181
1/2
✓ Branch 0 taken 28 times.
✗ Branch 1 not taken.
 28 delta_hh[ii] = (delta_hhh >= 0 ? delta_hh[ii] : -delta_hh[ii]);
1182 28 delta_hh[ii] = y[ii] + delta_hh[ii] - y[ii]; // Due to floating-point arithmetic rounding errors can result in: delta_hh[ii] != y[ii] + delta_hh[ii] - y[ii]
1183
1184 28 ysave[ii] = y[ii];
1185 28 y[ii] += delta_hh[ii];
1186
1187 28 delta_hh[ii] = 1. / delta_hh[ii];
1188 }
1189 }
1190 14 (*dasslData->residualFunction)(t, y, yprime, cj, dasslData->newdelta, &ires, rpar, ipar);
1191
1192 14 increaseJacContext(data);
1193
1194
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 14 times.
 42 for(ii = 0; ii < jacobian->sizeCols; ii++)
1195 {
1196
1/2
✓ Branch 0 taken 28 times.
✗ Branch 1 not taken.
 28 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
1197 {
1198 28 j = jacobian->sparsePattern->leadindex[ii];
1199
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 28 times.
 56 while(j < jacobian->sparsePattern->leadindex[ii+1])
1200 {
1201 28 l = jacobian->sparsePattern->index[j];
1202 28 k = l + ii*jacobian->sizeRows;
1203 28 matrixA[k] = (dasslData->newdelta[l] - delta[l]) * delta_hh[ii];
1204 // -I*cj will be added in callJacobian()
1205 28 j++;
1206 };
1207 28 y[ii] = ysave[ii];
1208 }
1209 }
1210 }
1211
1212 14 return 0;
1213 }
1214
1215 /**
1216 * @brief Calculate Jacobian.
1217 *
1218 * @param t Independent variable (time).
1219 * @param y Array with state variables, size dasslData->N.
1220 * @param yprime Array with state derivatives, size dasslData->N.
1221 * @param deltaD Previous F(t,y,y') := f(t,y) - y'
1222 * Array of size dasslData->N.
1223 * @param pd On output contains values of Jacobian matrix
1224 * J = (∂F)/(∂y) + cj * (∂F)/(∂y').
1225 * Array of size dasslData->N*dasslData->N, storing matrix in row-major order.
1226 * @param cj Coefficient from BDF method, cj = 1/alpha = h_n / alpha_n,0
1227 * @param h Step size of DASSL solver.
1228 * @param wt Array with error weights, size dasslData->N.
1229 * @param rpar Struct storing user data.
1230 * @param ipar Type: {DATA*, DASSL_DATA*, threadData_t*}
1231 * @return int Specified by library interface. Given to residualFunction, which ignores it.
1232 */
1233 14 static int callJacobian(double *t, double *y, double *yprime, double *deltaD,
1234 double *pd, double *cj, double *h, double *wt,
1235 double *rpar, int* ipar)
1236 {
1237 14 DATA* data = (DATA*)(void*)((double**)rpar)[0];
1238 14 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
1239 14 threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
1240 int i;
1241
1242 /* profiling */
1243
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
 14 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1244 14 rt_tick(SIM_TIMER_JACOBIAN);
1245
1246 /* Initialize dense Jacobian buffer. */
1247 14 memset(pd, 0, dasslData->N * dasslData->N * sizeof(double));
1248
1249 /* Compute J = (∂F)/(∂y) */
1250
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 14 times.
 14 if(dasslData->jacobianFunction(t, y, yprime, deltaD, pd, cj, h, wt, rpar, ipar))
1251 {
1252 ✗ throwStreamPrint(threadData, "Error, can not get Matrix A ");
1253 return 1;
1254 }
1255
1256 /* Compute J += cj * (∂F)/(∂y') = cj*(-I) */
1257
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 14 times.
 42 for(i = 0; i < dasslData->N*dasslData->N; i += dasslData->N + 1)
1258 {
1259 28 pd[i] -= *cj;
1260 }
1261
1262 /* debug */
1263 /* Compare evaluated Jacobian against a numerical reference.
1264 * Only meaningful when the configured method is not already numerical. */
1265
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
 14 if (OMC_ACTIVE_STREAM(OMC_LOG_JAC)
1266 ✗ && dasslData->dasslJacobian != COLOREDNUMJAC
1267 ✗ && dasslData->dasslJacobian != NUMJAC)
1268 {
1269 // print the analytical Jacobian for debugging
1270 ✗ printJacobianMatrix(OMC_LOG_JAC, "DASSL-Solver: analytical Jacobian pd (column-major)", pd,
1271 data, dasslData->N, *t);
1272
1273 // and print comparison to numerical Jacobian
1274 ✗ double* pdNumerical = (double*) calloc(dasslData->N * dasslData->N, sizeof(double));
1275 ✗ if (pdNumerical != NULL)
1276 {
1277 int row, col, k;
1278 double absDiff, relDiff;
1279 double maxAbsDiff = 0.0, maxRelDiff = 0.0;
1280 int maxAbsRow = 0, maxAbsCol = 0, maxRelRow = 0, maxRelCol = 0;
1281
1282 /* Compute numerical Jacobian ∂F/∂y using finite differences */
1283 ✗ jacA_num(t, y, yprime, deltaD, pdNumerical, cj, h, wt, rpar, ipar);
1284
1285 /* Apply the same cj * ∂F/∂y' = -cj*I correction */
1286 ✗ for (k = 0; k < dasslData->N * dasslData->N; k += dasslData->N + 1)
1287 {
1288 ✗ pdNumerical[k] -= *cj;
1289 }
1290
1291 /* Find maximum absolute and relative element-wise differences */
1292 ✗ for(col = 0; col < dasslData->N; col++)
1293 {
1294 ✗ for(row = 0; row < dasslData->N; row++)
1295 {
1296 ✗ int idx = col * dasslData->N + row;
1297 ✗ absDiff = fabs(pd[idx] - pdNumerical[idx]);
1298 ✗ relDiff = absDiff / fmax(fabs(pdNumerical[idx]), 1e-15);
1299 ✗ if(absDiff > maxAbsDiff) { maxAbsDiff = absDiff; maxAbsRow = row; maxAbsCol = col; }
1300 ✗ if(relDiff > maxRelDiff) { maxRelDiff = relDiff; maxRelRow = row; maxRelCol = col; }
1301 }
1302 }
1303
1304 ✗ infoStreamPrint(OMC_LOG_JAC, 1, "Jacobian verification: analytical vs. numerical");
1305 ✗ infoStreamPrint(OMC_LOG_JAC, 0,
1306 "Max absolute difference: %g at (row=%d:'%s', col=%d:'%s')",
1307 ✗ maxAbsDiff, maxAbsRow, data->modelData->realVarsData[maxAbsRow].info.name,
1308 ✗ maxAbsCol, data->modelData->realVarsData[maxAbsCol].info.name);
1309 ✗ infoStreamPrint(OMC_LOG_JAC, 0,
1310 "Max relative difference: %g at (row=%d:'%s', col=%d:'%s')",
1311 ✗ maxRelDiff, maxRelRow, data->modelData->realVarsData[maxRelRow].info.name,
1312 ✗ maxRelCol, data->modelData->realVarsData[maxRelCol].info.name);
1313 ✗ messageClose(OMC_LOG_JAC);
1314 ✗ free(pdNumerical);
1315 }
1316 }
1317
1318 /* set context for the start values extrapolation of non-linear algebraic loops */
1319 14 unsetContext(data);
1320
1321 /* profiling */
1322 14 rt_accumulate(SIM_TIMER_JACOBIAN);
1323
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
 14 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
1324
1325 return 0;
1326 }
1327
1328 #ifdef __cplusplus
1329 }
1330 #endif
1331