Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/radau.c
Warning:line 529, column 16
Dereference of undefined pointer value

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
31/*! \file radau.c
32 * author: Ru(n)ge, wbraun
33 * description: This file contains implicit Runge-Kutta methods of different order,
34 * based on Radau IIA and Lobatto IIA methods. The method with order one is
35 * corresponding to the implicit euler method. Further orders 2 to 6.
36 */
37
38#include <string.h>
39#include <math.h>
40
41#include "radau.h"
42#include "external_input.h"
43
44#include "simulation/options.h"
45#ifdef WITH_SUNDIALS
46
47/* adrpo: on mingw link with static sundials */
48#if defined(__MINGW32__)
49#define LINK_SUNDIALS_STATIC
50#endif
51
52#include <kinsol/kinsol.h>
53#include <kinsol/kinsol_dense.h>
54#include <kinsol/kinsol_spgmr.h>
55#include <kinsol/kinsol_sptfqmr.h>
56#include <sundials/sundials_types.h>
57#include <sundials/sundials_math.h>
58
59#ifdef __cplusplus /* wrapper to enable C++ usage */
60extern "C" {
61#endif
62int KINSpbcg(void *kinmem, int maxl);
63extern int KINSetErrHandlerFn(void * kmem, KINErrHandlerFn kinsol_errorHandler, void *);
64extern int KINSetInfoHandlerFn(void *kinmem, KINInfoHandlerFn ihfun, void *ih_data);
65#ifdef __cplusplus
66}
67#endif
68
69static int allocateNlpOde(KINODE *kinOde, int order);
70static int allocateKINSOLODE(KINODE *kinOde);
71
72static void kinsol_errorHandler(int error_code, const char* module, const char* function, char* msg, void* user_data);
73
74static void kinsol_infoHandler(const char* module, const char* function, char* msg, void* user_data);
75
76static int freeImOde(void *nlpode, int N);
77static int freeKinsol(void * kOde);
78
79static int boundsVars(KINODE *kinOde);
80
81static int radau1Coeff(KINODE *kinOd);
82static int radau3Coeff(KINODE *kinOde);
83static int radau5Coeff(KINODE *kinOd);
84static int lobatto4Coeff(KINODE *kinOd);
85static int lobatto6Coeff(KINODE *kinOd);
86
87static int radau1Res(N_Vector z, N_Vector f, void* user_data);
88static int radau3Res(N_Vector z, N_Vector f, void* user_data);
89static int radau5Res(N_Vector z, N_Vector f, void* user_data);
90static int lobatto2Res(N_Vector z, N_Vector f, void* user_data);
91static int lobatto4Res(N_Vector z, N_Vector f, void* user_data);
92static int lobatto6Res(N_Vector z, N_Vector f, void* user_data);
93
94int allocateKinOde(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, int order)
95{
96 int i;
97 KINODE *kinOde = (KINODE*) solverInfo->solverData;
98 kinOde->kData = (KDATAODE*) malloc(sizeof(KDATAODE));
99 kinOde->nlp = (NLPODE*) malloc(sizeof(NLPODE));
100 kinOde->order = order;
101 kinOde->N = ceil((double)order/2.0);
102 kinOde->data = data;
103 kinOde->threadData = threadData;
104 allocateNlpOde(kinOde, order);
1
Calling 'allocateNlpOde'
105 allocateKINSOLODE(kinOde);
106 kinOde->solverInfo = solverInfo;
107 kinOde->kData->mset = 50;
108
109 kinOde->kData->fnormtol = kinOde->data->simulationInfo->tolerance;
110 kinOde->kData->scsteptol = kinOde->data->simulationInfo->tolerance;
111
112 KINSetFuncNormTol(kinOde->kData->kmem, kinOde->kData->fnormtol);
113 KINSetScaledStepTol(kinOde->kData->kmem, kinOde->kData->scsteptol);
114 KINSetNumMaxIters(kinOde->kData->kmem, 10000);
115 if (ACTIVE_STREAM(LOG_SOLVER)(useStream[LOG_SOLVER])) {
116 KINSetPrintLevel(kinOde->kData->kmem,2);
117 }
118 /* KINSetEtaForm(kinOde->kData->kmem, KIN_ETACHOICE2); */
119 KINSetMaxSetupCalls(kinOde->kData->kmem, kinOde->kData->mset);
120
121 KINSetErrHandlerFn(kinOde->kData->kmem, kinsol_errorHandler, NULL((void*)0));
122 KINSetInfoHandlerFn(kinOde->kData->kmem, kinsol_infoHandler, NULL((void*)0));
123 switch(kinOde->order)
124 {
125 case 5:
126 KINInit(kinOde->kData->kmem, radau5Res, kinOde->kData->x);
127 break;
128 case 3:
129 KINInit(kinOde->kData->kmem, radau3Res, kinOde->kData->x);
130 break;
131 case 1:
132 KINInit(kinOde->kData->kmem, radau1Res, kinOde->kData->x);
133 break;
134 case 2:
135 KINInit(kinOde->kData->kmem, lobatto2Res, kinOde->kData->x);
136 break;
137 case 4:
138 KINInit(kinOde->kData->kmem, lobatto4Res, kinOde->kData->x);
139 break;
140 case 6:
141 KINInit(kinOde->kData->kmem, lobatto6Res, kinOde->kData->x);
142 break;
143
144 default:
145 assert(0)((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else __assert_fail
("0", "./simulation/solver/radau.c", 145, __extension__ __PRETTY_FUNCTION__
); }))
;
146 }
147 /* if FLAG_IMPRK_LS is set, choose ida linear solver method */
148 if (omc_flag[FLAG_IMPRK_LS])
149 {
150 for(i=1; i< IMPRK_LS_MAX;i++)
151 {
152 if (!strcmp((const char*)omc_flagValue[FLAG_IMPRK_LS], IMPRK_LS_METHOD[i])){
153 kinOde->lsMethod = (int)i;
154 break;
155 }
156 }
157 if (kinOde->lsMethod == IMPRK_LS_UNKNOWN)
158 {
159 if (ACTIVE_WARNING_STREAM(LOG_SOLVER)(showAllWarnings || useStream[LOG_SOLVER]))
160 {
161 warningStreamPrint(LOG_SOLVER, 1, "unrecognized linear solver method %s, current options are:", (const char*)omc_flagValue[FLAG_IMPRK_LS]);
162 for(i=1; i < IMPRK_LS_MAX; ++i)
163 {
164 warningStreamPrint(LOG_SOLVER, 0, "%-15s [%s]", IMPRK_LS_METHOD[i], IMPRK_LS_METHOD_DESC[i]);
165 }
166 messageClose(LOG_SOLVER);
167 }
168 throwStreamPrint(threadData,"unrecognized linear solver method %s", (const char*)omc_flagValue[FLAG_IMPRK_LS]);
169 }
170 }
171 else
172 {
173 kinOde->lsMethod = IMPRK_LS_ITERATIVE;
174 }
175
176 kinOde->kData->glstr = KIN_LINESEARCH1;
177
178 switch (kinOde->lsMethod){
179 case IMPRK_LS_ITERATIVE:
180 /*KINSpbcg*/
181 if(kinOde->nlp->nStates < 10)
182 KINSpgmr(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates+1);
183 else
184 KINSpbcg(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates+1);
185 break;
186 case IMPRK_LS_DENSE:
187 KINDense(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates);
188 break;
189 default:
190 throwStreamPrint(threadData,"unrecognized linear solver method %s", (const char*)omc_flagValue[FLAG_IMPRK_LS]);
191 break;
192 }
193 KINSetNoInitSetup(kinOde->kData->kmem, FALSE0);
194
195 return 0;
196}
197
198static int allocateNlpOde(KINODE *kinOde, int order)
199{
200 NLPODE * nlp = (NLPODE*) kinOde->nlp;
201 nlp->nStates = kinOde->data->modelData->nStates;
202
203 switch(order)
2
Control jumps to 'case 3:' at line 208
204 {
205 case 5:
206 radau5Coeff(kinOde);
207 break;
208 case 3:
209 radau3Coeff(kinOde);
3
Calling 'radau3Coeff'
210 break;
211 case 1:
212 radau1Coeff(kinOde);
213 break;
214 case 2:
215 radau1Coeff(kinOde); /* TODO: Is this right? */
216 break;
217 case 4:
218 lobatto4Coeff(kinOde);
219 break;
220 case 6:
221 lobatto6Coeff(kinOde);
222 break;
223 default:
224 assert(0)((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else __assert_fail
("0", "./simulation/solver/radau.c", 224, __extension__ __PRETTY_FUNCTION__
); }))
;
225 }
226
227 boundsVars(kinOde);
228 return 0;
229}
230
231static void kinsol_errorHandler(int error_code, const char* module, const char* function, char* msg, void* user_data)
232{
233 warningStreamPrint(LOG_SOLVER, 0, "[module] %s | [function] %s | [error_code] %d", module, function, error_code);
234 if (msg) warningStreamPrint(LOG_SOLVER, 0, "%s", msg);
235}
236
237static void kinsol_infoHandler(const char* module, const char* function, char* msg, void* user_data)
238{
239 infoStreamPrint(LOG_SOLVER, 0, " %s: %s ", module, function);
240 if (msg) infoStreamPrint(LOG_SOLVER, 0, "%s", msg);
241}
242
243static int allocateKINSOLODE(KINODE *kinOde)
244{
245 int m;
246 DATA *data = kinOde->data;
247 int n = data->modelData->nStates;
248 int i, j, k;
249 double* c;
250 KDATAODE * kData = (kinOde)->kData;
251 m = kinOde->N*n;
252 kData->x = N_VNew_Serial(m);
253 kData->sVars = N_VNew_Serial(m);
254 kData->sEqns = N_VNew_Serial(m);
255 kData->c = N_VNew_Serial(m);
256 kData->kmem = KINCreate();
257 c = NV_DATA_S(kData->c)( ( (N_VectorContent_Serial)(kData->c->content) )->data
)
;
258
259 for(j=0, k=0; j< kinOde->N; ++j)
260 for(i=0; i<n; ++i, ++k)
261 c[k] = 0;
262
263 KINSetUserData(kinOde->kData->kmem, (void*) kinOde);
264 KINSetConstraints(kinOde->kData->kmem, kData->c);
265 return 0;
266}
267
268int freeKinOde(DATA* data, SOLVER_INFO* solverInfo)
269{
270 KINODE *kinOde = (KINODE*) solverInfo->solverData;
271 freeImOde((void*) kinOde->nlp, kinOde->N);
272 freeKinsol((void*) kinOde->kData);
273 free(kinOde);
274 return 0;
275}
276
277static int freeImOde(void *nlpode, int N)
278{
279 int i;
280 NLPODE *nlp = (NLPODE*) nlpode;
281 free(nlp->min);
282 free(nlp->max);
283 free(nlp->s);
284
285 for(i=0; i<N; i++) {
286 free(nlp->c[i]);
287 }
288 free(nlp->c);
289
290 free(nlp->a);
291 return 0;
292}
293
294static int freeKinsol(void * kOde)
295{
296 KDATAODE *kData = (KDATAODE*) kOde;
297 N_VDestroy_Serial(kData->x);
298 N_VDestroy_Serial(kData->sVars);
299 N_VDestroy_Serial(kData->sEqns);
300 N_VDestroy_Serial(kData->c);
301 KINFree(&kData->kmem);
302 return 0;
303}
304
305static int initKinsol(KINODE *kinOde)
306{
307 int i,j,k,n;
308 double *scal_eq, *scal_var, *x, *f2;
309 long double tmp, h, hf, hf_min;
310 DATA *data;
311 NLPODE *nlp;
312 KDATAODE * kData;
313
314 n = kinOde->nlp->nStates;
315 data= kinOde->data;
316 kData = kinOde->kData;
317 x = NV_DATA_S(kData->x)( ( (N_VectorContent_Serial)(kData->x->content) )->data
)
;
318 nlp = kinOde->nlp;
319
320 nlp->currentStep = kinOde->solverInfo->currentStepSize;
321 f2 = data->localData[2]->realVars + n;
322 nlp->dt = nlp->currentStep;
323 nlp->derx = data->localData[0]->realVars + n;
324 nlp->x0 = data->localData[1]->realVars;
325 nlp->f0 = data->localData[1]->realVars + n;
326 nlp->t0 = data->localData[1]->timeValue;
327
328
329 scal_var = NV_DATA_S(kData->sVars)( ( (N_VectorContent_Serial)(kData->sVars->content) )->
data )
;
330 scal_eq = NV_DATA_S(kData->sEqns)( ( (N_VectorContent_Serial)(kData->sEqns->content) )->
data )
;
331
332 hf_min = 1e-6;
333 for(j=0, k=0; j<kinOde->N; ++j)
334 {
335 for(i=0;i<n;++i,++k)
336 {
337 hf = 0.5*nlp->dt*nlp->a[j]*(3*nlp->f0[i]-f2[i]);
338 if(fabsl(hf) < hf_min) hf_min = fabsl(hf);
339 x[k] = (nlp->x0[i] + hf);
340 tmp = fabs(x[k] + nlp->x0[i]) + 1e-12;
341 tmp = (tmp < 1e-9) ? nlp->s[i] : 2.0/tmp;
342 scal_var[k] = tmp + 1e-9;
343 scal_eq[k] = 1.0/(scal_var[k])+ 1e-12;
344 }
345 }
346 KINSetMaxNewtonStep(kinOde->kData->kmem, hf_min);
347 return 0;
348}
349
350int kinsolOde(SOLVER_INFO* solverInfo)
351{
352 KINODE *kinOde = (KINODE*) solverInfo->solverData;
353 KDATAODE *kData = kinOde->kData;
354 int flag, dense=0, try_again = 1;
355 int retries = 0;
356
357 infoStreamPrint(LOG_SOLVER, 1, "##IMPRK## new step from %.15g to %.15g", solverInfo->currentTime, solverInfo->currentTime + solverInfo->currentStepSize);
358 initKinsol(kinOde);
359
360 do
361 {
362 kData->error_code = KINSol(kData->kmem, /* KINSol memory block */
363 kData->x, /* initial guess on input; solution vector */
364 kData->glstr, /* global strategy choice */
365 kData->sVars, /* scaling vector, for the variable cc */
366 kData->sEqns /* scaling vector, for the residual eqns */
367 );
368
369 switch(kinOde->lsMethod){
370 retries++;
371 case IMPRK_LS_ITERATIVE:
372 if(retries == 0)
373 {
374 KINDense(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates);
375 dense=1;
376 warningStreamPrint(LOG_SOLVER,0,"Restart Kinsol: change linear solver to KINDense.");
377 }
378 else if(retries == 1)
379 {
380 KINSptfqmr(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates);
381 dense=0;
382 warningStreamPrint(LOG_SOLVER,0,"Restart Kinsol: change linear solver to KINSptfqmr.");
383 }
384 else if(retries == 2)
385 {
386 KINSpbcg(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates);
387 warningStreamPrint(LOG_SOLVER,0,"Restart Kinsol: change linear solver to KINSpbcg.");
388 }
389 else
390 {
391 try_again= 0;
392 }
393 break;
394 case IMPRK_LS_DENSE:
395 dense=1;
396 if (retries== 1)
397 {
398 kinOde->kData->glstr = KIN_NONE0;
399 }
400 else
401 {
402 try_again= 0;
403 }
404 break;
405 default:
406 kData->error_code = -42;
407 try_again = 0;
408 break;
409 }
410
411 }while(kData->error_code < 0 && try_again);
412
413 /* save stats */
414 /* steps */
415 solverInfo->solverStatsTmp[0] += 1;
416 /* functionODE evaluations */
417 {
418 long int tmp = 0;
419 flag = KINGetNumFuncEvals(kData->kmem, &tmp);
420 if (flag == KIN_SUCCESS0)
421 {
422 solverInfo->solverStatsTmp[1] += tmp;
423 }
424 }
425
426 /* Jacobians evaluations */
427 {
428 long int tmp = 0;
429 if (dense)
430 {
431 flag = KINDlsGetNumJacEvals(kData->kmem, &tmp);
432 }
433 else
434 {
435 flag = KINSpilsGetNumJtimesEvals(kData->kmem, &tmp);
436 }
437 if (flag == KIN_SUCCESS0)
438 {
439 solverInfo->solverStatsTmp[2] += tmp;
440 }
441 }
442
443 /* beta-condition failures evaluations */
444 {
445 long int tmp = 0;
446 flag = KINGetNumBetaCondFails(kData->kmem, &tmp);
447 if (flag == KIN_SUCCESS0)
448 {
449 solverInfo->solverStatsTmp[4] += tmp;
450 }
451 }
452 infoStreamPrint(LOG_SOLVER, 0, "##IMPRK## Integration step finished .");
453 /* closing new step message */
454 messageClose(LOG_SOLVER);
455
456 return (kData->error_code<0) ? -1 : 0;
457}
458
459
460static int boundsVars(KINODE *kinOde)
461{
462 int i;
463 double tmp;
464 NLPODE * nlp = (NLPODE*) kinOde->nlp;
465 DATA * data = (DATA*) kinOde->data;
466
467 nlp->min = (double*) malloc(nlp->nStates* sizeof(double));
468 nlp->max = (double*) malloc(nlp->nStates* sizeof(double));
469 nlp->s = (double*) malloc(nlp->nStates* sizeof(double));
470
471 for(i =0;i<nlp->nStates;i++)
472 {
473 nlp->min[i] = data->modelData->realVarsData[i].attribute.min;
474 nlp->max[i] = data->modelData->realVarsData[i].attribute.max;
475 tmp = fabs(data->modelData->realVarsData[i].attribute.nominal);
476 tmp = tmp >= 0.0 ? tmp : 1.0;
477 nlp->s[i] = 1.0/tmp;
478 }
479 return 0;
480}
481
482static int radau5Coeff(KINODE *kinOde)
483{
484 int i;
485 NLPODE * nlp = (NLPODE*) kinOde->nlp;
486
487 nlp->c = (long double**) malloc(kinOde->N * sizeof(long double*));
488 for(i = 0; i < kinOde->N; i++)
489 nlp->c[i] = (long double*) calloc(kinOde->N+1, sizeof(long double));
490
491 nlp->a = (double*) malloc(kinOde->N * sizeof(double));
492
493 nlp->c[0][0] = 4.1393876913398137178367408896470696703591369767880;
494 nlp->c[0][1] = 3.2247448713915890490986420373529456959829737403284;
495 nlp->c[0][2] = 1.1678400846904054949240412722156950122337492313015;
496 nlp->c[0][3] = 0.25319726474218082618594241992157103785758599484179;
497
498 nlp->c[1][0] = 1.7393876913398137178367408896470696703591369767880;
499 nlp->c[1][1] = 3.5678400846904054949240412722156950122337492313015;
500 nlp->c[1][2] = 0.7752551286084109509013579626470543040170262596716;
501 nlp->c[1][3] = 1.0531972647421808261859424199215710378575859948418;
502
503 nlp->c[2][0] = 3.0;
504 nlp->c[2][1] = 5.5319726474218082618594241992157103785758599484179;
505 nlp->c[2][2] = 7.5319726474218082618594241992157103785758599484179;
506 nlp->c[2][3] = 5.0;
507
508 nlp->a[0] = 0.15505102572168219018027159252941086080340525193433;
509 nlp->a[1] = 0.64494897427831780981972840747058913919659474806567;
510 nlp->a[2] = 1.0;
511 return 0;
512}
513
514static int radau3Coeff(KINODE *kinOde)
515{
516 int i;
517 NLPODE * nlp = (NLPODE*) kinOde->nlp;
518
519 nlp->c = (long double**) malloc(kinOde->N * sizeof(long double*));
4
Storing uninitialized value
520 for(i = 0; i < kinOde->N; i++)
5
Assuming the condition is true
6
Loop condition is true. Entering loop body
7
Assuming the condition is false
8
Loop condition is false. Execution continues on line 523
521 nlp->c[i] = (long double*) calloc(kinOde->N+1, sizeof(long double));
522
523 nlp->a = (double*) malloc(kinOde->N * sizeof(double));
524
525 nlp->c[0][0] = 2.0;
526 nlp->c[0][1] = 1.50;
527 nlp->c[0][2] = 0.50;;
528
529 nlp->c[1][0] = 2.0;
9
Dereference of undefined pointer value
530 nlp->c[1][1] = 4.50;
531 nlp->c[1][2] = 2.50;
532
533 nlp->a[0] = 1.0/3.0;
534 nlp->a[1] = 1.0;
535 return 0;
536}
537
538static int radau1Coeff(KINODE *kinOde)
539{
540 NLPODE * nlp = (NLPODE*) kinOde->nlp;
541 nlp->c = (long double**) malloc(kinOde->N * sizeof(long double*));
542 nlp->c[0] = (long double*) malloc(kinOde->N * sizeof(long double));
543 nlp->a = (double*) malloc(kinOde->N * sizeof(double));
544 nlp->a[0] = 1.0;
545 return 0;
546}
547
548static int lobatto4Coeff(KINODE *kinOde)
549{
550 NLPODE * nlp = (NLPODE*) kinOde->nlp;
551 nlp->c = (long double**) malloc(kinOde->N * sizeof(long double*));
552 nlp->c[0] = (long double*) malloc(kinOde->N * sizeof(long double));
553 nlp->c[1] = (long double*) malloc(kinOde->N * sizeof(long double));
554 nlp->a = (double*) malloc(kinOde->N * sizeof(double));
555 nlp->a[0] = 0.5;
556 nlp->a[1] = 1.0;
557 return 0;
558}
559
560static int lobatto6Coeff(KINODE *kinOde)
561{
562 int i;
563 NLPODE * nlp = (NLPODE*) kinOde->nlp;
564
565
566 nlp->c = (long double**) malloc(kinOde->N * sizeof(long double*));
567 for(i = 0; i < kinOde->N; ++i)
568 nlp->c[i] = (long double*) malloc((kinOde->N+2)* sizeof(long double));
569
570 nlp->a = (double*) malloc(kinOde->N * sizeof(double));
571
572 nlp->c[0][0] = 4.3013155617496424838955952368431696002490512113396;
573 nlp->c[0][1] = 3.6180339887498948482045868343656381177203091798058;
574 nlp->c[0][2] = 0.8541019662496845446137605030969143531609275394172;
575 nlp->c[0][3] = 0.17082039324993690892275210061938287063218550788345;
576 nlp->c[0][4] = 0.44721359549995793928183473374625524708812367192230;
577
578 nlp->c[1][0] = 3.3013155617496424838955952368431696002490512113396;
579 nlp->c[1][1] = 5.8541019662496845446137605030969143531609275394172;
580 nlp->c[1][2] = 1.3819660112501051517954131656343618822796908201942;
581 nlp->c[1][3] = 1.1708203932499369089227521006193828706321855078834;
582 nlp->c[1][4] = 0.44721359549995793928183473374625524708812367192230;
583
584 nlp->c[2][0] = 7.0;
585 nlp->c[2][1] = 11.180339887498948482045868343656381177203091798058;
586 nlp->c[2][2] = 11.180339887498948482045868343656381177203091798058;
587 nlp->c[2][3] = 7.0;
588 nlp->c[2][4] = 1.0;
589
590 nlp->a[0] = 0.27639320225002103035908263312687237645593816403885;
591 nlp->a[1] = 0.72360679774997896964091736687312762354406183596115;
592 nlp->a[2] = 1.0;
593 return 0;
594}
595
596static int refreshModell(DATA* data, threadData_t *threadData, double* x, double time)
597{
598 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
599
600 memcpy(sData->realVars, x, sizeof(double)*data->modelData->nStates);
601 sData->timeValue = time;
602 /* read input vars */
603 externalInputUpdate(data);
604 data->callback->input_function(data, threadData);
605 data->callback->functionODE(data, threadData);
606
607 return 0;
608}
609
610static int radau5Res(N_Vector x, N_Vector f, void* user_data)
611{
612 int i,k;
613 KINODE* kinOde = (KINODE*)user_data;
614 NLPODE *nlp = kinOde->nlp;
615 threadData_t *threadData = kinOde->threadData;
616
617 double *x0,*x1,*x2,*x3;
618 double*derx = nlp->derx;
619 double*a = nlp->a;
620 double* feq = NV_DATA_S(f)( ( (N_VectorContent_Serial)(f->content) )->data );
621
622 x0 = nlp->x0;
623 x1 = NV_DATA_S(x)( ( (N_VectorContent_Serial)(x->content) )->data );
624 x2 = x1 + nlp->nStates;
625 x3 = x2 + nlp->nStates;
626
627 refreshModell(kinOde->data, threadData, x1, nlp->t0 + a[0]*nlp->dt);
628 for(i = 0;i<nlp->nStates;i++)
629 {
630 feq[i] = (nlp->c[0][0]*x0[i] + nlp->c[0][3]*x3[i] + nlp->dt*derx[i]) -
631 (nlp->c[0][1]*x1[i] + nlp->c[0][2]*x2[i]);
632 if(isnan(feq[i])(sizeof ((feq[i])) == sizeof (float) ? __isnanf (feq[i]) : sizeof
((feq[i])) == sizeof (double) ? __isnan (feq[i]) : __isnanl (
feq[i]))
) return -1;
633 }
634
635 refreshModell(kinOde->data, threadData, x2, nlp->t0 + a[1]*nlp->dt);
636 for(i = 0, k=nlp->nStates; i<nlp->nStates; i++, k++)
637 {
638 feq[k] = (nlp->c[1][1]*x1[i] + nlp->dt*derx[i]) -
639 (nlp->c[1][0]*x0[i] + nlp->c[1][2]*x2[i] + nlp->c[1][3]*x3[i]);
640 if(isnan(feq[k])(sizeof ((feq[k])) == sizeof (float) ? __isnanf (feq[k]) : sizeof
((feq[k])) == sizeof (double) ? __isnan (feq[k]) : __isnanl (
feq[k]))
) return -1;
641 }
642
643 refreshModell(kinOde->data, threadData, x3, nlp->t0 + nlp->dt);
644 for(i = 0;i<nlp->nStates;i++,k++)
645 {
646 feq[k] = (nlp->c[2][0]*x0[i] + nlp->c[2][2]*x2[i] + nlp->dt*derx[i]) -
647 (nlp->c[2][1]*x1[i] + nlp->c[2][3]*x3[i]);
648 if(isnan(feq[k])(sizeof ((feq[k])) == sizeof (float) ? __isnanf (feq[k]) : sizeof
((feq[k])) == sizeof (double) ? __isnan (feq[k]) : __isnanl (
feq[k]))
) return -1;
649 }
650
651 return 0;
652}
653
654static int radau3Res(N_Vector x, N_Vector f, void* user_data)
655{
656 int i,k;
657 KINODE* kinOde = (KINODE*)user_data;
658 NLPODE *nlp = kinOde->nlp;
659 DATA *data = kinOde->data;
660 threadData_t *threadData = kinOde->threadData;
661
662 double* feq = NV_DATA_S(f)( ( (N_VectorContent_Serial)(f->content) )->data );
663
664 double *x0,*x1,*x2;
665
666 double*derx = nlp->derx;
667
668 x0 = nlp->x0;
669 x1 = NV_DATA_S(x)( ( (N_VectorContent_Serial)(x->content) )->data );
670 x2 = x1 + nlp->nStates;
671
672 refreshModell(data, threadData, x1, nlp->t0 + 0.5*nlp->dt);
673 for(i = 0;i<nlp->nStates;i++)
674 {
675 feq[i] = (nlp->c[0][0]*x0[i] + nlp->dt*derx[i]) -
676 (nlp->c[0][1]*x1[i] + nlp->c[0][2]*x2[i]);
677 if(isnan(feq[i])(sizeof ((feq[i])) == sizeof (float) ? __isnanf (feq[i]) : sizeof
((feq[i])) == sizeof (double) ? __isnan (feq[i]) : __isnanl (
feq[i]))
) return -1;
678 }
679
680 refreshModell(data, threadData, x2, nlp->t0 + nlp->dt);
681 for(i = 0, k=nlp->nStates;i<nlp->nStates;i++,k++)
682 {
683 feq[k] = (nlp->c[1][1]*x1[i] + nlp->dt*derx[i]) -
684 (nlp->c[1][0]*x0[i] + nlp->c[1][2]*x2[i]);
685 if(isnan(feq[k])(sizeof ((feq[k])) == sizeof (float) ? __isnanf (feq[k]) : sizeof
((feq[k])) == sizeof (double) ? __isnan (feq[k]) : __isnanl (
feq[k]))
) return -1;
686 }
687
688 return 0;
689}
690
691
692static int radau1Res(N_Vector x, N_Vector f, void* user_data)
693{
694 int i;
695 KINODE* kinOde = (KINODE*)user_data;
696 NLPODE *nlp = kinOde->nlp;
697 threadData_t *threadData = kinOde->threadData;
698 double* feq = NV_DATA_S(f)( ( (N_VectorContent_Serial)(f->content) )->data );
699
700 double *x0,*x1;
701 double*derx = nlp->derx;
702
703 x0 = nlp->x0;
704 x1 = NV_DATA_S(x)( ( (N_VectorContent_Serial)(x->content) )->data );
705
706 refreshModell(kinOde->data, threadData, x1, nlp->t0 + nlp->dt);
707 for(i = 0; i<nlp->nStates; ++i)
708 {
709 feq[i] = x0[i] - x1[i] + nlp->dt*derx[i];
710 if(isnan(feq[i])(sizeof ((feq[i])) == sizeof (float) ? __isnanf (feq[i]) : sizeof
((feq[i])) == sizeof (double) ? __isnan (feq[i]) : __isnanl (
feq[i]))
) return -1;
711 }
712 return 0;
713}
714
715static int lobatto2Res(N_Vector x, N_Vector f, void* user_data)
716{
717 int i;
718 KINODE* kinOde = (KINODE*)user_data;
719 NLPODE *nlp = kinOde->nlp;
720 DATA *data = kinOde->data;
721 threadData_t *threadData = kinOde->threadData;
722
723 double *feq = NV_DATA_S(f)( ( (N_VectorContent_Serial)(f->content) )->data );
724 double *x0,*x1, *f0;
725 double *derx = nlp->derx;
726
727 x0 = nlp->x0;
728 f0 = nlp->f0;
729 x1 = NV_DATA_S(x)( ( (N_VectorContent_Serial)(x->content) )->data );
730
731 refreshModell(data, threadData, x1, nlp->t0 + nlp->dt);
732 for(i = 0;i<nlp->nStates;i++)
733 {
734 feq[i] = x0[i] - x1[i] + 0.5*nlp->dt*(f0[i]+derx[i]);
735 if(isnan(feq[i])(sizeof ((feq[i])) == sizeof (float) ? __isnanf (feq[i]) : sizeof
((feq[i])) == sizeof (double) ? __isnan (feq[i]) : __isnanl (
feq[i]))
) return -1;
736 }
737 return 0;
738}
739
740static int lobatto4Res(N_Vector x, N_Vector f, void* user_data)
741{
742 int i,k;
743 KINODE* kinOde = (KINODE*)user_data;
744 NLPODE *nlp = kinOde->nlp;
745 DATA *data = kinOde->data;
746 threadData_t *threadData = kinOde->threadData;
747
748 double* feq = NV_DATA_S(f)( ( (N_VectorContent_Serial)(f->content) )->data );
749
750 double *x0, *x1, *x2, *f0;
751 double*derx = nlp->derx;
752
753 f0 = nlp->f0;
754 x0 = nlp->x0;
755 x1 = NV_DATA_S(x)( ( (N_VectorContent_Serial)(x->content) )->data );
756 x2 = x1 + nlp->nStates;
757
758 refreshModell(data, threadData, x1,nlp->t0 + 0.5*nlp->dt);
759 for(i = 0;i<nlp->nStates;i++)
760 {
761 feq[i] = (nlp->dt*(2.0*derx[i] +f0[i]) + 5.0*x0[i]) - (4*x1[i] + x2[i]);
762 if(isnan(feq[i])(sizeof ((feq[i])) == sizeof (float) ? __isnanf (feq[i]) : sizeof
((feq[i])) == sizeof (double) ? __isnan (feq[i]) : __isnanl (
feq[i]))
) return -1;
763 }
764
765 refreshModell(data, threadData, x2,nlp->t0 + nlp->dt);
766 for(i = 0,k=nlp->nStates;i<nlp->nStates;i++,k++)
767 {
768 feq[k] = (2.0*nlp->dt*derx[i] + 16.0*x1[i]) - (8.0*(x0[i] + x2[i]) +2.0*nlp->dt*f0[i]);
769 if(isnan(feq[k])(sizeof ((feq[k])) == sizeof (float) ? __isnanf (feq[k]) : sizeof
((feq[k])) == sizeof (double) ? __isnan (feq[k]) : __isnanl (
feq[k]))
) return -1;
770 }
771 return 0;
772}
773
774static int lobatto6Res(N_Vector x, N_Vector f, void* user_data)
775{
776 int i, k;
777 KINODE* kinOde = (KINODE*)user_data;
778 NLPODE *nlp = kinOde->nlp;
779 DATA *data = kinOde->data;
780 threadData_t *threadData = kinOde->threadData;
781
782 double* feq = NV_DATA_S(f)( ( (N_VectorContent_Serial)(f->content) )->data );
783 double *x0, *x1, *x2, *x3, *f0;
784 double*derx = nlp->derx;
785
786 f0 = nlp->f0;
787 x0 = nlp->x0;
788 x1 = NV_DATA_S(x)( ( (N_VectorContent_Serial)(x->content) )->data );
789 x2 = x1 + nlp->nStates;
790 x3 = x2 + nlp->nStates;
791
792 refreshModell(data, threadData, x1, nlp->t0 + nlp->a[0]*nlp->dt);
793 for(i = 0;i<nlp->nStates;i++)
794 {
795 feq[i] = (nlp->dt*(derx[i] + nlp->c[0][4]*f0[i]) + nlp->c[0][0]*x0[i] + nlp->c[0][3]*x3[i]) - (nlp->c[0][1]*x1[i] + nlp->c[0][2]*x2[i]);
796 if(isnan(feq[i])(sizeof ((feq[i])) == sizeof (float) ? __isnanf (feq[i]) : sizeof
((feq[i])) == sizeof (double) ? __isnan (feq[i]) : __isnanl (
feq[i]))
) return -1;
797 }
798
799 refreshModell(data, threadData, x2,nlp->t0 + nlp->a[1]*nlp->dt);
800 for(i = 0,k=nlp->nStates;i<nlp->nStates;i++,k++)
801 {
802 feq[k] = (nlp->dt*derx[i] + nlp->c[1][1]*x1[i]) - (nlp->dt*nlp->c[1][4]*f0[i] + nlp->c[1][0]*x0[i] + nlp->c[1][2]*x2[i] + nlp->c[1][3]*x3[i]);
803 if(isnan(feq[k])(sizeof ((feq[k])) == sizeof (float) ? __isnanf (feq[k]) : sizeof
((feq[k])) == sizeof (double) ? __isnan (feq[k]) : __isnanl (
feq[k]))
) return -1;
804 }
805
806 refreshModell(data, threadData, x3,nlp->t0 + nlp->dt);
807 for(i = 0;i<nlp->nStates;i++,k++)
808 {
809 feq[k] = (nlp->dt*(f0[i] + derx[i]) + nlp->c[2][0]*x0[i] + nlp->c[2][2]*x2[i]) - (nlp->c[2][1]*x1[i] + nlp->c[2][3]*x3[i]);
810 if(isnan(feq[k])(sizeof ((feq[k])) == sizeof (float) ? __isnanf (feq[k]) : sizeof
((feq[k])) == sizeof (double) ? __isnan (feq[k]) : __isnanl (
feq[k]))
) return -1;
811 }
812
813 return 0;
814}
815#else
816
817int kinsolOde(SOLVER_INFO* solverInfo)
818{
819 assert(0)((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else __assert_fail
("0", "./simulation/solver/radau.c", 819, __extension__ __PRETTY_FUNCTION__
); }))
;
820 return -1;
821}
822#endif