Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/ida_solver.c
Warning:line 1028, column 9
Value stored to 'flag' is never read

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2016, Open Source Modelica Consortium (OSMC),
5 * c/o Linköpings universitet, Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
11 * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
14 * ACCORDING TO RECIPIENTS CHOICE.
15 *
16 * The OpenModelica software and the OSMC (Open Source Modelica Consortium)
17 * Public License (OSMC-PL) are obtained from OSMC, either from the above
18 * address, from the URLs: http://www.openmodelica.org or
19 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica
20 * distribution. GNU version 3 is obtained from:
21 * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from:
22 * http://www.opensource.org/licenses/BSD-3-Clause.
23 *
24 * This program is distributed WITHOUT ANY WARRANTY; without even the implied
25 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS
26 * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE
27 * CONDITIONS OF OSMC-PL.
28 *
29 */
30
31/*! \file ida_solver.c
32 */
33
34#ifdef USE_PARJAC
35 #include <omp.h>
36 #define GC_THREADS
37 #include <gc/omc_gc.h>
38#endif
39
40#include <string.h>
41#include <setjmp.h>
42
43#include "omc_config.h"
44#include "openmodelica.h"
45#include "openmodelica_func.h"
46#include "simulation_data.h"
47
48#include "util/omc_error.h"
49#include "util/parallel_helper.h"
50#include "gc/omc_gc.h"
51
52#include "simulation/options.h"
53#include "simulation/simulation_runtime.h"
54#include "simulation/results/simulation_result.h"
55#include "simulation/solver/solver_main.h"
56#include "simulation/solver/model_help.h"
57#include "simulation/solver/external_input.h"
58#include "simulation/solver/epsilon.h"
59#include "simulation/solver/omc_math.h"
60#include "simulation/solver/ida_solver.h"
61#include "simulation/solver/dassl.h"
62#include "simulation/solver/dae_mode.h"
63#include "simulation/solver/jacobianSymbolical.h"
64
65#ifdef WITH_SUNDIALS
66
67/* adrpo: on mingw link with static sundials */
68#if defined(__MINGW32__)
69#define LINK_SUNDIALS_STATIC
70#endif
71
72/* readability */
73#define SCALE_MODE0 0
74#define RESCALE_MODE1 1
75
76
77#define MINIMAL_SCALE_FACTOR1e-8 1e-8
78
79
80#include <sundials/sundials_nvector.h>
81#include <nvector/nvector_serial.h>
82#include <idas/idas.h>
83#include <idas/idas_dense.h>
84#include <idas/idas_klu.h>
85#include <idas/idas_spgmr.h>
86#include <idas/idas_spbcgs.h>
87#include <idas/idas_sptfqmr.h>
88
89/* Function prototypes */
90static int callDenseJacobian(long int Neq, double tt, double cj,
91 N_Vector yy, N_Vector yp, N_Vector rr,
92 DlsMat Jac, void *user_data,
93 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
94
95static int callSparseJacobian(realtype currentTime, realtype cj,
96 N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, void *user_data,
97 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
98
99static int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, void* userData);
100int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void* userData);
101
102static int getScalingFactors(DATA* data, IDA_SOLVER *idaData, SlsMat scaleMatrix);
103
104static int idaScaleData(IDA_SOLVER *idaData);
105static int idaReScaleData(IDA_SOLVER *idaData);
106static int idaScaleVector(N_Vector vec, double* factors, unsigned int size);
107static int idaReScaleVector(N_Vector vec, double* factors, unsigned int size);
108
109int ida_event_update(DATA* data, threadData_t *threadData);
110
111/* Static variables */
112static IDA_SOLVER *idaDataGlobal;
113static int initializedSolver = 0;
114
115
116int checkIDAflag(int flag)
117{
118 TRACE_PUSH
119 int retVal;
120 switch(flag)
121 {
122 case IDA_SUCCESS0:
123 case IDA_TSTOP_RETURN1:
124 retVal = 0;
125 break;
126 default:
127 retVal = 1;
128 break;
129 }
130 TRACE_POP
131 return retVal;
132}
133
134void errOutputIDA(int error_code, const char *module, const char *function,
135 char *msg, void *userData)
136{
137 TRACE_PUSH
138 DATA* data = (DATA*)(((IDA_USERDATA*)((IDA_SOLVER*)userData)->simData)->data);
139 infoStreamPrint(LOG_SOLVER, 1, "#### IDA error message #####");
140 infoStreamPrint(LOG_SOLVER, 0, " -> error code %d\n -> module %s\n -> function %s", error_code, module, function);
141 infoStreamPrint(LOG_SOLVER, 0, " Message: %s", msg);
142 messageClose(LOG_SOLVER);
143 TRACE_POP
144}
145
146/* initialize main ida data */
147int ida_solver_initial(DATA* data, threadData_t *threadData,
148 SOLVER_INFO* solverInfo, IDA_SOLVER *idaData)
149{
150 TRACE_PUSH
151
152 int flag;
153 long int i;
154 double* tmp;
155
156 /* default max order */
157 int maxOrder = 5;
158
159 /* sim data */
160 idaData->simData = (IDA_USERDATA*)malloc(sizeof(IDA_USERDATA));
161 idaData->simData->data = data;
162 idaData->simData->threadData = threadData;
163
164 /* initialize constants */
165 idaData->setInitialSolution = 0;
166
167 /* start initialization routines of sundials */
168 idaData->ida_mem = IDACreate();
169 if (idaData->ida_mem == NULL((void*)0)){
170 throwStreamPrint(threadData, "##IDA## Initialization of IDA solver failed!");
171 }
172
173 idaData->daeMode = 0;
174 idaData->residualFunction = residualFunctionIDA;
175 idaData->N = (long int)data->modelData->nStates;
176
177 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
178
179 /* change parameter for DAE mode */
180 if (compiledInDAEMode)
181 {
182 idaData->daeMode = 1;
183 idaData->N = (long int) data->modelData->nStates + data->simulationInfo->daeModeData->nAlgebraicDAEVars;
184 }
185 infoStreamPrint(LOG_SOLVER, 1, "## IDA ## Initializing solver of size %ld %s.", idaData->N, idaData->daeMode?"in DAE mode":"");
186
187 /* initialize states and der(states) */
188 if (idaData->daeMode)
189 {
190 idaData->states = (double*) malloc(idaData->N*sizeof(double));
191 idaData->statesDer = (double*) calloc(idaData->N,sizeof(double));
192
193 memcpy(idaData->states, data->localData[0]->realVars, sizeof(double)*data->modelData->nStates);
194 // and also algebraic vars
195 getAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
196 memcpy(idaData->statesDer, data->localData[0]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
197
198 idaData->y = N_VMake_Serial(idaData->N, idaData->states);
199 idaData->yp = N_VMake_Serial(idaData->N, idaData->statesDer);
200 }
201 else
202 {
203 idaData->y = N_VMake_Serial(idaData->N, data->localData[0]->realVars);
204 idaData->yp = N_VMake_Serial(idaData->N, data->localData[0]->realVars + data->modelData->nStates);
205 }
206
207
208 flag = IDAInit(idaData->ida_mem,
209 idaData->residualFunction,
210 data->simulationInfo->startTime,
211 idaData->y,
212 idaData->yp);
213
214 /* allocate memory for jacobians calculation */
215 idaData->ysave = (double*) malloc(idaData->N*sizeof(double));
216 idaData->ypsave = (double*) malloc(idaData->N*sizeof(double));
217 idaData->delta_hh = (double*) malloc(idaData->N*sizeof(double));
218 idaData->errwgt = N_VNew_Serial(idaData->N);
219 idaData->newdelta = N_VNew_Serial(idaData->N);
220
221 /* allocate memory for initialization process */
222 tmp = (double*) malloc(idaData->N*sizeof(double));
223
224 if (checkIDAflag(flag)){
225 throwStreamPrint(threadData, "##IDA## Something goes wrong while initialize IDA solver!");
226 }
227
228 flag = IDASetUserData(idaData->ida_mem, idaData);
229 if (checkIDAflag(flag)){
230 throwStreamPrint(threadData, "##IDA## Something goes wrong while initialize IDA solver!");
231 }
232
233 flag = IDASetErrHandlerFn(idaData->ida_mem, errOutputIDA, idaData);
234 if (checkIDAflag(flag)){
235 throwStreamPrint(threadData, "##IDA## Something goes wrong while set error handler!");
236 }
237
238
239 /* set nominal values of the states for absolute tolerances */
240 infoStreamPrint(LOG_SOLVER, 1, "The relative tolerance is %g. Following absolute tolerances are used for the states: ", data->simulationInfo->tolerance);
241
242 for(i=0; i < data->modelData->nStates; ++i)
243 {
244 tmp[i] = fmax(fabs(data->modelData->realVarsData[i].attribute.nominal), 1e-32);
245 infoStreamPrint(LOG_SOLVER_V, 0, "%ld. %s -> %g", i+1, data->modelData->realVarsData[i].info.name, tmp[i]);
246 }
247
248 /* daeMode: set nominal values for algebraic variables */
249 if (idaData->daeMode)
250 {
251 getAlgebraicDAEVarNominals(data, tmp + data->modelData->nStates);
252 }
253 /* multiply by tolerance to obtain a relative tolerace */
254 for(i=0; i < idaData->N; ++i)
255 {
256 tmp[i] *= data->simulationInfo->tolerance;
257 }
258 messageClose(LOG_SOLVER);
259 flag = IDASVtolerances(idaData->ida_mem,
260 data->simulationInfo->tolerance,
261 N_VMake_Serial(idaData->N,tmp));
262 if (checkIDAflag(flag)){
263 throwStreamPrint(threadData, "##IDA## Setting tolerances fails while initialize IDA solver!");
264 }
265
266
267 if (omc_flag[FLAG_IDA_SCALING]) // idaNoScaling
268 {
269 /* allocate memory for scaling */
270 idaData->yScale = (double*) malloc(idaData->N*sizeof(double));
271 idaData->ypScale = (double*) malloc(idaData->N*sizeof(double));
272 idaData->resScale = (double*) malloc(idaData->N*sizeof(double));
273
274 /* set yScale from nominal values */
275 for(i=0; i < data->modelData->nStates; ++i)
276 {
277 idaData->yScale[i] = fabs(data->modelData->realVarsData[i].attribute.nominal);
278 idaData->ypScale[i] = 1.0;
279 }
280 /* daeMode: set nominal values for algebraic variables */
281 if (idaData->daeMode)
282 {
283 getAlgebraicDAEVarNominals(data, idaData->yScale + data->modelData->nStates);
284 for(i=data->modelData->nStates; i < idaData->N; ++i)
285 {
286 idaData->ypScale[i] = 1.0;
287 }
288 }
289 infoStreamPrint(LOG_SOLVER_V, 1, "The scale factors for all ida states: ");
290 for(i=0; i < idaData->N; ++i)
291 {
292 infoStreamPrint(LOG_SOLVER_V, 0, "%ld. scaleFactor: %g", i+1, idaData->yScale[i]);
293 }
294 messageClose(LOG_SOLVER_V);
295 }
296 /* initialize */
297 idaData->disableScaling = 0;
298
299 /* set root function */
300 /* if FLAG_NO_ROOTFINDING is set, solve perform without internal root finding*/
301 if(!omc_flag[FLAG_NO_ROOTFINDING])
302 {
303 solverInfo->solverRootFinding = 1;
304 flag = IDARootInit(idaData->ida_mem, data->modelData->nZeroCrossings, rootsFunctionIDA);
305 if (checkIDAflag(flag)){
306 throwStreamPrint(threadData, "##IDA## Setting root function fails while initialize IDA solver!");
307 }
308 }
309 infoStreamPrint(LOG_SOLVER, 0, "ida uses internal root finding method %s", solverInfo->solverRootFinding?"YES":"NO");
310
311 /* define maximum integration order of dassl */
312 if (omc_flag[FLAG_MAX_ORDER])
313 {
314 maxOrder = atoi(omc_flagValue[FLAG_MAX_ORDER]);
315
316 flag = IDASetMaxOrd(idaData->ida_mem, maxOrder);
317 if (checkIDAflag(flag)){
318 throwStreamPrint(threadData, "##IDA## Failed to set max integration order!");
319 }
320 }
321 infoStreamPrint(LOG_SOLVER, 0, "maximum integration order %d", maxOrder);
322
323
324 /* if FLAG_NOEQUIDISTANT_GRID is set, choose ida step method */
325 if (omc_flag[FLAG_NOEQUIDISTANT_GRID])
326 {
327 idaData->internalSteps = 1; /* TRUE */
328 solverInfo->solverNoEquidistantGrid = 1;
329 }
330 else
331 {
332 idaData->internalSteps = 0; /* FALSE */
333 }
334 infoStreamPrint(LOG_SOLVER, 0, "use equidistant time grid %s", idaData->internalSteps?"NO":"YES");
335
336 /* check if Flags FLAG_NOEQUIDISTANT_OUT_FREQ or FLAG_NOEQUIDISTANT_OUT_TIME are set */
337 if (idaData->internalSteps){
338 if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ])
339 {
340 idaData->stepsFreq = atoi(omc_flagValue[FLAG_NOEQUIDISTANT_OUT_FREQ]);
341 }
342 else if (omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME])
343 {
344 idaData->stepsTime = atof(omc_flagValue[FLAG_NOEQUIDISTANT_OUT_TIME]);
345 flag = IDASetMaxStep(idaData->ida_mem, idaData->stepsTime);
346 if (checkIDAflag(flag)){
347 throwStreamPrint(threadData, "##IDA## Setting max steps of the IDA solver!");
348 }
349 infoStreamPrint(LOG_SOLVER, 0, "maximum step size %g", idaData->stepsTime);
350 } else {
351 idaData->stepsFreq = 1;
352 idaData->stepsTime = 0.0;
353 }
354
355 if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ] && omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]){
356 warningStreamPrint(LOG_STDOUT, 0, "The flags are \"noEquidistantOutputFrequency\" "
357 "and \"noEquidistantOutputTime\" are in opposition "
358 "to each other. The flag \"noEquidistantOutputFrequency\" superiors.");
359 }
360 infoStreamPrint(LOG_SOLVER, 0, "as the output frequency control is used: %d", idaData->stepsFreq);
361 infoStreamPrint(LOG_SOLVER, 0, "as the output frequency time step control is used: %f", idaData->stepsTime);
362 }
363
364 /* if FLAG_IDA_LS is set, choose ida linear solver method */
365 if (omc_flag[FLAG_IDA_LS])
366 {
367 for(i=1; i< IDA_LS_MAX;i++)
368 {
369 if(!strcmp((const char*)omc_flagValue[FLAG_IDA_LS], IDA_LS_METHOD[i])){
370 idaData->linearSolverMethod = (int)i;
371 break;
372 }
373 }
374 if(idaData->linearSolverMethod == IDA_LS_UNKNOWN)
375 {
376 if (ACTIVE_WARNING_STREAM(LOG_SOLVER)(showAllWarnings || useStream[LOG_SOLVER]))
377 {
378 warningStreamPrint(LOG_SOLVER, 1, "unrecognized ida linear solver method %s, current options are:", (const char*)omc_flagValue[FLAG_IDA_LS]);
379 for(i=1; i < IDA_LS_MAX; ++i)
380 {
381 warningStreamPrint(LOG_SOLVER, 0, "%-15s [%s]", IDA_LS_METHOD[i], IDA_LS_METHOD_DESC[i]);
382 }
383 messageClose(LOG_SOLVER);
384 }
385 throwStreamPrint(threadData,"unrecognized ida linear solver method %s", (const char*)omc_flagValue[FLAG_IDA_LS]);
386 }
387 }
388 else
389 {
390 idaData->linearSolverMethod = IDA_LS_KLU;
391 }
392
393 switch (idaData->linearSolverMethod){
394 case IDA_LS_SPGMR:
395 flag = IDASpgmr(idaData->ida_mem, 0);
396 break;
397 case IDA_LS_SPBCG:
398 flag = IDASpbcg(idaData->ida_mem, 0);
399 break;
400 case IDA_LS_SPTFQMR:
401 flag = IDASptfqmr(idaData->ida_mem, 0);
402 break;
403 case IDA_LS_DENSE:
404 flag = IDADense(idaData->ida_mem, idaData->N);
405 break;
406 case IDA_LS_KLU:
407 /* Set KLU after initialized sparse pattern of the jacobian for nnz */
408 break;
409 default:
410 throwStreamPrint(threadData,"unrecognized linear solver method %s", (const char*)omc_flagValue[FLAG_IDA_LS]);
411 break;
412 }
413 if (checkIDAflag(flag)){
414 throwStreamPrint(threadData, "##IDA## Setting linear solver method fails while initialize IDA solver!");
415 } else {
416 infoStreamPrint(LOG_SOLVER, 0, "IDA linear solver method selected %s", IDA_LS_METHOD_DESC[idaData->linearSolverMethod]);
417 }
418
419
420 /* if FLAG_JACOBIAN is set, choose jacobian calculation method */
421 if (omc_flag[FLAG_JACOBIAN])
422 {
423 for(i=1; i< JAC_MAX;i++)
424 {
425 if(!strcmp((const char*)omc_flagValue[FLAG_JACOBIAN], JACOBIAN_METHOD[i])){
426 idaData->jacobianMethod = (int)i;
427 break;
428 }
429 }
430 if(idaData->daeMode && ( idaData->jacobianMethod == COLOREDSYMJAC || idaData->jacobianMethod == SYMJAC ))
431 {
432 warningStreamPrint(LOG_STDOUT, 1, "Symbolic Jacobians are currently not supported by DAE Mode. Switch back to numerical Jacobian!");
433 idaData->jacobianMethod = COLOREDNUMJAC;
434 }
435 if(idaData->jacobianMethod == SYMJAC)
436 {
437 warningStreamPrint(LOG_STDOUT, 1, "Symbolic Jacobians without coloring are currently not supported by IDA. Colored symbolical Jacobian will be used.");
438 idaData->jacobianMethod = COLOREDSYMJAC;
439 }
440 if(idaData->jacobianMethod == NUMJAC)
441 {
442 warningStreamPrint(LOG_STDOUT, 1, "Numerical Jacobians without coloring are currently not supported by IDA. Colored numerical Jacobian will be used.");
443 idaData->jacobianMethod = COLOREDNUMJAC;
444 }
445 if(idaData->jacobianMethod == JAC_UNKNOWN)
446 {
447 if (ACTIVE_WARNING_STREAM(LOG_SOLVER)(showAllWarnings || useStream[LOG_SOLVER]))
448 {
449 warningStreamPrint(LOG_SOLVER, 1, "unrecognized jacobian calculation method %s, current options are:", (const char*)omc_flagValue[FLAG_JACOBIAN]);
450 for(i=1; i < JAC_MAX; ++i)
451 {
452 warningStreamPrint(LOG_SOLVER, 0, "%-15s [%s]", JACOBIAN_METHOD[i], JACOBIAN_METHOD_DESC[i]);
453 }
454 messageClose(LOG_SOLVER);
455 }
456 throwStreamPrint(threadData,"unrecognized jacobian calculation method %s", (const char*)omc_flagValue[FLAG_JACOBIAN]);
457 }
458 /* default case colored numerical jacobian */
459 }
460 else
461 {
462 idaData->jacobianMethod = COLOREDNUMJAC;
463 }
464
465 /* selects the calculation method of the jacobian */
466 /* in daeMode sparse pattern is already initialized in DAEMODE_DATA */
467 if(!idaData->daeMode && (idaData->jacobianMethod == COLOREDNUMJAC ||
468 idaData->jacobianMethod == COLOREDSYMJAC ||
469 idaData->jacobianMethod == SYMJAC))
470 {
471 ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
472 if (data->callback->initialAnalyticJacobianA(data, threadData, jacobian))
473 {
474 infoStreamPrint(LOG_STDOUT, 0, "Jacobian or SparsePattern is not generated or failed to initialize! Switch back to normal.");
475 idaData->jacobianMethod = INTERNALNUMJAC;
476 if (idaData->linearSolverMethod == IDA_LS_KLU)
477 {
478 idaData->linearSolverMethod = IDA_LS_DENSE;
479 flag = IDADense(idaData->ida_mem, idaData->N);
480 warningStreamPrint(LOG_STDOUT, 0, "IDA linear solver method also switched back to %s", IDA_LS_METHOD_DESC[idaData->linearSolverMethod]);
481 }
482 } else {
483 ANALYTIC_JACOBIAN* jac = &data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A];
484 infoStreamPrint(LOG_SIMULATION, 1, "Initialized colored Jacobian:");
485 infoStreamPrint(LOG_SIMULATION, 0, "columns: %d rows: %d", jac->sizeCols, jac->sizeRows);
486 infoStreamPrint(LOG_SIMULATION, 0, "NNZ: %d colors: %d", jac->sparsePattern->numberOfNoneZeros, jac->sparsePattern->maxColors);
487 messageClose(LOG_SIMULATION);
488 }
489
490 }
491 /* set up the appropriate function pointer */
492 if (idaData->linearSolverMethod == IDA_LS_KLU)
493 {
494 if (idaData->daeMode)
495 {
496 idaData->NNZ = data->simulationInfo->daeModeData->sparsePattern->numberOfNoneZeros;
497 flag = IDAKLU(idaData->ida_mem, idaData->N, idaData->NNZ);
498 idaData->jacobianMethod = COLOREDNUMJAC; /* In DAE mode only numerical matrix is supported */
499 }
500 else
501 {
502 idaData->NNZ = data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern->numberOfNoneZeros;
503 flag = IDAKLU(idaData->ida_mem, idaData->N, idaData->NNZ);
504 /* to add a cj identety matrix */
505 idaData->tmpJac = NewSparseMat(idaData->N, idaData->N, idaData->N);
506 }
507 if (checkIDAflag(flag)){
508 throwStreamPrint(threadData, "##IDA## Setup of linear solver KLU failed!");
509 }
510
511 idaData->allocatedParMem = 0; /* false */
512
513 switch (idaData->jacobianMethod){
514 case SYMJAC:
515 case NUMJAC:
516 case COLOREDSYMJAC:
517 case COLOREDNUMJAC:
518 flag = IDASlsSetSparseJacFn(idaData->ida_mem, callSparseJacobian);
519#ifdef USE_PARJAC
520 allocateThreadLocalJacobians(data, &(idaData->jacColumns));
521 idaData->allocatedParMem = 1; /* true */
522#endif
523 break;
524 default:
525 throwStreamPrint(threadData,"For the klu solver jacobian calculation method has to be %s or %s", JACOBIAN_METHOD[COLOREDSYMJAC], JACOBIAN_METHOD[COLOREDNUMJAC]);
526 break;
527 }
528 }
529 else
530 {
531 switch (idaData->jacobianMethod){
532 case SYMJAC:
533 case COLOREDSYMJAC:
534 case COLOREDNUMJAC:
535 case NUMJAC:
536 /* set jacobian function */
537 flag = IDADlsSetDenseJacFn(idaData->ida_mem, callDenseJacobian);
538#ifdef USE_PARJAC
539 allocateThreadLocalJacobians(data, &(idaData->jacColumns));
540 idaData->allocatedParMem = 1; /* true */
541#endif
542 break;
543 case INTERNALNUMJAC:
544 break;
545 default:
546 throwStreamPrint(threadData,"unrecognized jacobian calculation method %s", (const char*)omc_flagValue[FLAG_JACOBIAN]);
547 break;
548 }
549 }
550
551 if (checkIDAflag(flag)){
552 throwStreamPrint(threadData, "##IDA## Setting jacobian function fails while initialize IDA solver!");
553 } else {
554 infoStreamPrint(LOG_SOLVER, 0, "jacobian is calculated by %s", JACOBIAN_METHOD_DESC[idaData->jacobianMethod]);
555 }
556
557 /* set max error test fails */
558 if (omc_flag[FLAG_IDA_MAXERRORTESTFAIL])
559 {
560 int maxErrorTestFails = atoi(omc_flagValue[FLAG_IDA_MAXERRORTESTFAIL]);
561 flag = IDASetMaxErrTestFails(idaData->ida_mem, maxErrorTestFails);
562 if (checkIDAflag(flag)){
563 throwStreamPrint(threadData, "##IDA## set IDASetMaxErrTestFails failed!");
564 }
565 }
566
567 /* set maximum number of nonlinear solver iterations at one step */
568 if (omc_flag[FLAG_IDA_MAXNONLINITERS])
569 {
570 int maxNonlinIters = atoi(omc_flagValue[FLAG_IDA_MAXNONLINITERS]);
571 flag = IDASetMaxNonlinIters(idaData->ida_mem, maxNonlinIters);
572 if (checkIDAflag(flag)){
573 throwStreamPrint(threadData, "##IDA## set IDASetMaxNonlinIters failed!");
574 }
575 }
576
577 /* maximum number of nonlinear solver convergence failures at one step */
578 if (omc_flag[FLAG_IDA_MAXCONVFAILS])
579 {
580 int maxConvFails = atoi(omc_flagValue[FLAG_IDA_MAXCONVFAILS]);
581 flag = IDASetMaxConvFails(idaData->ida_mem, maxConvFails);
582 if (checkIDAflag(flag)){
583 throwStreamPrint(threadData, "##IDA## set IDASetMaxConvFails failed!");
584 }
585 }
586
587 /* safety factor in the nonlinear convergence test */
588 if (omc_flag[FLAG_IDA_NONLINCONVCOEF])
589 {
590 double nonlinConvCoef = atof(omc_flagValue[FLAG_IDA_NONLINCONVCOEF]);
591 flag = IDASetNonlinConvCoef(idaData->ida_mem, nonlinConvCoef);
592 if (checkIDAflag(flag)){
593 throwStreamPrint(threadData, "##IDA## set IDASetNonlinConvCoef failed!");
594 }
595 }
596
597 /* configure algebraic variables as such */
598 if (idaData->daeMode)
599 {
600 if (omc_flag[FLAG_NO_SUPPRESS_ALG])
601 {
602 flag = IDASetSuppressAlg(idaData->ida_mem, TRUE1);
603 if (checkIDAflag(flag)){
604 throwStreamPrint(threadData, "##IDA## Suppress algebraic variables in the local error test failed");
605 }
606 }
607 for(i=0; i<idaData->N; ++i)
608 {
609 tmp[i] = (i<data->modelData->nStates)? 1.0: 0.0;
610 }
611
612 flag = IDASetId(idaData->ida_mem, N_VMake_Serial(idaData->N,tmp));
613 if (checkIDAflag(flag)){
614 throwStreamPrint(threadData, "##IDA## Mark algebraic variables as such failed!");
615 }
616 }
617
618 /* define initial step size */
619 if (omc_flag[FLAG_INITIAL_STEP_SIZE])
620 {
621 double initialStepSize = atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]);
622
623 assertStreamPrint(threadData, initialStepSize >= DASSL_STEP_EPS, "Selected initial step size %e is too small.", initialStepSize)if (!(initialStepSize >= DASSL_STEP_EPS)) {throwStreamPrint
((threadData), "Selected initial step size %e is too small.",
initialStepSize); ((void) sizeof ((0) ? 1 : 0), __extension__
({ if (0) ; else __assert_fail ("0", "./simulation/solver/ida_solver.c"
, 623, __extension__ __PRETTY_FUNCTION__); }));}
;
624
625 flag = IDASetInitStep(idaData->ida_mem, initialStepSize);
626 if (checkIDAflag(flag)){
627 throwStreamPrint(threadData, "##IDA## Set initial step size failed!");
628 }
629 infoStreamPrint(LOG_SOLVER, 0, "initial step size: %g", initialStepSize);
630 }
631 else
632 {
633 infoStreamPrint(LOG_SOLVER, 0, "initial step size is set automatically.");
634 }
635
636 /* configure the sensitivities part */
637 idaData->idaSmode = omc_flag[FLAG_IDAS] ? 1 : 0;
638
639 if (idaData->idaSmode)
640 {
641 idaData->Np = data->modelData->nSensitivityParamVars;
642 idaData->yS = N_VCloneVectorArray_Serial(idaData->Np, idaData->y);
643 idaData->ySp = N_VCloneVectorArray_Serial(idaData->Np, idaData->yp);
644
645 for(i=0; i<idaData->Np; ++i)
646 {
647 int j;
648 for(j=0; j<idaData->N; ++j)
649 {
650 NV_Ith_S(idaData->yS[i],j)( ( ( (N_VectorContent_Serial)(idaData->yS[i]->content)
)->data )[j] )
= 0;
651 NV_Ith_S(idaData->ySp[i],j)( ( ( (N_VectorContent_Serial)(idaData->ySp[i]->content
) )->data )[j] )
= 0;
652 }
653 }
654
655 flag = IDASensInit(idaData->ida_mem, idaData->Np, IDA_SIMULTANEOUS1, NULL((void*)0), idaData->yS, idaData->ySp);
656 if (checkIDAflag(flag)){
657 throwStreamPrint(threadData, "##IDA## set IDASensInit failed!");
658 }
659
660 flag = IDASetSensParams(idaData->ida_mem, data->simulationInfo->realParameter, NULL((void*)0), data->simulationInfo->sensitivityParList);
661 if (checkIDAflag(flag)){
662 throwStreamPrint(threadData, "##IDA## set IDASetSensParams failed!");
663 }
664
665 flag = IDASetSensDQMethod(idaData->ida_mem, IDA_FORWARD2, 0);
666 if (checkIDAflag(flag)){
667 throwStreamPrint(threadData, "##IDA## set IDASetSensDQMethod failed!");
668 }
669
670 flag = IDASensEEtolerances(idaData->ida_mem);
671 if (checkIDAflag(flag)){
672 throwStreamPrint(threadData, "##IDA## set IDASensEEtolerances failed!");
673 }
674/*
675 flag = IDASetSensErrCon(idaData->ida_mem, TRUE);
676 if (checkIDAflag(flag)){
677 throwStreamPrint(threadData, "##IDA## set IDASensEEtolerances failed!");
678 }
679*/
680 /* allocate result workspace */
681 idaData->ySResult = N_VCloneVectorArrayEmpty_Serial(idaData->Np, idaData->y);
682 for(i = 0; i < idaData->Np; ++i)
683 {
684 N_VSetArrayPointer_Serial((data->simulationInfo->sensitivityMatrix + i*idaData->N), idaData->ySResult[i]);
685 }
686 }
687 if (compiledInDAEMode){
688 idaDataGlobal = idaData;
689 initializedSolver = 1;
690 data->callback->functionDAE = ida_event_update;
691 }
692 messageClose(LOG_SOLVER);
693
694 if (measure_time_flag) rt_clear(SIM_TIMER_SOLVER12); /* Initialization should not add this timer... */
695
696 free(tmp);
697 TRACE_POP
698 return 0;
699}
700
701/* deinitialize ida data */
702int ida_solver_deinitial(IDA_SOLVER *idaData)
703{
704 TRACE_PUSH
705
706 free(idaData->simData);
707 free(idaData->ysave);
708 free(idaData->ypsave);
709 free(idaData->delta_hh);
710 if (!idaData->daeMode && idaData->linearSolverMethod == IDA_LS_KLU){
711 DestroySparseMat(idaData->tmpJac);
712 }
713
714 if (idaData->daeMode)
715 {
716 free(idaData->states);
717 free(idaData->statesDer);
718 }
719
720 if (idaData->idaSmode)
721 {
722 N_VDestroyVectorArray_Serial(idaData->yS, idaData->Np);
723 N_VDestroyVectorArray_Serial(idaData->ySp, idaData->Np);
724 N_VDestroyVectorArray_Serial(idaData->ySResult, idaData->Np);
725 }
726
727 N_VDestroy_Serial(idaData->errwgt);
728 N_VDestroy_Serial(idaData->newdelta);
729
730#ifdef USE_PARJAC
731 if (idaData->allocatedParMem) {
732 freeAnalyticalJacobian(&(idaData->jacColumns));
733 idaData->allocatedParMem = 0;
734 }
735#endif
736
737 IDAFree(&idaData->ida_mem);
738
739 TRACE_POP
740 return 0;
741}
742
743int ida_event_update(DATA* data, threadData_t *threadData)
744{
745 IDA_SOLVER *idaData = idaDataGlobal;
746 int flag;
747 long nonLinIters;
748 double init_h;
749 if (compiledInDAEMode){
750 if (initializedSolver){
751
752 data->simulationInfo->needToIterate = 0;
753
754 memcpy(idaData->states, data->localData[0]->realVars, sizeof(double)*data->modelData->nStates);
755 getAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
756 memcpy(idaData->statesDer, data->localData[0]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
757
758 /* update inner algebraic get new values from data */
759 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
760 evaluateDAEResiduals_wrapperEventUpdate(data, threadData);
761 getAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
762 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
763
764 infoStreamPrint(LOG_SOLVER, 0, "##IDA## do event update at %.15g", data->localData[0]->timeValue);
765 flag = IDAReInit(idaData->ida_mem,
766 data->localData[0]->timeValue,
767 idaData->y,
768 idaData->yp);
769
770 /* get initial step to provide a direction of the solution */
771 IDAGetActualInitStep(idaData->ida_mem, &init_h);
772 /* provide a feasible step-size if it's too small */
773 if (init_h < DBL_EPSILON2.2204460492503131e-16){
774 init_h = DBL_EPSILON2.2204460492503131e-16;
775 IDASetInitStep(idaData->ida_mem, init_h);
776 infoStreamPrint(LOG_SOLVER, 0, "##IDA## corrected step-size at %.15g", init_h);
777 }
778
779 /* increase limits of the non-linear solver */
780 IDASetMaxNumStepsIC(idaData->ida_mem, 2*idaData->N*10);
781 IDASetMaxNumJacsIC(idaData->ida_mem, 2*idaData->N*10);
782 IDASetMaxNumItersIC(idaData->ida_mem, 2*idaData->N*10);
783 /* Calc Consistent y_algebraic and y_prime with current y */
784 flag = IDACalcIC(idaData->ida_mem, IDA_YA_YDP_INIT1, data->localData[0]->timeValue+init_h);
785
786 /* debug */
787 IDAGetNumNonlinSolvIters(idaData->ida_mem, &nonLinIters);
788 infoStreamPrint(LOG_SOLVER, 0, "##IDA## IDACalcIC run status %d.\nIterations : %ld\n", flag, nonLinIters);
789
790 /* try again without line search if first try fails */
791 if (checkIDAflag(flag)){
792 infoStreamPrint(LOG_SOLVER, 0, "##IDA## first event iteration failed. Start next try without line search!");
793 IDASetLineSearchOffIC(idaData->ida_mem, TRUE1);
794 flag = IDACalcIC(idaData->ida_mem, IDA_YA_YDP_INIT1, data->localData[0]->timeValue+data->simulationInfo->tolerance);
795 IDAGetNumNonlinSolvIters(idaData->ida_mem, &nonLinIters);
796 infoStreamPrint(LOG_SOLVER, 0, "##IDA## IDACalcIC run status %d.\nIterations : %ld\n", flag, nonLinIters);
797 if (checkIDAflag(flag)){
798 throwStreamPrint(threadData, "##IDA## discrete update failed flag %d!", flag);
799 }
800 }
801 /* obtain consistent values of y_algebraic and y_prime */
802 IDAGetConsistentIC(idaData->ida_mem, idaData->y, idaData->yp);
803
804 /* update inner algebraic variables */
805 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
806 evaluateDAEResiduals_wrapperEventUpdate(data, threadData);
807
808 memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
809 // and also algebraic vars
810 setAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
811 memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
812 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
813
814 /* reset initial step size again to default */
815 IDASetInitStep(idaData->ida_mem, 0.0);
816 }
817 }
818 else{
819 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
820 data->callback->functionDAE(data, threadData);
821 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
822 }
823
824 return 0;
825}
826
827/* main ida function to make a step */
828int ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
829{
830 TRACE_PUSH
831 double tout = 0;
832 int i = 0, flag;
833 int retVal = 0, finished = FALSE0;
834 int saveJumpState;
835 long int tmp;
836 static unsigned int stepsOutputCounter = 1;
837 int stepsMode;
838 int restartAfterLSFail = 0;
839
840 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
841
842 IDA_SOLVER *idaData = (IDA_SOLVER*) solverInfo->solverData;
843
844 SIMULATION_DATA *sData = data->localData[0];
845 SIMULATION_DATA *sDataOld = data->localData[1];
846 MODEL_DATA *mData = (MODEL_DATA*) data->modelData;
847
848
849 /* alloc all work arrays */
850 if (!idaData->daeMode)
851 {
852 N_VSetArrayPointer(data->localData[0]->realVars, idaData->y);
853 N_VSetArrayPointer(data->localData[1]->realVars + data->modelData->nStates, idaData->yp);
854 }
855
856 if (solverInfo->didEventStep)
857 {
858 idaData->setInitialSolution = 0;
859 }
860
861 /* reinit solver */
862 if (!idaData->setInitialSolution)
863 {
864 debugStreamPrint(LOG_SOLVER, 0, "Re-initialized IDA Solver");
865
866 /* initialize states and der(states) */
867 if (idaData->daeMode)
868 {
869 memcpy(idaData->states, data->localData[0]->realVars, sizeof(double)*data->modelData->nStates);
870 /* and also algebraic vars */
871 getAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
872 memcpy(idaData->statesDer, data->localData[0]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
873 }
874
875 /* calculate matrix for residual scaling */
876 if (omc_flag[FLAG_IDA_SCALING])
877 {
878 getScalingFactors(data, idaData, NULL((void*)0));
879
880 /* scale idaData->y and idaData->yp */
881 infoStreamPrint(LOG_SOLVER_V, 1, "Scale y and yp");
882 idaScaleData(idaData);
883 messageClose(LOG_SOLVER_V);
884 }
885
886 flag = IDAReInit(idaData->ida_mem,
887 solverInfo->currentTime,
888 idaData->y,
889 idaData->yp);
890
891 if (checkIDAflag(flag)){
892 throwStreamPrint(threadData, "##IDA## Something goes wrong while reinit IDA solver after event!");
893 }
894
895 /* calculate matrix for residual scaling */
896 if (omc_flag[FLAG_IDA_SCALING])
897 {
898 /* scale idaData->y and idaData->yp */
899 idaReScaleData(idaData);
900 }
901
902 if (idaData->idaSmode)
903 {
904 for(i=0; i<idaData->Np; ++i)
905 {
906 int j;
907 for(j=0; j<idaData->N; ++j)
908 {
909 NV_Ith_S(idaData->yS[i],j)( ( ( (N_VectorContent_Serial)(idaData->yS[i]->content)
)->data )[j] )
= 0;
910 NV_Ith_S(idaData->ySp[i],j)( ( ( (N_VectorContent_Serial)(idaData->ySp[i]->content
) )->data )[j] )
= 0;
911 }
912 }
913 flag = IDASensReInit(idaData->ida_mem, IDA_SIMULTANEOUS1, idaData->yS, idaData->ySp);
914 if (checkIDAflag(flag)){
915 throwStreamPrint(threadData, "##IDA## set IDASensInit failed!");
916 }
917 }
918 idaData->setInitialSolution = 1;
919 }
920
921 saveJumpState = threadData->currentErrorStage;
922 threadData->currentErrorStage = ERROR_INTEGRATOR;
923
924 /* try */
925#if !defined(OMC_EMCC)
926 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
927#endif
928
929
930 /* Check that tout is not less than timeValue otherwise the solver
931 * will come in trouble.
932 * If that is the case we skip the current step. */
933 if (solverInfo->currentStepSize < DASSL_STEP_EPS)
934 {
935 infoStreamPrint(LOG_SOLVER, 0, "Desired step to small try next one");
936 infoStreamPrint(LOG_SOLVER, 0, "Interpolate linear");
937
938 /* linear extrapolation */
939 for(i = 0; i < idaData->N; i++)
940 {
941 NV_Ith_S(idaData->y, i)( ( ( (N_VectorContent_Serial)(idaData->y->content) )->
data )[i] )
= NV_Ith_S(idaData->y, i)( ( ( (N_VectorContent_Serial)(idaData->y->content) )->
data )[i] )
+ NV_Ith_S(idaData->yp, i)( ( ( (N_VectorContent_Serial)(idaData->yp->content) )->
data )[i] )
* solverInfo->currentStepSize;
942 }
943 sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
944 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
945 data->callback->functionODE(data, threadData);
946 solverInfo->currentTime = sData->timeValue;
947
948 TRACE_POP
949 return 0;
950 }
951
952
953 /* Calculate steps until TOUT is reached */
954 if (idaData->internalSteps)
955 {
956 /* If internalSteps are selected, let IDA run to stopTime or next sample event */
957 if (data->simulationInfo->nextSampleEvent < data->simulationInfo->stopTime)
958 {
959 tout = data->simulationInfo->nextSampleEvent;
960 }
961 else
962 {
963 tout = data->simulationInfo->stopTime;
964 }
965 stepsMode = IDA_ONE_STEP2;
966 flag = IDASetStopTime(idaData->ida_mem, tout);
967 if (checkIDAflag(flag)){
968 throwStreamPrint(threadData, "##IDA## Something goes wrong while set stopTime!");
969 }
970 }
971 else
972 {
973 tout = solverInfo->currentTime + solverInfo->currentStepSize;
974 stepsMode = IDA_NORMAL1;
975 }
976
977
978 do
979 {
980 infoStreamPrint(LOG_SOLVER, 1, "##IDA## new step from %.15g to %.15g", solverInfo->currentTime, tout);
981
982 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
983 /* read input vars */
984 externalInputUpdate(data);
985 data->callback->input_function(data, threadData);
986 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
987
988 if (omc_flag[FLAG_IDA_SCALING])
989 {
990 /* scale idaData->y and idaData->yp */
991 idaScaleData(idaData);
992 }
993
994 flag = IDASolve(idaData->ida_mem, tout, &solverInfo->currentTime, idaData->y, idaData->yp, stepsMode);
995
996 if (omc_flag[FLAG_IDA_SCALING])
997 {
998 /* rescale idaData->y and idaData->yp */
999 idaReScaleData(idaData);
1000 }
1001
1002 /* set time to current time */
1003 sData->timeValue = solverInfo->currentTime;
1004
1005 /* error handling */
1006 if ( !checkIDAflag(flag) && solverInfo->currentTime >= tout)
1007 {
1008 infoStreamPrint(LOG_SOLVER, 0, "##IDA## step done to time = %.15g", solverInfo->currentTime);
1009 finished = TRUE1;
1010 }
1011 else
1012 {
1013 if (!checkIDAflag(flag))
1014 {
1015 infoStreamPrint(LOG_SOLVER, 0, "##IDA## continue integration time = %.15g", solverInfo->currentTime);
1016 }
1017 else if (flag == IDA_ROOT_RETURN2)
1018 {
1019 infoStreamPrint(LOG_SOLVER, 0, "##IDA## root found at time = %.15g", solverInfo->currentTime);
1020 finished = TRUE1;
1021 }
1022 else if (flag == IDA_TOO_MUCH_WORK-1)
1023 {
1024 warningStreamPrint(LOG_SOLVER, 0, "##IDA## has done too much work with small steps at time = %.15g", solverInfo->currentTime);
1025 }
1026 else if (flag == IDA_LSETUP_FAIL-6 && !restartAfterLSFail )
1027 {
1028 flag = IDAReInit(idaData->ida_mem,
Value stored to 'flag' is never read
1029 solverInfo->currentTime,
1030 idaData->y,
1031 idaData->yp);
1032 restartAfterLSFail = 1;
1033 warningStreamPrint(LOG_SOLVER, 0, "##IDA## linear solver failed try once again = %.15g", solverInfo->currentTime);
1034 }
1035 else
1036 {
1037 infoStreamPrint(LOG_STDOUT, 0, "##IDA## %d error occurred at time = %.15g", flag, solverInfo->currentTime);
1038 finished = TRUE1;
1039 retVal = flag;
1040 }
1041 }
1042
1043 /* closing new step message */
1044 messageClose(LOG_SOLVER);
1045
1046 /* emit step, if step mode is selected */
1047 if (idaData->internalSteps)
1048 {
1049 infoStreamPrint(LOG_SOLVER, 0, "##IDA## noEquadistant stepsOutputCounter %d by freq %d at time = %.15g", stepsOutputCounter, idaData->stepsFreq, solverInfo->currentTime);
1050 if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ]){
1051 /* output every n-th time step */
1052 if (stepsOutputCounter >= idaData->stepsFreq){
1053 stepsOutputCounter = 1; /* next line set it to one */
1054 infoStreamPrint(LOG_SOLVER, 0, "##IDA## noEquadistant output %d by freq at time = %.15g", stepsOutputCounter, solverInfo->currentTime);
1055 break;
1056 }
1057 stepsOutputCounter++;
1058 } else if (omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]){
1059 /* output when time>=k*timeValue */
1060 if (solverInfo->currentTime > stepsOutputCounter * idaData->stepsTime){
1061 stepsOutputCounter++;
1062 infoStreamPrint(LOG_SOLVER, 0, "##IDA## noEquadistant output %d by time freq at time = %.15g", stepsOutputCounter, solverInfo->currentTime);
1063 break;
1064 }
1065 } else {
1066 break;
1067 }
1068 }
1069
1070 } while(!finished);
1071
1072#if !defined(OMC_EMCC)
1073 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1074#endif
1075 threadData->currentErrorStage = saveJumpState;
1076
1077 /* if a state event occurs than no sample event does need to be activated */
1078 if (data->simulationInfo->sampleActivated && solverInfo->currentTime < data->simulationInfo->nextSampleEvent)
1079 {
1080 data->simulationInfo->sampleActivated = 0;
1081 }
1082
1083 /* initialize states and der(states) */
1084 if (idaData->daeMode)
1085 {
1086 memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
1087 // and also algebraic vars
1088 setAlgebraicDAEVars(data, idaData->states + data->modelData->nStates);
1089 memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
1090 sData->timeValue = solverInfo->currentTime;
1091 }
1092
1093 /* sensitivity mode */
1094 if (idaData->idaSmode)
1095 {
1096 flag = IDAGetSens(idaData->ida_mem, &solverInfo->currentTime, idaData->ySResult);
1097 if (checkIDAflag(flag)){
1098 throwStreamPrint(threadData, "##IDA## Something goes wrong while obtain results for parameter sensitivities!");
1099 }
1100 }
1101
1102 /* save stats */
1103 /* steps */
1104 tmp = 0;
1105 flag = IDAGetNumSteps(idaData->ida_mem, &tmp);
1106 if (flag == IDA_SUCCESS0)
1107 {
1108 solverInfo->solverStatsTmp[0] = tmp;
1109 }
1110
1111 /* functionODE evaluations */
1112 tmp = 0;
1113 flag = IDAGetNumResEvals(idaData->ida_mem, &tmp);
1114 if (flag == IDA_SUCCESS0)
1115 {
1116 solverInfo->solverStatsTmp[1] = tmp;
1117 }
1118
1119 /* Jacobians evaluations */
1120 tmp = 0;
1121 if (idaData->linearSolverMethod == IDA_LS_KLU)
1122 {
1123 flag = IDASlsGetNumJacEvals(idaData->ida_mem, &tmp);
1124 }
1125 else
1126 {
1127 flag = IDADlsGetNumJacEvals(idaData->ida_mem, &tmp);
1128 }
1129
1130 if (flag == IDA_SUCCESS0)
1131 {
1132 solverInfo->solverStatsTmp[2] = tmp;
1133 }
1134
1135 /* local error test failures */
1136 tmp = 0;
1137 flag = IDAGetNumErrTestFails(idaData->ida_mem, &tmp);
1138 if (flag == IDA_SUCCESS0)
1139 {
1140 solverInfo->solverStatsTmp[3] = tmp;
1141 }
1142
1143 /* local error test failures */
1144 tmp = 0;
1145 flag = IDAGetNumNonlinSolvConvFails(idaData->ida_mem, &tmp);
1146 if (flag == IDA_SUCCESS0)
1147 {
1148 solverInfo->solverStatsTmp[4] = tmp;
1149 }
1150
1151 /* get more statistics */
1152 if (useStream[LOG_SOLVER_V])
1153 {
1154 long int tmp1,tmp2;
1155 double dtmp;
1156
1157 infoStreamPrint(LOG_SOLVER_V, 1, "### IDAStats ###");
1158 /* nonlinear stats */
1159 tmp1 = tmp2 = 0;
1160 flag = IDAGetNonlinSolvStats(idaData->ida_mem, &tmp1, &tmp2);
1161 infoStreamPrint(LOG_SOLVER_V, 0, " ## Cumulative number of nonlinear iterations performed: %ld", tmp1);
1162 infoStreamPrint(LOG_SOLVER_V, 0, " ## Cumulative number of nonlinear convergence failures that have occurred: %ld", tmp2);
1163
1164 /* others */
1165 flag = IDAGetTolScaleFactor(idaData->ida_mem, &dtmp);
1166 infoStreamPrint(LOG_SOLVER_V, 0, " ## Suggested scaling factor for user tolerances: %g", dtmp);
1167
1168 flag = IDAGetNumLinSolvSetups(idaData->ida_mem, &tmp1);
1169 infoStreamPrint(LOG_SOLVER_V, 0, " ## Number of calls made to the linear solver setup function: %ld", tmp1);
1170
1171 messageClose(LOG_SOLVER_V);
1172 }
1173
1174 infoStreamPrint(LOG_SOLVER, 0, "##IDA## Finished Integrator step.");
1175 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1176
1177 TRACE_POP
1178 return retVal;
1179}
1180
1181int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, void* userData)
1182{
1183 TRACE_PUSH
1184 IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
1185 DATA* data = (DATA*)(((IDA_USERDATA*)((IDA_SOLVER*)userData)->simData)->data);
1186 threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)userData)->simData)->threadData);
1187
1188 double timeBackup;
1189 long int i;
1190 int saveJumpState;
1191 int success = 0, retVal = 0;
1192 double *states = N_VGetArrayPointer(yy);
1193 double *statesDer = N_VGetArrayPointer(yp);
1194 double *delta = N_VGetArrayPointer(res);
1195
1196 infoStreamPrint(LOG_SOLVER_V, 1, "### eval residualFunctionIDA ###");
1197 /* rescale idaData->y and idaData->yp */
1198 if ((omc_flag[FLAG_IDA_SCALING] && !idaData->disableScaling))
1199 {
1200 idaReScaleData(idaData);
1201 }
1202
1203 if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
1204 {
1205 setContext(data, &time, CONTEXT_ODE);
1206 }
1207 data->localData[0]->timeValue = time;
1208
1209 saveJumpState = threadData->currentErrorStage;
1210 threadData->currentErrorStage = ERROR_INTEGRATOR;
1211
1212 /* try */
1213#if !defined(OMC_EMCC)
1214 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1215#endif
1216
1217 /* if sensitivity mode update also bound parameters*/
1218 if (idaData->idaSmode)
1219 {
1220 data->callback->updateBoundParameters(data, threadData);
1221 }
1222 /* if daeMode update also all dynamic algebraic equations */
1223 if (idaData->daeMode)
1224 {
1225 /* set state, state derivative and dynamic algebraic
1226 * variables for evaluateDAEResiduals evaluation
1227 */
1228 memcpy(data->localData[0]->realVars, states, sizeof(double)*data->modelData->nStates);
1229 memcpy(data->localData[0]->realVars + data->modelData->nStates, statesDer, sizeof(double)*data->modelData->nStates);
1230 setAlgebraicDAEVars(data, states + data->modelData->nStates);
1231 }
1232
1233 /* debug */
1234 if (ACTIVE_STREAM(LOG_DASSL_STATES)(useStream[LOG_DASSL_STATES])){
1235 printCurrentStatesVector(LOG_DASSL_STATES, data->localData[0]->realVars, data, time);
1236 printVector(LOG_DASSL_STATES, "yprime", data->localData[0]->realVars + data->modelData->nStates, data->modelData->nStates, time);
1237 if (idaData->daeMode)
1238 {
1239 printVector(LOG_DASSL_STATES, "yalg", states + data->modelData->nStates, data->simulationInfo->daeModeData->nAlgebraicDAEVars, time);
1240 }
1241 }
1242
1243 /* read input vars */
1244 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1245 externalInputUpdate(data);
1246 data->callback->input_function(data, threadData);
1247 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1248
1249 if (idaData->daeMode)
1250 {
1251 /* eval residual vars */
1252 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1253 data->simulationInfo->daeModeData->evaluateDAEResiduals(data, threadData, EVAL_DYNAMIC);
1254 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1255 /* get residual variables */
1256 for(i=0; i < idaData->N; i++)
1257 {
1258 NV_Ith_S(res, i)( ( ( (N_VectorContent_Serial)(res->content) )->data )[
i] )
= data->simulationInfo->daeModeData->residualVars[i];
1259 infoStreamPrint(LOG_SOLVER_V, 0, "%ld. residual = %e", i, NV_Ith_S(res, i)( ( ( (N_VectorContent_Serial)(res->content) )->data )[
i] )
);
1260 }
1261 }
1262 else
1263 {
1264 /* eval function ODE */
1265 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1266 data->callback->functionODE(data, threadData);
1267 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1268 for(i=0; i < idaData->N; i++)
1269 {
1270 NV_Ith_S(res, i)( ( ( (N_VectorContent_Serial)(res->content) )->data )[
i] )
= data->localData[0]->realVars[data->modelData->nStates + i] - NV_Ith_S(yp, i)( ( ( (N_VectorContent_Serial)(yp->content) )->data )[i
] )
;
1271 infoStreamPrint(LOG_SOLVER_V, 0, "%ld. residual = %e", i, NV_Ith_S(res, i)( ( ( (N_VectorContent_Serial)(res->content) )->data )[
i] )
);
1272 }
1273 }
1274
1275 /* scale res */
1276 if ((omc_flag[FLAG_IDA_SCALING] && !idaData->disableScaling))
1277 {
1278 infoStreamPrint(LOG_SOLVER_V, 1, "scale residuals");
1279 idaScaleVector(res, idaData->resScale, idaData->N);
1280 messageClose(LOG_SOLVER_V);
1281 idaScaleData(idaData);
1282 }
1283
1284 printVector(LOG_DASSL_STATES, "delta", delta, idaData->N, time);
1285 success = 1;
1286#if !defined(OMC_EMCC)
1287 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1288#endif
1289
1290 if (!success) {
1291 retVal = -1;
1292 }
1293
1294 threadData->currentErrorStage = saveJumpState;
1295
1296 if (data->simulationInfo->currentContext == CONTEXT_ODE){
1297 unsetContext(data);
1298 }
1299 messageClose(LOG_SOLVER_V);
1300 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1301
1302 TRACE_POP
1303 return retVal;
1304}
1305
1306int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void* userData)
1307{
1308 TRACE_PUSH
1309 IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
1310 DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
1311 threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)userData)->simData)->threadData);
1312 double *states = N_VGetArrayPointer(yy);
1313 double *statesDer = N_VGetArrayPointer(yp);
1314
1315 int saveJumpState;
1316
1317 infoStreamPrint(LOG_SOLVER_V, 1, "### eval rootsFunctionIDA ###");
1318
1319 if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
1320 {
1321 setContext(data, &time, CONTEXT_EVENTS);
1322 }
1323
1324 /* re-scale idaData->y and idaData->yp to evaluate the equations*/
1325 if (omc_flag[FLAG_IDA_SCALING])
1326 {
1327 idaReScaleData(idaData);
1328 }
1329
1330 saveJumpState = threadData->currentErrorStage;
1331 threadData->currentErrorStage = ERROR_EVENTSEARCH;
1332
1333 if (idaData->daeMode)
1334 {
1335 memcpy(data->localData[0]->realVars, states, sizeof(double)*data->modelData->nStates);
1336 setAlgebraicDAEVars(data, states + data->modelData->nStates);
1337 memcpy(data->localData[0]->realVars + data->modelData->nStates, statesDer, sizeof(double)*data->modelData->nStates);
1338 }
1339
1340 data->localData[0]->timeValue = time;
1341
1342 /* read input vars */
1343 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1344 externalInputUpdate(data);
1345 data->callback->input_function(data, threadData);
1346 /* eval needed equations*/
1347 if (idaData->daeMode){
1348 data->simulationInfo->daeModeData->evaluateDAEResiduals(data, threadData, EVAL_ZEROCROSS);
1349 }
1350 else
1351 {
1352 data->callback->function_ZeroCrossingsEquations(data, threadData);
1353 }
1354
1355 data->callback->function_ZeroCrossings(data, threadData, gout);
1356 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1357
1358 threadData->currentErrorStage = saveJumpState;
1359
1360 /* scale data again */
1361 if (omc_flag[FLAG_IDA_SCALING])
1362 {
1363 idaScaleData(idaData);
1364 }
1365
1366 if (data->simulationInfo->currentContext == CONTEXT_EVENTS){
1367 unsetContext(data);
1368 }
1369 messageClose(LOG_SOLVER_V);
1370 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1371
1372 TRACE_POP
1373 return 0;
1374}
1375
1376/*
1377 * function calculates a jacobian matrix by
1378 * numerical method finite differences with coloring
1379 * into a dense DlsMat matrix
1380 */
1381static
1382int jacColoredNumericalDense(double currentTime, N_Vector yy, N_Vector yp,
1383 N_Vector rr, DlsMat Jac, double cj, void *userData)
1384{
1385 TRACE_PUSH
1386 IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
1387 DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
1388 void* ida_mem = idaData->ida_mem;
1389 const int index = data->callback->INDEX_JAC_A;
1390
1391 /* prepare variables */
1392 double *states = N_VGetArrayPointer(yy);
1393 double *yprime = N_VGetArrayPointer(yp);
1394 double *delta = N_VGetArrayPointer(rr);
1395 double *newdelta = N_VGetArrayPointer(idaData->newdelta);
1396 double *errwgt = N_VGetArrayPointer(idaData->errwgt);
1397
1398 double *delta_hh = idaData->delta_hh;
1399 double *ysave = idaData->ysave;
1400 double *ypsave = idaData->ypsave;
1401
1402 double delta_h = numericalDifferentiationDeltaXsolver;
1403 double delta_hhh;
1404 long int i,j,l,ii;
1405
1406 double currentStep;
1407
1408 /* set values */
1409 IDAGetCurrentStep(ida_mem, &currentStep);
1410 if (!idaData->disableScaling){
1411 IDAGetErrWeights(ida_mem, idaData->errwgt);
1412 }
1413
1414 SPARSE_PATTERN* sparsePattern;
1415
1416 /* set sparse pattern */
1417 if (idaData->daeMode)
1418 {
1419 sparsePattern = data->simulationInfo->daeModeData->sparsePattern;
1420 }
1421 else
1422 {
1423 sparsePattern = data->simulationInfo->analyticJacobians[index].sparsePattern;
1424 }
1425
1426 setContext(data, &currentTime, CONTEXT_JACOBIAN);
1427
1428 for(i = 0; i < sparsePattern->maxColors; i++)
1429 {
1430 for(ii=0; ii < idaData->N; ii++)
1431 {
1432 if(sparsePattern->colorCols[ii]-1 == i)
1433 {
1434 delta_hhh = currentStep * yprime[ii];
1435 delta_hh[ii] = delta_h * fmax(fmax(fabs(states[ii]),fabs(delta_hhh)),fabs(1./errwgt[ii]));
1436 delta_hh[ii] = (delta_hhh >= 0 ? delta_hh[ii] : -delta_hh[ii]);
1437 delta_hh[ii] = (states[ii] + delta_hh[ii]) - states[ii];
1438 ysave[ii] = states[ii];
1439 states[ii] += delta_hh[ii];
1440
1441 if (idaData->daeMode){
1442 ypsave[ii] = yprime[ii];
1443 yprime[ii] += cj * delta_hh[ii];
1444 }
1445
1446 delta_hh[ii] = 1. / delta_hh[ii];
1447 }
1448 }
1449
1450 (*idaData->residualFunction)(currentTime, yy, yp, idaData->newdelta, userData);
1451
1452 increaseJacContext(data);
1453
1454 for(ii = 0; ii < idaData->N; ii++)
1455 {
1456 if(sparsePattern->colorCols[ii]-1 == i)
1457 {
1458 j = sparsePattern->leadindex[ii];
1459 while(j < sparsePattern->leadindex[ii+1])
1460 {
1461 l = sparsePattern->index[j];
1462 DENSE_ELEM(Jac, l, ii)((Jac->cols)[ii][l]) = (newdelta[l] - delta[l]) * delta_hh[ii];
1463 j++;
1464 };
1465 states[ii] = ysave[ii];
1466 if (idaData->daeMode)
1467 {
1468 yprime[ii] = ypsave[ii];
1469 }
1470 }
1471 }
1472 }
1473 unsetContext(data);
1474
1475 TRACE_POP
1476 return 0;
1477}
1478
1479/*
1480 * function calculates the Jacobian matrix symbolically with considering also
1481 * the coloring and pass it in a dense DlsMat matrix.
1482 */
1483static int jacColoredSymbolicalDense(double currentTime, N_Vector yy,
1484 N_Vector yp, N_Vector rr, DlsMat Jac,
1485 double cj, void *userData)
1486{
1487 TRACE_PUSH
1488 IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
1489 DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
1490 threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)idaData->simData)->threadData);
1491 void* ida_mem = idaData->ida_mem;
1492 long int N = idaData->N;
1493 const int index = data->callback->INDEX_JAC_A;
1494 unsigned int i,ii,j, nth;
1495 SPARSE_PATTERN* sparsePattern = data->simulationInfo->analyticJacobians[index].sparsePattern;
1496 ANALYTIC_JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[index]);
1497
1498 /* prepare variables */
1499 double *states = N_VGetArrayPointer(yy);
1500 double *yprime = N_VGetArrayPointer(yp);
1501
1502 setContext(data, &currentTime, CONTEXT_SYM_JACOBIAN); /* Reuse jacobian matrix in KLU solver */
1503
1504 /* Evaluate constant equations if available */
1505 if (jac->constantEqns != NULL((void*)0)) {
1506 jac->constantEqns(data, threadData, jac, NULL((void*)0));
1507 }
1508
1509#ifdef USE_PARJAC
1510 GC_allow_register_threads();
1511#endif
1512
1513#pragma omp parallel default(none) firstprivate(N) shared(i, sparsePattern, idaData, data, threadData, Jac) private(ii, j, nth)
1514{
1515#ifdef USE_PARJAC
1516 /* Register omp-thread in GC */
1517 if(!GC_thread_is_registered()) {
1518 struct GC_stack_base sb;
1519 memset (&sb, 0, sizeof(sb));
1520 GC_get_stack_base(&sb);
1521 GC_register_my_thread (&sb);
1522 }
1523 // ToDo Use always a thread local analyticJacobians (replace simulationInfo->analyticaJacobians)
1524 // These are not the Jacobians of the linear systems! (SimulationInfo->linearSystemData[idx].jacobian)
1525 ANALYTIC_JACOBIAN* t_jac = &(idaData->jacColumns[omc_get_thread_num()]);
1526#else
1527 ANALYTIC_JACOBIAN* t_jac = jac;
1528#endif
1529
1530#pragma omp for
1531 for(i = 0; i < sparsePattern->maxColors; i++)
1532 {
1533 for(ii=0; ii < N; ii++)
1534 {
1535 if(sparsePattern->colorCols[ii]-1 == i)
1536 {
1537 t_jac->seedVars[ii] = 1;
1538 }
1539 }
1540
1541 data->callback->functionJacA_column(data, threadData, t_jac, NULL((void*)0));
1542 increaseJacContext(data);
1543
1544 for(ii = 0; ii < N; ii++)
1545 {
1546 if(sparsePattern->colorCols[ii]-1 == i)
1547 {
1548 nth = sparsePattern->leadindex[ii];
1549 while(nth < sparsePattern->leadindex[ii+1])
1550 {
1551 j = sparsePattern->index[nth];
1552 infoStreamPrint(LOG_JAC, 0, "### symbolical jacobian at [%d,%d] = %f ###", j, ii, t_jac->resultVars[j]);
1553 DENSE_ELEM(Jac, j, ii)((Jac->cols)[ii][j]) = t_jac->resultVars[j];
1554 nth++;
1555 };
1556 }
1557 }
1558
1559 for(ii=0; ii < idaData->N; ii++)
1560 {
1561 t_jac->seedVars[ii] = 0;
1562 }
1563 } // for column
1564} // omp parallel
1565
1566 unsetContext(data);
1567
1568 TRACE_POP
1569 return 0;
1570}
1571
1572/*
1573 * wrapper function to call dense Jacobian
1574 */
1575static int callDenseJacobian(long int Neq, double tt, double cj,
1576 N_Vector yy, N_Vector yp, N_Vector rr,
1577 DlsMat Jac, void *user_data,
1578 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
1579{
1580 TRACE_PUSH
1581 IDA_SOLVER* idaData = (IDA_SOLVER*)user_data;
1582 threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)user_data)->simData)->threadData);
1583 int retVal;
1584
1585 /* profiling */
1586 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1587 rt_tick(SIM_TIMER_JACOBIAN5);
1588
1589 if (idaData->jacobianMethod == COLOREDNUMJAC || idaData->jacobianMethod == NUMJAC)
1590 {
1591 retVal = jacColoredNumericalDense(tt, yy, yp, rr, Jac, cj, user_data);
1592 }
1593 else if (idaData->jacobianMethod == COLOREDSYMJAC || idaData->jacobianMethod == SYMJAC)
1594 {
1595 retVal = jacColoredSymbolicalDense(tt, yy, yp, rr, Jac, cj, user_data);
1596 }
1597 else
1598 {
1599 throwStreamPrint(threadData, "##IDA## Something goes wrong while obtain jacobian matrix!");
1600 }
1601
1602 /* debug */
1603 if (ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC])){
1604 _omc_matrix* dumpJac = _omc_createMatrix(idaData->N, idaData->N, Jac->data);
1605 _omc_printMatrix(dumpJac, "IDA-Solver: Matrix A", LOG_JAC);
1606 _omc_destroyMatrix(dumpJac);
1607 }
1608
1609 /* add cj to diagonal elements and store in Jac */
1610 if (!idaData->daeMode)
1611 {
1612 long int i;
1613 for(i = 0; i < Neq; i++)
1614 {
1615 DENSE_ELEM(Jac, i, i)((Jac->cols)[i][i]) -= (double) cj;
1616 }
1617 }
1618
1619 /* profiling */
1620 rt_accumulate(SIM_TIMER_JACOBIAN5);
1621 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1622
1623 TRACE_POP
1624 return retVal;
1625}
1626
1627/* Element function for sparse matrix set */
1628static void setJacElementKluSparse(int row, int col, int nth, double value, void* spJac, int rows)
1629{
1630 (void) rows; // Unused, needed to match genericColoredSymbolicJacobianEvaluation
1631 SlsMat mat = (SlsMat)spJac;
1632 if (col > 0 && mat->colptrs[col] == 0){
1633 mat->colptrs[col] = nth;
1634 }
1635 mat->rowvals[nth] = row;
1636 mat->data[nth] = value;
1637}
1638
1639/* finish sparse matrix, by fixing colprts */
1640static void finishSparseColPtr(SlsMat mat, int nnz)
1641{
1642 int i;
1643 for(i=1; i<mat->N+1; ++i){
1644 if (mat->colptrs[i] == 0){
1645 mat->colptrs[i] = mat->colptrs[i-1];
1646 }
1647 }
1648 /* finish matrix colptrs */
1649 mat->colptrs[mat->N] = nnz;
1650}
1651
1652/*
1653 * function calculates a jacobian matrix by
1654 * numerical method finite differences with coloring
1655 * into a sparse SlsMat matrix
1656 */
1657static int jacoColoredNumericalSparse(double currentTime, N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, double cj, void *userData)
1658{
1659 TRACE_PUSH
1660 IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
1661 DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
1662 void* ida_mem = idaData->ida_mem;
1663 const int index = data->callback->INDEX_JAC_A;
1664
1665 /* prepare variables */
1666 double *states = N_VGetArrayPointer(yy);
1667 double *yprime = N_VGetArrayPointer(yp);
1668 double *delta = N_VGetArrayPointer(rr);
1669 double *newdelta = N_VGetArrayPointer(idaData->newdelta);
1670 double *errwgt = N_VGetArrayPointer(idaData->errwgt);
1671
1672 SPARSE_PATTERN* sparsePattern;
1673
1674 double *ysave = idaData->ysave;
1675 double *ypsave = idaData->ypsave;
1676
1677 double delta_h = numericalDifferentiationDeltaXsolver;
1678 double *delta_hh = idaData->delta_hh;
1679 double delta_hhh;
1680 double deltaInv;
1681
1682 long int i,j,ii;
1683 int nth = 0;
1684 int disBackup = idaData->disableScaling;
1685
1686 double currentStep;
1687
1688 infoStreamPrint(LOG_SOLVER_V, 1, "### eval jacobianSparseNumIDA ###");
1689 /* set values */
1690 IDAGetCurrentStep(ida_mem, &currentStep);
1691 if (!idaData->disableScaling){
1692 IDAGetErrWeights(ida_mem, idaData->errwgt);
1693 }
1694
1695 /* set sparse pattern */
1696 if (idaData->daeMode)
1697 {
1698 sparsePattern = data->simulationInfo->daeModeData->sparsePattern;
1699 }
1700 else
1701 {
1702 sparsePattern = data->simulationInfo->analyticJacobians[index].sparsePattern;
1703 }
1704
1705 /* it's needed to clear the matrix */
1706 SlsSetToZero(Jac);
1707
1708 setContext(data, &currentTime, CONTEXT_JACOBIAN);
1709
1710 /* rescale idaData->y and idaData->yp
1711 * the evaluation of the residual function
1712 * needs to be performed on unscaled values
1713 */
1714 if ((omc_flag[FLAG_IDA_SCALING] && !idaData->disableScaling))
1715 {
1716 idaReScaleVector(rr, idaData->resScale, idaData->N);
1717 idaReScaleData(idaData);
1718 }
1719
1720 for(i = 0; i < sparsePattern->maxColors; i++)
1721 {
1722 for(ii=0; ii < idaData->N; ii++)
1723 {
1724 if(sparsePattern->colorCols[ii]-1 == i)
1725 {
1726 delta_hhh = currentStep * yprime[ii];
1727 delta_hh[ii] = delta_h * fmax(fmax(fabs(states[ii]),fabs(delta_hhh)),fabs(1./errwgt[ii]));
1728 delta_hh[ii] = (delta_hhh >= 0 ? delta_hh[ii] : -delta_hh[ii]);
1729 delta_hh[ii] = (states[ii] + delta_hh[ii]) - states[ii];
1730 ysave[ii] = states[ii];
1731 states[ii] += delta_hh[ii];
1732
1733 if (idaData->daeMode){
1734 ypsave[ii] = yprime[ii];
1735 yprime[ii] += cj * delta_hh[ii];
1736 }
1737
1738 delta_hh[ii] = 1. / delta_hh[ii];
1739 }
1740 }
1741 idaData->disableScaling = 1;
1742 (*idaData->residualFunction)(currentTime, yy, yp, idaData->newdelta, userData);
1743 idaData->disableScaling = disBackup;
1744
1745 increaseJacContext(data);
1746
1747 for(ii = 0; ii < idaData->N; ii++)
1748 {
1749 if(sparsePattern->colorCols[ii]-1 == i)
1750 {
1751 nth = sparsePattern->leadindex[ii];
1752 while(nth < sparsePattern->leadindex[ii+1])
1753 {
1754 j = sparsePattern->index[nth];
1755 //setJacElementKluSparse(j, ii, nth, (newdelta[j] - delta[j]) * delta_hh[ii], Jac, -1);
1756 /* use row scaling for jacobian elements */
1757 if (idaData->disableScaling == 1 || !omc_flag[FLAG_IDA_SCALING]){
1758 setJacElementKluSparse(j, ii, nth, (newdelta[j] - delta[j]) * delta_hh[ii], Jac, -1);
1759 }else{
1760 setJacElementKluSparse(j, ii, nth, ((newdelta[j] - delta[j]) * delta_hh[ii]) / idaData->resScale[j] * idaData->yScale[ii], Jac, -1);
1761 }
1762 nth++;
1763 };
1764 states[ii] = ysave[ii];
1765 if (idaData->daeMode)
1766 {
1767 yprime[ii] = ypsave[ii];
1768 }
1769 }
1770 }
1771 }
1772 finishSparseColPtr(Jac, sparsePattern->numberOfNoneZeros);
1773
1774 /* scale idaData->y and idaData->yp again */
1775 if ((omc_flag[FLAG_IDA_SCALING] && !idaData->disableScaling))
1776 {
1777 idaScaleVector(rr, idaData->resScale, idaData->N);
1778 idaScaleData(idaData);
1779 }
1780
1781 unsetContext(data);
1782 messageClose(LOG_SOLVER_V);
1783
1784 TRACE_POP
1785 return 0;
1786}
1787
1788/*
1789 * This function calculates the jacobian matrix symbolically while exploiting coloring.
1790 */
1791int jacColoredSymbolicalSparse(double currentTime, N_Vector yy, N_Vector yp,
1792 N_Vector rr, SlsMat Jac, double cj,
1793 void *userData)
1794{
1795 TRACE_PUSH
1796 IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
1797 DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
1798 threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)idaData->simData)->threadData);
1799 const int index = data->callback->INDEX_JAC_A;
1800 ANALYTIC_JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[index]);
1801
1802 /* prepare variables */
1803 double *states = N_VGetArrayPointer(yy);
1804 double *yprime = N_VGetArrayPointer(yp);
1805
1806#ifdef USE_PARJAC
1807 ANALYTIC_JACOBIAN* t_jac = (idaData->jacColumns);
1808#else
1809 ANALYTIC_JACOBIAN* t_jac = jac;
1810#endif
1811 unsigned int columns = jac->sizeCols;
1812 unsigned int rows = jac->sizeRows;
1813 SPARSE_PATTERN* sparsePattern = jac->sparsePattern;
1814 int maxColors = sparsePattern->maxColors;
1815
1816 /* it's needed to clear the matrix */
1817 SlsSetToZero(Jac);
1818
1819 setContext(data, &currentTime, CONTEXT_SYM_JACOBIAN); /* Reuse jacobian matrix in KLU solver */
1820
1821 /* Evaluate constant equations if available */
1822 if (jac->constantEqns != NULL((void*)0)) {
1823 jac->constantEqns(data, threadData, jac, NULL((void*)0));
1824 }
1825
1826 genericColoredSymbolicJacobianEvaluation(rows, columns, sparsePattern, Jac, t_jac,
1827 data, threadData, &setJacElementKluSparse);
1828
1829 finishSparseColPtr(Jac, sparsePattern->numberOfNoneZeros);
1830 unsetContext(data);
1831
1832 TRACE_POP
1833 return 0;
1834}
1835
1836/*
1837 * Wrapper function to call numerical or symbolical jacobian matrix
1838 */
1839static int callSparseJacobian(double currentTime, double cj,
1840 N_Vector yy, N_Vector yp, N_Vector rr,
1841 SlsMat Jac, void *user_data,
1842 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
1843{
1844 TRACE_PUSH
1845
1846 IDA_SOLVER* idaData = (IDA_SOLVER*)user_data;
1847 DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
1848 threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)user_data)->simData)->threadData);
1849
1850 /* profiling */
1851 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1852 rt_tick(SIM_TIMER_JACOBIAN5);
1853
1854 if (idaData->jacobianMethod == COLOREDSYMJAC || idaData->jacobianMethod == SYMJAC)
1855 {
1856 jacColoredSymbolicalSparse(currentTime, yy, yp, rr, Jac, cj, user_data);
1857 }
1858 else if (idaData->jacobianMethod == COLOREDNUMJAC || idaData->jacobianMethod == NUMJAC)
1859 {
1860 jacoColoredNumericalSparse(currentTime, yy, yp, rr, Jac, cj, user_data);
1861 }
1862
1863 /* debug */
1864 if (ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC])){
1865 infoStreamPrint(LOG_JAC, 0, "##IDA## Sparse Matrix A.");
1866 PrintSparseMat(Jac);
1867 }
1868
1869 /* add cj to diagonal elements and store in Jac */
1870 if (!idaData->daeMode)
1871 {
1872 int i;
1873 for (i=0; i < idaData->N; ++i){
1874 idaData->tmpJac->colptrs[i] = i;
1875 idaData->tmpJac->rowvals[i] = i;
1876 idaData->tmpJac->data[i] = -cj;
1877 }
1878 idaData->tmpJac->colptrs[idaData->N] = idaData->N;
1879 SlsAddMat(Jac, idaData->tmpJac);
1880 }
1881
1882 /* profiling */
1883 rt_accumulate(SIM_TIMER_JACOBIAN5);
1884 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1885
1886 TRACE_POP
1887 return 0;
1888}
1889
1890static int getScalingFactors(DATA* data, IDA_SOLVER *idaData, SlsMat inScaleMatrix)
1891{
1892 int i;
1893
1894 N_Vector tmp1 = N_VNew_Serial(idaData->N);
1895 N_Vector tmp2 = N_VNew_Serial(idaData->N);
1896 N_Vector tmp3 = N_VNew_Serial(idaData->N);
1897
1898 N_Vector rres = N_VNew_Serial(idaData->N);
1899
1900 /* fill errwgt, since it is needed by the Jacobian calculation function */
1901 double *errwgt = N_VGetArrayPointer(idaData->errwgt);
1902 _omc_fillVector(_omc_createVector(idaData->N, errwgt), 1.);
1903
1904 SlsMat scaleMatrix;
1905
1906
1907 if (inScaleMatrix == NULL((void*)0)){
1908
1909 infoStreamPrint(LOG_SOLVER_V, 1, "##IDA## get new scaling matrix.");
1910
1911 /* use y scale to scale jacobian, but y and yp are not scaled */
1912 idaData->disableScaling = 1;
1913
1914 /* eval residual function first */
1915 residualFunctionIDA(data->localData[0]->timeValue, idaData->y, idaData->yp, rres, (void*) idaData);
1916
1917 /* choose the jacobian sparse vs. dense */
1918 if (idaData->linearSolverMethod == IDA_LS_KLU)
1919 {
1920 scaleMatrix = NewSparseMat(idaData->N, idaData->N, idaData->NNZ);
1921 callSparseJacobian(data->localData[0]->timeValue, 1.0, idaData->y, idaData->yp, rres,
1922 scaleMatrix, idaData, tmp1, tmp2, tmp3);
1923 }
1924 else
1925 {
1926 DlsMat denseMatrix = NewDenseMat(idaData->N, idaData->N);
1927 callDenseJacobian(idaData->N, data->localData[0]->timeValue, 1.0, idaData->y, idaData->yp, rres,
1928 denseMatrix, idaData, tmp1, tmp2, tmp3);
1929 scaleMatrix = SlsConvertDls(denseMatrix);
1930 }
1931 /* enable scaled jacobian again */
1932 idaData->disableScaling = 0;
1933 }
1934 else
1935 {
1936 infoStreamPrint(LOG_SOLVER_V, 1, "##IDA## use given scaling matrix.");
1937 scaleMatrix = inScaleMatrix;
1938 }
1939
1940 /* set resScale factors */
1941 _omc_fillVector(_omc_createVector(idaData->N,idaData->resScale), MINIMAL_SCALE_FACTOR1e-8);
1942 for(i=0; i<scaleMatrix->NNZ; ++i){
1943 if (idaData->resScale[scaleMatrix->rowvals[i]] < fabs(scaleMatrix->data[i])){
1944 idaData->resScale[scaleMatrix->rowvals[i]] = fabs(scaleMatrix->data[i]);
1945 }
1946 }
1947
1948 printVector(LOG_SOLVER_V, "Prime scale factors", idaData->ypScale, idaData->N, 0.0);
1949 printVector(LOG_SOLVER_V, "Residual scale factors", idaData->resScale, idaData->N, 0.0);
1950
1951 messageClose(LOG_SOLVER_V);
1952 N_VDestroy_Serial(tmp1);
1953 N_VDestroy_Serial(tmp2);
1954 N_VDestroy_Serial(tmp3);
1955 N_VDestroy_Serial(rres);
1956
1957 return 0;
1958}
1959
1960static int idaScaleVector(N_Vector vec, double* factors, unsigned int size)
1961{
1962 int i;
1963 double *data = N_VGetArrayPointer(vec);
1964 printVector(LOG_SOLVER_V, "un-scaled", data, size, 0.0);
1965 for(i=0; i < size; ++i)
1966 {
1967 data[i] = data[i] / factors[i];
1968 }
1969 printVector(LOG_SOLVER_V, "scaled", data, size, 0.0);
1970 return 0;
1971}
1972
1973static int idaReScaleVector(N_Vector vec, double* factors, unsigned int size)
1974{
1975 int i;
1976 double *data = N_VGetArrayPointer(vec);
1977
1978 printVector(LOG_SOLVER_V, "scaled", data, size, 0.0);
1979 for(i=0; i < size; ++i)
1980 {
1981 data[i] = data[i] * factors[i];
1982 }
1983 printVector(LOG_SOLVER_V, "un-scaled", data, size, 0.0);
1984 return 0;
1985}
1986
1987static int idaScaleData(IDA_SOLVER *idaData)
1988{
1989 infoStreamPrint(LOG_SOLVER_V, 1, "Scale y");
1990 idaScaleVector(idaData->y, idaData->yScale, idaData->N);
1991 messageClose(LOG_SOLVER_V);
1992 infoStreamPrint(LOG_SOLVER_V, 1, "Scale yp");
1993 idaScaleVector(idaData->yp, idaData->ypScale, idaData->N);
1994 messageClose(LOG_SOLVER_V);
1995
1996 return 0;
1997}
1998
1999static int idaReScaleData(IDA_SOLVER *idaData)
2000{
2001 infoStreamPrint(LOG_SOLVER_V, 1, "Re-Scale y");
2002 idaReScaleVector(idaData->y, idaData->yScale, idaData->N);
2003 messageClose(LOG_SOLVER_V);
2004 infoStreamPrint(LOG_SOLVER_V, 1, "Re-Scale yp");
2005 idaReScaleVector(idaData->yp, idaData->ypScale, idaData->N);
2006 messageClose(LOG_SOLVER_V);
2007
2008 return 0;
2009}
2010
2011
2012#endif