Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/dassl.c
Warning:line 173, column 19
Value stored to 'tmpSimData' during its initialization 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-2014, 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#ifdef USE_PARJAC
31 #include <omp.h>
32 #define GC_THREADS
33 #include <gc/omc_gc.h>
34#endif
35
36#include <string.h>
37#include <setjmp.h>
38
39#include "openmodelica.h"
40#include "openmodelica_func.h"
41#include "simulation_data.h"
42
43#include "util/omc_error.h"
44#include "util/parallel_helper.h"
45#include "gc/omc_gc.h"
46
47#include "simulation/options.h"
48#include "simulation/simulation_runtime.h"
49#include "simulation/results/simulation_result.h"
50#include "simulation/solver/solver_main.h"
51#include "simulation/solver/model_help.h"
52#include "simulation/solver/external_input.h"
53#include "simulation/solver/epsilon.h"
54#include "simulation/solver/omc_math.h"
55#include "simulation/solver/jacobianSymbolical.h"
56#include "simulation/solver/dassl.h"
57#include "meta/meta_modelica.h"
58
59#ifdef __cplusplus
60extern "C" {
61#endif
62
63/* experimental flag for SKF TLM Master Solver Interface
64 * - it's used with -noEquidistantTimeGrid flag.
65 * - it's set to 1 if the continuous system is evaluated
66 * when dassl finished a step, otherwise it's 0.
67 */
68int RHSFinalFlag;
69
70/* provides a dummy Jacobian to be used with DASSL */
71static int dummy_Jacobian(double *t, double *y, double *yprime,double *deltaD,
72 double *delta, double *cj, double *h, double *wt,
73 double *rpar, int* ipar) {
74 return 0;
75}
76
77/* provides a dummy zero crossing function to be used with DASSL */
78static int dummy_zeroCrossing(int *neqm, double *t, double *y, double *yp,
79 int *ng, double *gout, double *rpar, int* ipar) {
80 return 0;
81}
82
83/* provides a dumm precondition function to be used with DASSL */
84static int dummy_precondition(int *neq, double *t, double *y, double *yprime,
85 double *savr, double *pwk, double *cj,
86 double *wt, double *wp, int *iwp, double *b,
87 double eplin, int* ires, double *rpar, int* ipar){
88 return 0;
89}
90
91/* Function prototypes */
92static int callJacobian(double *t, double *y, double *yprime, double *deltaD,
93 double *pd, double *cj, double *h, double *wt,
94 double *rpar, int* ipar);
95
96static int jacA_num(double *t, double *y, double *yprime, double *deltaD,
97 double *pd, double *cj, double *h, double *wt,
98 double *rpar, int* ipar);
99
100static int jacA_numColored(double *t, double *y, double *yprime,
101 double *deltaD, double *pd, double *cj, double *h,
102 double *wt, double *rpar, int* ipar);
103
104static int jacA_sym(double *t, double *y, double *yprime, double *deltaD,
105 double *pd, double *cj, double *h, double *wt,
106 double *rpar, int* ipar);
107
108static int jacA_symColored(double *t, double *y, double *yprime,
109 double *deltaD, double *pd, double *cj, double *h,
110 double *wt, double *rpar, int* ipar);
111
112static void setJacElementDasslSparse(int l, int k, int nth, double val,
113 void* matrixA, int rows);
114
115void DDASKR_daskr_ddaskr_(
116 int (*res) (double *t, double *y, double *yprime, double* cj, double *delta, int *ires, double *rpar, int* ipar),
117 int *neq,
118 double *t,
119 double *y,
120 double *yprime,
121 double *tout,
122 int *info,
123 double *rtol,
124 double *atol,
125 int *idid,
126 double *rwork,
127 int *lrw,
128 int *iwork,
129 int *liw,
130 double *rpar,
131 int *ipar,
132 int (*jac) (double *t, double *y, double *yprime, double *deltaD, double *delta, double *cj, double *h, double *wt, double *rpar, int* ipar),
133 int (*psol) (int *neq, double *t, double *y, double *yprime, double *savr, double *pwk, double *cj, double *wt, double *wp, int *iwp, double *b, double eplin, int* ires, double *rpar, int* ipar),
134 int (*g) (int *neqm, double *t, double *y, double *yp, int *ng, double *gout, double *rpar, int* ipar),
135 int *ng,
136 int *jroot
137);
138
139static int continue_DASSL(int* idid, double* tolarence);
140
141/* function for calculating state values on residual form */
142static int functionODE_residual(double *t, double *x, double *xprime,
143 double *cj, double *delta, int *ires,
144 double *rpar, int* ipar);
145
146/* function for calculating zeroCrossings */
147static int function_ZeroCrossingsDASSL(int *neqm, double *t, double *y,
148 double *yp, int *ng, double *gout,
149 double *rpar, int* ipar);
150
151
152/*
153 * \brief Configure DASSL solver
154 *
155 * Allocate memory for intern data of `dasslData`.
156 * Configures DASSL:
157 * - Set relative and absolute tolerance
158 * - Set maximum step size, initial step size
159 * - Set maximum integration order
160 * - Set time grid
161 * - Set method for jacobian computation
162 * - Set root finding method
163 * - Set event handling and restart option
164 *
165 */
166int dassl_initial(DATA* data, threadData_t *threadData,
167 SOLVER_INFO* solverInfo, DASSL_DATA *dasslData)
168{
169 TRACE_PUSH
170 /* work arrays for DASSL */
171 unsigned int i;
172 long N;
173 SIMULATION_DATA tmpSimData = {0};
Value stored to 'tmpSimData' during its initialization is never read
174
175 dasslData->residualFunction = functionODE_residual;
176 N = data->modelData->nStates;
177
178 dasslData->N = N;
179
180 RHSFinalFlag = 0;
181
182 dasslData->liw = 40 + N;
183 dasslData->lrw = 60 + ((maxOrder + 4) * N) + (N * N) + (3*data->modelData->nZeroCrossings);
184 dasslData->rwork = (double*) calloc(dasslData->lrw, sizeof(double));
185 assertStreamPrint(threadData, 0 != dasslData->rwork,"out of memory")if (!(0 != dasslData->rwork)) {throwStreamPrint((threadData
), "out of memory"); ((void) sizeof ((0) ? 1 : 0), __extension__
({ if (0) ; else __assert_fail ("0", "./simulation/solver/dassl.c"
, 185, __extension__ __PRETTY_FUNCTION__); }));}
;
186 dasslData->iwork = (int*) calloc(dasslData->liw, sizeof(int));
187 assertStreamPrint(threadData, 0 != dasslData->iwork,"out of memory")if (!(0 != dasslData->iwork)) {throwStreamPrint((threadData
), "out of memory"); ((void) sizeof ((0) ? 1 : 0), __extension__
({ if (0) ; else __assert_fail ("0", "./simulation/solver/dassl.c"
, 187, __extension__ __PRETTY_FUNCTION__); }));}
;
188 dasslData->ng = (int) data->modelData->nZeroCrossings;
189 dasslData->jroot = (int*) calloc(data->modelData->nZeroCrossings, sizeof(int));
190 dasslData->rpar = (double**) malloc(3*sizeof(double*));
191 dasslData->ipar = (int*) malloc(sizeof(int));
192 dasslData->ipar[0] = ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC]);
193 assertStreamPrint(threadData, 0 != dasslData->ipar,"out of memory")if (!(0 != dasslData->ipar)) {throwStreamPrint((threadData
), "out of memory"); ((void) sizeof ((0) ? 1 : 0), __extension__
({ if (0) ; else __assert_fail ("0", "./simulation/solver/dassl.c"
, 193, __extension__ __PRETTY_FUNCTION__); }));}
;
194 dasslData->atol = (double*) malloc(N*sizeof(double));
195 dasslData->rtol = (double*) malloc(N*sizeof(double));
196 dasslData->info = (int*) calloc(infoLength, sizeof(int));
197 assertStreamPrint(threadData, 0 != dasslData->info,"out of memory")if (!(0 != dasslData->info)) {throwStreamPrint((threadData
), "out of memory"); ((void) sizeof ((0) ? 1 : 0), __extension__
({ if (0) ; else __assert_fail ("0", "./simulation/solver/dassl.c"
, 197, __extension__ __PRETTY_FUNCTION__); }));}
;
198
199 dasslData->idid = 0;
200
201 dasslData->ysave = (double*) malloc(N*sizeof(double));
202 dasslData->delta_hh = (double*) malloc(N*sizeof(double));
203 dasslData->newdelta = (double*) malloc(N*sizeof(double));
204 dasslData->stateDer = (double*) calloc(N, sizeof(double));
205 dasslData->states = (double*) malloc(N*sizeof(double));
206 dasslData->allocatedParMem = 0; /* false */
207
208 data->simulationInfo->currentContext = CONTEXT_ALGEBRAIC;
209
210 /* ### start configuration of dassl ### */
211 infoStreamPrint(LOG_SOLVER, 1, "Configuration of the dassl code:");
212
213
214
215 /* set nominal values of the states for absolute tolerances */
216 dasslData->info[1] = 1;
217 infoStreamPrint(LOG_SOLVER, 1, "The relative tolerance is %g. Following absolute tolerances are used for the states: ", data->simulationInfo->tolerance);
218 for(i=0; i<dasslData->N; ++i)
219 {
220 dasslData->rtol[i] = data->simulationInfo->tolerance;
221 dasslData->atol[i] = data->simulationInfo->tolerance * fmax(fabs(data->modelData->realVarsData[i].attribute.nominal), 1e-32);
222 infoStreamPrint(LOG_SOLVER_V, 0, "%d. %s -> %g", i+1, data->modelData->realVarsData[i].info.name, dasslData->atol[i]);
223 }
224 messageClose(LOG_SOLVER);
225
226
227
228 /* let dassl return at every internal step */
229 dasslData->info[2] = 1;
230
231
232 /* define maximum step size, which is dassl is allowed to go */
233 if (omc_flag[FLAG_MAX_STEP_SIZE])
234 {
235 double maxStepSize = atof(omc_flagValue[FLAG_MAX_STEP_SIZE]);
236
237 assertStreamPrint(threadData, maxStepSize >= DASSL_STEP_EPS, "Selected maximum step size %e is too small.", maxStepSize)if (!(maxStepSize >= DASSL_STEP_EPS)) {throwStreamPrint((threadData
), "Selected maximum step size %e is too small.", maxStepSize
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "./simulation/solver/dassl.c", 237, __extension__
__PRETTY_FUNCTION__); }));}
;
238
239 dasslData->rwork[1] = maxStepSize;
240 dasslData->info[6] = 1;
241 infoStreamPrint(LOG_SOLVER, 0, "maximum step size %g", dasslData->rwork[1]);
242 }
243 else
244 {
245 infoStreamPrint(LOG_SOLVER, 0, "maximum step size not set");
246 }
247
248
249 /* define initial step size, which is dassl is used every time it restarts */
250 if (omc_flag[FLAG_INITIAL_STEP_SIZE])
251 {
252 double initialStepSize = atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]);
253
254 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/dassl.c"
, 254, __extension__ __PRETTY_FUNCTION__); }));}
;
255
256 dasslData->rwork[2] = initialStepSize;
257 dasslData->info[7] = 1;
258 infoStreamPrint(LOG_SOLVER, 0, "initial step size %g", dasslData->rwork[2]);
259 }
260 else
261 {
262 infoStreamPrint(LOG_SOLVER, 0, "initial step size not set");
263 }
264
265
266 /* define maximum integration order of dassl */
267 if (omc_flag[FLAG_MAX_ORDER])
268 {
269 int maxOrder = atoi(omc_flagValue[FLAG_MAX_ORDER]);
270
271 assertStreamPrint(threadData, maxOrder >= 1 && maxOrder <= 5, "Selected maximum order %d is out of range (1-5).", maxOrder)if (!(maxOrder >= 1 && maxOrder <= 5)) {throwStreamPrint
((threadData), "Selected maximum order %d is out of range (1-5)."
, maxOrder); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if
(0) ; else __assert_fail ("0", "./simulation/solver/dassl.c"
, 271, __extension__ __PRETTY_FUNCTION__); }));}
;
272
273 dasslData->iwork[2] = maxOrder;
274 dasslData->info[8] = 1;
275 }
276 infoStreamPrint(LOG_SOLVER, 0, "maximum integration order %d", dasslData->info[8]?dasslData->iwork[2]:maxOrder);
277
278
279 /* if FLAG_NOEQUIDISTANT_GRID is set, choose dassl step method */
280 if (omc_flag[FLAG_NOEQUIDISTANT_GRID])
281 {
282 dasslData->dasslSteps = 1; /* TRUE */
283 solverInfo->solverNoEquidistantGrid = 1;
284 }
285 else
286 {
287 dasslData->dasslSteps = 0; /* FALSE */
288 }
289 infoStreamPrint(LOG_SOLVER, 0, "use equidistant time grid %s", dasslData->dasslSteps?"NO":"YES");
290
291 /* check if Flags FLAG_NOEQUIDISTANT_OUT_FREQ or FLAG_NOEQUIDISTANT_OUT_TIME are set */
292 if (dasslData->dasslSteps){
293 if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ])
294 {
295 dasslData->dasslStepsFreq = atoi(omc_flagValue[FLAG_NOEQUIDISTANT_OUT_FREQ]);
296 }
297 else if (omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME])
298 {
299 dasslData->dasslStepsTime = atof(omc_flagValue[FLAG_NOEQUIDISTANT_OUT_TIME]);
300 dasslData->rwork[1] = dasslData->dasslStepsTime;
301 dasslData->info[6] = 1;
302 infoStreamPrint(LOG_SOLVER, 0, "maximum step size %g", dasslData->rwork[1]);
303 } else {
304 dasslData->dasslStepsFreq = 1;
305 dasslData->dasslStepsTime = 0.0;
306 }
307
308 if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ] && omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]){
309 warningStreamPrint(LOG_STDOUT, 0, "The flags are \"noEquidistantOutputFrequency\" "
310 "and \"noEquidistantOutputTime\" are in opposition "
311 "to each other. The flag \"noEquidistantOutputFrequency\" superiors.");
312 }
313 infoStreamPrint(LOG_SOLVER, 0, "as the output frequency control is used: %d", dasslData->dasslStepsFreq);
314 infoStreamPrint(LOG_SOLVER, 0, "as the output frequency time step control is used: %f", dasslData->dasslStepsTime);
315 }
316
317 /* if FLAG_JACOBIAN is set, choose dassl jacobian calculation method */
318 if (omc_flag[FLAG_JACOBIAN])
319 {
320 for(i=1; i< JAC_MAX;i++)
321 {
322 if(!strcmp((const char*)omc_flagValue[FLAG_JACOBIAN], JACOBIAN_METHOD[i])){
323 dasslData->dasslJacobian = (int)i;
324 break;
325 }
326 }
327 if(dasslData->dasslJacobian == JAC_UNKNOWN)
328 {
329 if (ACTIVE_WARNING_STREAM(LOG_SOLVER)(showAllWarnings || useStream[LOG_SOLVER]))
330 {
331 warningStreamPrint(LOG_SOLVER, 1, "unrecognized jacobian calculation method %s, current options are:", (const char*)omc_flagValue[FLAG_JACOBIAN]);
332 for(i=1; i < JAC_MAX; ++i)
333 {
334 warningStreamPrint(LOG_SOLVER, 0, "%-15s [%s]", JACOBIAN_METHOD[i], JACOBIAN_METHOD_DESC[i]);
335 }
336 messageClose(LOG_SOLVER);
337 }
338 throwStreamPrint(threadData,"unrecognized jacobian calculation method %s", (const char*)omc_flagValue[FLAG_JACOBIAN]);
339 }
340 /* default case colored numerical jacobian */
341 }
342 else
343 {
344 dasslData->dasslJacobian = COLOREDNUMJAC;
345 }
346
347 /* selects the calculation method of the jacobian */
348 if(dasslData->dasslJacobian == COLOREDNUMJAC ||
349 dasslData->dasslJacobian == COLOREDSYMJAC ||
350 dasslData->dasslJacobian == SYMJAC)
351 {
352 ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
353 if (data->callback->initialAnalyticJacobianA(data, threadData, jacobian))
354 {
355 infoStreamPrint(LOG_STDOUT, 0, "Jacobian or SparsePattern is not generated or failed to initialize! Switch back to normal.");
356 dasslData->dasslJacobian = INTERNALNUMJAC;
357 } else {
358 ANALYTIC_JACOBIAN* jac = &data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A];
359 infoStreamPrint(LOG_SIMULATION, 1, "Initialized colored Jacobian:");
360 infoStreamPrint(LOG_SIMULATION, 0, "columns: %d rows: %d", jac->sizeCols, jac->sizeRows);
361 infoStreamPrint(LOG_SIMULATION, 0, "NNZ: %d colors: %d", jac->sparsePattern->numberOfNoneZeros, jac->sparsePattern->maxColors);
362 messageClose(LOG_SIMULATION);
363 }
364 }
365 /* default use a user sub-routine for JAC */
366 dasslData->info[4] = 1;
367
368 /* set up the appropriate function pointer */
369 switch (dasslData->dasslJacobian){
370 case COLOREDNUMJAC:
371 data->simulationInfo->jacobianEvals = data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern->maxColors;
372 dasslData->jacobianFunction = jacA_numColored;
373 break;
374 case COLOREDSYMJAC:
375 data->simulationInfo->jacobianEvals = data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern->maxColors;
376 dasslData->jacobianFunction = jacA_symColored;
377#ifdef USE_PARJAC
378 allocateThreadLocalJacobians(data, &(dasslData->jacColumns));
379 dasslData->allocatedParMem = 1; /* true */
380#endif
381 break;
382 case SYMJAC:
383 dasslData->jacobianFunction = jacA_sym;
384#ifdef USE_PARJAC
385 allocateThreadLocalJacobians(data, &(dasslData->jacColumns));
386 dasslData->allocatedParMem = 1; /* true */
387#endif
388 break;
389 case NUMJAC:
390 dasslData->jacobianFunction = jacA_num;
391 break;
392 case INTERNALNUMJAC:
393 dasslData->jacobianFunction = dummy_Jacobian;
394 /* no user sub-routine for JAC */
395 dasslData->info[4] = 0;
396 break;
397 default:
398 throwStreamPrint(threadData,"unrecognized jacobian calculation method %s", (const char*)omc_flagValue[FLAG_JACOBIAN]);
399 break;
400 }
401 infoStreamPrint(LOG_SOLVER, 0, "jacobian is calculated by %s", JACOBIAN_METHOD_DESC[dasslData->dasslJacobian]);
402
403 /* if FLAG_NO_ROOTFINDING is set, choose dassl with out internal root finding */
404 if(omc_flag[FLAG_NO_ROOTFINDING])
405 {
406 dasslData->dasslRootFinding = 0;
407 dasslData->zeroCrossingFunction = dummy_zeroCrossing;
408 dasslData->ng = 0;
409 }
410 else
411 {
412 solverInfo->solverRootFinding = 1;
413 dasslData->dasslRootFinding = 1;
414 dasslData->zeroCrossingFunction = function_ZeroCrossingsDASSL;
415 }
416 infoStreamPrint(LOG_SOLVER, 0, "dassl uses internal root finding method %s", dasslData->dasslRootFinding?"YES":"NO");
417
418
419 /* if FLAG_NO_RESTART is set, choose dassl step method */
420 if (omc_flag[FLAG_NO_RESTART])
421 {
422 dasslData->dasslAvoidEventRestart = 1; /* TRUE */
423 }
424 else
425 {
426 dasslData->dasslAvoidEventRestart = 0; /* FALSE */
427 }
428 infoStreamPrint(LOG_SOLVER, 0, "dassl performs an restart after an event occurs %s", dasslData->dasslAvoidEventRestart?"NO":"YES");
429
430 /* ### end configuration of dassl ### */
431
432 messageClose(LOG_SOLVER);
433 TRACE_POP
434 return 0;
435}
436
437
438/*
439 * \brief Deallocates `DASSL_DATA`
440 */
441int dassl_deinitial(DASSL_DATA *dasslData)
442{
443 TRACE_PUSH
444 unsigned int i;
445
446 /* free work arrays for DASSL */
447 free(dasslData->rwork);
448 free(dasslData->iwork);
449 free(dasslData->jroot);
450 free(dasslData->rpar);
451 free(dasslData->ipar);
452 free(dasslData->atol);
453 free(dasslData->rtol);
454 free(dasslData->info);
455 free(dasslData->ysave);
456 free(dasslData->delta_hh);
457 free(dasslData->newdelta);
458 free(dasslData->stateDer);
459 free(dasslData->states);
460
461
462#ifdef USE_PARJAC
463 if (dasslData->allocatedParMem) {
464 freeAnalyticalJacobian(&(dasslData->jacColumns));
465 dasslData->allocatedParMem = 0;
466 }
467#endif
468
469 free(dasslData);
470
471 TRACE_POP
472 return 0;
473}
474
475/* \fn printCurrentStatesVector(int logLevel, double* y, DATA* data, double time)
476 *
477 * \param [in] [logLevel]
478 * \param [in] [states]
479 * \param [in] [data]
480 * \param [in] [time]
481 *
482 * This function outputs states vector.
483 *
484 */
485int printCurrentStatesVector(int logLevel, double* states, DATA* data, double time)
486{
487 int i;
488 infoStreamPrint(logLevel, 1, "states at time=%g", time);
489 for(i=0;i<data->modelData->nStates;++i)
490 {
491 infoStreamPrint(logLevel, 0, "%d. %s = %g", i+1, data->modelData->realVarsData[i].info.name, states[i]);
492 }
493 messageClose(logLevel);
494
495 return 0;
496}
497
498/* \fn printVector(int logLevel, double* y, DATA* data, double time)
499 *
500 * \param [in] [logLevel]
501 * \param [in] [name]
502 * \param [in] [vec]
503 * \param [in] [size]
504 * \param [in] [time]
505 *
506 * This function outputs a vector of size
507 *
508 */
509int printVector(int logLevel, const char* name, double* vec, int n, double time)
510{
511 int i;
512 infoStreamPrint(logLevel, 1, "%s at time=%g", name, time);
513 for(i=0; i<n; ++i)
514 {
515 infoStreamPrint(logLevel, 0, "%d. %g", i+1, vec[i]);
516 }
517 messageClose(logLevel);
518
519 return 0;
520}
521
522
523/**********************************************************************************************
524 * DASSL with synchronous treating of when equation
525 * - without integrated ZeroCrossing method.
526 * + ZeroCrossing are handled outside DASSL.
527 * + if no event occurs outside DASSL performs a warm-start
528 **********************************************************************************************/
529int dassl_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
530{
531 TRACE_PUSH
532 double tout = 0;
533 int i = 0;
534 unsigned int ui = 0;
535 int retVal = 0;
536 int saveJumpState;
537 static unsigned int dasslStepsOutputCounter = 1;
538
539 DASSL_DATA *dasslData = (DASSL_DATA*) solverInfo->solverData;
540
541 SIMULATION_DATA *sData = data->localData[0];
542 SIMULATION_DATA *sDataOld = data->localData[1];
543
544 modelica_real* states = sData->realVars;
545 modelica_real* stateDer = dasslData->stateDer;
546
547
548 MODEL_DATA *mData = (MODEL_DATA*) data->modelData;
549
550 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
551
552 memcpy(stateDer, data->localData[1]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
553
554 dasslData->rpar[0] = (double*) (void*) data;
555 dasslData->rpar[1] = (double*) (void*) dasslData;
556 dasslData->rpar[2] = (double*) (void*) threadData;
557
558 saveJumpState = threadData->currentErrorStage;
559 threadData->currentErrorStage = ERROR_INTEGRATOR;
560
561 /* try */
562#if !defined(OMC_EMCC)
563 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) {
564#endif
565
566 assertStreamPrint(threadData, 0 != dasslData->rpar, "could not passed to DDASKR")if (!(0 != dasslData->rpar)) {throwStreamPrint((threadData
), "could not passed to DDASKR"); ((void) sizeof ((0) ? 1 : 0
), __extension__ ({ if (0) ; else __assert_fail ("0", "./simulation/solver/dassl.c"
, 566, __extension__ __PRETTY_FUNCTION__); }));}
;
567
568 /* If an event is triggered and processed restart dassl. */
569 if(!dasslData->dasslAvoidEventRestart && (solverInfo->didEventStep || 0 == dasslData->idid))
570 {
571 debugStreamPrint(LOG_EVENTS_V, 0, "Event-management forced reset of DDASKR");
572 /* obtain reset */
573 dasslData->info[0] = 0;
574 dasslData->idid = 0;
575
576 }
577
578 /* Calculate steps until TOUT is reached */
579 if (dasslData->dasslSteps)
580 {
581 /* If dasslsteps is selected, the dassl run to stopTime or next sample event */
582 if (data->simulationInfo->nextSampleEvent < data->simulationInfo->stopTime)
583 {
584 tout = data->simulationInfo->nextSampleEvent;
585 }
586 else
587 {
588 tout = data->simulationInfo->stopTime;
589 }
590 }
591 else
592 {
593 tout = solverInfo->currentTime + solverInfo->currentStepSize;
594 }
595
596 /* Check that tout is not less than timeValue
597 * else will dassl get in trouble. If that is the case we skip the current step. */
598 if (solverInfo->currentStepSize < DASSL_STEP_EPS)
599 {
600 infoStreamPrint(LOG_DASSL, 0, "Desired step to small try next one");
601 infoStreamPrint(LOG_DASSL, 0, "Interpolate linear");
602
603 /*euler step*/
604 for(i = 0; i < data->modelData->nStates; i++)
605 {
606 sData->realVars[i] = sDataOld->realVars[i] + stateDer[i] * solverInfo->currentStepSize;
607 }
608 sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
609 data->callback->functionODE(data, threadData);
610 solverInfo->currentTime = sData->timeValue;
611
612 TRACE_POP
613 return 0;
614 }
615
616 do
617 {
618 infoStreamPrint(LOG_DASSL, 1, "new step at time = %.15g", solverInfo->currentTime);
619
620 /* rhs final flag is FALSE during for dassl evaluation */
621 RHSFinalFlag = 0;
622
623 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
624 /* read input vars */
625 externalInputUpdate(data);
626 data->callback->input_function(data, threadData);
627 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
628
629 DDASKR_daskr_ddaskr_(dasslData->residualFunction, (int*) &dasslData->N,
630 &solverInfo->currentTime, states, stateDer, &tout,
631 dasslData->info, dasslData->rtol, dasslData->atol, &dasslData->idid,
632 dasslData->rwork, &dasslData->lrw, dasslData->iwork, &dasslData->liw,
633 (double*) (void*) dasslData->rpar, dasslData->ipar, callJacobian, dummy_precondition,
634 dasslData->zeroCrossingFunction, (int*) &dasslData->ng, dasslData->jroot);
635
636 /* closing new step message */
637 messageClose(LOG_DASSL);
638
639 /* set ringbuffer time to current time */
640 sData->timeValue = solverInfo->currentTime;
641
642 /* rhs final flag is TRUE during for output evaluation */
643 RHSFinalFlag = 1;
644
645 if(dasslData->idid == -1)
646 {
647 fflush(stderrstderr);
648 fflush(stdoutstdout);
649 warningStreamPrint(LOG_DASSL, 0, "A large amount of work has been expended.(About 500 steps). Trying to continue ...");
650 infoStreamPrint(LOG_DASSL, 0, "DASSL will try again...");
651 dasslData->info[0] = 1; /* try again */
652 if (solverInfo->currentTime <= data->simulationInfo->stopTime)
653 continue;
654 }
655 else if(dasslData->idid < 0)
656 {
657 fflush(stderrstderr);
658 fflush(stdoutstdout);
659 retVal = continue_DASSL(&dasslData->idid, &data->simulationInfo->tolerance);
660 warningStreamPrint(LOG_STDOUT, 0, "can't continue. time = %f", sData->timeValue);
661 TRACE_POP
662 break;
663 }
664 else if(dasslData->idid == 5)
665 {
666 threadData->currentErrorStage = ERROR_EVENTSEARCH;
667 }
668
669 /* emit step, if dasslsteps is selected */
670 if (dasslData->dasslSteps)
671 {
672 if (omc_flag[FLAG_NOEQUIDISTANT_OUT_FREQ]){
673 /* output every n-th time step */
674 if (dasslStepsOutputCounter >= dasslData->dasslStepsFreq){
675 dasslStepsOutputCounter = 1; /* next line set it to one */
676 break;
677 }
678 dasslStepsOutputCounter++;
679 } else if (omc_flag[FLAG_NOEQUIDISTANT_OUT_TIME]){
680 /* output when time>=k*timeValue */
681 if (solverInfo->currentTime > dasslStepsOutputCounter * dasslData->dasslStepsTime){
682 dasslStepsOutputCounter++;
683 break;
684 }
685 } else {
686 break;
687 }
688 }
689
690 } while(dasslData->idid == 1);
691
692 states = dasslData->states;
693
694#if !defined(OMC_EMCC)
695 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
696#endif
697 threadData->currentErrorStage = saveJumpState;
698
699 /* if a state event occurs than no sample event does need to be activated */
700 if (data->simulationInfo->sampleActivated && solverInfo->currentTime < data->simulationInfo->nextSampleEvent)
701 {
702 data->simulationInfo->sampleActivated = 0;
703 }
704
705
706 if(ACTIVE_STREAM(LOG_DASSL)(useStream[LOG_DASSL]))
707 {
708 infoStreamPrint(LOG_DASSL, 1, "dassl call statistics: ");
709 infoStreamPrint(LOG_DASSL, 0, "value of idid: %d", (int)dasslData->idid);
710 infoStreamPrint(LOG_DASSL, 0, "current time value: %0.4g", solverInfo->currentTime);
711 infoStreamPrint(LOG_DASSL, 0, "current integration time value: %0.4g", dasslData->rwork[3]);
712 infoStreamPrint(LOG_DASSL, 0, "step size H to be attempted on next step: %0.4g", dasslData->rwork[2]);
713 infoStreamPrint(LOG_DASSL, 0, "step size used on last successful step: %0.4g", dasslData->rwork[6]);
714 infoStreamPrint(LOG_DASSL, 0, "the order of the method used on the last step: %d", dasslData->iwork[7]);
715 infoStreamPrint(LOG_DASSL, 0, "the order of the method to be attempted on the next step: %d", dasslData->iwork[8]);
716 infoStreamPrint(LOG_DASSL, 0, "number of steps taken so far: %d", (int)dasslData->iwork[10]);
717 infoStreamPrint(LOG_DASSL, 0, "number of calls of functionODE() : %d", (int)dasslData->iwork[11]);
718 infoStreamPrint(LOG_DASSL, 0, "number of calculation of jacobian : %d", (int)dasslData->iwork[12]);
719 infoStreamPrint(LOG_DASSL, 0, "total number of convergence test failures: %d", (int)dasslData->iwork[13]);
720 infoStreamPrint(LOG_DASSL, 0, "total number of error test failures: %d", (int)dasslData->iwork[14]);
721 messageClose(LOG_DASSL);
722 }
723
724 /* save dassl stats */
725 for(ui = 0; ui < numStatistics; ui++)
726 {
727 assert(10 + ui < dasslData->liw)((void) sizeof ((10 + ui < dasslData->liw) ? 1 : 0), __extension__
({ if (10 + ui < dasslData->liw) ; else __assert_fail (
"10 + ui < dasslData->liw", "./simulation/solver/dassl.c"
, 727, __extension__ __PRETTY_FUNCTION__); }))
;
728 solverInfo->solverStatsTmp[ui] = dasslData->iwork[10 + ui];
729 }
730
731 infoStreamPrint(LOG_DASSL, 0, "Finished DASSL step.");
732 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
733
734 TRACE_POP
735 return retVal;
736}
737
738static int continue_DASSL(int* idid, double* atol)
739{
740 TRACE_PUSH
741 int retValue = -1;
742
743 switch(*idid)
744 {
745 case 1:
746 case 2:
747 case 3:
748 /* 1-4 means success */
749 break;
750 case -1:
751 warningStreamPrint(LOG_DASSL, 0, "A large amount of work has been expended.(About 500 steps). Trying to continue ...");
752 retValue = 1; /* adrpo: try to continue */
753 break;
754 case -2:
755 warningStreamPrint(LOG_STDOUT, 0, "The error tolerances are too stringent");
756 retValue = -2;
757 break;
758 case -3:
759 /* wbraun: don't throw at this point let the solver handle it */
760 /* throwStreamPrint("DDASKR: THE LAST STEP TERMINATED WITH A NEGATIVE IDID value"); */
761 retValue = -3;
762 break;
763 case -6:
764 warningStreamPrint(LOG_STDOUT, 0, "DDASSL had repeated error test failures on the last attempted step.");
765 retValue = -6;
766 break;
767 case -7:
768 warningStreamPrint(LOG_STDOUT, 0, "The corrector could not converge.");
769 retValue = -7;
770 break;
771 case -8:
772 warningStreamPrint(LOG_STDOUT, 0, "The matrix of partial derivatives is singular.");
773 retValue = -8;
774 break;
775 case -9:
776 warningStreamPrint(LOG_STDOUT, 0, "The corrector could not converge. There were repeated error test failures in this step.");
777 retValue = -9;
778 break;
779 case -10:
780 warningStreamPrint(LOG_STDOUT, 0, "A Modelica assert prevents the integrator to continue. For more information use -lv LOG_SOLVER");
781 retValue = -10;
782 break;
783 case -11:
784 warningStreamPrint(LOG_STDOUT, 0, "IRES equal to -2 was encountered and control is being returned to the calling program.");
785 retValue = -11;
786 break;
787 case -12:
788 warningStreamPrint(LOG_STDOUT, 0, "DDASSL failed to compute the initial YPRIME.");
789 retValue = -12;
790 break;
791 case -33:
792 warningStreamPrint(LOG_STDOUT, 0, "The code has encountered trouble from which it cannot recover.");
793 retValue = -33;
794 break;
795 }
796
797 TRACE_POP
798 return retValue;
799}
800
801int functionODE_residual(double *t, double *y, double *yd, double* cj,
802 double *delta, int *ires, double *rpar, int *ipar)
803{
804 TRACE_PUSH
805 DATA* data = (DATA*)((double**)rpar)[0];
806 DASSL_DATA* dasslData = (DASSL_DATA*)((double**)rpar)[1];
807 threadData_t *threadData = (threadData_t*)((double**)rpar)[2];
808
809 double timeBackup;
810 long i;
811 int saveJumpState;
812 int success = 0;
813
814 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
815 if (measure_time_flag) rt_tick(SIM_TIMER_RESIDUALS9);
816
817 if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
818 {
819 setContext(data, t, CONTEXT_ODE);
820 }
821 printCurrentStatesVector(LOG_DASSL_STATES, y, data, *t);
822 printVector(LOG_DASSL_STATES, "yd", yd, data->modelData->nStates, *t);
823
824 timeBackup = data->localData[0]->timeValue;
825 data->localData[0]->timeValue = *t;
826
827 saveJumpState = threadData->currentErrorStage;
828 threadData->currentErrorStage = ERROR_INTEGRATOR;
829
830 /* try */
831#if !defined(OMC_EMCC)
832 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) {
833#endif
834
835 /* read input vars */
836 externalInputUpdate(data);
837 data->callback->input_function(data, threadData);
838
839 /* eval input vars */
840 data->callback->functionODE(data, threadData);
841
842 /* get the difference between the temp_xd(=localData->statesDerivatives)
843 and xd(=statesDerivativesBackup) */
844 for(i=0; i < data->modelData->nStates; i++)
845 {
846 delta[i] = data->localData[0]->realVars[data->modelData->nStates + i] - yd[i];
847 }
848 printVector(LOG_DASSL_STATES, "dd", delta, data->modelData->nStates, *t);
849 success = 1;
850#if !defined(OMC_EMCC)
851 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
852#endif
853
854 if (!success) {
855 *ires = -1;
856 }
857
858 threadData->currentErrorStage = saveJumpState;
859
860 data->localData[0]->timeValue = timeBackup;
861
862 if (data->simulationInfo->currentContext == CONTEXT_ODE){
863 unsetContext(data);
864 }
865
866 if (measure_time_flag) rt_accumulate(SIM_TIMER_RESIDUALS9);
867 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
868
869 TRACE_POP
870 return 0;
871}
872
873int function_ZeroCrossingsDASSL(int *neqm, double *t, double *y, double *yp,
874 int *ng, double *gout, double *rpar, int* ipar)
875{
876 TRACE_PUSH
877 DATA* data = (DATA*)(void*)((double**)rpar)[0];
878 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
879 threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
880
881 double timeBackup;
882 int saveJumpState;
883
884 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
885 if (measure_time_flag) rt_tick(SIM_TIMER_EVENT4);
886
887 if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
888 {
889 setContext(data, t, CONTEXT_EVENTS);
890 }
891
892 saveJumpState = threadData->currentErrorStage;
893 threadData->currentErrorStage = ERROR_EVENTSEARCH;
894
895 timeBackup = data->localData[0]->timeValue;
896 data->localData[0]->timeValue = *t;
897
898 /* read input vars */
899 externalInputUpdate(data);
900 data->callback->input_function(data, threadData);
901 /* eval needed equations*/
902 data->callback->function_ZeroCrossingsEquations(data, threadData);
903
904 data->callback->function_ZeroCrossings(data, threadData, gout);
905
906 threadData->currentErrorStage = saveJumpState;
907 data->localData[0]->timeValue = timeBackup;
908
909 if (data->simulationInfo->currentContext == CONTEXT_EVENTS){
910 unsetContext(data);
911 }
912
913 if (measure_time_flag) rt_accumulate(SIM_TIMER_EVENT4);
914 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
915
916 TRACE_POP
917 return 0;
918}
919
920/*
921 * Sets element (l,j) in matrixA to given value val.
922 */
923void setJacElementDasslSparse(int l, int j, int nth, double val, void* matrixA,
924 int rows)
925{
926 int k = j*rows + l;
927 ((double*) matrixA)[k]=val;
928}
929
930/* \fn jacA_symColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
931 double *rpar, int* ipar)
932 *
933 *
934 * This function calculates symbolically the jacobian matrix and exploiting the coloring.
935 */
936int jacA_symColored(double *t, double *y, double *yprime, double *delta,
937 double *matrixA, double *cj, double *h, double *wt,
938 double *rpar, int *ipar)
939{
940 TRACE_PUSH
941 DATA* data = (DATA*)(void*)((double**)rpar)[0];
942 threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
943 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
944 const int index = data->callback->INDEX_JAC_A;
945 ANALYTIC_JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[index]);
946
947#ifdef USE_PARJAC
948 ANALYTIC_JACOBIAN* t_jac = (dasslData->jacColumns);
949#else
950 ANALYTIC_JACOBIAN* t_jac = jac;
951#endif
952
953 unsigned int columns = jac->sizeCols;
954 unsigned int rows = jac->sizeRows;
955 unsigned int sizeTmpVars = jac->sizeTmpVars;
956 SPARSE_PATTERN* spp = jac->sparsePattern;
957
958 /* Evaluate constant equations if available */
959 if (jac->constantEqns != NULL((void*)0)) {
960 jac->constantEqns(data, threadData, jac, NULL((void*)0));
961 }
962
963 genericColoredSymbolicJacobianEvaluation(rows, columns, spp, matrixA, t_jac,
964 data, threadData, &setJacElementDasslSparse);
965
966 TRACE_POP
967 return 0;
968}
969
970/* \fn jacA_sym(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
971 double *rpar, int* ipar)
972 *
973 *
974 * This function calculates symbolically the jacobian matrix.
975 * Can calculate the jacobian in parallel.
976 */
977int jacA_sym(double *t, double *y, double *yprime, double *delta,
978 double *matrixA, double *cj, double *h, double *wt, double *rpar,
979 int *ipar)
980{
981 TRACE_PUSH
982
983 DATA* data = (DATA*)(void*)((double**)rpar)[0];
984 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
985 threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
986
987 const int index = data->callback->INDEX_JAC_A;
988 ANALYTIC_JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[index]);
989 unsigned int columns = jac->sizeCols;
990 unsigned int rows = jac->sizeRows;
991 unsigned int sizeTmpVars = jac->sizeTmpVars;
992 unsigned int i;
993
994 /* Evaluate constant equations if available */
995 if (jac->constantEqns != NULL((void*)0)) {
996 jac->constantEqns(data, threadData, jac, NULL((void*)0));
997 }
998
999#ifdef USE_PARJAC
1000 GC_allow_register_threads();
1001#endif
1002
1003#pragma omp parallel default(none) firstprivate(columns, rows, sizeTmpVars) shared(i, matrixA, data, threadData, dasslData)
1004{
1005#ifdef USE_PARJAC
1006 /* Register omp-thread in GC */
1007 if(!GC_thread_is_registered()) {
1008 struct GC_stack_base sb;
1009 memset (&sb, 0, sizeof(sb));
1010 GC_get_stack_base(&sb);
1011 GC_register_my_thread (&sb);
1012 }
1013 // Use thread local analytic Jacobians
1014 ANALYTIC_JACOBIAN* t_jac = &(dasslData->jacColumns[omc_get_thread_num()]);
1015 //printf("index= %d, t_jac->sizeCols= %d, t_jac->sizeRows = %d, t_jac->sizeTmpVars = %d \n",index, t_jac->sizeCols , t_jac->sizeRows, t_jac->sizeTmpVars);
1016#else
1017 ANALYTIC_JACOBIAN* t_jac = jac;
1018#endif
1019 unsigned int j;
1020
1021#pragma omp for schedule(runtime)
1022 for(i=0; i < columns; i++)
1023 {
1024 t_jac->seedVars[i] = 1.0;
1025 data->callback->functionJacA_column(data, threadData, t_jac, NULL((void*)0));
1026
1027 for(j = 0; j < rows; j++)
1028 matrixA[i*columns+j] = t_jac->resultVars[j];
1029
1030 t_jac->seedVars[i] = 0.0;
1031 } // for loop
1032} // omp parallel
1033
1034 TRACE_POP
1035 return 0;
1036}
1037
1038/* \fn jacA_num(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
1039 double *rpar, int* ipar)
1040 *
1041 *
1042 * This function calculates a jacobian matrix by
1043 * numerical with forward finite differences.
1044 */
1045int jacA_num(double *t, double *y, double *yprime, double *delta,
1046 double *matrixA, double *cj, double *h, double *wt, double *rpar,
1047 int *ipar)
1048{
1049 TRACE_PUSH
1050 DATA* data = (DATA*)(void*)((double**)rpar)[0];
1051 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
1052 threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
1053
1054 double delta_h = numericalDifferentiationDeltaXsolver;
1055 double delta_hh,delta_hhh, deltaInv;
1056 double ysave;
1057 int ires;
1058 int i,j;
1059
1060 /* set context for the start values extrapolation of non-linear algebraic loops */
1061 setContext(data, t, CONTEXT_JACOBIAN);
1062
1063 for(i=dasslData->N-1; i >= 0; i--)
1064 {
1065 delta_hhh = *h * yprime[i];
1066 delta_hh = delta_h * fmax(fmax(fabs(y[i]),fabs(delta_hhh)),fabs(1. / wt[i]));
1067 delta_hh = (delta_hhh >= 0 ? delta_hh : -delta_hh);
1068 delta_hh = y[i] + delta_hh - y[i];
1069 deltaInv = 1. / delta_hh;
1070 ysave = y[i];
1071 y[i] += delta_hh;
1072
1073 /* internal dassl numerical jacobian is
1074 * calculated by adding cj to yprime.
1075 * This lead to numerical cancellations.
1076 */
1077 /*yprime[i] += *cj * delta_hh;*/
1078
1079 (*dasslData->residualFunction)(t, y, yprime, cj, dasslData->newdelta, &ires, rpar, ipar);
1080
1081 increaseJacContext(data);
1082
1083 for(j = dasslData->N-1; j >= 0 ; j--)
1084 {
1085 matrixA[i*dasslData->N+j] = (dasslData->newdelta[j] - delta[j]) * deltaInv;
1086 }
1087 y[i] = ysave;
1088 }
1089
1090 TRACE_POP
1091 return 0;
1092}
1093
1094/* \fn jacA_numColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
1095 double *rpar, int* ipar)
1096 *
1097 *
1098 * This function calculates a jacobian matrix by
1099 * numerical with forward finite differences and exploiting the coloring.
1100 */
1101int jacA_numColored(double *t, double *y, double *yprime, double *delta,
1102 double *matrixA, double *cj, double *h, double *wt,
1103 double *rpar, int *ipar)
1104{
1105 TRACE_PUSH
1106
1107 DATA* data = (DATA*)(void*)((double**)rpar)[0];
1108 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
1109 threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
1110
1111 const int index = data->callback->INDEX_JAC_A;
1112 ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
1113
1114 double delta_h = numericalDifferentiationDeltaXsolver;
1115 double delta_hhh;
1116 int ires;
1117 double* delta_hh = dasslData->delta_hh;
1118 double* ysave = dasslData->ysave;
1119
1120 unsigned int i,j,l,k,ii;
1121
1122 /* set context for the start values extrapolation of non-linear algebraic loops */
1123 setContext(data, t, CONTEXT_JACOBIAN);
1124
1125 for(i = 0; i < jacobian->sparsePattern->maxColors; i++)
1126 {
1127 for(ii=0; ii < jacobian->sizeCols; ii++)
1128 {
1129 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
1130 {
1131 delta_hhh = *h * yprime[ii];
1132 delta_hh[ii] = delta_h * fmax(fmax(fabs(y[ii]),fabs(delta_hhh)),fabs(1./wt[ii]));
1133 delta_hh[ii] = (delta_hhh >= 0 ? delta_hh[ii] : -delta_hh[ii]);
1134 delta_hh[ii] = y[ii] + delta_hh[ii] - y[ii];
1135
1136 ysave[ii] = y[ii];
1137 y[ii] += delta_hh[ii];
1138
1139 delta_hh[ii] = 1. / delta_hh[ii];
1140 }
1141 }
1142 (*dasslData->residualFunction)(t, y, yprime, cj, dasslData->newdelta, &ires, rpar, ipar);
1143
1144 increaseJacContext(data);
1145
1146 for(ii = 0; ii < jacobian->sizeCols; ii++)
1147 {
1148 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
1149 {
1150 j = jacobian->sparsePattern->leadindex[ii];
1151 while(j < jacobian->sparsePattern->leadindex[ii+1])
1152 {
1153 l = jacobian->sparsePattern->index[j];
1154 k = l + ii*jacobian->sizeRows;
1155 matrixA[k] = (dasslData->newdelta[l] - delta[l]) * delta_hh[ii];
1156 j++;
1157 };
1158 y[ii] = ysave[ii];
1159 }
1160 }
1161 }
1162
1163 TRACE_POP
1164 return 0;
1165}
1166
1167/* \fn callJacobian(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
1168 double *rpar, int* ipar)
1169 *
1170 *
1171 * This function is called by dassl to calculate the jacobian matrix.
1172 *
1173 */
1174static int callJacobian(double *t, double *y, double *yprime, double *deltaD,
1175 double *pd, double *cj, double *h, double *wt,
1176 double *rpar, int* ipar)
1177{
1178 TRACE_PUSH
1179 DATA* data = (DATA*)(void*)((double**)rpar)[0];
1180 DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
1181 threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
1182
1183 /* profiling */
1184 if (measure_time_flag) rt_accumulate(SIM_TIMER_SOLVER12);
1185 rt_tick(SIM_TIMER_JACOBIAN5);
1186
1187 if(dasslData->jacobianFunction(t, y, yprime, deltaD, pd, cj, h, wt, rpar, ipar))
1188 {
1189 throwStreamPrint(threadData, "Error, can not get Matrix A ");
1190 TRACE_POP
1191 return 1;
1192 }
1193
1194 /* debug */
1195 if (ACTIVE_STREAM(LOG_JAC)(useStream[LOG_JAC])){
1196 _omc_matrix* dumpJac = _omc_createMatrix(dasslData->N, dasslData->N, pd);
1197 _omc_printMatrix(dumpJac, "DASSL-Solver: Matrix A", LOG_JAC);
1198 _omc_destroyMatrix(dumpJac);
1199 }
1200
1201 int i,j = 0;
1202 for(i = 0; i < dasslData->N; i++)
1203 {
1204 pd[j] -= (double) *cj;
1205 j += dasslData->N + 1;
1206 }
1207 /* set context for the start values extrapolation of non-linear algebraic loops */
1208 unsetContext(data);
1209
1210 /* profiling */
1211 rt_accumulate(SIM_TIMER_JACOBIAN5);
1212 if (measure_time_flag) rt_tick(SIM_TIMER_SOLVER12);
1213
1214 TRACE_POP
1215 return 0;
1216}
1217
1218#ifdef __cplusplus
1219}
1220#endif