/[ascend]/trunk/base/generic/integrator/ida.c
ViewVC logotype

Contents of /trunk/base/generic/integrator/ida.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1218 - (show annotations) (download) (as text)
Tue Jan 23 13:42:49 2007 UTC (17 years, 4 months ago) by johnpye
File MIME type: text/x-csrc
File size: 62544 byte(s)
Fixed integrator_ida_check_lists.
1 /* ASCEND modelling environment
2 Copyright (C) 2006-2007 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//**
19 @file
20 Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
21 as part of the GPL-licensed SUNDIALS solver package from LLNL.
22
23 IDA provides the non-linear parts, as well as a number of pluggable linear
24 solvers: dense, banded and krylov types.
25
26 We also implement here an EXPERIMENTAL direct sparse linear solver for IDA
27 using the ASCEND linsolqr routines.
28
29 @see http://www.llnl.gov/casc/sundials/
30 *//*
31 by John Pye, May 2006
32 */
33
34 #include "ida.h"
35
36 #include <signal.h>
37
38 /* SUNDIALS includes */
39 #ifdef ASC_WITH_IDA
40 # include <sundials/sundials_config.h>
41 # include <sundials/sundials_dense.h>
42 # include <ida/ida.h>
43 # include <nvector/nvector_serial.h>
44 # include <ida/ida_spgmr.h>
45 # include <ida/ida_spbcgs.h>
46 # include <ida/ida_sptfqmr.h>
47 # include <ida/ida_dense.h>
48 # ifndef IDA_SUCCESS
49 # error "Failed to include SUNDIALS IDA header file"
50 # endif
51 #else
52 # error "If you're building this file, you should have ASC_WITH_IDA"
53 #endif
54
55 #ifdef ASC_WITH_MMIO
56 # include <mmio.h>
57 #endif
58
59 #include <utilities/ascConfig.h>
60 #include <utilities/error.h>
61 #include <utilities/ascSignal.h>
62 #include <utilities/ascPanic.h>
63 #include <compiler/instance_enum.h>
64
65 #include <solver/slv_client.h>
66 #include <solver/block.h>
67 #include <solver/slvDOF.h>
68 #include <solver/relman.h>
69 #ifdef ASC_IDA_NEW_ANALYSE
70 # include <solver/diffvars.h>
71 # include <solver/slv_stdcalls.h>
72 #endif
73
74 #include "idalinear.h"
75
76 /*
77 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
78 */
79 #ifndef SUNDIALS_VERSION_MINOR
80 # ifdef __GNUC__
81 # warning "GUESSING SUNDIALS VERSION 2.2"
82 # endif
83 # define SUNDIALS_VERSION_MINOR 2
84 #endif
85 #ifndef SUNDIALS_VERSION_MAJOR
86 # define SUNDIALS_VERSION_MAJOR 2
87 #endif
88
89 /* check that we've got what we expect now */
90 #ifndef ASC_IDA_H
91 # error "Failed to include ASCEND IDA header file"
92 #endif
93
94 /* #define FEX_DEBUG */
95 #define JEX_DEBUG
96 #define DJEX_DEBUG
97 #define SOLVE_DEBUG
98 #define STATS_DEBUG
99 #define PREC_DEBUG
100 /* #define ANALYSE_DEBUG */
101 /**
102 Everthing that the outside world needs to know about IDA
103 */
104 const IntegratorInternals integrator_ida_internals = {
105 integrator_ida_create
106 ,integrator_ida_params_default
107 #ifdef ASC_IDA_NEW_ANALYSE
108 ,integrator_ida_analyse
109 #else
110 ,integrator_analyse_dae /* note, this routine is back in integrator.c */
111 #endif
112 ,integrator_ida_solve
113 ,NULL /* writematrixfn */
114 ,integrator_ida_debug
115 ,integrator_ida_free
116 ,INTEG_IDA
117 ,"IDA"
118 };
119
120 /*-------------------------------------------------------------
121 FORWARD DECLS
122 */
123
124 /* forward dec needed for IntegratorIdaPrecFreeFn */
125 struct IntegratorIdaDataStruct;
126
127 /* functions for allocating storage for and freeing preconditioner data */
128 typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *blsys);
129 typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);
130
131 /**
132 Struct containing any stuff that IDA needs that doesn't fit into the
133 common IntegratorSystem struct.
134 */
135 typedef struct IntegratorIdaDataStruct{
136 struct rel_relation **rellist; /**< NULL terminated list of rels */
137 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
138 int nrels;
139 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
140 var_filter_t vfilter;
141 rel_filter_t rfilter; /**< Used to filter relations from rellist (@TODO needs work) */
142 void *precdata; /**< For use by the preconditioner */
143 IntegratorIdaPrecFreeFn *pfree; /**< Store instructions here on how to free precdata */
144 } IntegratorIdaData;
145
146 typedef struct IntegratorIdaPrecDJStruct{
147 N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */
148 } IntegratorIdaPrecDataJacobi;
149
150 /**
151 Hold all the function pointers associated with a particular preconditioner
152 We don't need to store the 'pfree' function here as it is allocated to the enginedata struct
153 by the pcreate function (ensures that corresponding 'free' and 'create' are always used)
154
155 @note IDA uses a different convention for function pointer types, so no '*'.
156 */
157 typedef struct IntegratorIdaPrecStruct{
158 IntegratorIdaPrecCreateFn *pcreate;
159 IDASpilsPrecSetupFn psetup;
160 IDASpilsPrecSolveFn psolve;
161 } IntegratorIdaPrec;
162
163 /* residual function forward declaration */
164 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
165
166 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
167 , N_Vector v, N_Vector Jv, realtype c_j
168 , void *jac_data, N_Vector tmp1, N_Vector tmp2
169 );
170
171 /* error handler forward declaration */
172 void integrator_ida_error(int error_code
173 , const char *module, const char *function
174 , char *msg, void *eh_data
175 );
176
177 /* dense jacobian evaluation for IDADense dense direct linear solver */
178 int integrator_ida_djex(long int Neq, realtype tt
179 , N_Vector yy, N_Vector yp, N_Vector rr
180 , realtype c_j, void *jac_data, DenseMat Jac
181 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
182 );
183
184 /* sparse jacobian evaluation for ASCEND's sparse direct solver */
185 IntegratorSparseJacFn integrator_ida_sjex;
186
187 typedef struct IntegratorIdaStatsStruct{
188 long nsteps;
189 long nrevals;
190 long nlinsetups;
191 long netfails;
192 int qlast, qcur;
193 realtype hinused, hlast, hcur;
194 realtype tcur;
195 } IntegratorIdaStats;
196
197 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
198 void integrator_ida_write_stats(IntegratorIdaStats *stats);
199 void integrator_ida_write_incidence(IntegratorSystem *blsys);
200 static int integrator_ida_check_lists(const IntegratorSystem *sys);
201
202 /*------
203 Jacobi preconditioner -- experimental
204 */
205
206 int integrator_ida_psetup_jacobi(realtype tt,
207 N_Vector yy, N_Vector yp, N_Vector rr,
208 realtype c_j, void *prec_data,
209 N_Vector tmp1, N_Vector tmp2,
210 N_Vector tmp3
211 );
212
213 int integrator_ida_psolve_jacobi(realtype tt,
214 N_Vector yy, N_Vector yp, N_Vector rr,
215 N_Vector rvec, N_Vector zvec,
216 realtype c_j, realtype delta, void *prec_data,
217 N_Vector tmp
218 );
219
220 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys);
221
222 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata);
223
224 static const IntegratorIdaPrec prec_jacobi = {
225 integrator_ida_pcreate_jacobi
226 , integrator_ida_psetup_jacobi
227 , integrator_ida_psolve_jacobi
228 };
229
230 /*-------------------------------------------------------------
231 SETUP/TEARDOWN ROUTINES
232 */
233 void integrator_ida_create(IntegratorSystem *blsys){
234 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
235 IntegratorIdaData *enginedata;
236 enginedata = ASC_NEW(IntegratorIdaData);
237 enginedata->rellist = NULL;
238 enginedata->safeeval = 0;
239 enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED;
240 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE;
241 enginedata->pfree = NULL;
242
243 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
244 enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
245
246 blsys->enginedata = (void *)enginedata;
247
248 integrator_ida_params_default(blsys);
249 }
250
251 void integrator_ida_free(void *enginedata){
252 CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
253 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
254 if(d->pfree){
255 /* free the preconditioner data, whatever it happens to be */
256 (d->pfree)(enginedata);
257 }
258 /* note, we don't own the rellist, so don't need to free it */
259 ASC_FREE(d);
260 }
261
262 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
263 IntegratorIdaData *d;
264 assert(blsys!=NULL);
265 assert(blsys->enginedata!=NULL);
266 assert(blsys->engine==INTEG_IDA);
267 d = ((IntegratorIdaData *)(blsys->enginedata));
268 return d;
269 }
270
271 /*-------------------------------------------------------------
272 PARAMETERS FOR IDA
273 */
274
275 enum ida_parameters{
276 IDA_PARAM_LINSOLVER
277 ,IDA_PARAM_MAXL
278 ,IDA_PARAM_AUTODIFF
279 ,IDA_PARAM_CALCIC
280 ,IDA_PARAM_SAFEEVAL
281 ,IDA_PARAM_RTOL
282 ,IDA_PARAM_ATOL
283 ,IDA_PARAM_ATOLVECT
284 ,IDA_PARAM_GSMODIFIED
285 ,IDA_PARAM_MAXNCF
286 ,IDA_PARAM_PREC
287 ,IDA_PARAMS_SIZE
288 };
289
290 /**
291 Here the full set of parameters is defined, along with upper/lower bounds,
292 etc. The values are stuck into the blsys->params structure.
293
294 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
295 list. Then add a slv_param_*(...) statement below to define the type, description and range
296 for the new parameter.
297
298 @return 0 on success
299 */
300 int integrator_ida_params_default(IntegratorSystem *blsys){
301 asc_assert(blsys!=NULL);
302 asc_assert(blsys->engine==INTEG_IDA);
303 slv_parameters_t *p;
304 p = &(blsys->params);
305
306 slv_destroy_parms(p);
307
308 if(p->parms==NULL){
309 CONSOLE_DEBUG("params NULL");
310 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
311 if(p->parms==NULL)return -1;
312 p->dynamic_parms = 1;
313 }else{
314 CONSOLE_DEBUG("params not NULL");
315 }
316
317 /* reset the number of parameters to zero so that we can check it at the end */
318 p->num_parms = 0;
319
320 slv_param_bool(p,IDA_PARAM_AUTODIFF
321 ,(SlvParameterInitBool){{"autodiff"
322 ,"Use auto-diff?",1
323 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
324 }, TRUE}
325 );
326
327 slv_param_char(p,IDA_PARAM_CALCIC
328 ,(SlvParameterInitChar){{"calcic"
329 ,"Initial conditions calcuation",1
330 ,"Use specified values of ydot to solve for inital y (Y),"
331 " or use the the values of the differential variables (yd) to solve"
332 " for the pure algebraic variables (ya) along with the derivatives"
333 " of the differential variables (yddot) (YA_YDP), or else don't solve"
334 " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"
335 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE",NULL}
336 );
337
338 slv_param_bool(p,IDA_PARAM_SAFEEVAL
339 ,(SlvParameterInitBool){{"safeeval"
340 ,"Use safe evaluation?",1
341 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
342 "throw SIGFPE errors which will then halt integration."
343 }, FALSE}
344 );
345
346
347 slv_param_bool(p,IDA_PARAM_ATOLVECT
348 ,(SlvParameterInitBool){{"atolvect"
349 ,"Use 'ode_atol' values as specified?",1
350 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
351 " in the integration. If FALSE, a scalar absolute tolerance value"
352 " is shared by all variables. See IDA manual, section 5.5.1"
353 }, TRUE }
354 );
355
356 slv_param_real(p,IDA_PARAM_ATOL
357 ,(SlvParameterInitReal){{"atol"
358 ,"Scalar absolute error tolerance",1
359 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
360 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
361 }, 1e-5, 1e-15, 1.0e15 }
362 );
363
364 slv_param_real(p,IDA_PARAM_RTOL
365 ,(SlvParameterInitReal){{"rtol"
366 ,"Scalar relative error tolerance",1
367 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
368 " it's not possible to set per-variable relative tolerances as it is"
369 " with LSODE)."
370 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
371 }, 1e-4, 0, 1 }
372 );
373
374 slv_param_char(p,IDA_PARAM_LINSOLVER
375 ,(SlvParameterInitChar){{"linsolver"
376 ,"Linear solver",1
377 ,"See IDA manual, section 5.5.3. Choose 'ASCEND' to use the linsolqr"
378 " direct linear solver bundled with ASCEND, 'DENSE' to use the dense"
379 " solver bundled with IDA, or one of the Krylov solvers SPGMR, SPBCG"
380 " or SPTFQMR (which still need preconditioners to be implemented"
381 " before they can be very useful."
382 }, "SPGMR"}, (char *[]){"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
383 );
384
385 slv_param_int(p,IDA_PARAM_MAXL
386 ,(SlvParameterInitInt){{"maxl"
387 ,"Maximum Krylov dimension",0
388 ,"The maximum dimension of Krylov space used by the linear solver"
389 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
390 " The default of 0 results in IDA using its internal default, which"
391 " is currently a value of 5."
392 }, 0, 0, 20 }
393 );
394
395 slv_param_bool(p,IDA_PARAM_GSMODIFIED
396 ,(SlvParameterInitBool){{"gsmodified"
397 ,"Gram-Schmidt Orthogonalisation Scheme", 2
398 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
399 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
400 }, TRUE}
401 );
402
403 slv_param_int(p,IDA_PARAM_MAXNCF
404 ,(SlvParameterInitInt){{"maxncf"
405 ,"Max nonlinear solver convergence failures per step", 2
406 ,"Maximum number of allowable nonlinear solver convergence failures"
407 " on one step. See IDA manual section 5.5.6.1."
408 }, 10,0,1000 }
409 );
410
411 slv_param_char(p,IDA_PARAM_PREC
412 ,(SlvParameterInitChar){{"prec"
413 ,"Preconditioner",1
414 ,"See IDA manual, section section 5.6.8."
415 },"NONE"}, (char *[]){"NONE","DIAG",NULL}
416 );
417
418 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
419
420 CONSOLE_DEBUG("Created %d params", p->num_parms);
421
422 return 0;
423 }
424
425 /*------------------------------------------------------------------------------
426 ANALYSIS ROUTINE (new implementation)
427 */
428
429 static int integrator_ida_block_check(IntegratorSystem *sys){
430 int res;
431 int dof;
432 #ifdef ANALYSE_DEBUG
433 long *vlist, *vp, i, nv, nv_ok;
434 char *varname;
435 struct var_variable **solversvars;
436
437 var_filter_t vfilt = {
438 VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED,
439 VAR_ACTIVE | VAR_INCIDENT | 0
440 };
441
442 nv = slv_get_num_solvers_vars(sys->system);
443 solversvars = slv_get_solvers_var_list(sys->system);
444 CONSOLE_DEBUG("-------------- nv = %ld -------------",nv);
445 for(nv_ok=0, i=0; i < nv; ++i){
446 if(var_apply_filter(solversvars[i],&vfilt)){
447 varname = var_make_name(sys->system,solversvars[i]);
448 fprintf(stderr,"%s\n",varname);
449 ASC_FREE(varname);
450 nv_ok++;
451 }
452 }
453 CONSOLE_DEBUG("----------- got %ld ok -------------",nv_ok);
454 #endif
455
456 if(!slvDOF_status(sys->system, &res, &dof)){
457 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
458 return -1;
459 }
460 switch(res){
461 case 1: CONSOLE_DEBUG("With derivatives fixed, system is underspecified (%d degrees of freedom)",dof);break;
462 case 2: CONSOLE_DEBUG("With derivatives fixed, system is square"); return 0; /* all OK */
463 case 3: CONSOLE_DEBUG("With derivatives fixed, system is structurally singular"); break;
464 case 4: CONSOLE_DEBUG("With derivatives fixed, system is overspecified"); break;
465 default:
466 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
467 return -2;
468 }
469
470 #ifdef ANALYSE_DEBUG
471 /* if it was underspecified, what vars could be fixed? */
472 if(res==1){
473 CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
474 solversvars = slv_get_solvers_var_list(sys->system);
475 if(!slvDOF_eligible(sys->system, &vlist)){
476 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
477 return -3;
478 }
479 for(vp=vlist;*vp!=-1;++vp){
480 varname = var_make_name(sys->system, solversvars[*vp]);
481 CONSOLE_DEBUG("Fixable var: %s",varname);
482 ASC_FREE(varname);
483 }
484 CONSOLE_DEBUG("(Found %ld fixable vars)",(int)(vp-vlist));
485 return 1;
486 }
487 #endif
488
489 return res;
490 }
491
492
493 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
494 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
495 IntegratorVarVisitorFn integrator_dae_classify_var;
496 void integrator_dae_show_var(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
497
498 #ifdef ASC_IDA_NEW_ANALYSE
499 /**
500 Perform additional problem analysis to prepare problem for integration with
501 IDA.
502
503 We assume that the analyse_generate_diffvars routine has been called, so
504 that we just need to call slv_get_diffvars to access the derivative
505 chains.
506
507 We can also assume that the independent variable has been found.
508
509 See mailing list, ~Jan 2007.
510
511 Note, the stuff for identifying the static and output sub-problems should
512 be part of integrator.c, not this file. We will assume this is handled
513
514 @return 0 on success
515 @see integrator_analyse
516 */
517 int integrator_ida_analyse(IntegratorSystem *sys){
518 /* struct var_variable **solversvars;
519 unsigned long nsolversvars;
520 struct rel_relation **solversrels;
521 unsigned long nsolversrels;*/
522 const SolverDiffVarCollection *diffvars;
523 SolverDiffVarSequence seq;
524 long i, j, n_y, n_ydot, n_dyn, n_skipped_diff, n_skipped_alg, n_skipped_deriv;
525 int res;
526
527 struct var_variable *v;
528 char *varname;
529
530 CONSOLE_DEBUG("NEW integrator_ida_analyse------------------>");
531
532 asc_assert(sys->engine==INTEG_IDA);
533
534 /* partition into static, dynamic and output problems */
535 CONSOLE_DEBUG("Block partitioning system...");
536 if(slv_block_partition(sys->system)){
537 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to block-partition system");
538 return 1;
539 }
540
541 diffvars = slv_get_diffvars(sys->system);
542 if(diffvars==NULL){
543 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
544 return 6;
545 }
546
547 CONSOLE_DEBUG("Got %ld chains, maxorder = %ld",diffvars->nseqs,diffvars->maxorder);
548
549 if(diffvars->maxorder > 2){
550 ERROR_REPORTER_HERE(ASC_USER_ERROR
551 ,"System higher-order derivatives. You must manually reduce the"
552 " system to a first-order system of DAEs. (maxorder=%d)"
553 ,diffvars->maxorder
554 );
555 return 1;
556 };
557
558 n_dyn = diffvars->nalg + 2 * diffvars->ndiff;
559
560 /* set up the dynamic problem */
561 CONSOLE_DEBUG("Setting up the dynamic problem");
562
563 asc_assert(sys->y == NULL);
564
565 sys->y = ASC_NEW_ARRAY(struct var_variable *,diffvars->nalg + diffvars->ndiff);
566 sys->ydot = ASC_NEW_ARRAY(struct var_variable *,diffvars->nalg + diffvars->ndiff);
567 sys->y_id = ASC_NEW_ARRAY(long,n_dyn);
568 n_y = 0;
569 n_skipped_alg = 0; n_skipped_diff = 0; n_skipped_deriv = 0;
570
571 /* initialise y_id to n_dyn (i.e. off limits) */
572 for(i=0;i<n_dyn;++i){
573 sys->y_id[i] = n_dyn;
574 }
575
576 /* add the variables from the derivative chains */
577 CONSOLE_DEBUG("Counting vars (%ld)",n_y);
578 for(i=0; i<diffvars->nseqs; ++i){
579 seq = diffvars->seqs[i];
580 asc_assert(seq.n >= 1);
581 if(seq.n == 1){
582 v = seq.vars[0];
583 varname = var_make_name(sys->system,v);
584 if(var_fixed(v)){
585 CONSOLE_DEBUG("'%s' is fixed",varname);
586 ASC_FREE(varname);
587 n_skipped_alg++;
588 continue;
589 }else if(!var_incident(v)){
590 CONSOLE_DEBUG("'%s' is not incident",varname);
591 ASC_FREE(varname);
592 n_skipped_alg++;
593 continue;
594 }
595 CONSOLE_DEBUG("'%s' is algebraic (%ld)",varname,n_y);
596 }else{
597 v = seq.vars[0];
598 varname = var_make_name(sys->system,v);
599 asc_assert(var_active(v));
600 if(var_fixed(v)){
601 CONSOLE_DEBUG("Differential var '%s' is fixed",varname);
602 n_skipped_diff++;
603 ASC_FREE(varname);
604 for(j=1; j<seq.n; ++j){
605 v = seq.vars[j];
606 varname = var_make_name(sys->system,v);
607 var_set_active(v,FALSE);
608 var_set_value(v,0);
609 CONSOLE_DEBUG("Derivative '%s' SET INACTIVE",varname);
610 ASC_FREE(varname);
611 n_skipped_deriv++;
612 }
613 continue;
614 }else if(var_fixed(seq.vars[1])){
615 /* diff var with fixed derivative */
616 CONSOLE_DEBUG("Derivative of var '%s' is fixed; CONVERTING TO ALGEBRAIC",varname);
617 n_skipped_deriv++;
618 }else{
619 /* seq.n > 1, var is not fixed, deriv is not fixed */
620 asc_assert(var_active(seq.vars[1]));
621 asc_assert(seq.n == 2);
622 sys->y[n_y] = v;
623 sys->y_id[var_sindex(v)] = n_y;
624 CONSOLE_DEBUG("'%s' is differential (%ld)",varname,n_y);
625 ASC_FREE(varname);
626 v = seq.vars[1];
627 sys->ydot[n_y] = v;
628 sys->y_id[var_sindex(v)] = -n_y-1;
629 varname = var_make_name(sys->system,v);
630 CONSOLE_DEBUG("'%s' is derivative (%ld, %p, %p)",varname,n_y,v,sys->ydot[-sys->y_id[var_sindex(v)]-1]);
631 ASC_FREE(varname);
632 n_y++;
633 continue;
634 }
635 }
636 /* fall through: v is algebraic */
637 ASC_FREE(varname);
638 asc_assert(var_active(v));
639 sys->y[n_y] = v;
640 sys->ydot[n_y] = NULL;
641 sys->y_id[var_sindex(v)]=n_y;
642 n_y++;
643 continue;
644 }
645 sys->n_y = n_y;
646 n_ydot = n_y;
647
648 /* FIX the derivatives */
649 CONSOLE_DEBUG("Checking system with derivatives fixed...");
650 for(i=0;i<n_y;++i){
651 if(sys->ydot[i])var_set_fixed(sys->ydot[i],1);
652 }
653
654 slv_block_partition(sys->system);
655 res = integrator_ida_block_check(sys);
656 if(res)return 100 + res;
657
658 /* FREE the derivatives, FIX the diffvars */
659 for(i=0;i<n_y;++i){
660 if(sys->ydot[i]){
661 var_set_fixed(sys->ydot[i],0);
662 var_set_fixed(sys->y[i],1);
663 }
664 }
665
666 CONSOLE_DEBUG("Checking system with differential variables fixed...");
667 slv_block_partition(sys->system);
668 res = integrator_ida_block_check(sys);
669 if(res)return 200 + res;
670
671 /* FREE the diffvars */
672 for(i=0;i<n_y;++i){
673 if(sys->ydot[i]){
674 var_set_fixed(sys->y[i],0);
675 }
676 }
677
678 /*if(sys->n_y != nsolversrels){
679 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square: n_y = %d, n_rels = %d"
680 " (n_alg = %d - %d = %d"
681 ", n_diff = 2 * %d - %d - %d)"
682 ,sys->n_y, nsolversrels
683 , diffvars->nalg, n_skipped_alg, diffvars->nalg - n_skipped_alg
684 , diffvars->ndiff, n_skipped_diff, n_skipped_deriv
685 );
686 return 2;
687 }*/
688
689 /* check the indep var is same as was located elsewhere */
690 if(diffvars->nindep>1){
691 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
692 return 3;
693 }else if(diffvars->nindep<1){
694 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
695 return 4;
696 }else if(diffvars->indep[0]!=sys->x){
697 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
698 return 5;
699 }
700
701 /* get the observations */
702 sys->n_obs = diffvars->nobs;
703 sys->obs = ASC_NEW_ARRAY(struct var_variable *,sys->n_obs);
704 for(i=0;i<sys->n_obs;++i){
705 sys->obs[i] = diffvars->obs[i];
706 varname = var_make_name(sys->system,sys->obs[i]);
707 CONSOLE_DEBUG("'%s' is observation",varname);
708 ASC_FREE(varname);
709 }
710
711 /* - 'y' list as [ya|yd] */
712 /* - sparsity pattern for dF/dy and dF/dy' */
713 /* - sparsity pattern for union of above */
714 /* - block decomposition based on above */
715 /* - block decomposition results in reordering of y and y' */
716 /* - boundaries (optional) */
717 /* ERROR_REPORTER_HERE(ASC_PROG_ERR,"Implementation incomplete");
718 return -1; */
719
720 if(integrator_ida_check_lists(sys)){
721 integrator_debug(sys,stderr);
722 }
723 return 0;
724 }
725 #endif
726
727 /*-------------------------------------------------------------
728 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
729 */
730
731 /*static double div1(double a, double b){
732 return a/b;
733 }*/
734
735 typedef int IdaFlagFn(void *,int *);
736 typedef char *IdaFlagNameFn(int);
737
738 /* return 0 on success */
739 int integrator_ida_solve(
740 IntegratorSystem *blsys
741 , unsigned long start_index
742 , unsigned long finish_index
743 ){
744 void *ida_mem;
745 int size, flag, flag1, t_index;
746 realtype t0, reltol, abstol, t, tret, tout1;
747 N_Vector y0, yp0, abstolvect, ypret, yret, id;
748 IntegratorIdaData *enginedata;
749 char *linsolver;
750 int maxl;
751 IdaFlagFn *flagfn;
752 IdaFlagNameFn *flagnamefn;
753 const char *flagfntype;
754 char *pname = NULL;
755 char *varname;
756 int i;
757 const IntegratorIdaPrec *prec = NULL;
758 int icopt; /* initial conditions strategy */
759
760 CONSOLE_DEBUG("STARTING IDA...");
761
762 enginedata = integrator_ida_enginedata(blsys);
763
764 enginedata->safeeval = SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_SAFEEVAL);
765 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
766
767 /* store reference to list of relations (in enginedata) */
768 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
769 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
770 enginedata->bndlist = slv_get_solvers_bnd_list(blsys->system);
771
772 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
773 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
774 size = blsys->n_y;
775
776 if(enginedata->nrels!=size){
777 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
778 return 1; /* failure */
779 }
780
781 if(integrator_ida_check_lists(blsys)){
782 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Inconsistent integrator lists");
783 return 1;
784 }
785
786 /* retrieve initial values from the system */
787
788 /** @TODO fix this, the starting time != first sample */
789 t0 = integrator_get_t(blsys);
790 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
791
792 CONSOLE_DEBUG("RETRIEVING y0");
793
794 y0 = N_VNew_Serial(size);
795 integrator_get_y(blsys,NV_DATA_S(y0));
796
797 #ifdef SOLVE_DEBUG
798 CONSOLE_DEBUG("RETRIEVING yp0");
799 #endif
800
801 yp0 = N_VNew_Serial(size);
802 integrator_get_ydot(blsys,NV_DATA_S(yp0));
803
804 #ifdef SOLVE_DEBUG
805 N_VPrint_Serial(yp0);
806 CONSOLE_DEBUG("yp0 is at %p",&yp0);
807 #endif
808
809 /* create IDA object */
810 ida_mem = IDACreate();
811
812 /* relative error tolerance */
813 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
814 CONSOLE_DEBUG("rtol = %8.2e",reltol);
815
816 /* allocate internal memory */
817 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
818 /* vector of absolute tolerances */
819 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
820 abstolvect = N_VNew_Serial(size);
821 integrator_get_atol(blsys,NV_DATA_S(abstolvect));
822
823 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
824
825 N_VDestroy_Serial(abstolvect);
826 }else{
827 /* scalar absolute tolerance (one value for all) */
828 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
829 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
830 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
831 }
832
833 if(flag==IDA_MEM_NULL){
834 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
835 return 2;
836 }else if(flag==IDA_MEM_FAIL){
837 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
838 return 3;
839 }else if(flag==IDA_ILL_INPUT){
840 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
841 return 4;
842 }/* else success */
843
844 /* set optional inputs... */
845 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
846 IDASetRdata(ida_mem, (void *)blsys);
847 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
848 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
849 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
850 if(integrator_get_minstep(blsys)>0){
851 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
852 }
853
854 CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
855 IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
856
857 /* there's no capability for setting *minimum* step size in IDA */
858
859
860 /* attach linear solver module, using the default value of maxl */
861 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
862 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
863 if(strcmp(linsolver,"ASCEND")==0){
864 CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",size);
865 IDAASCEND(ida_mem,size);
866 IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)blsys);
867
868 flagfntype = "IDAASCEND";
869 flagfn = &IDAASCENDGetLastFlag;
870 flagnamefn = &IDAASCENDGetReturnFlagName;
871
872 }else if(strcmp(linsolver,"DENSE")==0){
873 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
874 flag = IDADense(ida_mem, size);
875 switch(flag){
876 case IDADENSE_SUCCESS: break;
877 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 5;
878 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 5;
879 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 5;
880 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;
881 }
882
883 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
884 CONSOLE_DEBUG("USING AUTODIFF");
885 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
886 switch(flag){
887 case IDADENSE_SUCCESS: break;
888 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
889 }
890 }else{
891 CONSOLE_DEBUG("USING NUMERICAL DIFF");
892 }
893
894 flagfntype = "IDADENSE";
895 flagfn = &IDADenseGetLastFlag;
896 flagnamefn = &IDADenseGetReturnFlagName;
897 }else{
898 /* remaining methods are all SPILS */
899 CONSOLE_DEBUG("IDA SPILS");
900
901 maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);
902 CONSOLE_DEBUG("maxl = %d",maxl);
903
904 /* what preconditioner? */
905 pname = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_PREC);
906 if(strcmp(pname,"NONE")==0){
907 prec = NULL;
908 }else if(strcmp(pname,"JACOBI")==0){
909 prec = &prec_jacobi;
910 }else{
911 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
912 return 7;
913 }
914
915 /* which SPILS linear solver? */
916 if(strcmp(linsolver,"SPGMR")==0){
917 CONSOLE_DEBUG("IDA SPGMR");
918 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
919 }else if(strcmp(linsolver,"SPBCG")==0){
920 CONSOLE_DEBUG("IDA SPBCG");
921 flag = IDASpbcg(ida_mem, maxl);
922 }else if(strcmp(linsolver,"SPTFQMR")==0){
923 CONSOLE_DEBUG("IDA SPTFQMR");
924 flag = IDASptfqmr(ida_mem,maxl);
925 }else{
926 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
927 return 8;
928 }
929
930 if(prec){
931 /* assign the preconditioner to the linear solver */
932 (prec->pcreate)(blsys);
933 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)blsys);
934 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
935 }else{
936 CONSOLE_DEBUG("No preconditioner");
937 }
938
939 flagfntype = "IDASPILS";
940 flagfn = &IDASpilsGetLastFlag;
941 flagnamefn = &IDASpilsGetReturnFlagName;
942
943 if(flag==IDASPILS_MEM_NULL){
944 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
945 return 9;
946 }else if(flag==IDASPILS_MEM_FAIL){
947 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
948 return 9;
949 }/* else success */
950
951 /* assign the J*v function */
952 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
953 CONSOLE_DEBUG("USING AUTODIFF");
954 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
955 if(flag==IDASPILS_MEM_NULL){
956 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
957 return 10;
958 }else if(flag==IDASPILS_LMEM_NULL){
959 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
960 return 10;
961 }/* else success */
962 }else{
963 CONSOLE_DEBUG("USING NUMERICAL DIFF");
964 }
965
966 if(strcmp(linsolver,"SPGMR")==0){
967 /* select Gram-Schmidt orthogonalisation */
968 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_GSMODIFIED)){
969 CONSOLE_DEBUG("USING MODIFIED GS");
970 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
971 if(flag!=IDASPILS_SUCCESS){
972 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
973 return 11;
974 }
975 }else{
976 CONSOLE_DEBUG("USING CLASSICAL GS");
977 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
978 if(flag!=IDASPILS_SUCCESS){
979 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
980 return 11;
981 }
982 }
983 }
984 }
985
986 /* set linear solver optional inputs...
987 ...nothing here at the moment...
988 */
989
990 /* calculate initial conditions */
991 icopt = 0;
992 if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"Y")==0){
993 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
994 icopt = IDA_Y_INIT;
995 asc_assert(icopt!=0);
996 }else if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
997 CONSOLE_DEBUG("Solving initial conditions using values of yd");
998 icopt = IDA_YA_YDP_INIT;
999 asc_assert(icopt!=0);
1000 id = N_VNew_Serial(blsys->n_y);
1001 for(i=0; i < blsys->n_y; ++i){
1002 if(blsys->ydot[i] == NULL){
1003 NV_Ith_S(id,i) = 0.0;
1004 varname = var_make_name(blsys->system,blsys->y[i]);
1005 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
1006 ASC_FREE(varname);
1007 }else{
1008 CONSOLE_DEBUG("y[%d] is differential",i);
1009 NV_Ith_S(id,i) = 1.0;
1010 }
1011 }
1012 IDASetId(ida_mem, id);
1013 N_VDestroy_Serial(id);
1014 }else{
1015 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
1016 }
1017
1018 if(icopt){
1019 blsys->currentstep=0;
1020 t_index=start_index + 1;
1021 tout1 = samplelist_get(blsys->samples, t_index);
1022
1023 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
1024
1025 #ifdef ASC_SIGNAL_TRAPS
1026 /* catch SIGFPE if desired to */
1027 if(enginedata->safeeval){
1028 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
1029 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1030 }else{
1031 #ifdef FEX_DEBUG
1032 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
1033 #endif
1034 Asc_SignalHandlerPushDefault(SIGFPE);
1035 }
1036 if (setjmp(g_fpe_env)==0) {
1037 #endif
1038
1039 //CONSOLE_DEBUG("Raising signal...");
1040 //CONSOLE_DEBUG("1/0 = %f", div1(1.0,0.0));
1041 //CONSOLE_DEBUG("Still here...");
1042
1043 /* correct initial values, given derivatives */
1044 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
1045 /* note the new API from version 2.3 and onwards */
1046 flag = IDACalcIC(ida_mem, icopt, tout1);
1047 # else
1048 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
1049 # endif
1050
1051 switch(flag){
1052 case IDA_SUCCESS:
1053 CONSOLE_DEBUG("Initial conditions solved OK");
1054 break;
1055
1056 case IDA_LSETUP_FAIL:
1057 case IDA_LINIT_FAIL:
1058 case IDA_LSOLVE_FAIL:
1059 case IDA_NO_RECOVERY:
1060 flag1 = -999;
1061 flag = (flagfn)(ida_mem,&flag1);
1062 if(flag){
1063 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to retrieve error code from %s (err %d)",flagfntype,flag);
1064 return 12;
1065 }
1066 ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);
1067 return 12;
1068
1069 default:
1070 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
1071 return 12;
1072 }
1073 #ifdef ASC_SIGNAL_TRAPS
1074 }else{
1075 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
1076 return 13;
1077 }
1078
1079 if(enginedata->safeeval){
1080 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1081 }else{
1082 CONSOLE_DEBUG("pop...");
1083 Asc_SignalHandlerPopDefault(SIGFPE);
1084 CONSOLE_DEBUG("...pop");
1085 }
1086 #endif
1087
1088 }
1089
1090 /* optionally, specify ROO-FINDING PROBLEM */
1091
1092 /* -- set up the IntegratorReporter */
1093 integrator_output_init(blsys);
1094
1095 /* -- store the initial values of all the stuff */
1096 integrator_output_write(blsys);
1097 integrator_output_write_obs(blsys);
1098
1099 /* specify where the returned values should be stored */
1100 yret = y0;
1101 ypret = yp0;
1102
1103 /* advance solution in time, return values as yret and derivatives as ypret */
1104 blsys->currentstep=1;
1105 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
1106 t = samplelist_get(blsys->samples, t_index);
1107 t0 = integrator_get_t(blsys);
1108 asc_assert(t > t0);
1109
1110 #ifdef SOLVE_DEBUG
1111 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
1112 #endif
1113
1114 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
1115
1116 /* pass the values of everything back to the compiler */
1117 integrator_set_t(blsys, (double)tret);
1118 integrator_set_y(blsys, NV_DATA_S(yret));
1119 integrator_set_ydot(blsys, NV_DATA_S(ypret));
1120
1121 if(flag<0){
1122 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
1123 break;
1124 }
1125
1126 /* -- do something so that blsys knows the values of tret, yret and ypret */
1127
1128 /* -- store the current values of all the stuff */
1129 integrator_output_write(blsys);
1130 integrator_output_write_obs(blsys);
1131 }
1132
1133 /* -- close the IntegratorReporter */
1134 integrator_output_close(blsys);
1135
1136 /* get optional outputs */
1137 #ifdef STATS_DEBUG
1138 IntegratorIdaStats stats;
1139 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
1140 integrator_ida_write_stats(&stats);
1141 }
1142 #endif
1143
1144 /* free solution memory */
1145 N_VDestroy_Serial(yret);
1146 N_VDestroy_Serial(ypret);
1147
1148 /* free solver memory */
1149 IDAFree(ida_mem);
1150
1151 if(flag < 0){
1152 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1153 return 14;
1154 }
1155
1156 /* all done, success */
1157 return 0;
1158 }
1159
1160 /*--------------------------------------------------
1161 RESIDUALS AND JACOBIAN
1162 */
1163 /**
1164 Function to evaluate system residuals, in the form required for IDA.
1165
1166 Given tt, yy and yp, we need to evaluate and return rr.
1167
1168 @param tt current value of indep variable (time)
1169 @param yy current values of dependent variable vector
1170 @param yp current values of derivatives of dependent variables
1171 @param rr the output residual vector (is we're returning data to)
1172 @param res_data pointer to our stuff (blsys in this case).
1173
1174 @return 0 on success, positive on recoverable error, and
1175 negative on unrecoverable error.
1176 */
1177 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
1178 IntegratorSystem *blsys;
1179 IntegratorIdaData *enginedata;
1180 int i, calc_ok, is_error;
1181 struct rel_relation** relptr;
1182 double resid;
1183 char *relname;
1184 #ifdef FEX_DEBUG
1185 char *varname;
1186 char diffname[30];
1187 #endif
1188
1189 blsys = (IntegratorSystem *)res_data;
1190 enginedata = integrator_ida_enginedata(blsys);
1191
1192 #ifdef FEX_DEBUG
1193 /* fprintf(stderr,"\n\n"); */
1194 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
1195 #endif
1196
1197 /* pass the values of everything back to the compiler */
1198 integrator_set_t(blsys, (double)tt);
1199 integrator_set_y(blsys, NV_DATA_S(yy));
1200 integrator_set_ydot(blsys, NV_DATA_S(yp));
1201
1202 if(NV_LENGTH_S(rr)!=enginedata->nrels){
1203 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
1204 return -1; /* unrecoverable */
1205 }
1206
1207 /**
1208 @TODO does this function (fex) do bounds checking already?
1209 */
1210
1211 /* evaluate each residual in the rellist */
1212 is_error = 0;
1213 relptr = enginedata->rellist;
1214
1215 #ifdef ASC_SIGNAL_TRAPS
1216 if(enginedata->safeeval){
1217 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1218 }else{
1219 # ifdef FEX_DEBUG
1220 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
1221 # endif
1222 Asc_SignalHandlerPushDefault(SIGFPE);
1223 }
1224
1225 if (SETJMP(g_fpe_env)==0) {
1226 #endif
1227 for(i=0, relptr = enginedata->rellist;
1228 i< enginedata->nrels && relptr != NULL;
1229 ++i, ++relptr
1230 ){
1231 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
1232
1233 NV_Ith_S(rr,i) = resid;
1234 if(!calc_ok){
1235 relname = rel_make_name(blsys->system, *relptr);
1236 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1237 ASC_FREE(relname);
1238 /* presumable some output already made? */
1239 is_error = 1;
1240 }/*else{
1241 CONSOLE_DEBUG("Calc OK");
1242 }*/
1243 }
1244
1245 #ifdef ASC_SIGNAL_TRAPS
1246 }else{
1247 relname = rel_make_name(blsys->system, *relptr);
1248 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1249 ASC_FREE(relname);
1250 is_error = 1;
1251 }
1252
1253 if(enginedata->safeeval){
1254 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1255 }else{
1256 Asc_SignalHandlerPopDefault(SIGFPE);
1257 }
1258 #endif
1259
1260 #ifdef FEX_DEBUG
1261 /* output residuals to console */
1262 CONSOLE_DEBUG("RESIDUAL OUTPUT");
1263 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1264 for(i=0; i<blsys->n_y; ++i){
1265 varname = var_make_name(blsys->system,blsys->y[i]);
1266 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1267 if(blsys->ydot[i]){
1268 varname = var_make_name(blsys->system,blsys->ydot[i]);
1269 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1270 }else{
1271 snprintf(diffname,99,"diff(%s)",varname);
1272 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1273 }
1274 ASC_FREE(varname);
1275 relname = rel_make_name(blsys->system,enginedata->rellist[i]);
1276 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1277 }
1278 #endif
1279
1280 if(is_error){
1281 return 1;
1282 }
1283
1284 #ifdef FEX_DEBUG
1285 CONSOLE_DEBUG("RESIDUAL OK");
1286 #endif
1287 return 0;
1288 }
1289
1290 /**
1291 Dense Jacobian evaluation. Only suitable for small problems!
1292 */
1293 int integrator_ida_djex(long int Neq, realtype tt
1294 , N_Vector yy, N_Vector yp, N_Vector rr
1295 , realtype c_j, void *jac_data, DenseMat Jac
1296 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1297 ){
1298 IntegratorSystem *blsys;
1299 IntegratorIdaData *enginedata;
1300 char *relname;
1301 #ifdef DJEX_DEBUG
1302 struct var_variable **varlist;
1303 char *varname;
1304 #endif
1305 int status;
1306 struct rel_relation **relptr;
1307 int i;
1308 double *derivatives;
1309 int *variables;
1310 int count, j;
1311 long var_yindex;
1312
1313 blsys = (IntegratorSystem *)jac_data;
1314 enginedata = integrator_ida_enginedata(blsys);
1315
1316 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1317 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
1318 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1319
1320 /* pass the values of everything back to the compiler */
1321 integrator_set_t(blsys, (double)tt);
1322 integrator_set_y(blsys, NV_DATA_S(yy));
1323 integrator_set_ydot(blsys, NV_DATA_S(yp));
1324
1325 #ifdef DJEX_DEBUG
1326 varlist = slv_get_solvers_var_list(blsys->system);
1327
1328 /* print vars */
1329 for(i=0; i < blsys->n_y; ++i){
1330 varname = var_make_name(blsys->system, blsys->y[i]);
1331 CONSOLE_DEBUG("%s = %f = %f",varname,NV_Ith_S(yy,i),var_value(blsys->y[i]));
1332 ASC_FREE(varname);
1333 }
1334
1335 /* print derivatives */
1336 for(i=0; i < blsys->n_y; ++i){
1337 if(blsys->ydot[i]){
1338 varname = var_make_name(blsys->system, blsys->ydot[i]);
1339 CONSOLE_DEBUG("%s = %f =%f",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
1340 ASC_FREE(varname);
1341 }else{
1342 varname = var_make_name(blsys->system, blsys->y[i]);
1343 CONSOLE_DEBUG("diff(%s) = %f",varname,NV_Ith_S(yp,i));
1344 ASC_FREE(varname);
1345 }
1346 }
1347
1348 /* print step size */
1349 CONSOLE_DEBUG("<c_j> = %f",c_j);
1350 #endif
1351
1352 /* build up the dense jacobian matrix... */
1353 status = 0;
1354 for(i=0, relptr = enginedata->rellist;
1355 i< enginedata->nrels && relptr != NULL;
1356 ++i, ++relptr
1357 ){
1358 /* get derivatives for this particular relation */
1359 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1360
1361 if(status){
1362 relname = rel_make_name(blsys->system, *relptr);
1363 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1364 ASC_FREE(relname);
1365 break;
1366 }
1367
1368 /* output what's going on here ... */
1369 #ifdef DJEX_DEBUG
1370 relname = rel_make_name(blsys->system, *relptr);
1371 CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
1372 fprintf(stderr,"%d: '%s': ",i,relname);
1373 ASC_FREE(relname);
1374 for(j=0;j<count;++j){
1375 varname = var_make_name(blsys->system, varlist[variables[j]]);
1376 var_yindex = blsys->y_id[variables[j]];
1377 if(var_yindex >=0){
1378 fprintf(stderr," var[%d]='%s'=y[%ld]",variables[j],varname,var_yindex);
1379 }else{
1380 fprintf(stderr," var[%d]='%s'=ydot[%ld]",variables[j],varname,-var_yindex-1);
1381 }
1382 ASC_FREE(varname);
1383 }
1384 fprintf(stderr,"\n");
1385 #endif
1386 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1387 asc_assert(blsys->y_id);
1388 for(j=0; j < count; ++j){
1389 var_yindex = blsys->y_id[variables[j]];
1390 /* the SUNDIALS headers seem not to store 'N' on Windows */
1391 ASC_ASSERT_RANGE(var_yindex, -blsys->n_y, blsys->n_y);
1392
1393 if(var_yindex >= 0){
1394 #ifdef DJEX_DEBUG
1395 varname = var_make_name(blsys->system, blsys->y[var_yindex]);
1396 CONSOLE_DEBUG("blsys->y[var_yindex] = '%s'",varname);
1397 ASC_FREE(varname);
1398 varname = var_make_name(blsys->system, varlist[variables[j]]);
1399 CONSOLE_DEBUG("varlist[variables[j]] = '%s'",varname);
1400 ASC_FREE(varname);
1401 ASC_ASSERT_EQ_PTR(blsys->y[var_yindex],varlist[variables[j]]);
1402 #endif
1403 DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
1404 }else{
1405 #ifdef DJEX_DEBUG
1406 asc_assert(blsys->ydot[-var_yindex-1]==varlist[variables[j]]);
1407 #endif
1408 DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
1409 }
1410 }
1411 }
1412
1413 #ifdef DJEX_DEBUG
1414 CONSOLE_DEBUG("PRINTING JAC");
1415 fprintf(stderr,"\t");
1416 for(j=0; j < blsys->n_y; ++j){
1417 if(j)fprintf(stderr,"\t");
1418 varname = var_make_name(blsys->system,blsys->y[j]);
1419 fprintf(stderr,"%11s",varname);
1420 ASC_FREE(varname);
1421 }
1422 fprintf(stderr,"\n");
1423 for(i=0; i < enginedata->nrels; ++i){
1424 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
1425 fprintf(stderr,"%s\t",relname);
1426 ASC_FREE(relname);
1427
1428 for(j=0; j < blsys->n_y; ++j){
1429 if(j)fprintf(stderr,"\t");
1430 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1431 }
1432 fprintf(stderr,"\n");
1433 }
1434 #endif
1435
1436 if(status){
1437 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1438 return 1;
1439 }
1440
1441 #ifdef DJEX_DEBUG
1442 CONSOLE_DEBUG("DJEX RETURNING 0");
1443 #endif
1444 return 0;
1445 }
1446
1447 /**
1448 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1449
1450 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1451
1452 @param tt current value of the independent variable (time, t)
1453 @param yy current value of the dependent variable vector, y(t).
1454 @param yp current value of y'(t).
1455 @param rr current value of the residual vector F(t, y, y').
1456 @param v the vector by which the Jacobian must be multiplied to the right.
1457 @param Jv the output vector computed
1458 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1459 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1460 @param tmp1 @see tmp2
1461 @param tmp2 (as well as tmp1) pointers to memory allocated for variables of type N_Vector for use here as temporary storage or work space.
1462 @return 0 on success
1463 */
1464 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1465 , N_Vector v, N_Vector Jv, realtype c_j
1466 , void *jac_data, N_Vector tmp1, N_Vector tmp2
1467 ){
1468 IntegratorSystem *blsys;
1469 IntegratorIdaData *enginedata;
1470 int i, j, is_error=0;
1471 struct rel_relation** relptr;
1472 char *relname;
1473 int status;
1474 double Jv_i;
1475 long var_yindex;
1476
1477 int *variables;
1478 double *derivatives;
1479 int count;
1480 struct var_variable **varlist;
1481 #ifdef JEX_DEBUG
1482
1483 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1484 #endif
1485
1486 blsys = (IntegratorSystem *)jac_data;
1487 enginedata = integrator_ida_enginedata(blsys);
1488 varlist = slv_get_solvers_var_list(blsys->system);
1489
1490 /* pass the values of everything back to the compiler */
1491 integrator_set_t(blsys, (double)tt);
1492 integrator_set_y(blsys, NV_DATA_S(yy));
1493 integrator_set_ydot(blsys, NV_DATA_S(yp));
1494 /* no real use for residuals (rr) here, I don't think? */
1495
1496 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1497
1498 i = NV_LENGTH_S(yy) * 2;
1499 #ifdef JEX_DEBUG
1500 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1501 #endif
1502 variables = ASC_NEW_ARRAY(int, i);
1503 derivatives = ASC_NEW_ARRAY(double, i);
1504
1505 /* evaluate the derivatives... */
1506 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1507
1508 #ifdef ASC_SIGNAL_TRAPS
1509 Asc_SignalHandlerPushDefault(SIGFPE);
1510 if (SETJMP(g_fpe_env)==0) {
1511 #endif
1512 for(i=0, relptr = enginedata->rellist;
1513 i< enginedata->nrels && relptr != NULL;
1514 ++i, ++relptr
1515 ){
1516 /* get derivatives for this particular relation */
1517 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1518 #ifdef JEX_DEBUG
1519 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1520 #endif
1521
1522 if(status){
1523 relname = rel_make_name(blsys->system, *relptr);
1524 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1525 ASC_FREE(relname);
1526 is_error = 1;
1527 break;
1528 }
1529
1530 /*
1531 Now we have the derivatives wrt each alg/diff variable in the
1532 present equation. variables[] points into the varlist. need
1533 a mapping from the varlist to the y and ydot lists.
1534 */
1535
1536 Jv_i = 0;
1537 asc_assert(blsys->y_id);
1538 for(j=0; j < count; ++j){
1539 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
1540 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1541 if(varname){
1542 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1543 ASC_FREE(varname);
1544 }else{
1545 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1546 }
1547 */
1548
1549 /* we don't calculate derivatives wrt indep var */
1550 asc_assert(variables[j]>=0);
1551 if(varlist[variables[j]] == blsys->x) continue;
1552
1553 var_yindex = blsys->y_id[variables[j]];
1554 #ifdef JEX_DEBUG
1555 CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %ld",j,variables[j],var_yindex);
1556 #endif
1557
1558 ASC_ASSERT_RANGE(-var_yindex-1, -NV_LENGTH_S(v),NV_LENGTH_S(v));
1559
1560 if(var_yindex >= 0){
1561 #ifdef JEX_DEBUG
1562 asc_assert(blsys->y[var_yindex]==varlist[variables[j]]);
1563 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%ld] = %f, v[%ld] = %f)\n", i
1564 , derivatives[j] * NV_Ith_S(v,var_yindex)
1565 , i, var_yindex, derivatives[j]
1566 , var_yindex, NV_Ith_S(v,var_yindex)
1567 );
1568 #endif
1569 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
1570 }else{
1571 #ifdef JEX_DEBUG
1572 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%ld] = %f, v[%ld] = %f)\n", i
1573 , derivatives[j] * NV_Ith_S(v,-var_yindex-1)
1574 , i, -var_yindex-1, derivatives[j]
1575 , -var_yindex-1, NV_Ith_S(v,-var_yindex-1)
1576 );
1577 #endif
1578 asc_assert(blsys->ydot[-var_yindex-1]==varlist[variables[j]]);
1579 Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;
1580 }
1581 }
1582
1583 NV_Ith_S(Jv,i) = Jv_i;
1584 #ifdef JEX_DEBUG
1585 CONSOLE_DEBUG("rel = %p",*relptr);
1586 relname = rel_make_name(blsys->system, *relptr);
1587 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1588 //ASC_FREE(relname);
1589 return 1;
1590 #endif
1591 }
1592 #ifdef ASC_SIGNAL_TRAPS
1593 }else{
1594 relname = rel_make_name(blsys->system, *relptr);
1595 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1596 ASC_FREE(relname);
1597 is_error = 1;
1598 }
1599 Asc_SignalHandlerPopDefault(SIGFPE);
1600 #endif
1601
1602 if(is_error){
1603 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1604 return 1;
1605 }
1606 return 0;
1607 }
1608
1609 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
1610 int integrator_ida_sjex(long int Neq, realtype tt
1611 , N_Vector yy, N_Vector yp, N_Vector rr
1612 , realtype c_j, void *jac_data, mtx_matrix_t Jac
1613 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1614 ){
1615 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1616 return -1;
1617 }
1618
1619 /*----------------------------------------------
1620 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1621 */
1622
1623 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys){
1624 IntegratorIdaData * enginedata =blsys->enginedata;
1625 IntegratorIdaPrecDataJacobi *precdata;
1626 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1627
1628 asc_assert(blsys->n_y);
1629 precdata->PIii = N_VNew_Serial(blsys->n_y);
1630
1631 enginedata->pfree = &integrator_ida_pfree_jacobi;
1632 enginedata->precdata = precdata;
1633 CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1634 }
1635
1636 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1637 if(enginedata->precdata){
1638 IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1639 N_VDestroy_Serial(precdata->PIii);
1640
1641 ASC_FREE(precdata);
1642 enginedata->precdata = NULL;
1643 CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1644 }
1645 enginedata->pfree = NULL;
1646 }
1647
1648 /**
1649 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1650
1651 'setup' function.
1652 */
1653 int integrator_ida_psetup_jacobi(realtype tt,
1654 N_Vector yy, N_Vector yp, N_Vector rr,
1655 realtype c_j, void *p_data,
1656 N_Vector tmp1, N_Vector tmp2,
1657 N_Vector tmp3
1658 ){
1659 int i, j;
1660 IntegratorSystem *blsys;
1661 IntegratorIdaData *enginedata;
1662 IntegratorIdaPrecDataJacobi *precdata;
1663 struct rel_relation **relptr;
1664
1665 blsys = (IntegratorSystem *)p_data;
1666 enginedata = blsys->enginedata;
1667 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1668 double *derivatives;
1669 int *variables;
1670 int count, status;
1671 char *relname;
1672 int var_yindex;
1673
1674 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1675
1676 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
1677 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1678
1679 /**
1680 @TODO FIXME here we are using the very inefficient and contorted approach
1681 of calculating the whole jacobian, then extracting just the diagonal elements.
1682 */
1683
1684 for(i=0, relptr = enginedata->rellist;
1685 i< enginedata->nrels && relptr != NULL;
1686 ++i, ++relptr
1687 ){
1688
1689 /* get derivatives for this particular relation */
1690 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1691 if(status){
1692 relname = rel_make_name(blsys->system, *relptr);
1693 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1694 ASC_FREE(relname);
1695 break;
1696 }
1697 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1698 /* find the diagonal elements */
1699 asc_assert(blsys->y_id);
1700 for(j=0; j<count; ++j){
1701 if(variables[j]==i){
1702 var_yindex = blsys->y_id[variables[j]];
1703 if(var_yindex >= 0){
1704 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1705 }else{
1706 NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
1707 }
1708 }
1709 }
1710 #ifdef PREC_DEBUG
1711 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1712 #endif
1713 }
1714
1715 if(status){
1716 CONSOLE_DEBUG("Error found when evaluating derivatives");
1717 return 1; /* recoverable */
1718 }
1719
1720 integrator_ida_write_incidence(blsys);
1721
1722 ASC_FREE(variables);
1723 ASC_FREE(derivatives);
1724
1725 return 0;
1726 };
1727
1728 /**
1729 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1730
1731 'solve' function.
1732 */
1733 int integrator_ida_psolve_jacobi(realtype tt,
1734 N_Vector yy, N_Vector yp, N_Vector rr,
1735 N_Vector rvec, N_Vector zvec,
1736 realtype c_j, realtype delta, void *p_data,
1737 N_Vector tmp
1738 ){
1739 IntegratorSystem *blsys;
1740 IntegratorIdaData *data;
1741 IntegratorIdaPrecDataJacobi *precdata;
1742 blsys = (IntegratorSystem *)p_data;
1743 data = blsys->enginedata;
1744 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1745
1746 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
1747 N_VProd(precdata->PIii, rvec, zvec);
1748 return 0;
1749 };
1750
1751 /*----------------------------------------------
1752 STATS
1753 */
1754
1755 /**
1756 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1757 status in a struct instead of separately.
1758
1759 @return IDA_SUCCESS on sucess.
1760 */
1761 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
1762 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1763 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1764 ,&s->hlast, &s->hcur, &s->tcur
1765 );
1766 }
1767
1768 /**
1769 This routine just outputs the stats to the CONSOLE_DEBUG routine.
1770
1771 @TODO provide a GUI way of stats reporting from IDA.
1772 */
1773 void integrator_ida_write_stats(IntegratorIdaStats *stats){
1774 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
1775 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
1776 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
1777 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
1778 SI(qlast); SI(qcur);
1779 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
1780 # undef SL
1781 # undef SI
1782 # undef SR
1783 }
1784
1785 /*------------------------------------------------------------------------------
1786 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
1787 */
1788
1789 enum integrator_ida_write_jac_enum{
1790 II_WRITE_Y
1791 , II_WRITE_YDOT
1792 };
1793
1794 /**
1795 @TODO COMPLETE THIS...
1796 */
1797 void integrator_ida_write_jacobian(IntegratorSystem *blsys, realtype c_j, FILE *f, enum integrator_ida_write_jac_enum type){
1798 IntegratorIdaData *enginedata;
1799 MM_typecode matcode;
1800 int nnz, rhomax;
1801 double *derivatives;
1802 int *variables;
1803 struct rel_relation **relptr;
1804 int i, j, status, count, var_yindex;
1805 char *relname;
1806
1807 var_filter_t vfiltery = {
1808 VAR_SVAR | VAR_FIXED | VAR_DERIV
1809 , VAR_SVAR
1810 };
1811 var_filter_t vfilteryd = {
1812 VAR_SVAR | VAR_FIXED | VAR_DERIV
1813 , VAR_SVAR | VAR_DERIV
1814 };
1815
1816 enginedata = (IntegratorIdaData *)blsys->enginedata;
1817
1818 /* number of non-zeros for all the non-FIXED solver_vars,
1819 in all the active included equality relations.
1820 */
1821 nnz = relman_jacobian_count(enginedata->rellist, enginedata->nrels
1822 , &enginedata->vfilter, &enginedata->rfilter
1823 , &rhomax
1824 );
1825
1826 /* we must have found the same number of relations */
1827 asc_assert(rhomax == enginedata->nrels);
1828
1829 /* output the mmio file header, now that we know our size*/
1830 /* note that we are asserting that our problem is square */
1831 mm_initialize_typecode(&matcode);
1832 mm_set_matrix(&matcode);
1833 mm_set_coordinate(&matcode);
1834 mm_set_real(&matcode);
1835 mm_write_banner(f, matcode);
1836 mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);
1837
1838 variables = ASC_NEW_ARRAY(int, blsys->n_y * 2);
1839 derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1840
1841 CONSOLE_DEBUG("Writing sparse Jacobian to file...");
1842
1843 for(i=0, relptr = enginedata->rellist;
1844 i< enginedata->nrels && relptr != NULL;
1845 ++i, ++relptr
1846 ){
1847 relname = rel_make_name(blsys->system, *relptr);
1848
1849 /* get derivatives of y */
1850 status = relman_diff2(*relptr, &vfiltery, derivatives, variables, &count, enginedata->safeeval);
1851 if(status){
1852 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1853 ASC_FREE(relname);
1854 break;
1855 }
1856
1857 /* get derivatives of y */
1858 status = relman_diff2(*relptr, &vfilteryd, derivatives, variables, &count, enginedata->safeeval);
1859 if(status){
1860 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1861 ASC_FREE(relname);
1862 break;
1863 }
1864
1865 asc_assert(blsys->y_id);
1866 for(j=0; j<count; ++j){
1867 var_yindex = blsys->y_id[variables[j]];
1868 if(var_yindex >= 0 && type == II_WRITE_Y){
1869 fprintf(f, "%d %d %10.3g\n", i, var_yindex, derivatives[j]);
1870 }else if(var_yindex < 0 && type == II_WRITE_YDOT){
1871 fprintf(f, "%d %d %10.3g\n", i, -var_yindex-1, derivatives[j]);
1872 }
1873 }
1874 }
1875 ASC_FREE(variables);
1876 ASC_FREE(derivatives);
1877 }
1878
1879 /**
1880 This routine outputs matrix structure in a crude text format, for the sake
1881 of debugging.
1882 */
1883 void integrator_ida_write_incidence(IntegratorSystem *blsys){
1884 int i, j;
1885 struct rel_relation **relptr;
1886 IntegratorIdaData *enginedata = blsys->enginedata;
1887 double *derivatives;
1888 int *variables;
1889 int count, status;
1890 char *relname;
1891 int var_yindex;
1892
1893 if(enginedata->nrels > 100){
1894 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
1895 return;
1896 }
1897
1898 variables = ASC_NEW_ARRAY(int, blsys->n_y * 2);
1899 derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1900
1901 CONSOLE_DEBUG("Outputting incidence information to console...");
1902
1903 for(i=0, relptr = enginedata->rellist;
1904 i< enginedata->nrels && relptr != NULL;
1905 ++i, ++relptr
1906 ){
1907 relname = rel_make_name(blsys->system, *relptr);
1908
1909 /* get derivatives for this particular relation */
1910 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1911 if(status){
1912 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1913 ASC_FREE(relname);
1914 break;
1915 }
1916
1917 fprintf(stderr,"%3d:%-15s:",i,relname);
1918 ASC_FREE(relname);
1919
1920 asc_assert(blsys->y_id);
1921 for(j=0; j<count; ++j){
1922 var_yindex = blsys->y_id[variables[j]];
1923 if(var_yindex >= 0){
1924 fprintf(stderr," %d:y[%d]",variables[j],var_yindex);
1925 }else{
1926 fprintf(stderr," %d:ydot[%d]",variables[j],-var_yindex-1);
1927 }
1928 }
1929 fprintf(stderr,"\n");
1930 }
1931 ASC_FREE(variables);
1932 ASC_FREE(derivatives);
1933 }
1934
1935 /* @return 0 on success */
1936 int integrator_ida_debug(const IntegratorSystem *sys, FILE *fp){
1937 char *varname, *relname;
1938 struct var_variable **vlist, *var;
1939 struct rel_relation **rlist, *rel;
1940 long vlen, rlen;
1941 long i, y_id;
1942
1943 fprintf(fp,"THERE ARE %ld VARIABLES IN THE INTEGRATION SYSTEM\n\n",sys->n_y);
1944
1945 /* if(integrator_sort_obs_vars(sys))return 10; */
1946
1947 fprintf(fp,"RESULTS OF ANALYSIS\n\n");
1948 fprintf(fp,"index\ty\tydot\n");
1949 fprintf(fp,"-----\t-----\t-----\n");
1950 for(i=0;i<sys->n_y;++i){
1951 varname = var_make_name(sys->system, sys->y[i]);
1952 fprintf(fp,"%ld\t%s\t",i,varname);
1953 if(sys->ydot[i]){
1954 ASC_FREE(varname);
1955 varname = var_make_name(sys->system, sys->ydot[i]);
1956 fprintf(fp,"%s\n",varname);
1957 ASC_FREE(varname);
1958 }else{
1959 fprintf(fp,".\n");
1960 ASC_FREE(varname);
1961 }
1962 }
1963
1964 fprintf(fp,"\n\nCORRESPONDENCE OF SOLVER VARS TO INTEGRATOR VARS\n\n");
1965 fprintf(fp,"sindex\t%-15s\ty_id\ty\tydot\n","name");
1966 fprintf(fp,"------\t%-15s\t-----\t-----\t-----\n","----");
1967
1968
1969 /* visit all the slv_system_t master var lists to collect vars */
1970 /* find the vars mostly in this one */
1971 vlist = slv_get_solvers_var_list(sys->system);
1972 vlen = slv_get_num_solvers_vars(sys->system);
1973 for(i=0;i<vlen;i++){
1974 var = vlist[i];
1975
1976 varname = var_make_name(sys->system, var);
1977 fprintf(fp,"%ld\t%-15s\t",i,varname);
1978
1979 if(var_fixed(var)){
1980 // it's fixed, so not really a DAE var
1981 fprintf(fp,"(fixed)\n");
1982 }else if(!var_active(var)){
1983 // inactive
1984 fprintf(fp,"(inactive)\n");
1985 }else{
1986 // a DAE var, otherwise
1987 y_id = sys->y_id[i];
1988
1989 fprintf(fp,"%ld", y_id);
1990 if(y_id >= 0){
1991 fprintf(fp,"\ty[%ld]\t.\n",y_id);
1992 }else{
1993 fprintf(fp,"\t.\tydot[%ld]\n",-y_id-1);
1994 }
1995
1996 ASC_ASSERT_LT(i,1e7L);
1997 ASC_ASSERT_LT(y_id, 9999999L);
1998 ASC_ASSERT_LT(-9999999L, y_id);
1999 }
2000 ASC_FREE(varname);
2001 }
2002
2003 /* let's write out the relations too */
2004 rlist = slv_get_solvers_rel_list(sys->system);
2005 rlen = slv_get_num_solvers_rels(sys->system);
2006
2007 fprintf(fp,"\nRELATIONS (%ld)\n\n",rlen);
2008 fprintf(fp,"index\tname\n");
2009 fprintf(fp,"-----\t----\n");
2010 for(i=0; i<rlen; ++i){
2011 rel = rlist[i];
2012 relname = rel_make_name(sys->system,rel);
2013 fprintf(fp,"%ld\t%s\n",i,relname);
2014 ASC_FREE(relname);
2015 }
2016
2017 /* and lets write block debug output */
2018 block_debug(sys->system, fp);
2019
2020 return 0; /* success */
2021 }
2022
2023 static int integrator_ida_check_lists(const IntegratorSystem *sys){
2024 long i, y_id;
2025 struct var_variable **svars;
2026 char *varname1, *varname2;
2027 long n;
2028 int err = 0;
2029 svars = slv_get_solvers_var_list(sys->system);
2030 n = slv_get_num_solvers_vars(sys->system);
2031 for(i=0; i<n; ++i){
2032 y_id = sys->y_id[i];
2033 if(y_id >= sys->n_y || y_id <= -sys->n_y - 1){
2034 CONSOLE_DEBUG("Skipping y_id[%ld] = %ld",i,y_id);
2035 continue;
2036 }
2037 varname1 = var_make_name(sys->system,svars[i]);
2038 if(y_id >= 0){
2039 if(sys->y[y_id]!=svars[i]){
2040 varname2 = var_make_name(sys->system,sys->y[y_id]);
2041 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error y: y_id[%ld]=%ld-->'%s', but svars[%ld]='%s'"
2042 ,i,y_id,varname2,i,varname1
2043 );
2044 ASC_FREE(varname2);
2045 err++;
2046 }
2047 }else if(sys->ydot[-y_id-1] != NULL){
2048 CONSOLE_DEBUG("y_id = %ld --> -y_id-1 = %ld --> sys->ydot[%ld] = %p"
2049 ,y_id, -y_id-1, -y_id-1, sys->ydot[-y_id-1]
2050 );
2051 varname2 = var_make_name(sys->system,sys->ydot[-y_id-1]);
2052 if(sys->ydot[-y_id-1]!=svars[i]){
2053 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error y: y_id[%ld]=%ld-->'%s', but svars[%ld]='%s'"
2054 ,i,-y_id-1,varname2,i,varname1
2055 );
2056 err++;
2057 }
2058 ASC_FREE(varname2);
2059 }
2060 ASC_FREE(varname1);
2061 }
2062 return err;
2063 }
2064
2065 /*----------------------------------------------
2066 ERROR REPORTING
2067 */
2068 /**
2069 Error message reporter function to be passed to IDA. All error messages
2070 will trigger a call to this function, so we should find everything
2071 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
2072 panel (in the case of PyGTK).
2073 */
2074 void integrator_ida_error(int error_code
2075 , const char *module, const char *function
2076 , char *msg, void *eh_data
2077 ){
2078 IntegratorSystem *blsys;
2079 error_severity_t sev;
2080
2081 /* cast back the IntegratorSystem, just in case we need it */
2082 blsys = (IntegratorSystem *)eh_data;
2083
2084 /* severity depends on the sign of the error_code value */
2085 if(error_code <= 0){
2086 sev = ASC_PROG_ERR;
2087 }else{
2088 sev = ASC_PROG_WARNING;
2089 }
2090
2091 /* use our all-purpose error reporting to get stuff back to the GUI */
2092 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
2093 }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22