Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 835
Functions: 0.0% 0 / 0 / 20
Branches: 0.0% 0 / 0 / 463
OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.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 /*! \file ida_solver.c
29 */
30
31 #include <float.h>
32 #include <math.h>
33 #include <string.h>
34 #include <setjmp.h>
35
36 #include "omc_config.h"
37 #include "openmodelica.h"
38 #include "openmodelica_func.h"
39 #include "simulation_data.h"
40
41 #include "gc/omc_gc.h"
42 #include "util/context.h"
43 #include "util/omc_error.h"
44
45 #include "ida_solver.h"
46
47 #include "sundials_error.h"
48 #include "sundials_util.h"
49
50 #include "../arrayIndex.h"
51 #include "dae_mode.h"
52 #include "dassl.h"
53 #include "epsilon.h"
54 #include "external_input.h"
55 #include "simulation/jacobian_util.h"
56 #include "model_help.h"
57 #include "omc_math.h"
58 #include "simulation/options.h"
59 #include "simulation/results/simulation_result.h"
60 #include "simulation/simulation_runtime.h"
61 #include "solver_main.h"
62
63 #ifdef WITH_SUNDIALS
64
65
66 /* Private function prototypes */
67 static int callDenseJacobian(sunrealtype tt, sunrealtype cj, N_Vector yy,
68 N_Vector yp, N_Vector rr, SUNMatrix Jac,
69 void *user_data, N_Vector tmp1, N_Vector tmp2,
70 N_Vector tmp3);
71
72 static int callSparseJacobian(sunrealtype currentTime, sunrealtype cj,
73 N_Vector yy, N_Vector yp, N_Vector rr, SUNMatrix Jac, void *user_data,
74 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
75
76 static int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, void* user_data);
77 static int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void* userData);
78
79 static int getScalingFactors(DATA* data, IDA_SOLVER *idaData, SUNMatrix scaleMatrix);
80
81 static void idaScaleData(IDA_SOLVER *idaData);
82 static void idaReScaleData(IDA_SOLVER *idaData);
83 static void idaScaleVector(N_Vector vec, double* factors, unsigned int size);
84 static void idaReScaleVector(N_Vector vec, double* factors, unsigned int size);
85
86 int ida_event_update(DATA* data, threadData_t *threadData);
87
88 /* Static variables */
89 /* TODO: Don't use global variables */
90 static IDA_SOLVER *idaDataGlobal;
91
92
93 /**
94 * @brief Return true if flag signals success.
95 *
96 * If flag is IDA_SUCCESS or IDA_TSTOP_RETURN return true.
97 * Warnings (IDA_WARNING) or encountering roots (IDA_ROOT_RETURN) doesn't count as success.
98 *
99 * @param flag Value of IDA/IDAS flag.
100 * @return modelica_boolean Return true if flag signals success, otherwise return false.
101 */
102 ✗ modelica_boolean IDAflagIsSuccess(int flag) {
103 ✗ switch (flag)
104 {
105 case IDA_SUCCESS:
106 case IDA_TSTOP_RETURN:
107 case IDA_ROOT_RETURN:
108 return TRUE;
109 ✗ default:
110 ✗ return FALSE;
111 }
112 }
113
114
115 /**
116 * @brief Read the unknowns' nominal values into the tolerances and the scaling.
117 *
118 * Re-read by updateSolverNominals once initialization has computed the nominals
119 * that are parameter expressions.
120 *
121 * @param data Runtime data struct.
122 * @param threadData Thread data for error handling
123 * @param idaData IDA solver data with its arrays already allocated.
124 * @return int Returns 0 on success.
125 */
126 ✗ int ida_solver_setNominals(DATA* data, threadData_t *threadData, IDA_SOLVER* idaData)
127 {
128 int flag;
129 long int i;
130 ✗ double* abstol = N_VGetArrayPointer_Serial(idaData->absoluteTolerance);
131
132 ✗ infoStreamPrint(OMC_LOG_SOLVER, 1, "The relative tolerance is %g. Following absolute tolerances are used for the states: ", data->simulationInfo->tolerance);
133 ✗ for(i=0; i < data->modelData->nStates; ++i) {
134 ✗ const modelica_real nominal = getNominalFromScalarIdx(data->simulationInfo, data->modelData, VAR_KIND_STATE, i);
135 ✗ idaData->nominal[i] = fmax(fabs(nominal), 1e-32);
136 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, "%ld. %s -> %g", i+1, data->modelData->realVarsData[data->simulationInfo->realVarsReverseIndex[i].array_idx].info.name, idaData->nominal[i]);
137 }
138
139 /* daeMode: set nominal values for algebraic variables */
140 ✗ if (idaData->daeMode) {
141 ✗ getAlgebraicDAEVarNominals(data, idaData->nominal + data->modelData->nStates);
142 ✗ for(i=data->modelData->nStates; i < idaData->N; ++i) {
143 ✗ idaData->nominal[i] = fmax(fabs(idaData->nominal[i]), 1e-32);
144 }
145 }
146 /* multiply by tolerance to obtain a relative tolerace */
147 ✗ for(i=0; i < idaData->N; ++i) {
148 ✗ abstol[i] = idaData->nominal[i] * data->simulationInfo->tolerance;
149 }
150 ✗ messageClose(OMC_LOG_SOLVER);
151 ✗ flag = IDASVtolerances(idaData->ida_mem, data->simulationInfo->tolerance, idaData->absoluteTolerance);
152 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASVtolerances");
153
154 ✗ if (idaData->yScale != NULL) {
155 ✗ for(i=0; i < idaData->N; ++i) {
156 ✗ idaData->yScale[i] = idaData->nominal[i];
157 }
158 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "The scale factors for all ida states: ");
159 ✗ for (i=0; i < idaData->N; ++i) {
160 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, "%ld. scaleFactor: %g", i+1, idaData->yScale[i]);
161 }
162 ✗ messageClose(OMC_LOG_SOLVER_V);
163 }
164
165 ✗ return 0;
166 }
167
168 /**
169 * @brief Initialize main IDA data.
170 *
171 * Allocate memory for IDA_SOLVER struct and initialize IDA solver.
172 *
173 * @param data Runtime data struct.
174 * @param threadData Thread data for error handling
175 * @param solverInfo Information about main solver. Unused at the moment.
176 * @param idaData IDA solver data.
177 * @return int Returns 0 on success, aborts with throw on failure.
178 */
179 ✗ int ida_solver_initial(DATA* data, threadData_t *threadData,
180 SOLVER_INFO* solverInfo, IDA_SOLVER* idaData)
181 {
182 /* Variables */
183 int flag;
184 long int i;
185 int maxOrder;
186
187 /* Initialize constants */
188 ✗ idaData->setInitialSolution = FALSE;
189 ✗ idaData->homotopyRampActive = 0;
190 ✗ idaData->homotopyTramp = -1.0;
191
192 /* Instantiate IDA solver object */
193 /* Create the SUNDIALS context every other SUNDIALS object is created with */
194 ✗ flag = SUNContext_Create(SUN_COMM_NULL, &idaData->sunctx);
195 ✗ assertStreamPrint(threadData, flag == SUN_SUCCESS, "SUNDIALS_ERROR: SUNContext_Create failed.");
196 ✗ sundialsSilenceLogger(idaData->sunctx);
197
198 /* Set error handler */
199 ✗ flag = SUNContext_PushErrHandler(idaData->sunctx, sundialsErrorHandlerFunction, idaData);
200 ✗ assertStreamPrint(threadData, flag == SUN_SUCCESS, "SUNDIALS_ERROR: SUNContext_PushErrHandler failed.");
201
202 ✗ idaData->ida_mem = IDACreate(idaData->sunctx);
203 ✗ if (idaData->ida_mem == NULL) {
204 ✗ throwStreamPrint(threadData, "##IDA## Initialization of IDA solver failed!");
205 }
206
207 ✗ idaData->residualFunction = residualFunctionIDA;
208
209 /* Start measuring time */ /* TODO: Why start here? */
210 ✗ if (measure_time_flag) {
211 ✗ rt_tick(SIM_TIMER_SOLVER);
212 }
213
214 /* change parameter for DAE mode */
215 ✗ if (compiledInDAEMode) {
216 ✗ idaData->daeMode = TRUE;
217 ✗ idaData->N = (long int) data->modelData->nStates + data->simulationInfo->daeModeData->nAlgebraicDAEVars;
218 }
219 else {
220 ✗ idaData->daeMode = FALSE;
221 ✗ idaData->N = (long int) data->modelData->nStates;
222 }
223 ✗ infoStreamPrint(OMC_LOG_SOLVER, 1, "## IDA ## Initializing solver of size %ld %s.", idaData->N, idaData->daeMode?"in DAE mode":"");
224 ✗ idaData->NNZ = -1;
225
226 /* initialize states and der(states) */
227 ✗ if (idaData->daeMode)
228 {
229 ✗ idaData->states = (double*) malloc(idaData->N*sizeof(double));
230 ✗ idaData->statesDer = (double*) calloc(idaData->N,sizeof(double));
231
232 ✗ memcpy(idaData->states, data->localData[0]->realVars, sizeof(double)*data->modelData->nStates);
233 // and also algebraic vars
234 ✗ getAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
235 ✗ memcpy(idaData->statesDer, data->localData[0]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
236
237 ✗ idaData->y = N_VMake_Serial(idaData->N, idaData->states, idaData->sunctx);
238 ✗ idaData->yp = N_VMake_Serial(idaData->N, idaData->statesDer, idaData->sunctx);
239 }
240 else {
241 ✗ idaData->states = NULL;
242 ✗ idaData->statesDer = NULL;
243 ✗ idaData->y = N_VMake_Serial(idaData->N, data->localData[0]->realVars, idaData->sunctx);
244 ✗ idaData->yp = N_VMake_Serial(idaData->N, data->localData[0]->realVars + data->modelData->nStates, idaData->sunctx);
245 }
246
247 ✗ flag = IDAInit(idaData->ida_mem, idaData->residualFunction,
248 ✗ data->simulationInfo->startTime, idaData->y, idaData->yp);
249 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAInit");
250
251 /* Allocate memory for jacobians calculation */
252 ✗ idaData->ysave = (double*) malloc(idaData->N*sizeof(double));
253 ✗ idaData->ypsave = (double*) malloc(idaData->N*sizeof(double));
254 ✗ idaData->delta_hh = (double*) malloc(idaData->N*sizeof(double));
255 ✗ idaData->nominal = (double*) malloc(idaData->N*sizeof(double));
256 ✗ idaData->newdelta = N_VNew_Serial(idaData->N, idaData->sunctx);
257
258 /* Allocate memory for linear solver */
259 ✗ idaData->y_linSol = N_VNew_Serial(idaData->N, idaData->sunctx);
260
261 /* Set user data */
262 ✗ idaData->userData = (IDA_USERDATA*) malloc(sizeof(IDA_USERDATA));
263 ✗ idaData->userData->data = data;
264 ✗ idaData->userData->threadData = threadData;
265 ✗ flag = IDASetUserData(idaData->ida_mem, idaData);
266 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetUserData");
267
268 ✗ idaData->jacNominalFactor = omc_flag[FLAG_JACOBIAN_NOMINAL_FACTOR]
269 ✗ ? atof(omc_flagValue[FLAG_JACOBIAN_NOMINAL_FACTOR]) : 1.0;
270
271 ✗ idaData->absoluteTolerance = N_VNew_Serial(idaData->N, idaData->sunctx);
272 ✗ idaData->id = NULL;
273
274 ✗ if (omc_flag[FLAG_IDA_SCALING]) { /* idaNoScaling */
275 /* allocate memory for scaling */
276 ✗ idaData->yScale = (double*) malloc(idaData->N*sizeof(double));
277 ✗ idaData->ypScale = (double*) malloc(idaData->N*sizeof(double));
278 ✗ idaData->resScale = (double*) malloc(idaData->N*sizeof(double));
279
280 ✗ for(i=0; i < idaData->N; ++i) {
281 ✗ idaData->ypScale[i] = 1.0; // TODO: 1 is not a good scaling value. Use something like nominal value / number of intervals
282 }
283 } else {
284 ✗ idaData->yScale = NULL;
285 ✗ idaData->ypScale = NULL;
286 ✗ idaData->resScale = NULL;
287 }
288
289 ✗ ida_solver_setNominals(data, threadData, idaData);
290 /* initialize */
291 ✗ idaData->useScaling = TRUE;
292
293 /* Set root functions unless flag FLAG_NO_ROOTFINDING is set */
294 ✗ if (!omc_flag[FLAG_NO_ROOTFINDING]) {
295 ✗ solverInfo->solverRootFinding = 1;
296 ✗ flag = IDARootInit(idaData->ida_mem, data->modelData->nZeroCrossings, rootsFunctionIDA);
297 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDARootInit");
298 }
299 else {
300 ✗ solverInfo->solverRootFinding = 0;
301 }
302 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "IDA uses internal root finding method %s", solverInfo->solverRootFinding?"YES":"NO");
303
304 /* Define maximum integration order of IDA */
305 ✗ if (omc_flag[FLAG_MAX_ORDER]) {
306 ✗ maxOrder = atoi(omc_flagValue[FLAG_MAX_ORDER]);
307
308 ✗ flag = IDASetMaxOrd(idaData->ida_mem, maxOrder);
309 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetMaxOrd");
310 } else {
311 maxOrder = 5; /* Default max order for IDA */
312 }
313 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Maximum integration order %d", maxOrder);
314
315 /* if FLAG_NOEQUIDISTANT_GRID is set, choose ida step method */
316 ✗ if (omc_flag[FLAG_NOEQUIDISTANT_GRID]) {
317 ✗ idaData->internalSteps = 1; /* TRUE */
318 ✗ solverInfo->solverNoEquidistantGrid = TRUE;
319 } else {
320 ✗ idaData->internalSteps = 0; /* FALSE */
321 }
322 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "use equidistant time grid %s", idaData->internalSteps?"NO":"YES");
323
324 /* check if Flags FLAG_NOEQUIDISTANT_OUT_FREQ or FLAG_NOEQUIDISTANT_OUT_TIME are set */
325 ✗ if (idaData->internalSteps) {
326 ✗ if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ]) {
327 ✗ idaData->stepsFreq = atoi(omc_flagValue[FLAG_NOEQUIDISTANT_OUT_FREQ]);
328 ✗ } else if (omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]) {
329 ✗ idaData->stepsTime = atof(omc_flagValue[FLAG_NOEQUIDISTANT_OUT_TIME]);
330 ✗ flag = IDASetMaxStep(idaData->ida_mem, idaData->stepsTime);
331 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetMaxStep");
332 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "maximum step size %g", idaData->stepsTime);
333 } else {
334 ✗ idaData->stepsFreq = 1;
335 ✗ idaData->stepsTime = 0.0;
336 }
337
338 ✗ if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ] && omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]) {
339 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The flags are \"noEquidistantOutputFrequency\" "
340 "and \"noEquidistantOutputTime\" are in opposition "
341 "to each other. The flag \"noEquidistantOutputFrequency\" superiors.");
342 }
343 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "as the output frequency control is used: %d", idaData->stepsFreq);
344 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "as the output frequency time step control is used: %f", idaData->stepsTime);
345 }
346
347 /* if FLAG_IDA_LS is set, choose ida linear solver method */
348 ✗ if (omc_flag[FLAG_IDA_LS]) {
349 ✗ for (i=1; i< IDA_LS_MAX; i++) {
350 ✗ if (!strcmp((const char*)omc_flagValue[FLAG_IDA_LS], IDA_LS_METHOD_NAME[i])) {
351 ✗ idaData->linearSolverMethod = (enum IDA_LS)i;
352 ✗ break;
353 }
354 }
355 ✗ if (idaData->linearSolverMethod == IDA_LS_UNKNOWN) {
356 ✗ if (OMC_ACTIVE_WARNING_STREAM(OMC_LOG_SOLVER)) {
357 ✗ warningStreamPrint(OMC_LOG_SOLVER, 1, "unrecognized ida linear solver method %s, current options are:", (const char*)omc_flagValue[FLAG_IDA_LS]);
358 ✗ for(i=1; i < IDA_LS_MAX; ++i) {
359 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "%-15s [%s]", IDA_LS_METHOD_NAME[i], IDA_LS_METHOD_DESC[i]);
360 }
361 ✗ messageCloseWarning(OMC_LOG_SOLVER);
362 }
363 ✗ throwStreamPrint(threadData,"unrecognized ida linear solver method %s", (const char*)omc_flagValue[FLAG_IDA_LS]);
364 }
365 } else {
366 ✗ idaData->linearSolverMethod = IDA_LS_KLU;
367 }
368
369 /* Choose and initialize the ODE Jacobian. The mapping from the `-jacobian` flag to
370 * the forward / adjoint / bidirectional Jacobian is shared with DASSL and GBODE. */
371 ✗ idaData->jacobianMethod = getRequestedJacobianMethod(threadData);
372 ✗ JACOBIAN* jacobian = initSymbolicOdeJacobian(data, threadData, &idaData->jacobianMethod, FALSE);
373 /* IDA always needs a column oriented pattern. */
374 ✗ const SPARSE_PATTERN* cscPattern = getJacobianCscPattern(jacobian);
375
376 // change IDA specific jacobian method
377 ✗ if(idaData->jacobianMethod == SYMJAC) {
378 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Symbolic Jacobians without coloring are currently not supported by IDA."
379 " Colored symbolical Jacobian will be used.");
380 ✗ idaData->jacobianMethod = COLOREDSYMJAC;
381 ✗ }else if(idaData->jacobianMethod == NUMJAC) {
382 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Numerical Jacobians without coloring are currently not supported by IDA."
383 " Colored numerical Jacobian will be used.");
384 ✗ idaData->jacobianMethod = COLOREDNUMJAC;
385 ✗ }else if(idaData->jacobianMethod == INTERNALNUMJAC && idaData->linearSolverMethod == IDA_LS_KLU) {
386 ✗ if ((!idaData->daeMode && jacobian->sparsePattern == NULL) || (idaData->daeMode && data->simulationInfo->daeModeData->sparsePattern == NULL)) {
387 ✗ throwStreamPrint(threadData, "##IDA## Internal Numerical Jacobians require a sparse pattern for the jacobian but no sparse pattern is generated.");
388 } else {
389 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Internal Numerical Jacobians without coloring are currently not supported by IDA with KLU."
390 " Colored numerical Jacobian will be used.");
391 ✗ idaData->jacobianMethod = COLOREDNUMJAC;
392 }
393 }
394
395 /* Set NNZ */
396 ✗ if (idaData->daeMode) {
397 ✗ idaData->NNZ = data->simulationInfo->daeModeData->sparsePattern->nnz;
398 } else {
399 // could also be jacobian->sparsePattern->nnz but the cscPattern is the one that is actually used.
400 ✗ idaData->NNZ = cscPattern->nnz;
401 }
402
403 ✗ switch (idaData->linearSolverMethod){
404 ✗ case IDA_LS_SPGMR:
405 ✗ idaData->J = NULL;
406 ✗ idaData->linSol = SUNLinSol_SPGMR(idaData->y_linSol, SUN_PREC_NONE, idaData->N, idaData->sunctx);
407 ✗ if (idaData->linSol == NULL) {
408 ✗ throwStreamPrint(threadData, "##IDA## In function SUNLinSol_SPGMR: Input incompatible.");
409 }
410 ✗ idaData->jacobianMethod = INTERNALNUMJAC;
411 ✗ break;
412 ✗ case IDA_LS_SPBCG:
413 ✗ idaData->J = NULL;
414 ✗ idaData->linSol = SUNLinSol_SPBCGS(idaData->y_linSol, SUN_PREC_NONE, idaData->N, idaData->sunctx);
415 ✗ if (idaData->linSol == NULL) {
416 ✗ throwStreamPrint(threadData, "##IDA## In function SUNLinSol_SPBCGS: Input incompatible.");
417 }
418 ✗ idaData->jacobianMethod = INTERNALNUMJAC;
419 ✗ break;
420 ✗ case IDA_LS_SPTFQMR:
421 ✗ idaData->J = NULL;
422 ✗ idaData->linSol = SUNLinSol_SPTFQMR(idaData->y_linSol, SUN_PREC_NONE, idaData->N, idaData->sunctx);
423 ✗ if (idaData->linSol == NULL) {
424 ✗ throwStreamPrint(threadData, "##IDA## In function SUNLinSol_SPTFQMR: Input incompatible.");
425 }
426 ✗ idaData->jacobianMethod = INTERNALNUMJAC;
427 ✗ break;
428 ✗ case IDA_LS_DENSE:
429 ✗ idaData->J = SUNDenseMatrix(idaData->N, idaData->N, idaData->sunctx);
430 ✗ idaData->linSol = SUNLinSol_Dense(idaData->y_linSol, idaData->J, idaData->sunctx);
431 ✗ if (idaData->linSol == NULL) {
432 ✗ throwStreamPrint(threadData, "##IDA## In function SUNLinSol_Dense: Input incompatible.");
433 }
434 break;
435 ✗ case IDA_LS_KLU:
436 /* Set KLU after initialized sparse pattern of the jacobian for nnz */
437 ✗ if (idaData->NNZ < 0) {
438 ✗ throwStreamPrint(threadData, "##IDA## idaData->NNZ not set.");
439 }
440 ✗ idaData->J = SUNSparseMatrix(idaData->N, idaData->N, idaData->NNZ + idaData->N, SUN_CSC_MAT, idaData->sunctx);
441 ✗ idaData->linSol = SUNLinSol_KLU(idaData->y_linSol, idaData->J, idaData->sunctx);
442 ✗ if (idaData->linSol == NULL) {
443 ✗ throwStreamPrint(threadData, "##IDA## In function SUNLinSol_KLU: Input incompatible.");
444 }
445 break;
446 ✗ default:
447 ✗ throwStreamPrint(threadData,"unrecognized linear solver method %s", (const char*)omc_flagValue[FLAG_IDA_LS]);
448 break;
449 }
450
451 ✗ flag = IDASetLinearSolver(idaData->ida_mem, idaData->linSol, idaData->J);
452 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDALS_FLAG, "IDASetLinearSolver");
453 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "IDA linear solver method selected %s", IDA_LS_METHOD_DESC[idaData->linearSolverMethod]);
454
455 /* Set Jacobian function */
456 ✗ idaData->scaleMatrix = NULL; /* allocated on demand, see getScalingFactors */
457
458 /* Use sparse jacobian evaluation */
459 ✗ if (idaData->linearSolverMethod == IDA_LS_KLU) {
460 /* Set Jacobian function for matrix based linear solvers */
461 ✗ switch (idaData->jacobianMethod){
462 ✗ case SYMJAC:
463 case NUMJAC:
464 case COLOREDSYMJAC:
465 case COLOREDSYMJACADJ:
466 case BICOLOREDSYMJAC:
467 case COLOREDNUMJAC:
468 ✗ flag = IDASetJacFn(idaData->ida_mem, callSparseJacobian);
469
470 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDALS_FLAG, "IDASetJacFn");
471 ✗ if (omc_flag[FLAG_IDA_SCALING]) {
472 ✗ idaData->scaleMatrix = SUNSparseMatrix(idaData->N, idaData->N, idaData->NNZ + idaData->N, SUN_CSC_MAT, idaData->sunctx);
473 }
474 break;
475 ✗ default:
476 ✗ throwStreamPrint(threadData,"For the klu solver jacobian calculation method has to be one of %s, %s, %s or %s",
477 JACOBIAN_METHOD_NAME[COLOREDSYMJAC], JACOBIAN_METHOD_NAME[COLOREDSYMJACADJ],
478 JACOBIAN_METHOD_NAME[BICOLOREDSYMJAC], JACOBIAN_METHOD_NAME[COLOREDNUMJAC]);
479 break;
480 }
481 /* Use dense jacobian evaluation */
482 } else {
483 ✗ switch (idaData->jacobianMethod){
484 ✗ case SYMJAC:
485 case NUMJAC:
486 case COLOREDSYMJAC:
487 case COLOREDSYMJACADJ:
488 case BICOLOREDSYMJAC:
489 case COLOREDNUMJAC:
490 ✗ flag = IDASetJacFn(idaData->ida_mem, callDenseJacobian);
491 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDALS_FLAG, "IDASetJacFn");
492 ✗ break;
493 case INTERNALNUMJAC:
494 /* TODO: Set a preconditioner if possible */
495 break;
496 ✗ default:
497 ✗ throwStreamPrint(threadData,"unrecognized jacobian calculation method %s", (const char*)omc_flagValue[FLAG_JACOBIAN]);
498 break;
499 }
500 }
501 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Jacobian is calculated by \"%s\"", JACOBIAN_METHOD_DESC[idaData->jacobianMethod]);
502
503 /* Set max error test fails */
504 ✗ if (omc_flag[FLAG_IDA_MAXERRORTESTFAIL])
505 {
506 ✗ int maxErrorTestFails = atoi(omc_flagValue[FLAG_IDA_MAXERRORTESTFAIL]);
507 ✗ flag = IDASetMaxErrTestFails(idaData->ida_mem, maxErrorTestFails);
508 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetMaxErrTestFails");
509 }
510
511 /* set maximum number of nonlinear solver iterations at one step */
512 ✗ if (omc_flag[FLAG_IDA_MAXNONLINITERS])
513 {
514 ✗ int maxNonlinIters = atoi(omc_flagValue[FLAG_IDA_MAXNONLINITERS]);
515 ✗ flag = IDASetMaxNonlinIters(idaData->ida_mem, maxNonlinIters);
516 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetMaxNonlinIters");
517 }
518
519 /* maximum number of nonlinear solver convergence failures at one step */
520 ✗ if (omc_flag[FLAG_IDA_MAXCONVFAILS])
521 {
522 ✗ int maxConvFails = atoi(omc_flagValue[FLAG_IDA_MAXCONVFAILS]);
523 ✗ flag = IDASetMaxConvFails(idaData->ida_mem, maxConvFails);
524 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetMaxConvFails");
525 }
526
527 /* safety factor in the nonlinear convergence test */
528 ✗ if (omc_flag[FLAG_IDA_NONLINCONVCOEF])
529 {
530 ✗ double nonlinConvCoef = atof(omc_flagValue[FLAG_IDA_NONLINCONVCOEF]);
531 ✗ flag = IDASetNonlinConvCoef(idaData->ida_mem, nonlinConvCoef);
532 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetNonlinConvCoef");
533 }
534
535 /* configure algebraic variables as such */
536 ✗ if (idaData->daeMode) {
537 ✗ if (omc_flag[FLAG_NO_SUPPRESS_ALG]) {
538 ✗ flag = IDASetSuppressAlg(idaData->ida_mem, 1 /* TRUE */);
539 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetSuppressAlg");
540 }
541 ✗ idaData->id = N_VNew_Serial(idaData->N, idaData->sunctx);
542 ✗ for (i=0; i<idaData->N; ++i) {
543 ✗ NV_Ith_S(idaData->id, i) = (i<data->modelData->nStates)? 1.0: 0.0;
544 }
545
546 ✗ flag = IDASetId(idaData->ida_mem, idaData->id);
547 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetId");
548 }
549
550 /* define initial step size */
551 ✗ if (omc_flag[FLAG_INITIAL_STEP_SIZE]) {
552 ✗ double initialStepSize = atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]);
553
554 ✗ assertStreamPrint(threadData, initialStepSize >= DASSL_STEP_EPS, "Selected initial step size %e is too small.", initialStepSize);
555
556 ✗ flag = IDASetInitStep(idaData->ida_mem, initialStepSize);
557 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetInitStep");
558 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "initial step size: %g", initialStepSize);
559 } else {
560 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "initial step size is set automatically.");
561 }
562
563 /* Initialize sensitivities analysis */
564 ✗ idaData->idaSmode = omc_flag[FLAG_IDAS] ? 1 : 0;
565
566 ✗ if (idaData->idaSmode) {
567 ✗ idaData->Ns = data->modelData->nSensitivityParamVars;
568 ✗ idaData->yS = N_VCloneVectorArray(idaData->Ns, idaData->y);
569 ✗ idaData->ySp = N_VCloneVectorArray(idaData->Ns, idaData->yp);
570
571 ✗ for (i=0; i<idaData->Ns; ++i) {
572 ✗ N_VConst_Serial(0.0, idaData->yS[i]);
573 ✗ N_VConst_Serial(0.0, idaData->ySp[i]);
574 }
575
576 ✗ flag = IDASensInit(idaData->ida_mem, idaData->Ns, IDA_SIMULTANEOUS, NULL, idaData->yS, idaData->ySp);
577 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASensInit");
578
579 ✗ flag = IDASetSensParams(idaData->ida_mem, data->simulationInfo->realParameter, NULL, data->simulationInfo->sensitivityParList);
580 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetSensParams");
581 ✗ flag = IDASetSensDQMethod(idaData->ida_mem, IDA_FORWARD, 0);
582 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetSensDQMethod");
583
584 ✗ flag = IDASensEEtolerances(idaData->ida_mem);
585 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASensEEtolerances");
586 /*
587 flag = IDASetSensErrCon(idaData->ida_mem, TRUE);
588 checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetSensErrCon");
589 */
590 /* allocate result workspace */
591 ✗ idaData->ySResult = N_VNewVectorArray(idaData->Ns, idaData->sunctx);
592 ✗ for(i = 0; i < idaData->Ns; ++i)
593 {
594 ✗ idaData->ySResult[i] = N_VMake_Serial(idaData->N,
595 ✗ data->simulationInfo->sensitivityMatrix + i*idaData->N,
596 idaData->sunctx);
597 }
598 }
599 ✗ if (compiledInDAEMode){
600 ✗ idaDataGlobal = idaData;
601 ✗ data->callback->functionDAE = ida_event_update;
602 }
603 ✗ messageClose(OMC_LOG_SOLVER);
604
605 ✗ if (measure_time_flag) rt_clear(SIM_TIMER_SOLVER); /* TODO Initialization should not add to this timer... */
606
607 ✗ return 0;
608 }
609
610 /**
611 * @brief Deinitialize IDA data.
612 *
613 * @param idaData Pointer to IDA solver data struct.
614 */
615 ✗ void ida_solver_deinitial(IDA_SOLVER *idaData)
616 {
617 ✗ if (omc_flag[FLAG_IDA_SCALING]) {
618 /* free scaling data */
619 ✗ free(idaData->yScale);
620 ✗ free(idaData->ypScale);
621 ✗ free(idaData->resScale);
622 ✗ SUNMatDestroy(idaData->scaleMatrix);
623 }
624
625 /* free work arrays */
626 ✗ free(idaData->userData);
627 ✗ free(idaData->ysave);
628 ✗ free(idaData->ypsave);
629 ✗ free(idaData->delta_hh);
630
631 /* Free linear solver data */
632 ✗ N_VDestroy_Serial(idaData->y_linSol);
633 ✗ SUNMatDestroy(idaData->J);
634 ✗ SUNLinSolFree(idaData->linSol);
635
636 /* Free dae-mode data */
637 ✗ if (idaData->daeMode) {
638 ✗ free(idaData->states);
639 ✗ free(idaData->statesDer);
640 }
641
642 /* Free sensitivity-mode data */
643 ✗ if (idaData->idaSmode) {
644 ✗ N_VDestroyVectorArray(idaData->yS, idaData->Ns);
645 ✗ N_VDestroyVectorArray(idaData->ySp, idaData->Ns);
646 ✗ N_VDestroyVectorArray(idaData->ySResult, idaData->Ns);
647 }
648
649 ✗ free(idaData->nominal);
650 ✗ N_VDestroy_Serial(idaData->absoluteTolerance);
651 ✗ if (idaData->id != NULL) {
652 ✗ N_VDestroy_Serial(idaData->id);
653 }
654 ✗ N_VDestroy_Serial(idaData->newdelta);
655
656 ✗ IDAFree(&idaData->ida_mem);
657
658 ✗ SUNContext_Free(&idaData->sunctx);
659 ✗ }
660
661
662 /**
663 * @brief EventHandle for DAE mode.
664 *
665 * Handles events by reinitialize main IDA solver, initializing next step and
666 * evaluate DAE residual equations.
667 *
668 * @param data Runtime data struct.
669 * @param threadData Thread data for error handling.
670 * @return int Return 0.
671 */
672 ✗ int ida_event_update(DATA* data, threadData_t *threadData)
673 {
674 ✗ IDA_SOLVER *idaData = idaDataGlobal;
675 int flag;
676 long nonLinIters;
677 double init_h;
678
679 ✗ if (!compiledInDAEMode){
680 ✗ throwStreamPrint(threadData, "Function ida_event_update only callable in DAE mode");
681 }
682
683 ✗ data->simulationInfo->needToIterate = 0 /* FALSE */;
684
685 ✗ memcpy(idaData->states, data->localData[0]->realVars, sizeof(double)*data->modelData->nStates);
686 ✗ getAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
687 ✗ memcpy(idaData->statesDer, data->localData[0]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
688
689 /* update inner algebraic get new values from data */
690 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
691 ✗ evaluateDAEResiduals_wrapperEventUpdate(data, threadData);
692 ✗ getAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
693 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
694
695 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## do event update at %.15g", data->localData[0]->timeValue);
696 ✗ memcpy(idaData->states, data->localData[0]->realVars, sizeof(double)*data->modelData->nStates);
697 ✗ memcpy(idaData->statesDer, data->localData[0]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
698 ✗ memcpy(NV_DATA_S(idaData->y), idaData->states, idaData->N);
699 ✗ memcpy(NV_DATA_S(idaData->yp), idaData->statesDer, idaData->N);
700 ✗ flag = IDAReInit(idaData->ida_mem,
701 ✗ data->localData[0]->timeValue,
702 idaData->y,
703 idaData->yp);
704 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAReInit");
705
706 /* get initial step to provide a direction of the solution */
707 ✗ flag = IDAGetActualInitStep(idaData->ida_mem, &init_h);
708 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAGetActualInitStep");
709 /* provide a feasible step-size if it's too small */
710 ✗ if (init_h < DBL_EPSILON){
711 ✗ init_h = DBL_EPSILON;
712 ✗ flag = IDASetInitStep(idaData->ida_mem, init_h);
713 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetInitStep");
714 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## corrected step-size at %.15g", init_h);
715 }
716
717 /* increase limits of the non-linear solver */
718 ✗ IDASetMaxNumStepsIC(idaData->ida_mem, 2*idaData->N*10);
719 ✗ IDASetMaxNumJacsIC(idaData->ida_mem, 2*idaData->N*10);
720 ✗ IDASetMaxNumItersIC(idaData->ida_mem, 2*idaData->N*10);
721 /* Calc Consistent y_algebraic and y_prime with current y */
722 ✗ flag = IDACalcIC(idaData->ida_mem, IDA_YA_YDP_INIT, data->localData[0]->timeValue+init_h);
723
724 /* debug */
725 ✗ IDAGetNumNonlinSolvIters(idaData->ida_mem, &nonLinIters);
726 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## IDACalcIC run status %d.\nIterations : %ld\n", flag, nonLinIters);
727
728 /* try again without line search if first try fails */
729 ✗ if (!IDAflagIsSuccess(flag)){
730 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## first event iteration failed. Start next try without line search!");
731 ✗ IDASetLineSearchOffIC(idaData->ida_mem, 1 /* TRUE */);
732 ✗ flag = IDACalcIC(idaData->ida_mem, IDA_YA_YDP_INIT, data->localData[0]->timeValue+data->simulationInfo->tolerance);
733 ✗ IDAGetNumNonlinSolvIters(idaData->ida_mem, &nonLinIters);
734 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## IDACalcIC run status %d.\nIterations : %ld\n", flag, nonLinIters);
735 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAGetNumNonlinSolvIters");
736 }
737 /* obtain consistent values of y_algebraic and y_prime */
738 ✗ IDAGetConsistentIC(idaData->ida_mem, idaData->y, idaData->yp);
739
740 /* update inner algebraic variables */
741 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
742 ✗ evaluateDAEResiduals_wrapperEventUpdate(data, threadData);
743
744 ✗ memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
745 // and also algebraic vars
746 ✗ setAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
747 ✗ memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
748 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
749
750 /* reset initial step size again to default */
751 ✗ IDASetInitStep(idaData->ida_mem, 0.0);
752
753 ✗ return 0;
754 }
755
756 /**
757 * @brief Main IDA solver step.
758 *
759 * Call IDASolve to make solver steps.
760 *
761 * @param data Runtime data struct.
762 * @param threadData Thread data for error handling.
763 * @param solverInfo Main ODE/DAE solver info.
764 * @return int Return 0 on success or IDA flag on failure.
765 */
766 /* Activate the daeMode homotopy ramp after a singular initial DAE Jacobian was
767 detected (IDA_LSETUP_FAIL). Sets the ramp window and caps the integrator step
768 so the lambda 0->1 transition is resolved gradually regardless of the
769 requested output interval count; the cap is lifted again in ida_solver_step
770 once the ramp is complete. */
771 ✗ static void idaActivateHomotopyRamp(IDA_SOLVER *idaData, DATA *data)
772 {
773 ✗ const char *e = getenv("OMC_DAE_HOMOTOPY_TRAMP");
774 ✗ idaData->homotopyTramp = (e != NULL) ? atof(e)
775 ✗ : 0.1 * (data->simulationInfo->stopTime - data->simulationInfo->startTime);
776 ✗ idaData->homotopyRampActive = 1;
777 ✗ if (idaData->homotopyTramp > 0.0)
778 ✗ IDASetMaxStep(idaData->ida_mem, idaData->homotopyTramp / 50.0);
779 ✗ }
780
781 /* Whether a failed IDASolve is the degenerate initial operating point the ramp
782 above recovers from. The Jacobian there is singular, but a numerical one is
783 only singular to the accuracy of its finite difference: a step that reaches
784 past the regularized region gives a matrix that factorizes, and the corrector
785 fails to converge instead. lambda is ramped over [startTime, startTime +
786 t_ramp], so this only applies while still at startTime. */
787 static modelica_boolean idaHomotopyRampRecovers(IDA_SOLVER *idaData, DATA *data,
788 SOLVER_INFO *solverInfo, int flag)
789 {
790 ✗ if (!idaData->daeMode || idaData->homotopyRampActive ||
791 ✗ data->simulationInfo->homotopySteps <= 0)
792 return FALSE;
793 ✗ if (flag != IDA_LSETUP_FAIL && flag != IDA_CONV_FAIL && flag != IDA_ERR_FAIL)
794 return FALSE;
795 ✗ return solverInfo->currentTime <= data->simulationInfo->startTime;
796 }
797
798 ✗ int ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
799 {
800 double tout = 0;
801 int i = 0, flag;
802 ✗ int retVal = 0, finished = 0 /* FALSE */;
803 int saveJumpState;
804 long int tmp;
805 static unsigned int stepsOutputCounter = 1;
806 int stepsMode; /* Has to be IDA_NORMAL (1) or IDA_ONE_STEP (2) */
807 ✗ int restartAfterLSFail = 0;
808 modelica_boolean rampRecovery;
809
810 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
811
812 ✗ IDA_SOLVER *idaData = (IDA_SOLVER*) solverInfo->solverData;
813
814 ✗ SIMULATION_DATA *sData = data->localData[0];
815 SIMULATION_DATA *sDataOld = data->localData[1];
816 MODEL_DATA *mData = (MODEL_DATA*) data->modelData;
817
818 /* DAE-mode homotopy ramp: lambda is set smoothly as a function of time in the
819 residual callback (residualFunctionIDA). Once the ramp window has elapsed,
820 lift the step-size cap that was applied while ramping and pin lambda to 1
821 (the actual model) for the rest of the integration. */
822 ✗ if (idaData->homotopyRampActive && idaData->homotopyTramp > 0.0 &&
823 ✗ solverInfo->currentTime >= data->simulationInfo->startTime + idaData->homotopyTramp)
824 {
825 ✗ IDASetMaxStep(idaData->ida_mem, 0.0); /* 0 = no limit */
826 ✗ data->simulationInfo->lambda = 1.0;
827 ✗ idaData->homotopyRampActive = 0;
828 }
829
830
831 /* alloc all work arrays */
832 ✗ if (!idaData->daeMode)
833 {
834 ✗ N_VSetArrayPointer_Serial(data->localData[0]->realVars, idaData->y);
835 ✗ N_VSetArrayPointer_Serial(data->localData[1]->realVars + data->modelData->nStates, idaData->yp);
836 }
837
838 ✗ if (solverInfo->didEventStep)
839 {
840 ✗ idaData->setInitialSolution = 0;
841 }
842
843 /* reinit solver */
844 ✗ if (!idaData->setInitialSolution)
845 {
846 /* initialize states and der(states) */
847 ✗ if (idaData->daeMode)
848 {
849 ✗ memcpy(idaData->states, data->localData[0]->realVars, sizeof(double)*data->modelData->nStates);
850 /* and also algebraic vars */
851 ✗ getAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
852 ✗ memcpy(idaData->statesDer, data->localData[0]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
853 }
854
855 /* calculate matrix for residual scaling */
856 ✗ if (omc_flag[FLAG_IDA_SCALING])
857 {
858 ✗ getScalingFactors(data, idaData, NULL);
859
860 /* scale idaData->y and idaData->yp */
861 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "Scale y and yp");
862 ✗ idaScaleData(idaData);
863 ✗ messageClose(OMC_LOG_SOLVER_V);
864 }
865
866 ✗ flag = IDAReInit(idaData->ida_mem,
867 solverInfo->currentTime,
868 idaData->y,
869 idaData->yp);
870 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAReInit");
871
872 /* calculate matrix for residual scaling */
873 ✗ if (omc_flag[FLAG_IDA_SCALING])
874 {
875 /* scale idaData->y and idaData->yp */
876 ✗ idaReScaleData(idaData);
877 }
878
879 ✗ if (idaData->idaSmode)
880 {
881 ✗ for(i=0; i<idaData->Ns; ++i)
882 {
883 int j;
884 ✗ for(j=0; j<idaData->N; ++j)
885 {
886 ✗ NV_Ith_S(idaData->yS[i],j) = 0;
887 ✗ NV_Ith_S(idaData->ySp[i],j) = 0;
888 }
889 }
890 ✗ flag = IDASensReInit(idaData->ida_mem, IDA_SIMULTANEOUS, idaData->yS, idaData->ySp);
891 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASensReInit");
892 }
893 ✗ idaData->setInitialSolution = 1;
894 }
895
896 ✗ saveJumpState = threadData->currentErrorStage;
897 ✗ threadData->currentErrorStage = ERROR_INTEGRATOR;
898
899 /* try */
900 #if !defined(OMC_EMCC)
901 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
902 #endif
903
904
905 /* Check that tout is not less than timeValue otherwise the solver
906 * will come in trouble.
907 * If that is the case we skip the current step. */
908 ✗ if (solverInfo->currentStepSize < DASSL_STEP_EPS)
909 {
910 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Desired step to small try next one");
911 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Interpolate linear");
912
913 /* linear extrapolation */
914 ✗ for(i = 0; i < idaData->N; i++)
915 {
916 ✗ NV_Ith_S(idaData->y, i) = NV_Ith_S(idaData->y, i) + NV_Ith_S(idaData->yp, i) * solverInfo->currentStepSize;
917 }
918 ✗ sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
919 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
920 ✗ data->callback->functionODE(data, threadData);
921 ✗ solverInfo->currentTime = sData->timeValue;
922
923 ✗ return 0;
924 }
925
926
927 /* Calculate steps until TOUT is reached */
928 ✗ if (idaData->internalSteps)
929 {
930 /* If internalSteps are selected, let IDA run to stopTime or next sample event */
931 ✗ if (data->simulationInfo->nextSampleEvent < data->simulationInfo->stopTime)
932 {
933 tout = data->simulationInfo->nextSampleEvent;
934 }
935 else
936 {
937 tout = data->simulationInfo->stopTime;
938 }
939 stepsMode = IDA_ONE_STEP;
940 ✗ flag = IDASetStopTime(idaData->ida_mem, tout);
941 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDASetStopTime");
942 }
943 else
944 {
945 ✗ tout = solverInfo->currentTime + solverInfo->currentStepSize;
946 stepsMode = IDA_NORMAL;
947 /* Never step past the next time event */
948 ✗ if (data->simulationInfo->nextSampleEvent < DBL_MAX)
949 {
950 ✗ IDASetStopTime(idaData->ida_mem, fmax(data->simulationInfo->nextSampleEvent, tout));
951 }
952 else
953 {
954 ✗ IDAClearStopTime(idaData->ida_mem);
955 }
956 }
957
958
959 do
960 {
961 ✗ infoStreamPrint(OMC_LOG_SOLVER, 1, "##IDA## new step from %.15g to %.15g", solverInfo->currentTime, tout);
962
963 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
964 /* read input vars */
965 ✗ externalInputUpdate(data);
966 ✗ data->callback->input_function(data, threadData);
967 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
968
969 ✗ if (omc_flag[FLAG_IDA_SCALING])
970 {
971 /* scale idaData->y and idaData->yp */
972 ✗ idaScaleData(idaData);
973 }
974
975 ✗ flag = IDASolve(idaData->ida_mem, tout, &solverInfo->currentTime, idaData->y, idaData->yp, stepsMode);
976
977 ✗ if (omc_flag[FLAG_IDA_SCALING])
978 {
979 /* rescale idaData->y and idaData->yp */
980 ✗ idaReScaleData(idaData);
981 }
982
983 /* set time to current time */
984 ✗ sData->timeValue = solverInfo->currentTime;
985
986 rampRecovery = idaHomotopyRampRecovers(idaData, data, solverInfo, flag);
987
988 /* error handling */
989 ✗ if (IDAflagIsSuccess(flag) && solverInfo->currentTime >= tout)
990 {
991 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## step done to time = %.15g", solverInfo->currentTime);
992 finished = 1 /* TRUE */;
993 }
994 ✗ else if (flag == IDA_ROOT_RETURN)
995 {
996 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## root found at time = %.15g", solverInfo->currentTime);
997 finished = 1 /* TRUE */;
998 }
999 ✗ else if (flag == IDA_SUCCESS || flag == IDA_TSTOP_RETURN)
1000 {
1001 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## continue integration time = %.15g", solverInfo->currentTime);
1002 }
1003 ✗ else if (flag == IDA_TOO_MUCH_WORK)
1004 {
1005 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## has done too much work with small steps at time = %.15g", solverInfo->currentTime);
1006 }
1007 ✗ else if ((flag == IDA_LSETUP_FAIL || rampRecovery) && !restartAfterLSFail)
1008 {
1009 ✗ if (rampRecovery)
1010 {
1011 ✗ idaActivateHomotopyRamp(idaData, data);
1012 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## degenerate DAE operating point at t = %.15g (flag %d); activating homotopy ramp", solverInfo->currentTime, flag);
1013 }
1014 ✗ flag = IDAReInit(idaData->ida_mem,
1015 solverInfo->currentTime,
1016 idaData->y,
1017 idaData->yp);
1018 ✗ restartAfterLSFail = 1;
1019 ✗ warningStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## solver failed, try once again at time = %.15g", solverInfo->currentTime);
1020 }
1021 else
1022 {
1023 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "##IDA## %d error occurred at time = %.15g", flag, solverInfo->currentTime);
1024 finished = 1 /* TRUE */;
1025 retVal = flag;
1026 }
1027
1028 /* closing new step message */
1029 ✗ messageClose(OMC_LOG_SOLVER);
1030
1031 /* emit step, if step mode is selected */
1032 ✗ if (idaData->internalSteps)
1033 {
1034 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## noEquadistant stepsOutputCounter %d by freq %d at time = %.15g", stepsOutputCounter, idaData->stepsFreq, solverInfo->currentTime);
1035 ✗ if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ]){
1036 /* output every n-th time step */
1037 ✗ if (stepsOutputCounter >= idaData->stepsFreq){
1038 ✗ stepsOutputCounter = 1; /* next line set it to one */
1039 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## noEquadistant output %d by freq at time = %.15g", stepsOutputCounter, solverInfo->currentTime);
1040 break;
1041 }
1042 ✗ stepsOutputCounter++;
1043 ✗ } else if (omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]){
1044 /* output when time>=k*timeValue */
1045 ✗ if (solverInfo->currentTime > stepsOutputCounter * idaData->stepsTime){
1046 ✗ stepsOutputCounter++;
1047 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## noEquadistant output %d by time freq at time = %.15g", stepsOutputCounter, solverInfo->currentTime);
1048 break;
1049 }
1050 } else {
1051 break;
1052 }
1053 }
1054
1055 ✗ } while(!finished && !OMC_ERROR_RAISED());
1056 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
1057
1058 #if !defined(OMC_EMCC)
1059 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1060 #endif
1061 ✗ threadData->currentErrorStage = saveJumpState;
1062
1063 /* if a state event occurs than no sample event does need to be activated */
1064 ✗ if (data->simulationInfo->sampleActivated && solverInfo->currentTime < data->simulationInfo->nextSampleEvent)
1065 {
1066 ✗ data->simulationInfo->sampleActivated = 0;
1067 }
1068
1069 /* initialize states and der(states) */
1070 ✗ if (idaData->daeMode)
1071 {
1072 ✗ memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
1073 // and also algebraic vars
1074 ✗ setAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
1075 ✗ memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
1076 ✗ sData->timeValue = solverInfo->currentTime;
1077 }
1078
1079 /* sensitivity mode */
1080 ✗ if (idaData->idaSmode)
1081 {
1082 ✗ flag = IDAGetSens(idaData->ida_mem, &solverInfo->currentTime, idaData->ySResult);
1083 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAGetSens");
1084 }
1085
1086 /* save stats */
1087 /* steps */
1088 ✗ tmp = 0;
1089 ✗ flag = IDAGetNumSteps(idaData->ida_mem, &tmp);
1090 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAGetNumSteps");
1091 ✗ solverInfo->solverStatsTmp.nStepsTaken = tmp;
1092
1093 /* functionODE evaluations */
1094 ✗ tmp = 0;
1095 ✗ flag = IDAGetNumResEvals(idaData->ida_mem, &tmp);
1096 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAGetNumResEvals");
1097 ✗ solverInfo->solverStatsTmp.nCallsODE = tmp;
1098
1099 /* Jacobians evaluations */
1100 ✗ tmp = 0;
1101 ✗ flag = IDAGetNumJacEvals(idaData->ida_mem, &tmp);
1102 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAGetNumJacEvals");
1103 ✗ solverInfo->solverStatsTmp.nCallsJacobian = tmp;
1104
1105 /* local error test failures */
1106 ✗ tmp = 0;
1107 ✗ flag = IDAGetNumErrTestFails(idaData->ida_mem, &tmp);
1108 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAGetNumErrTestFails");
1109 ✗ solverInfo->solverStatsTmp.nErrorTestFailures = tmp;
1110
1111 /* local error test failures */
1112 ✗ tmp = 0;
1113 ✗ flag = IDAGetNumNonlinSolvConvFails(idaData->ida_mem, &tmp);
1114 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_IDA_FLAG, "IDAGetNumNonlinSolvConvFails");
1115 ✗ solverInfo->solverStatsTmp.nConvergenceTestFailures = tmp;
1116
1117 /* get more statistics */
1118 ✗ if (omc_useStream[OMC_LOG_SOLVER_V])
1119 {
1120 long int tmp1,tmp2;
1121 double dtmp;
1122
1123 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "### IDAStats ###");
1124 /* nonlinear stats */
1125 ✗ tmp1 = tmp2 = 0;
1126 ✗ flag = IDAGetNonlinSolvStats(idaData->ida_mem, &tmp1, &tmp2);
1127 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, " ## Cumulative number of nonlinear iterations performed: %ld", tmp1);
1128 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, " ## Cumulative number of nonlinear convergence failures that have occurred: %ld", tmp2);
1129
1130 /* others */
1131 ✗ flag = IDAGetTolScaleFactor(idaData->ida_mem, &dtmp);
1132 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, " ## Suggested scaling factor for user tolerances: %g", dtmp);
1133
1134 ✗ flag = IDAGetNumLinSolvSetups(idaData->ida_mem, &tmp1);
1135 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, " ## Number of calls made to the linear solver setup function: %ld", tmp1);
1136
1137 ✗ messageClose(OMC_LOG_SOLVER_V);
1138 }
1139
1140 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "##IDA## Finished Integrator step.");
1141 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1142
1143 ✗ return retVal;
1144 }
1145
1146 /**
1147 * @brief Compute residual F(t, y, y').
1148 *
1149 * This function has to be of type IDAResFn.
1150 * See section 4.6.1 Residual function of SUNDIALS v5.4.0 IDA documentation.
1151 *
1152 * @param time Value of the independent variable (time).
1153 * @param yy Vector of state variables y(t).
1154 * @param yp Vector of derivatives y'(t).
1155 * @param rr Output residual vector F(t, y, y').
1156 * @param user_data Pointer to user data of type IDA_SOLVER*, set with IDASetUserDat.
1157 * @return int Return 0 on success, positive value on recoverable error and negative value otherwise.
1158 */
1159 ✗ static int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data)
1160 {
1161 IDA_SOLVER* idaData = (IDA_SOLVER*) user_data;
1162 ✗ DATA* data = idaData->userData->data;
1163 ✗ threadData_t* threadData = idaData->userData->threadData;
1164
1165 double timeBackup;
1166 long int i;
1167 int saveJumpState;
1168 ✗ int success = 0, retVal = 0;
1169 ✗ double *states = N_VGetArrayPointer_Serial(yy);
1170 ✗ double *statesDer = N_VGetArrayPointer_Serial(yp);
1171 ✗ double *delta = N_VGetArrayPointer_Serial(rr);
1172
1173 /* DAE-mode homotopy ramp for a degenerate initial operating point: a
1174 homotopy()-regularized characteristic (e.g. a pump at zero flow/speed) makes
1175 the actual (lambda=1) DAE Jacobian singular near t=0, so IDA's LU
1176 factorization fails. The simplified (lambda<1) branch is non-singular, so
1177 once such a failure is detected (homotopyRampActive set in ida_solver_step)
1178 we ramp lambda 0->1 smoothly over [startTime, startTime+t_ramp]. This is
1179 activated only on failure, so models that integrate normally are unaffected.
1180 t_ramp defaults to 10% of the simulation interval (tunable via env var). */
1181 ✗ if (idaData->daeMode && idaData->homotopyRampActive && data->simulationInfo->homotopySteps > 0) {
1182 ✗ double tRamp = idaData->homotopyTramp;
1183 ✗ data->simulationInfo->lambda = (tRamp > 0.0 && time < data->simulationInfo->startTime + tRamp)
1184 ✗ ? ((time - data->simulationInfo->startTime) / tRamp) : 1.0;
1185 }
1186
1187 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "### eval residualFunctionIDA ###");
1188 /* rescale idaData->y and idaData->yp */
1189 ✗ if ((omc_flag[FLAG_IDA_SCALING] && idaData->useScaling))
1190 {
1191 ✗ idaReScaleData(idaData);
1192 }
1193
1194 ✗ if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
1195 {
1196 ✗ setContext(data, time, CONTEXT_ODE);
1197 }
1198 ✗ data->localData[0]->timeValue = time;
1199
1200 ✗ saveJumpState = threadData->currentErrorStage;
1201 ✗ threadData->currentErrorStage = ERROR_INTEGRATOR;
1202
1203 /* try */
1204 #if !defined(OMC_EMCC)
1205 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1206 #endif
1207
1208 /* if sensitivity mode update also bound parameters*/
1209 ✗ if (idaData->idaSmode)
1210 {
1211 ✗ data->callback->updateBoundParameters(data, threadData);
1212 }
1213 /* if daeMode update also all dynamic algebraic equations */
1214 ✗ if (idaData->daeMode)
1215 {
1216 /* set state, state derivative and dynamic algebraic
1217 * variables for evaluateDAEResiduals evaluation
1218 */
1219 ✗ memcpy(data->localData[0]->realVars, states, sizeof(double)*data->modelData->nStates);
1220 ✗ memcpy(data->localData[0]->realVars + data->modelData->nStates, statesDer, sizeof(double)*data->modelData->nStates);
1221 ✗ setAlgebraicDAEVars(data, states + data->modelData->nStates);
1222 }
1223
1224 /* debug */
1225 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER_V)){
1226 ✗ printCurrentStatesVector(OMC_LOG_SOLVER_V, data->localData[0]->realVars, data, time);
1227 ✗ printVector(OMC_LOG_SOLVER_V, "yprime", data->localData[0]->realVars + data->modelData->nStates, data->modelData->nStates, time);
1228 ✗ if (idaData->daeMode)
1229 {
1230 ✗ printVector(OMC_LOG_SOLVER_V, "yalg", states + data->modelData->nStates, data->simulationInfo->daeModeData->nAlgebraicDAEVars, time);
1231 }
1232 }
1233
1234 /* read input vars */
1235 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1236 ✗ externalInputUpdate(data);
1237 ✗ data->callback->input_function(data, threadData);
1238 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
1239
1240 ✗ if (idaData->daeMode)
1241 {
1242 /* eval residual vars */
1243 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1244 ✗ data->simulationInfo->daeModeData->evaluateDAEResiduals(data, threadData, EVAL_DYNAMIC);
1245 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
1246 /* get residual variables */
1247 ✗ for(i=0; i < idaData->N; i++)
1248 {
1249 ✗ NV_Ith_S(rr, i) = data->simulationInfo->daeModeData->residualVars[i];
1250 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, "%ld. residual = %e", i, NV_Ith_S(rr, i));
1251 }
1252 }
1253 else
1254 {
1255 /* In sensitivity mode IDAS evaluates the residual on perturbed copies of the
1256 * state vector (yy != idaData->y, which otherwise shares storage with realVars).
1257 * Copy the perturbed states into realVars so functionODE sees them; otherwise the
1258 * dF/dy term is lost from the sensitivity difference quotient. The base states are
1259 * backed up and restored afterwards so idaData->y is not corrupted for the
1260 * remaining sensitivity equations of the same difference-quotient sweep. */
1261 ✗ const int perturbedStates = (idaData->idaSmode && data->localData[0]->realVars != states);
1262 if (perturbedStates)
1263 {
1264 ✗ memcpy(idaData->ysave, data->localData[0]->realVars, sizeof(double)*data->modelData->nStates);
1265 ✗ memcpy(data->localData[0]->realVars, states, sizeof(double)*data->modelData->nStates);
1266 }
1267 /* eval function ODE */
1268 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1269 ✗ data->callback->functionODE(data, threadData);
1270 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
1271 ✗ for(i=0; i < idaData->N; i++)
1272 {
1273 ✗ NV_Ith_S(rr, i) = data->localData[0]->realVars[data->modelData->nStates + i] - NV_Ith_S(yp, i);
1274 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 0, "%ld. residual = %e", i, NV_Ith_S(rr, i));
1275 }
1276 ✗ if (perturbedStates)
1277 ✗ memcpy(data->localData[0]->realVars, idaData->ysave, sizeof(double)*data->modelData->nStates);
1278 }
1279
1280 /* scale ressidual rr */
1281 ✗ if ((omc_flag[FLAG_IDA_SCALING] && idaData->useScaling))
1282 {
1283 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "scale residuals");
1284 ✗ idaScaleVector(rr, idaData->resScale, idaData->N);
1285 ✗ messageClose(OMC_LOG_SOLVER_V);
1286 ✗ idaScaleData(idaData);
1287 }
1288
1289 ✗ printVector(OMC_LOG_DASSL_STATES, "delta", delta, idaData->N, time);
1290 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { success = 1; }
1291 #if !defined(OMC_EMCC)
1292 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1293 #endif
1294
1295 ✗ if (!success) {
1296 retVal = 1; /* Recoverable error, reduce step size and retry */
1297 }
1298
1299 ✗ threadData->currentErrorStage = saveJumpState;
1300
1301 ✗ if (data->simulationInfo->currentContext == CONTEXT_ODE){
1302 ✗ unsetContext(data);
1303 }
1304 ✗ messageClose(OMC_LOG_SOLVER_V);
1305 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1306
1307 ✗ return retVal;
1308 }
1309
1310 /**
1311 * @brief Evaluate zero crossings for IDA root finding.
1312 *
1313 * Set by IDARootInit, has be of type IDARootFn.
1314 * See section 4.5.6 Rootfinding initialization function of SUNDIALS v5.4.0 IDA documentation.
1315 *
1316 * @param time Independent variable (time).
1317 * @param yy Vector of state variables y.
1318 * @param yp Vector of state derivatives y'.
1319 * @param gout Output array: ZeroCrossings g(t, y, y').
1320 * @param user_data Pointer to user data of type `IDA_SOLVER*`.
1321 * @return int Return 0 on success and otherwise error.
1322 */
1323 ✗ static int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void* user_data)
1324 {
1325 IDA_SOLVER* idaData = (IDA_SOLVER*) user_data;
1326 ✗ DATA* data = idaData->userData->data;
1327 ✗ threadData_t* threadData = idaData->userData->threadData;
1328 ✗ double *states = N_VGetArrayPointer_Serial(yy);
1329 ✗ double *statesDer = N_VGetArrayPointer_Serial(yp);
1330
1331 int saveJumpState;
1332
1333 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "### eval rootsFunctionIDA ###");
1334
1335 ✗ if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
1336 {
1337 ✗ setContext(data, time, CONTEXT_EVENTS);
1338 }
1339
1340 /* re-scale idaData->y and idaData->yp to evaluate the equations */
1341 ✗ if (omc_flag[FLAG_IDA_SCALING])
1342 {
1343 ✗ idaReScaleData(idaData);
1344 }
1345
1346 ✗ saveJumpState = threadData->currentErrorStage;
1347 ✗ threadData->currentErrorStage = ERROR_EVENTSEARCH;
1348
1349 ✗ if (idaData->daeMode)
1350 {
1351 ✗ memcpy(data->localData[0]->realVars, states, sizeof(double)*data->modelData->nStates);
1352 ✗ setAlgebraicDAEVars(data, states + data->modelData->nStates);
1353 ✗ memcpy(data->localData[0]->realVars + data->modelData->nStates, statesDer, sizeof(double)*data->modelData->nStates);
1354 }
1355
1356 ✗ data->localData[0]->timeValue = time;
1357
1358 /* Exlude zero-crossings eval from sim timer */
1359 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1360
1361 /* Update inputs and evaluate zero crossings */
1362 ✗ externalInputUpdate(data);
1363 ✗ data->callback->input_function(data, threadData);
1364 ✗ if (idaData->daeMode){
1365 ✗ data->simulationInfo->daeModeData->evaluateDAEResiduals(data, threadData, EVAL_ZEROCROSS);
1366 }
1367 else
1368 {
1369 ✗ data->callback->function_ZeroCrossingsEquations(data, threadData);
1370 }
1371 ✗ data->callback->function_ZeroCrossings(data, threadData, gout);
1372
1373 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
1374 ✗ threadData->currentErrorStage = saveJumpState;
1375
1376 /* scale data again */
1377 ✗ if (omc_flag[FLAG_IDA_SCALING])
1378 {
1379 ✗ idaScaleData(idaData);
1380 }
1381
1382 ✗ if (data->simulationInfo->currentContext == CONTEXT_EVENTS){
1383 ✗ unsetContext(data);
1384 }
1385 ✗ messageClose(OMC_LOG_SOLVER_V);
1386 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER); // TODO: Why do we have two rt_tick calls? Keep only this one?
1387
1388 ✗ return 0;
1389 }
1390
1391 /**
1392 * @brief Compute colored numeric Jacobian.
1393 *
1394 * Calculate Jacobian matrix using finite differences method
1395 * with coloring from sparsity pattern.
1396 *
1397 * @param currentTime Independent variable (time).
1398 * @param cj Scalar in the system Jacobian, proportional to the inverse of the step size.
1399 * @param yy Vector of state variables y.
1400 * @param yp Vector of derivatives y'.
1401 * @param rr Vector of residual vector F(y,y').
1402 * @param Jac Output Jacobian: J = (∂F)/(∂y).
1403 * @param idaData Pointer to IDA user data.
1404 * @return int Return 0 on success, positive value on recoverable error and negative value otherwise.
1405 */
1406 ✗ static int jacColoredNumericalDense(double currentTime, double cj, N_Vector yy, N_Vector yp,
1407 N_Vector rr, SUNMatrix Jac, IDA_SOLVER *idaData)
1408 {
1409 ✗ DATA* data = idaData->userData->data;
1410 ✗ void* ida_mem = idaData->ida_mem;
1411 ✗ const int index = data->callback->INDEX_JAC_A;
1412
1413 /* prepare variables */
1414 ✗ double *states = N_VGetArrayPointer_Serial(yy);
1415 ✗ double *yprime = N_VGetArrayPointer_Serial(yp);
1416 ✗ double *delta = N_VGetArrayPointer_Serial(rr);
1417 ✗ double *newdelta = N_VGetArrayPointer_Serial(idaData->newdelta);
1418
1419 ✗ double *delta_hh = idaData->delta_hh;
1420 ✗ double *ysave = idaData->ysave;
1421 ✗ double *ypsave = idaData->ypsave;
1422
1423 double delta_hhh;
1424 ✗ double *abstol = N_VGetArrayPointer_Serial(idaData->absoluteTolerance);
1425 ✗ double rtol = data->simulationInfo->tolerance;
1426 long int i,j,l,ii;
1427
1428 double currentStep;
1429
1430 /* set values */
1431 ✗ IDAGetCurrentStep(ida_mem, &currentStep);
1432
1433 SPARSE_PATTERN* sparsePattern;
1434
1435 /* set sparse pattern */
1436 ✗ if (idaData->daeMode)
1437 {
1438 ✗ sparsePattern = data->simulationInfo->daeModeData->sparsePattern;
1439 }
1440 else
1441 {
1442 ✗ sparsePattern = data->simulationInfo->analyticJacobians[index].sparsePattern;
1443 }
1444
1445 ✗ setContext(data, currentTime, CONTEXT_JACOBIAN);
1446
1447 ✗ for(i = 0; i < sparsePattern->maxColors; i++)
1448 {
1449 ✗ for(ii=0; ii < idaData->N; ii++)
1450 {
1451 ✗ if(sparsePattern->colorCols[ii]-1 == i)
1452 {
1453 ✗ delta_hhh = currentStep * yprime[ii];
1454 ✗ delta_hh[ii] = numericalJacobianStep(states[ii], delta_hhh, rtol*fabs(states[ii]) + abstol[ii],
1455 ✗ idaData->jacNominalFactor * idaData->nominal[ii]);
1456 ✗ delta_hh[ii] = (delta_hhh >= 0 ? delta_hh[ii] : -delta_hh[ii]);
1457 ✗ delta_hh[ii] = (states[ii] + delta_hh[ii]) - states[ii]; // Due to floating-point arithmetic rounding errors can result in: delta_hh[ii] != (states[ii] + delta_hh[ii]) - states[ii]
1458 ✗ ysave[ii] = states[ii];
1459 ✗ states[ii] += delta_hh[ii];
1460
1461 ✗ if (idaData->daeMode){
1462 ✗ ypsave[ii] = yprime[ii];
1463 ✗ yprime[ii] += cj * delta_hh[ii];
1464 }
1465
1466 ✗ delta_hh[ii] = 1. / delta_hh[ii];
1467 }
1468 }
1469
1470 ✗ idaData->residualFunction(currentTime, yy, yp, idaData->newdelta, (void*) idaData); /* Points to residualFunctionIDA */
1471
1472 ✗ increaseJacContext(data);
1473
1474 ✗ for(ii = 0; ii < idaData->N; ii++)
1475 {
1476 ✗ if(sparsePattern->colorCols[ii]-1 == i)
1477 {
1478 ✗ j = sparsePattern->leadindex[ii];
1479 ✗ while(j < sparsePattern->leadindex[ii+1])
1480 {
1481 ✗ l = sparsePattern->index[j];
1482 ✗ SM_ELEMENT_D(Jac, l, ii) = (newdelta[l] - delta[l]) * delta_hh[ii];
1483 ✗ j++;
1484 };
1485 ✗ states[ii] = ysave[ii];
1486 ✗ if (idaData->daeMode)
1487 {
1488 ✗ yprime[ii] = ypsave[ii];
1489 }
1490 }
1491 }
1492 }
1493 ✗ unsetContext(data);
1494
1495 ✗ return 0;
1496 }
1497
1498 /**
1499 * @brief Compute colored symbolic Jacobian.
1500 *
1501 * Calculate the Jacobian matrix with the shared symbolic Jacobian evaluation, which
1502 * transparently handles forward (coloredSymbolical), adjoint (coloredSymbolicalAdjoint)
1503 * and bidirectional (bicoloredSymbolical) evaluation.
1504 *
1505 * @param currentTime Independent variable (time).
1506 * @param cj Scalar in the system Jacobian, proportional to the inverse of the step size.
1507 * @param yy Vector of state variables y.
1508 * @param yp Vector of derivatives y'.
1509 * @param rr Vector of residual vector F(y,y').
1510 * @param Jac Output Jacobian: J = (∂F)/(∂y) + cj * (∂F)/(∂y').
1511 * @param idaData Pointer to IDA user data.
1512 * @return int Return 0 on success, positive value on recoverable error and negative value otherwise.
1513 */
1514 ✗ static int jacColoredSymbolicalDense(double currentTime, double cj, N_Vector yy,
1515 N_Vector yp, N_Vector rr, SUNMatrix Jac,
1516 IDA_SOLVER *idaData)
1517 {
1518 ✗ DATA* data = idaData->userData->data;
1519 ✗ threadData_t* threadData = idaData->userData->threadData;
1520 ✗ JACOBIAN* jac = getSymbolicOdeJacobian(data);
1521 ✗ jac->dae_cj = cj;
1522
1523 ✗ setContext(data, currentTime, CONTEXT_SYM_JACOBIAN); /* Reuse jacobian matrix in KLU solver */
1524 ✗ evalJacobian(data, threadData, jac, NULL, SM_DATA_D(Jac), TRUE);
1525
1526 ✗ unsetContext(data);
1527
1528 ✗ return 0;
1529 }
1530
1531 /**
1532 * @brief Compute colored Jacobian matrix of ODE/DAE system.
1533 *
1534 * Available methods:
1535 * - Colored Numeric Jacobian --> jacColoredNumericalDense
1536 * - Colored Symbolic Jacobian --> jacColoredSymbolicalDense
1537 *
1538 * See Section 4.6.5 in IDA documentation of SUNDIALS v5.4.0 for more details.
1539 *
1540 * @param tt Independent variable (time).
1541 * @param cj Scalar in the system Jacobian, proportional to the inverse of the step size.
1542 * @param yy Vector of state variables y.
1543 * @param yp Vector of state derivatives y'.
1544 * @param rr Vector of residual vector F(y,y').
1545 * @param Jac Output Jacobian: J = (∂F)/(∂y) + cj * (∂F)/(∂y').
1546 * @param user_data Pointer to user data of type `IDA_SOLVER*`.
1547 * @param tmp1 Work array that can be used by, currently unused.
1548 * @param tmp2 Work array that can be used by, currently unused.
1549 * @param tmp3 Work array that can be used by, currently unused.
1550 * @return int Return 0 on success, positive value on recoverable error and negative value otherwise.
1551 */
1552 ✗ static int callDenseJacobian(sunrealtype tt, sunrealtype cj, N_Vector yy,
1553 N_Vector yp, N_Vector rr, SUNMatrix Jac,
1554 void *user_data, N_Vector tmp1, N_Vector tmp2,
1555 N_Vector tmp3) {
1556 IDA_SOLVER* idaData = (IDA_SOLVER*) user_data;
1557 ✗ threadData_t* threadData = idaData->userData->threadData;
1558 int retVal;
1559
1560 /* profiling */
1561 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1562 ✗ rt_tick(SIM_TIMER_JACOBIAN);
1563
1564 ✗ if (idaData->jacobianMethod == COLOREDNUMJAC || idaData->jacobianMethod == NUMJAC)
1565 {
1566 ✗ retVal = jacColoredNumericalDense(tt, cj, yy, yp, rr, Jac, idaData);
1567 }
1568 ✗ else if (idaData->jacobianMethod == COLOREDSYMJAC || idaData->jacobianMethod == SYMJAC
1569 ✗ || idaData->jacobianMethod == COLOREDSYMJACADJ || idaData->jacobianMethod == BICOLOREDSYMJAC)
1570 {
1571 ✗ retVal = jacColoredSymbolicalDense(tt, cj, yy, yp, rr, Jac, idaData);
1572 }
1573 else
1574 {
1575 ✗ throwStreamPrint(threadData, "##IDA## Something went wrong while obtain Jacobian matrix.");
1576 }
1577
1578 /* debug */
1579 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC)){
1580 ✗ _omc_matrix* dumpJac = _omc_createMatrix(idaData->N, idaData->N, SM_DATA_D(Jac));
1581 ✗ _omc_printMatrix(dumpJac, "IDA-Solver: Matrix A", OMC_LOG_JAC);
1582 ✗ _omc_destroyMatrix(dumpJac);
1583 }
1584
1585 /* add cj to diagonal elements and store in Jac */
1586 ✗ if (!idaData->daeMode)
1587 {
1588 ✗ for(int i = 0; i < SM_COLUMNS_D(Jac); i++)
1589 {
1590 ✗ SM_ELEMENT_D(Jac, i, i) -= (double) cj;
1591 }
1592 }
1593
1594 /* profiling */
1595 ✗ rt_accumulate(SIM_TIMER_JACOBIAN);
1596 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
1597
1598 ✗ return retVal;
1599 }
1600
1601 /*
1602 * function calculates a jacobian matrix by
1603 * numerical method finite differences with coloring
1604 * into a sparse SlsMat matrix
1605 */
1606 ✗ static int jacoColoredNumericalSparse(double currentTime, N_Vector yy,
1607 N_Vector yp, N_Vector rr, SUNMatrix Jac,
1608 double cj, void *userData) {
1609 IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
1610 ✗ DATA* data = (DATA*)(((IDA_USERDATA*)idaData->userData)->data);
1611 ✗ void* ida_mem = idaData->ida_mem;
1612 ✗ const int index = data->callback->INDEX_JAC_A;
1613
1614 /* prepare variables */
1615 ✗ double *states = N_VGetArrayPointer_Serial(yy);
1616 ✗ double *yprime = N_VGetArrayPointer_Serial(yp);
1617 ✗ double *delta = N_VGetArrayPointer_Serial(rr);
1618 ✗ double *newdelta = N_VGetArrayPointer_Serial(idaData->newdelta);
1619
1620 SPARSE_PATTERN* sparsePattern;
1621
1622 ✗ double *ysave = idaData->ysave;
1623 ✗ double *ypsave = idaData->ypsave;
1624
1625 ✗ double *delta_hh = idaData->delta_hh;
1626 double delta_hhh;
1627 ✗ double *abstol = N_VGetArrayPointer_Serial(idaData->absoluteTolerance);
1628 ✗ double rtol = data->simulationInfo->tolerance;
1629
1630 long int i,j,ii;
1631 int nth = 0;
1632 ✗ int disBackup = idaData->useScaling;
1633
1634 double currentStep;
1635
1636 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "### eval jacobianSparseNumIDA ###");
1637 /* set values */
1638 ✗ IDAGetCurrentStep(ida_mem, &currentStep);
1639
1640 /* set sparse pattern */
1641 ✗ if (idaData->daeMode)
1642 {
1643 ✗ sparsePattern = data->simulationInfo->daeModeData->sparsePattern;
1644 }
1645 else
1646 {
1647 ✗ sparsePattern = data->simulationInfo->analyticJacobians[index].sparsePattern;
1648 }
1649
1650 /* Reset Jacobian matrix */
1651 ✗ SUNMatZero(Jac);
1652
1653 ✗ setContext(data, currentTime, CONTEXT_JACOBIAN);
1654
1655 /* rescale idaData->y and idaData->yp
1656 * the evaluation of the residual function
1657 * needs to be performed on unscaled values
1658 */
1659 ✗ if (omc_flag[FLAG_IDA_SCALING] && idaData->useScaling)
1660 {
1661 ✗ idaReScaleVector(rr, idaData->resScale, idaData->N);
1662 ✗ idaReScaleData(idaData);
1663 }
1664
1665 ✗ for(i = 0; i < sparsePattern->maxColors; i++)
1666 {
1667 ✗ for(ii=0; ii < idaData->N; ii++)
1668 {
1669 ✗ if(sparsePattern->colorCols[ii]-1 == i)
1670 {
1671 ✗ delta_hhh = currentStep * yprime[ii];
1672 ✗ delta_hh[ii] = numericalJacobianStep(states[ii], delta_hhh, rtol*fabs(states[ii]) + abstol[ii],
1673 ✗ idaData->jacNominalFactor * idaData->nominal[ii]);
1674 ✗ delta_hh[ii] = (delta_hhh >= 0 ? delta_hh[ii] : -delta_hh[ii]);
1675 ✗ delta_hh[ii] = (states[ii] + delta_hh[ii]) - states[ii]; // Due to floating-point arithmetic rounding errors can result in: delta_hh[ii] != (states[ii] + delta_hh[ii]) - states[ii]
1676 ✗ ysave[ii] = states[ii];
1677 ✗ states[ii] += delta_hh[ii];
1678
1679 ✗ if (idaData->daeMode){
1680 ✗ ypsave[ii] = yprime[ii];
1681 ✗ yprime[ii] += cj * delta_hh[ii];
1682 }
1683
1684 ✗ delta_hh[ii] = 1. / delta_hh[ii];
1685 }
1686 }
1687 ✗ idaData->useScaling = FALSE;
1688 ✗ idaData->residualFunction(currentTime, yy, yp, idaData->newdelta, userData); /* Points to residualFunctionIDA */
1689 ✗ idaData->useScaling = disBackup;
1690
1691 ✗ increaseJacContext(data);
1692
1693 ✗ for(ii = 0; ii < idaData->N; ii++)
1694 {
1695 ✗ if(sparsePattern->colorCols[ii]-1 == i)
1696 {
1697 ✗ nth = sparsePattern->leadindex[ii];
1698 ✗ while(nth < sparsePattern->leadindex[ii+1])
1699 {
1700 ✗ j = sparsePattern->index[nth];
1701 /* use row scaling for jacobian elements */
1702 ✗ if (!idaData->useScaling || !omc_flag[FLAG_IDA_SCALING]){
1703 ✗ setJacElementSundialsSparse(j, ii, nth, (newdelta[j] - delta[j]) * delta_hh[ii], Jac, SM_CONTENT_S(Jac)->M);
1704 } else {
1705 ✗ setJacElementSundialsSparse(j, ii, nth, ((newdelta[j] - delta[j]) * delta_hh[ii]) / idaData->resScale[j] * idaData->yScale[ii], Jac, SM_CONTENT_S(Jac)->M);
1706 }
1707 ✗ nth++;
1708 };
1709 ✗ states[ii] = ysave[ii];
1710 ✗ if (idaData->daeMode)
1711 {
1712 ✗ yprime[ii] = ypsave[ii];
1713 }
1714 }
1715 }
1716 }
1717 ✗ setSundialsSparseColPtrs(sparsePattern, Jac);
1718
1719 /* scale idaData->y and idaData->yp again */
1720 ✗ if ((omc_flag[FLAG_IDA_SCALING] && idaData->useScaling))
1721 {
1722 ✗ idaScaleVector(rr, idaData->resScale, idaData->N);
1723 ✗ idaScaleData(idaData);
1724 }
1725
1726 ✗ unsetContext(data);
1727 ✗ messageClose(OMC_LOG_SOLVER_V);
1728
1729 ✗ return 0;
1730 }
1731
1732 /*
1733 * This function calculates the jacobian matrix symbolically while exploiting coloring.
1734 * ToDo: backend: generate seeds for der(x)
1735 here: always set der(x) seeds to cj when setting seed for x
1736 */
1737 ✗ int jacColoredSymbolicalSparse(double currentTime, N_Vector yy, N_Vector yp,
1738 N_Vector rr, SUNMatrix Jac, double cj,
1739 void *userData)
1740 {
1741 IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
1742 ✗ DATA* data = (DATA*)(((IDA_USERDATA*)idaData->userData)->data);
1743 ✗ threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)idaData->userData)->threadData);
1744 ✗ JACOBIAN* jac = getSymbolicOdeJacobian(data);
1745 ✗ jac->dae_cj = cj;
1746
1747 /* Reset Jacobian matrix */
1748 ✗ SUNMatZero(Jac);
1749
1750 ✗ setContext(data, currentTime, CONTEXT_SYM_JACOBIAN); /* Reuse jacobian matrix in KLU solver */
1751
1752 ✗ setSundialsSparsePattern(jac, Jac);
1753 ✗ evalJacobian(data, threadData, jac, NULL, SM_DATA_S(Jac), FALSE);
1754
1755 ✗ unsetContext(data);
1756
1757 ✗ return 0;
1758 }
1759
1760 /*
1761 * Wrapper function to call numerical or symbolical jacobian matrix
1762 */
1763 ✗ static int callSparseJacobian(double currentTime, double cj,
1764 N_Vector yy, N_Vector yp, N_Vector rr,
1765 SUNMatrix Jac, void *user_data,
1766 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
1767 {
1768 IDA_SOLVER* idaData = (IDA_SOLVER*)user_data;
1769 DATA* data = (DATA*)(((IDA_USERDATA*)idaData->userData)->data);
1770 threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)user_data)->userData)->threadData);
1771 int i;
1772 SUNErrCode flag;
1773
1774 /* profiling */
1775 ✗ if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER);
1776 ✗ rt_tick(SIM_TIMER_JACOBIAN);
1777
1778 ✗ if (idaData->jacobianMethod == COLOREDSYMJAC || idaData->jacobianMethod == SYMJAC
1779 ✗ || idaData->jacobianMethod == COLOREDSYMJACADJ || idaData->jacobianMethod == BICOLOREDSYMJAC)
1780 {
1781 ✗ jacColoredSymbolicalSparse(currentTime, yy, yp, rr, Jac, cj, user_data);
1782 }
1783 ✗ else if (idaData->jacobianMethod == COLOREDNUMJAC || idaData->jacobianMethod == NUMJAC)
1784 {
1785 ✗ jacoColoredNumericalSparse(currentTime, yy, yp, rr, Jac, cj, user_data);
1786 }
1787
1788 /* debug */
1789 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_JAC)) {
1790 ✗ infoStreamPrint(OMC_LOG_JAC, 0, "##IDA## Sparse Matrix A.");
1791 ✗ SUNSparseMatrix_Print(Jac, stdout);
1792 }
1793 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_DEBUG)) {
1794 ✗ sundialsPrintSparseMatrix(Jac, "A", OMC_LOG_JAC);
1795 }
1796
1797 /* add cj to diagonal elements and store in Jac */
1798 ✗ if (!idaData->daeMode) {
1799 ✗ flag = _omc_SUNMatScaleIAdd_Sparse(-cj, Jac);
1800 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_MATRIX_FLAG, "_omc_SUNMatScaleIAdd_Sparse");
1801 }
1802
1803 /* profiling */
1804 ✗ rt_accumulate(SIM_TIMER_JACOBIAN);
1805 ✗ if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER);
1806
1807 ✗ return 0;
1808 }
1809
1810
1811 /* TODO: Unify with nlsKinsolFScaling from kinsolSolver.c? */
1812 ✗ static int getScalingFactors(DATA* data, IDA_SOLVER* idaData, SUNMatrix inScaleMatrix)
1813 {
1814 int i;
1815
1816 ✗ N_Vector tmp1 = N_VNew_Serial(idaData->N, idaData->sunctx);
1817 ✗ N_Vector tmp2 = N_VNew_Serial(idaData->N, idaData->sunctx);
1818 ✗ N_Vector tmp3 = N_VNew_Serial(idaData->N, idaData->sunctx);
1819
1820 ✗ N_Vector rres = N_VNew_Serial(idaData->N, idaData->sunctx);
1821
1822 SUNMatrix denseMatrix;
1823
1824 ✗ if (inScaleMatrix == NULL)
1825 {
1826 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "##IDA## get new scaling matrix.");
1827
1828 /* use y scale to scale jacobian, but y and yp are not scaled */
1829 ✗ idaData->useScaling = FALSE;
1830
1831 /* eval residual function first */
1832 ✗ idaData->residualFunction(data->localData[0]->timeValue, idaData->y, idaData->yp, rres, (void*) idaData); /* Points to residualFunctionIDA */
1833
1834 /* choose the jacobian sparse vs. dense */
1835 ✗ if (idaData->linearSolverMethod == IDA_LS_KLU)
1836 {
1837 ✗ if (idaData->NNZ < 0)
1838 {
1839 ✗ throwStreamPrint(NULL, "##IDA## idaData->NNZ not set.");
1840 }
1841 ✗ if (idaData->scaleMatrix == NULL) {
1842 ✗ idaData->scaleMatrix = SUNSparseMatrix(idaData->N, idaData->N, idaData->NNZ + idaData->N, SUN_CSC_MAT, idaData->sunctx);
1843 }
1844 ✗ callSparseJacobian(data->localData[0]->timeValue, 1.0, idaData->y, idaData->yp, rres,
1845 idaData->scaleMatrix, idaData, tmp1, tmp2, tmp3);
1846 }
1847 else
1848 {
1849 ✗ denseMatrix = SUNDenseMatrix(idaData->N, idaData->N, idaData->sunctx);
1850 ✗ callDenseJacobian(data->localData[0]->timeValue, 1.0, idaData->y,
1851 idaData->yp, rres, denseMatrix, idaData, tmp1, tmp2,
1852 tmp3);
1853 ✗ SUNMatDestroy(idaData->scaleMatrix);
1854 ✗ idaData->scaleMatrix = SUNSparseFromDenseMatrix(denseMatrix, DBL_MIN, SUN_CSC_MAT);
1855 ✗ if (idaData->scaleMatrix == NULL) {
1856 ✗ errorStreamPrint(
1857 OMC_LOG_STDOUT, 0,
1858 "##IDA## In function SUNSparseFromDenseMatrix: Requirements are "
1859 "violated, or matrix storage request cannot be satisfied.");
1860 }
1861 ✗ SUNMatDestroy(denseMatrix);
1862 }
1863 /* enable scaled jacobian again */
1864 ✗ idaData->useScaling = TRUE;
1865 }
1866 else
1867 {
1868 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "##IDA## use given scaling matrix.");
1869 ✗ idaData->scaleMatrix = inScaleMatrix;
1870 }
1871
1872 /* set resScale factors */
1873 ✗ _omc_fillVector(_omc_createVector(idaData->N,idaData->resScale), MINIMAL_SCALE_FACTOR);
1874 ✗ for (i=0; i<SM_INDEXPTRS_S(idaData->scaleMatrix)[idaData->N]; ++i) {
1875 ✗ if (idaData->resScale[SM_INDEXVALS_S(idaData->scaleMatrix)[i]] < fabs(SM_DATA_S(idaData->scaleMatrix)[i])) {
1876 ✗ idaData->resScale[SM_INDEXVALS_S(idaData->scaleMatrix)[i]] = fabs(SM_DATA_S(idaData->scaleMatrix)[i]);
1877 }
1878 }
1879
1880 ✗ printVector(OMC_LOG_SOLVER_V, "Prime scale factors", idaData->ypScale, idaData->N, 0.0);
1881 ✗ printVector(OMC_LOG_SOLVER_V, "Residual scale factors", idaData->resScale, idaData->N, 0.0);
1882
1883 /* Free memory */
1884 ✗ messageClose(OMC_LOG_SOLVER_V);
1885 ✗ N_VDestroy_Serial(tmp1);
1886 ✗ N_VDestroy_Serial(tmp2);
1887 ✗ N_VDestroy_Serial(tmp3);
1888 ✗ N_VDestroy_Serial(rres);
1889
1890 ✗ return 0;
1891 }
1892
1893 /**
1894 * @brief Scale NVector by factors.
1895 *
1896 * @param vec Vector to scale.
1897 * @param factors Array with scaling factors.
1898 * @param size Length of array factors and vector vec.
1899 */
1900 ✗ static void idaScaleVector(N_Vector vec, double* factors, unsigned int size)
1901 {
1902 int i;
1903 ✗ double *data = N_VGetArrayPointer_Serial(vec);
1904 ✗ printVector(OMC_LOG_SOLVER_V, "un-scaled", data, size, 0.0);
1905 ✗ for(i=0; i < size; ++i)
1906 {
1907 ✗ data[i] = data[i] / factors[i];
1908 }
1909 ✗ printVector(OMC_LOG_SOLVER_V, "scaled", data, size, 0.0);
1910 ✗ }
1911
1912 /**
1913 * @brief Scale state and state derivate vector.
1914 *
1915 * @param idaData Containing state vector y(t) and state derivate vector y'(t)
1916 * as well as scaling arrays yScale and ypScale.
1917 */
1918 ✗ static void idaScaleData(IDA_SOLVER *idaData)
1919 {
1920 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "Scale y");
1921 ✗ idaScaleVector(idaData->y, idaData->yScale, idaData->N);
1922 ✗ messageClose(OMC_LOG_SOLVER_V);
1923 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "Scale yp");
1924 ✗ idaScaleVector(idaData->yp, idaData->ypScale, idaData->N);
1925 ✗ messageClose(OMC_LOG_SOLVER_V);
1926 ✗ }
1927
1928 /**
1929 * @brief Rescale NVector by factors.
1930 *
1931 * @param vec Vector to rescale.
1932 * @param factors Array with scaling factors.
1933 * @param size Length of array factors and vector vec.
1934 */
1935 ✗ static void idaReScaleVector(N_Vector vec, double* factors, unsigned int size)
1936 {
1937 int i;
1938 ✗ double *data = N_VGetArrayPointer_Serial(vec);
1939
1940 ✗ printVector(OMC_LOG_SOLVER_V, "scaled", data, size, 0.0);
1941 ✗ for(i=0; i < size; ++i)
1942 {
1943 ✗ data[i] = data[i] * factors[i];
1944 }
1945 ✗ printVector(OMC_LOG_SOLVER_V, "un-scaled", data, size, 0.0);
1946 ✗ }
1947
1948 /**
1949 * @brief Rescale state and state derivate vector.
1950 * Undo scaling of function idaScaleData.
1951 *
1952 * @param idaData Containing state vector y(t) and state derivate vector y'(t)
1953 * as well as scaling arrays yScale and ypScale.
1954 */
1955 ✗ static void idaReScaleData(IDA_SOLVER *idaData)
1956 {
1957 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "Re-Scale y");
1958 ✗ idaReScaleVector(idaData->y, idaData->yScale, idaData->N);
1959 ✗ messageClose(OMC_LOG_SOLVER_V);
1960 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "Re-Scale yp");
1961 ✗ idaReScaleVector(idaData->yp, idaData->ypScale, idaData->N);
1962 ✗ messageClose(OMC_LOG_SOLVER_V);
1963 ✗ }
1964
1965 #endif /* #ifdef WITH_SUNDIALS */
1966