/[ascend]/branches/harry/solvers/ida/ida.c
ViewVC logotype

Diff of /branches/harry/solvers/ida/ida.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3031 by raymondj, Fri May 22 08:39:50 2015 UTC revision 3032 by raymondj, Wed Jul 29 03:30:04 2015 UTC
# Line 1  Line 1 
1  /*  ASCEND modelling environment  /*  ASCEND modelling environment
2   Copyright (C) 2006-2011 Carnegie Mellon University      Copyright (C) 2011 Carnegie Mellon University
3    
4   This program is free software; you can redistribute it and/or modify      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)      it under the terms of the GNU General Public License as published by
7   any later version.      the Free Software Foundation; either version 2, or (at your option)
8        any later version.
9   This program is distributed in the hope that it will be useful,  
10   but WITHOUT ANY WARRANTY; without even the implied warranty of      This program is distributed in the hope that it will be useful,
11   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      but WITHOUT ANY WARRANTY; without even the implied warranty of
12   GNU General Public License for more details.      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    
14   You should have received a copy of the GNU General Public License      GNU General Public License for more details.
15   along with this program.  If not, see <http://www.gnu.org/licenses/>.  
16  *//**      You should have received a copy of the GNU General Public License
17   @file      along with this program; if not, write to the Free Software
18   Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes      Foundation, Inc., 59 Temple Place - Suite 330,
19   as part of the GPL-licensed SUNDIALS solver package from LLNL.      Boston, MA 02111-1307, USA.
20    
21   IDA provides the non-linear parts, as well as a number of pluggable linear  *//*
22   solvers: dense, banded and krylov types.      by Harry, 4th June 2015
23        This file contains the ida_main_solve function.
24   We also implement here an EXPERIMENTAL direct sparse linear solver for IDA      This is the main solver function that calls other
25   using the ASCEND linsolqr routines.      sub-functions to complete integration.
26        The current outline of this function is as below:
27   @see http://www.llnl.gov/casc/sundials/      Boundary flags initialisation: all the condition equations are evaluated
28   *//*          and the values of corresponding boundary flags are set.
29   by John Pye, May 2006      Boolean variables initialisation: the system is solved using the logical
30   */          equations solver in order to get the correct values of Boolean variables.
31        Using the logical values found at step 2 the configuration of the system is chosen.
32  #define _GNU_SOURCE      Solver initialisation: setting solver parameters, allocating memory, et cetera.
33        Loop through all time steps:
34  #include "ida.h"          If a root has not been found at the previous iteration then then reset all
35  #include "idalinear.h"              the root direction settings. It means that IDA will be searching for
36  #include "idaanalyse.h"              roots of condition equations in both directions.
37  #include "idatypes.h"          If a root has been found and we are far from desired output then make a very
38  #include "idaprec.h"              small time step and solve the system with IDA. Else make a usual time
39  #include "idacalc.h"              step which is set by the user.
40  #include "idaio.h"          If we have made a small time step then call a logical solver in order to update
41  #include "idaboundary.h"              the values of logical variables.
42            Check if a root has been found at this step. If yes:
43  #include <signal.h>              Call a function for processing this boundary.
44  #include <setjmp.h>              If there have been any changes in the system, then reinitialise the solver.
45  #include <fenv.h>              Set the direction in which the roots of condition equations will be detected:
46  #include <math.h>                  if a boundary has been crossed, then at the next step (which is an
47                    auxiliary small step) it will not be crossed in the same direction.
48  #ifdef ASC_WITH_MMIO              If we are still far from desired output then repeat steps.
49  # include <mmio.h>  
50  #endif      So basically we have two modes here: the first is when we make an ordinary time step, check all
51        roots and if a boundary is found then call the function for processing it. The second mode is entered.
52  #include <ascend/general/platform.h>  */
53  #include <ascend/utilities/error.h>  #define _GNU_SOURCE                            
54  #include <ascend/utilities/ascSignal.h>  #include "ida.h"
55  #include <ascend/general/panic.h>  #include "idalinear.h"
56  #include <ascend/compiler/instance_enum.h>  #include "idaanalyse.h"
57    #include "idatypes.h"  
58    #include "prepare_integrator.h"                                 /* The list of includes needs cleaning up! Once basic files are written,  
59  #include <ascend/system/slv_client.h>  #include "idaprec.h"                                        all commonly called functions can be added in just one header file.
60  #include <ascend/system/relman.h>  #include "idacalc.h"                                        Once that's done, these includes must be revisited*/
61  #include <ascend/system/block.h>  #include "idaio.h"
62  #include <ascend/system/slv_stdcalls.h>  #include "idaboundary.h"
63  #include <ascend/system/jacobian.h>  #include <signal.h>                                     /*Check if all these includes are necessary*/
64  #include <ascend/system/bndman.h>  #include <setjmp.h>
65    #include <fenv.h>
66  #include <ascend/utilities/config.h>  #include <math.h>
67  #include <ascend/integrator/integrator.h>  #ifdef ASC_WITH_MMIO
68    # include <mmio.h>
69  /* #define FEX_DEBUG */  #endif
70  /* #define SOLVE_DEBUG */  #include <ascend/general/platform.h>
71  #define IDA_BND_DEBUG  #include <ascend/utilities/error.h>
72  /* #define STATS_DEBUG */  #include <ascend/utilities/ascSignal.h>
73  /* #define DESTROY_DEBUG */  #include <ascend/general/panic.h>
74    #include <ascend/compiler/instance_enum.h>
75  /*-------------------------------------------------------------  #include <ascend/system/slv_client.h>
76   SOLVER REGISTRATION  #include <ascend/system/relman.h>
77   */  #include <ascend/system/block.h>
78    #include <ascend/system/slv_stdcalls.h>
79  /*  #include <ascend/system/jacobian.h>
80      The following functions are declared static because we don't want them  #include <ascend/system/bndman.h>
81      ever to be called directly from outside code. Instead, they are call via  #include <ascend/utilities/config.h>
82      function pointers passed through the ida_register function.  #include <ascend/integrator/integrator.h>
83  */  
84  static IntegratorCreateFn integrator_ida_create;  
85  static IntegratorParamsDefaultFn integrator_ida_params_default;  int integrate_ida_solve(IntegratorSystem *integ, unsigned long start_index, unsigned long finish_index){
86  static IntegratorSolveFn integrator_ida_solve;      void *ida_mem;
87  static IntegratorFreeFn integrator_ida_free;      int t_index;
88        realtype t0, tout, tret, tol;
89  /**      realtype tmin;              /** < The length of a small additional step made after an event is triggered */
90   This data structure contains pointers to the various functions that are      N_Vector ypret, yret;
91   exposed to libascend in order for ASCEND to drive this solver.      IntegratorIdaData *enginedata;
92   */      int i, flag;
93  static const IntegratorInternals integrator_ida_internals = {      int *rootsfound;            /** < IDA rootfinder reports root index in here */
94          integrator_ida_create, integrator_ida_params_default,      int *rootdir = NULL;                /** < Used to tell IDA to ignore doulve crossings */
95          integrator_ida_analyse, integrator_ida_solve,      int *bnd_cond_states;       /** < Record of boundary states so that IDA can tell LRSlv how to evaluate a boundary crossing */
96          integrator_ida_write_matrix, integrator_ida_debug, integrator_ida_free,      int *bnd_not_set;
97          INTEG_IDA, "IDA" };      int all_bnds_set = 1;
98        int need_to_reconfigure;    /** < Flag to indicate system rebuild after crossing */
99  /**      int need_to_reinteg = 0;    /** < Flag for when crossings happen on or very close to timesteps */
100   This function is accessed by libascend when loading this solver. The      int skipping_output;        /** < Flag to skip output to reporter */
101   function will register the integrator such that it can then be applied      int qrslv_ind, lrslv_ind;
102   to solving problems.      int after_root = 0;
103   */      int peaw = 0;                  /*Flag returned by process_events_and_whens*/
104  extern ASC_EXPORT int ida_register(void) {      int subpeaw = 0;        /*Sub-flag used to control auxiliay integration steps*/
105      //CONSOLE_DEBUG("Registering IDA...");      int auxcount = 1;
106      return integrator_register(&integrator_ida_internals);      int first_run = 1;
107  }  #if defined(SOLVE_DEBUG) || defined(IDA_BND_DEBUG)
108        char *relname;
109  /*-------------------------------------------------------------      //CONSOLE_DEBUG("STARTING IDA...");
110   FORWARD DECLS  #endif
111   */  #ifdef SOLVE_DEBUG
112        integrator_ida_debug(integ, stderr);
113  typedef void ( IntegratorVarVisitorFn)(IntegratorSystem *integ,  #endif
114          struct var_variable *var, const int *varindx);      /* create IDA object */
115        ida_mem = IDACreate();
116  /*static IntegratorVarVisitorFn integrator_dae_classify_var;      /* solve the initial conditions, allocate memory, other stuff... */
117   static void integrator_visit_system_vars(IntegratorSystem *integ,IntegratorVarVisitorFn *visitor);      ida_prepare_integrator(integ);  /*Change call to new function*/
118   static void integrator_dae_show_var(IntegratorSystem *integ, struct var_variable *var, const int *varindx); */      /* store reference to list of relations (in enginedata) */
119        ida_load_rellist(integ);
120  static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);      
121        tol = 0.0001*(samplelist_get(integ->samples, finish_index)
122  /*-------------------------------------------------------------                      - samplelist_get(integ->samples, start_index))
123   SETUP/TEARDOWN ROUTINES                  /samplelist_length(integ->samples);
124   */      tmin = tol;
125        rootdir = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);                  
126  /**      for(i = 0; i < enginedata->nbnds; i++){
127   Allocate memory inside the IntegratorIdaData structure at          rootdir[i] = 0;
128   the time when IDA is assigned to a particular system as its integrator.  #ifdef IDA_BND_DEBUG
129   */          char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);
130  static void integrator_ida_create(IntegratorSystem *integ) {          CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d]=%d, bndman_calc_satisfied=%d; (trigger dirn=%s)"
131      //CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");          ,n, i, bnd_cond_states[i], bndman_calc_satisfied(enginedata->bndlist[i])
132      IntegratorIdaData *enginedata;          ,rootdir[i]==1?"UP":(rootdir[i]==-1?"DOWN":"both")
133      enginedata = ASC_NEW(IntegratorIdaData);          );
134      //CONSOLE_DEBUG("enginedata = %p",enginedata);          ASC_FREE(n);
135      enginedata->rellist = NULL;  #endif
136      enginedata->safeeval = 0;      }
137      enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE      
138              | VAR_FIXED;  
139      enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0;      for(t_index = start_index + 1; t_index <= finish_index; ++t_index, ++integ->currentstep){
140      enginedata->pfree = NULL;          tout = samplelist_get(integ->samples, t_index);
141            t0 = integrator_get_t(integ);
142      enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;          asc_assert(tout > t0);
143      enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;  
144    #ifdef SOLVE_DEBUG
145      enginedata->flagfn = NULL;          CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, tout);
146      enginedata->flagfntype = NULL;  #endif          
147      enginedata->flagnamefn = NULL;          if(enginedata->nbnds || first_run){
148                peaw = process_events_and_whens(integ, ida_mem, t0, tout, rootdir, first_run);
149      integ->enginedata = (void *) enginedata;              if(first_run == 1){
150                    first_run = 0;
151      integrator_ida_params_default(integ);              }
152  }          }
153        
154  /**  
155   Free any memory specifically allocated for use by our ASCEND wrapper for          if(peaw==143){                              /*Flag for rootsfound - system is reconfigured already!*/
156   the IDA integrator, other than that which the IDA code itself may be  
157   managing.              do{
158   */                  subpeaw = 0;                            /*Flag for taking auxiliary small steps after root*/
159  static void integrator_ida_free(void *enginedata) {                  flag = IDASolve(ida_mem, t0 + auxcount*tmin, &t0, yret, ypret, IDA_NORMAL);
160  #ifdef DESTROY_DEBUG                  /*Todo: Error care for integrator*/
161      CONSOLE_DEBUG("DESTROYING IDA engine data at %p",enginedata);                  /*Now, take small timesteps and check for further roots*/
162  #endif                  for(i = 0; i < enginedata->nbnds; i++) {
163      IntegratorIdaData *d = (IntegratorIdaData *) enginedata;                      rootdir[i] = -1*rootsfound[i];
164      asc_assert(d);  #ifdef IDA_BND_DEBUG
165      if (d->pfree) {                      char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);
166          CONSOLE_DEBUG("DESTROYING preconditioner data using fn at %p",d->pfree);                      CONSOLE_DEBUG("Set direction=%d for boundary '%s'",rootdir[i],n);
167          /* free the preconditioner data, whatever it happens to be */                      ASC_FREE(n);
168          (d->pfree)(enginedata);  #endif
169      }                  }
170                    subpeaw = process_events_and_whens(integ, ida_mem, t0, t0 + auxcount*tmin, rootdir,first_run);
171      ASC_FREE(d->rellist);              }while(subpeaw == 143 && auxcount < 20);        /*Arbitrary limit. Still better limit: auxcount < (tout-t0)/tmin ?*/
172                
173  #ifdef DESTROY_DEBUG              /*Important Todo ---- Time needs to be reset*/
174      CONSOLE_DEBUG("Now destroying the enginedata");                          
175  #endif  
176      ASC_FREE(d);          }
177  #ifdef DESTROY_DEBUG  
178      CONSOLE_DEBUG("enginedata freed");          
179  #endif          if(peaw==657){
180  }              flag = IDASolve(ida_mem, tout, &t0, yret, ypret, IDA_NORMAL);  
181            }
182  IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *integ) {  
183      IntegratorIdaData *d;  
184      assert(integ!=NULL);      }/*End of Main integration For Loop*/
185      assert(integ->enginedata!=NULL);  
186      assert(integ->engine==INTEG_IDA);  /* -- set up the IntegratorReporter */
187      d = ((IntegratorIdaData *) (integ->enginedata));      integrator_output_init(integ);
188      return d;      /* -- store the initial values of all the stuff */
189  }      integrator_output_write(integ);
190        integrator_output_write_obs(integ);
191  /*-------------------------------------------------------------      /* specify where the returned values should be stored */
192   PARAMETERS FOR IDA      yret    = ida_bnd_new_zero_NV(integ->n_y);
193   */      ypret   = ida_bnd_new_zero_NV(integ->n_y);
194        /* advance solution in time, return values as yret and derivatives as ypret */
195  enum ida_parameters {      integ->currentstep = 1;
196      IDA_PARAM_LINSOLVER,  return 0;  
197      IDA_PARAM_MAXL,  
198      IDA_PARAM_MAXORD,  }
     IDA_PARAM_AUTODIFF,  
     IDA_PARAM_CALCIC,  
     IDA_PARAM_SAFEEVAL,  
     IDA_PARAM_RTOL,  
     IDA_PARAM_ATOL,  
     IDA_PARAM_ATOLVECT,  
     IDA_PARAM_GSMODIFIED,  
     IDA_PARAM_MAXNCF,  
     IDA_PARAM_PREC,  
     IDA_PARAMS_SIZE  
 };  
   
 /**  
  Here the full set of parameters is defined, along with upper/lower bounds,  
  etc. The values are stuck into the integ->params structure.  
   
  To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters  
  list. Then add a slv_param_*(...) statement below to define the type, description and range  
  for the new parameter.  
   
  @return 0 on success  
  */  
 static int integrator_ida_params_default(IntegratorSystem *integ) {  
     asc_assert(integ!=NULL);  
     asc_assert(integ->engine==INTEG_IDA);  
     slv_parameters_t *p;  
     p = &(integ->params);  
   
     slv_destroy_parms(p);  
   
     if (p->parms == NULL) {  
         //CONSOLE_DEBUG("params NULL");  
         p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);  
         if (p->parms == NULL)  
             return -1;  
         p->dynamic_parms = 1;  
     } else {  
         //CONSOLE_DEBUG("params not NULL");  
     }  
   
     /* reset the number of parameters to zero so that we can check it at the end */  
     p->num_parms = 0;  
   
     slv_param_bool(  
             p,  
             IDA_PARAM_AUTODIFF,  
             (SlvParameterInitBool) { {"autodiff"  
                             ,"Use auto-diff?",1  
                             ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"  
                         }, TRUE}  
                     );  
   
                     slv_param_char(p,IDA_PARAM_CALCIC  
                             ,(SlvParameterInitChar) { {"calcic"  
                                     ,"Initial conditions calcuation",1  
                                     ,"Use specified values of ydot to solve for inital y (Y),"  
                                     " or use the the values of the differential variables (yd) to solve"  
                                     " for the pure algebraic variables (ya) along with the derivatives"  
                                     " of the differential variables (yddot) (YA_YDP), or else don't solve"  
                                     " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"  
                                 }, "YA_YDP"}, (char *[]) {"Y", "YA_YDP", "NONE",NULL}  
                     );  
   
                     slv_param_bool(p,IDA_PARAM_SAFEEVAL  
                             ,(SlvParameterInitBool) { {"safeeval"  
                                     ,"Use safe evaluation?",1  
                                     ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "  
                                     "throw SIGFPE errors which will then halt integration."  
                                 }, FALSE}  
                     );  
   
                     slv_param_bool(p,IDA_PARAM_ATOLVECT  
                             ,(SlvParameterInitBool) { {"atolvect"  
                                     ,"Use 'ode_atol' values as specified?",1  
                                     ,"If TRUE, values of 'ode_atol' are taken from your model and used "  
                                     " in the integration. If FALSE, a scalar absolute tolerance value"  
                                     " is shared by all variables. See IDA manual, section 5.5.1"  
                                 }, TRUE}  
                     );  
   
                     slv_param_real(p,IDA_PARAM_ATOL  
                             ,(SlvParameterInitReal) { {"atol"  
                                     ,"Scalar absolute error tolerance",1  
                                     ,"Value of the scalar absolute error tolerance. See also 'atolvect'."  
                                     " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"  
                                 }, 1e-5, 0, 1e10}  
                     );  
   
                     slv_param_real(p,IDA_PARAM_RTOL  
                             ,(SlvParameterInitReal) { {"rtol"  
                                     ,"Scalar relative error tolerance",1  
                                     ,"Value of the scalar relative error tolerance. (Note that for IDA,"  
                                     " it's not possible to set per-variable relative tolerances as it is"  
                                     " with LSODE)."  
                                     " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"  
                                 }, 1e-4, 0, 1}  
                     );  
   
                     slv_param_char(p,IDA_PARAM_LINSOLVER  
                             ,(SlvParameterInitChar) { {"linsolver"  
                                     ,"Linear solver",1  
                                     ,"See IDA manual, section 5.5.3. Choose 'ASCEND' to use the linsolqr"  
                                     " direct linear solver bundled with ASCEND, 'DENSE' to use the dense"  
                                     " solver bundled with IDA, or one of the Krylov solvers SPGMR, SPBCG"  
                                     " or SPTFQMR (which still need preconditioners to be implemented"  
                                     " before they can be very useful."  
                                 }, "DENSE"}, (char *[]) {"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}  
                     );  
   
                     slv_param_int(p,IDA_PARAM_MAXL  
                             ,(SlvParameterInitInt) { {"maxl"  
                                     ,"Maximum Krylov dimension",0  
                                     ,"The maximum dimension of Krylov space used by the linear solver"  
                                     " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."  
                                     " The default of 0 results in IDA using its internal default, which"  
                                     " is currently a value of 5."  
                                 }, 0, 0, 20}  
                     );  
   
                     slv_param_int(p,IDA_PARAM_MAXORD  
                             ,(SlvParameterInitInt) { {"maxord"  
                                     ,"Maximum order of linear multistep method",0  
                                     ,"The maximum order of the linear multistep method with IDA. See"  
                                     " IDA manual p 38."  
                                 }, 5, 1, 5}  
                     );  
   
                     slv_param_bool(p,IDA_PARAM_GSMODIFIED  
                             ,(SlvParameterInitBool) { {"gsmodified"  
                                     ,"Gram-Schmidt Orthogonalisation Scheme", 2  
                                     ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"  
                                     " 5.5.6.6. Only applies when linsolve=SPGMR is selected."  
                                 }, TRUE}  
                     );  
   
                     slv_param_int(p,IDA_PARAM_MAXNCF  
                             ,(SlvParameterInitInt) { {"maxncf"  
                                     ,"Max nonlinear solver convergence failures per step", 2  
                                     ,"Maximum number of allowable nonlinear solver convergence failures"  
                                     " on one step. See IDA manual section 5.5.6.1."  
                                 }, 10,0,1000}  
                     );  
   
                     slv_param_char(p,IDA_PARAM_PREC  
                             ,(SlvParameterInitChar) { {"prec"  
                                     ,"Preconditioner",1  
                                     ,"See IDA manual, section section 5.6.8."  
                                 },"NONE"}, (char *[]) {"NONE","DIAG",NULL}  
                     );  
   
                     asc_assert(p->num_parms == IDA_PARAMS_SIZE);  
   
                     //CONSOLE_DEBUG("Created %d params", p->num_parms);  
   
                     return 0;  
                 }  
   
   
   
 /*******************************************  
  * SOLVE SETUP FUNCTIONS  
  *******************************************/  
   
 int ida_load_rellist(IntegratorSystem *integ) {  
     IntegratorIdaData *enginedata;  
     struct rel_relation **rels;  
     int i, j, n_solverrels, n_active_rels;  
   
     enginedata = integrator_ida_enginedata(integ);  
   
     n_solverrels = slv_get_num_solvers_rels(integ->system);  
     n_active_rels = slv_count_solvers_rels(integ->system, &integrator_ida_rel);  
     rels = slv_get_solvers_rel_list(integ->system);  
   
     if (enginedata->rellist != NULL) {  
         ASC_FREE(enginedata->rellist);  
         enginedata->rellist = NULL;  
     }  
   
     enginedata->rellist  
             = ASC_NEW_ARRAY(struct rel_relation *, n_active_rels);  
   
 #ifdef SOLVE_DEBUG  
     CONSOLE_DEBUG("rels matchbits:  0x%x",integrator_ida_rel.matchbits);  
     CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);  
   
     CONSOLE_DEBUG("Number of relations: %d",n_solverrels);  
     CONSOLE_DEBUG("Number of active relations: %d",n_active_rels);  
     CONSOLE_DEBUG("Number of dependent vars: %d",integ->n_y);  
     CONSOLE_DEBUG("Number of boundaries: %d",enginedata->nbnds);  
 #endif  
   
   
     j = 0;  
     for (i = 0; i < n_solverrels; ++i) {  
         if (rel_apply_filter(rels[i], &integrator_ida_rel)) {  
 #ifdef SOLVE_DEBUG  
             {  
                 char *relname = rel_make_name(integ->system, rels[i]);  
                 CONSOLE_DEBUG("rel '%s': 0x%x", relname, rel_flags(rels[i]));  
                 ASC_FREE(relname);  
             }  
 #endif  
             enginedata->rellist[j++] = rels[i];  
         }  
     }  
   
     asc_assert(j == n_active_rels);  
     enginedata->nrels = n_active_rels;  
   
     if (enginedata->nrels != integ->n_y) {  
         ERROR_REPORTER_HERE(ASC_USER_ERROR  
                 ,"Integration problem is not square (%d active rels, %d vars)"  
                 ,n_active_rels, integ->n_y  
         );  
         return 1; /* failure */  
     }  
   
     return 0;  
 }  
   
 /*  
  * Retrieve initial values from the system  
  */  
 int ida_retrieve_IVs(IntegratorSystem *integ, realtype t0, N_Vector y0,  
         N_Vector yp0) {  
   
 #ifdef SOLVE_DEBUG  
     char *varname;  
     char diffname[30];  
     int i;  
     CONSOLE_DEBUG("RETRIEVING INITIAL VALUES:");  
     CONSOLE_DEBUG("t0 = %f",t0);  
 #endif  
   
     integrator_get_y(integ, NV_DATA_S(y0));  
     integrator_get_ydot(integ, NV_DATA_S(yp0));  
   
 #ifdef SOLVE_DEBUG  
     fprintf(stderr, "index\t%25s\t%25s\n", "y", "ydot");  
     for (i = 0; i < integ->n_y; ++i) {  
         varname = var_make_name(integ->system, integ->y[i]);  
         fprintf(stderr, "%d\t%15s=%10f\t", i, varname, NV_Ith_S(y0,i));  
   
         if (integ->ydot[i]) {  
             ASC_FREE(varname);  
             varname = var_make_name(integ->system, integ->ydot[i]);  
             fprintf(stderr, "%15s=%10f\t\n", varname, NV_Ith_S(yp0,i));  
         } else {  
             snprintf(diffname,99,"diff(%s)",varname);  
             fprintf(stderr,"%15s=%10f\t\n",diffname,NV_Ith_S(yp0,i));  
         }  
         ASC_FREE(varname);  
     }  
 #endif  
   
     return 0;  
 }  
   
 /*  
  * Assign internal memory for the appropriate sundials version and error  
  * tolerance parameter settings.  
  */  
 int ida_malloc(IntegratorSystem *integ, void *ida_mem, realtype t0,  
         N_Vector y0, N_Vector yp0) {  
     int flag;  
     N_Vector abstolvect;  
     realtype reltol, abstol;  
   
     /* relative error tolerance */  
     reltol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_RTOL);  
     //CONSOLE_DEBUG("rtol = %8.2e",reltol);  
   
   
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
     flag = IDAInit(ida_mem, &integrator_ida_fex, t0, y0 ,yp0);  
 #else  
     if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)) {  
         /* vector of absolute tolerances */  
         CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");  
         abstolvect = N_VNew_Serial(integ->n_y);  
         integrator_get_atol(integ, NV_DATA_S(abstolvect));  
   
         flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV,  
                 reltol, abstolvect);  
   
         N_VDestroy_Serial(abstolvect);  
     }else{  
         /* scalar absolute tolerance (one value for all) */  
         abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);  
         CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);  
         flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS,  
                 reltol, &abstol);  
     }  
 #endif  
   
     if(flag == IDA_MEM_NULL) {  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");  
         return 2;  
     }else if(flag == IDA_MEM_FAIL) {  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");  
         return 3;  
     }else if(flag == IDA_ILL_INPUT) {  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");  
         return 4;  
     }  
   
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
     //CONSOLE_DEBUG("Assigning tolerances...");  
     /* assign tolerances */  
     if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)) {  
         //CONSOLE_DEBUG("using vector of atol values");  
         abstolvect = N_VNew_Serial(integ->n_y);  
         integrator_get_atol(integ,NV_DATA_S(abstolvect));  
         IDASVtolerances(ida_mem, reltol, abstolvect);  
         N_VDestroy_Serial(abstolvect);  
     }else{  
         /* scalar tolerances */  
         abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);  
         CONSOLE_DEBUG("using scalar atol value = %8.2e",abstol);  
         IDASStolerances(ida_mem, reltol, abstol);  
     }  
 #endif  
   
     /* success */  
     return 0;  
 }  
   
 /*  
  * Set parameter inputs for step size, the linear solver module and preconditioner  
  */  
 int ida_set_optional_inputs(IntegratorSystem *integ, void *ida_mem) {  
     int flag;  
     char *linsolver;  
     char *pname = NULL;  
     int maxl;  
     const IntegratorIdaPrec *prec = NULL;  
   
     IntegratorIdaData *enginedata = integ->enginedata;  
   
     IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *) integ);  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
     IDASetUserData(ida_mem, (void *)integ);  
 #else  
     IDASetRdata(ida_mem, (void *) integ);  
 #endif  
     IDASetMaxStep(ida_mem, integrator_get_maxstep(integ));  
     IDASetInitStep(ida_mem, integrator_get_stepzero(integ));  
     IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(integ));  
     if (integrator_get_minstep(integ) > 0) {  
         ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");  
     }  
   
     //CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));  
     IDASetMaxConvFails(ida_mem, SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));  
   
     //CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));  
     IDASetMaxOrd(ida_mem, SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));  
   
     /* there's no capability for setting *minimum* step size in IDA */  
   
     /* attach linear solver module, using the default value of maxl */  
     linsolver = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_LINSOLVER);  
     //CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);  
     if(strcmp(linsolver, "ASCEND") == 0) {  
         CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",integ->n_y);  
         IDAASCEND(ida_mem, integ->n_y);  
         IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *) integ);  
   
         enginedata->flagfntype = "IDAASCEND";  
         enginedata->flagfn = &IDAASCENDGetLastFlag;  
         enginedata->flagnamefn = &IDAASCENDGetReturnFlagName;  
   
     }else if (strcmp(linsolver, "DENSE") == 0) {  
         //CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",integ->n_y);  
         flag = IDADense(ida_mem, integ->n_y);  
         switch(flag){  
         case IDADENSE_SUCCESS:  
             break;  
         case IDADENSE_MEM_NULL:  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");  
             return 5;  
         case IDADENSE_ILL_INPUT:  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module");  
             return 5;  
         case IDADENSE_MEM_FAIL:  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE");  
             return 5;  
         default:  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return");  
             return 5;  
         }  
   
         if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)) {  
             //CONSOLE_DEBUG("USING AUTODIFF");  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
             flag = IDADlsSetDenseJacFn(ida_mem, &integrator_ida_djex);  
 #else  
             flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex,  
                     (void *) integ);  
 #endif  
             switch (flag) {  
             case IDADENSE_SUCCESS:  
                 break;  
             default:  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn");  
                 return 6;  
             }  
         }else{  
             CONSOLE_DEBUG("USING NUMERICAL DIFF");  
         }  
   
         enginedata->flagfntype = "IDADENSE";  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
         enginedata->flagfn = &IDADlsGetLastFlag;  
         enginedata->flagnamefn = &IDADlsGetReturnFlagName;  
 #else  
         enginedata->flagfn = &IDADenseGetLastFlag;  
         enginedata->flagnamefn = &IDADenseGetReturnFlagName;  
 #endif  
     }else{  
         /* remaining methods are all SPILS */  
         CONSOLE_DEBUG("IDA SPILS");  
   
         maxl = SLV_PARAM_INT(&(integ->params),IDA_PARAM_MAXL);  
         CONSOLE_DEBUG("maxl = %d",maxl);  
   
         /* what preconditioner? */  
         pname = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_PREC);  
         if(strcmp(pname, "NONE") == 0){  
             prec = NULL;  
         }else if(strcmp(pname, "JACOBI") == 0){  
             prec = &prec_jacobi;  
         }else{  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);  
             return 7;  
         }  
   
         /* which SPILS linear solver? */  
         if(strcmp(linsolver, "SPGMR") == 0){  
             CONSOLE_DEBUG("IDA SPGMR");  
             flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */  
         }else if(strcmp(linsolver, "SPBCG") == 0){  
             CONSOLE_DEBUG("IDA SPBCG");  
             flag = IDASpbcg(ida_mem, maxl);  
         }else if(strcmp(linsolver, "SPTFQMR") == 0){  
             CONSOLE_DEBUG("IDA SPTFQMR");  
             flag = IDASptfqmr(ida_mem, maxl);  
         }else{  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);  
             return 8;  
         }  
   
         if(prec){  
             /* assign the preconditioner to the linear solver */  
             (prec->pcreate)(integ);  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
             IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve);  
 #else  
             IDASpilsSetPreconditioner(ida_mem, prec->psetup, prec->psolve,  
                     (void *) integ);  
 #endif  
             CONSOLE_DEBUG("PRECONDITIONER = %s",pname);  
         }else{  
             CONSOLE_DEBUG("No preconditioner");  
         }  
   
         enginedata->flagfntype = "IDASPILS";  
         enginedata->flagfn = &IDASpilsGetLastFlag;  
         enginedata->flagnamefn = &IDASpilsGetReturnFlagName;  
   
         if(flag == IDASPILS_MEM_NULL){  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");  
             return 9;  
         }else if (flag == IDASPILS_MEM_FAIL) {  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");  
             return 9;  
         }/* else success */  
   
   
         /* assign the J*v function */  
         if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)) {  
             CONSOLE_DEBUG("USING AUTODIFF");  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
             flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex);  
 #else  
             flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex,  
                     (void *) integ);  
 #endif  
             if(flag == IDASPILS_MEM_NULL) {  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");  
                 return 10;  
             }else if (flag == IDASPILS_LMEM_NULL) {  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");  
                 return 10;  
             }/* else success */  
         }else{  
             CONSOLE_DEBUG("USING NUMERICAL DIFF");  
         }  
   
         if(strcmp(linsolver, "SPGMR") == 0) {  
             /* select Gram-Schmidt orthogonalisation */  
             if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_GSMODIFIED)) {  
                 CONSOLE_DEBUG("USING MODIFIED GS");  
                 flag = IDASpilsSetGSType(ida_mem, MODIFIED_GS);  
                 if(flag != IDASPILS_SUCCESS) {  
                     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");  
                     return 11;  
                 }  
             }else{  
                 CONSOLE_DEBUG("USING CLASSICAL GS");  
                 flag = IDASpilsSetGSType(ida_mem, CLASSICAL_GS);  
                 if(flag != IDASPILS_SUCCESS) {  
                     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");  
                     return 11;  
                 }  
             }  
         }  
     }  
   
     /* set linear solver optional inputs...  
      ...nothing here at the moment...  
      */  
   
     return 0;  
 } /* ida_set_optional_inputs */  
   
 /**  
  * Calculate initial conditions using IDACalcIC  
  */  
   
 int ida_setup_IC(IntegratorSystem *integ, void *ida_mem,  
         realtype tout1, realtype t0, N_Vector y0, N_Vector yp0) {  
     int i, flag, flag1;  
     int icopt; /* initial conditions strategy */  
     N_Vector id;  
   
     IntegratorIdaData *enginedata = integ->enginedata;  
   
 #ifdef SOLVE_DEBUG  
     char *varname;  
 #endif  
   
 #ifdef IDA_BND_DEBUG  
     CONSOLE_DEBUG("Solving initial conditions...");  
 #endif  
   
     icopt = 0;  
     if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "Y") == 0) {  
 #ifdef SOLVE_DEBUG  
         CONSOLE_DEBUG("Solving initial conditions using values of yddot");  
 #endif  
         icopt = IDA_Y_INIT;  
         asc_assert(icopt!=0);  
     }else if(0==strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "YA_YDP")){  
 #ifdef SOLVE_DEBUG  
         CONSOLE_DEBUG("Solving initial conditions using values of yd");  
 #endif  
         icopt = IDA_YA_YDP_INIT;  
         asc_assert(icopt!=0);  
         id = N_VNew_Serial(integ->n_y);  
         for(i = 0; i < integ->n_y; ++i) {  
             if(integ->ydot[i] == NULL) {  
                 NV_Ith_S(id, i) = 0.0;  
 #ifdef SOLVE_DEBUG  
                 varname = var_make_name(integ->system, integ->y[i]);  
                 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);  
                 ASC_FREE(varname);  
 #endif  
             }else{  
 #ifdef SOLVE_DEBUG  
                 CONSOLE_DEBUG("y[%d] is differential",i);  
 #endif  
                 NV_Ith_S(id, i) = 1.0;  
             }  
         }  
         IDASetId(ida_mem, id);  
         N_VDestroy_Serial(id);  
     }else if (strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "NONE")  
             == 0) {  
         ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");  
     } else {  
         ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");  
     }  
   
     if (icopt) {  
 #ifdef SOLVE_DEBUG  
         CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);  
 #endif  
   
 #ifdef ASC_SIGNAL_TRAPS  
         /* catch SIGFPE if desired to */  
         if(enginedata->safeeval) {  
             CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");  
             Asc_SignalHandlerPush(SIGFPE, SIG_DFL);  
         }else {  
 # ifdef FEX_DEBUG  
             CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");  
 # endif  
             Asc_SignalHandlerPushDefault(SIGFPE);  
         }  
         if(setjmp(g_fpe_env) == 0) {  
 #endif  
   
 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=3  
             flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3  */  
 # else  
             flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);  
 # endif  
             /* check flags and output status */  
             switch(flag){  
             case IDA_SUCCESS:  
 #ifdef SOLVE_DEBUG  
                 CONSOLE_DEBUG("Initial conditions solved OK");  
 #endif  
                 break;  
   
             case IDA_LSETUP_FAIL:  
             case IDA_LINIT_FAIL:  
             case IDA_LSOLVE_FAIL:  
             case IDA_NO_RECOVERY:  
                 flag1 = -999;  
                 flag = (enginedata->flagfn)(ida_mem, &flag1);  
                 if (flag) {  
                     ERROR_REPORTER_HERE(ASC_PROG_ERR  
                             ,"Unable to retrieve error code from %s (err %d)"  
                             ,enginedata->flagfntype,flag  
                     );  
                     return 12;  
                 }  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR  
                         ,"%s returned flag '%s' (value = %d)"  
                         ,enginedata->flagfntype,(enginedata->flagnamefn)(flag1),flag1  
                 );  
                 return 12;  
   
             default:  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");  
                 return 12;  
             }  
 #ifdef ASC_SIGNAL_TRAPS  
         }else{  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");  
             return 13;  
         }  
   
         if(enginedata->safeeval){  
             Asc_SignalHandlerPop(SIGFPE, SIG_DFL);  
         }else{  
             /* CONSOLE_DEBUG("pop..."); */  
             Asc_SignalHandlerPopDefault(SIGFPE);  
             /* CONSOLE_DEBUG("...pop"); */  
         }  
 #endif  
     }/* icopt */  
   
     return 0;  
 } /* ida_setup_IC */  
   
   
 int ida_root_init(IntegratorSystem *integ, void *ida_mem) {  
     IntegratorIdaData *enginedata = integ->enginedata;  
   
     if(enginedata->nbnds) {  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
         IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn);  
 #else  
         IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn,  
                 (void *) integ);  
 #endif  
     }  
   
     return 0;  
   
 }  
   
 /**  
  * Allocate memory and prepare for integration.  
  *  
  * @param t the first point at which a solution is desired. Required by IdaCalcIC  
  *  
  * @return 0 on success  
  */  
 int ida_prepare_integrator(IntegratorSystem *integ, void *ida_mem,  
         realtype tout1) {  
     realtype t0;  
     N_Vector y0, yp0;  
   
     y0  = ida_bnd_new_zero_NV(integ->n_y);  
     yp0 = ida_bnd_new_zero_NV(integ->n_y);  
   
 #if 0  
     int i;  
     double val;  
     CONSOLE_DEBUG("Values of the derivatives present in the model");  
     for(i=0; i < integ->n_y; i++) {  
         if(integ->ydot[i]){  
             val = var_value(integ->ydot[i]);  
             CONSOLE_DEBUG("ydot[%d]= %g", i, val);  
         }  
     }  
 #endif  
   
     t0 = integrator_get_t(integ);  
     ida_retrieve_IVs(integ, t0, y0, yp0);  
   
 #ifdef IDA_BND_DEBUG  
     CONSOLE_DEBUG("Initial values BEFORE IDACalcIC: y0 =");  
     N_VPrint_Serial(y0);  
     CONSOLE_DEBUG("yp0 = ");  
     N_VPrint_Serial(yp0);  
 #endif  
   
     /* allocate internal memory  */  
     ida_malloc(integ, ida_mem, t0, y0, yp0);  
   
     /* set optional inputs... */  
     ida_set_optional_inputs(integ, ida_mem);  
   
     /* calculate initial conditions */  
     ida_setup_IC(integ, ida_mem, tout1, t0, y0, yp0);  
   
     /* specify ROOT-FINDING problem (if necessary) */  
     ida_root_init(integ, ida_mem);  
   
     /* Clean up */  
     N_VDestroy_Serial(y0);  
     N_VDestroy_Serial(yp0);  
     return 0;  
 }  
   
 /**  
  * Reinitialise IDA after boundary-crossing. We assume that IDAInit has been called.  
  * Works with Sundials 2.4.0 and later versions.  
  *  
  * @param t the first point at which a solution is desired. Required by IdaCalcIC  
  *  
  * @return 0 on success  
  */  
 int ida_reinit_integrator(IntegratorSystem *integ, void *ida_mem,  
         realtype tout1) {  
     realtype t0;  
     N_Vector y0, yp0;  
   
     y0  = ida_bnd_new_zero_NV(integ->n_y);  
     yp0 = ida_bnd_new_zero_NV(integ->n_y);  
   
     t0 = integrator_get_t(integ);  
     ida_retrieve_IVs(integ, t0, y0, yp0);  
   
     int flag;  
   
     flag = IDAReInit(ida_mem, t0, y0, yp0);  
     if(flag!=IDA_SUCCESS){  
         ERROR_REPORTER_HERE(ASC_PROG_ERR, "Reinitialisation failed.");  
     }  
   
     /* calculate initial conditions */  
     ida_setup_IC(integ, ida_mem, tout1, t0, y0, yp0);  
   
     /* Clean up */  
     N_VDestroy_Serial(y0);  
     N_VDestroy_Serial(yp0);  
     return 0;  
 }  
   
     /*-------------------------------------------------------------  
      MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.  
      */  
   
 /**  
  Main IDA solver routine. We can assume that the integrator_ida_analyse  
  function will already have been run.  
   
  The presence of 'start_index' and 'finish_index' is a hangover from the  
  Tcl/Tk GUI. We would like to get rid of them eventually.  
   
  @return 0 on success */  
 static int integrator_ida_solve(IntegratorSystem *integ,  
         unsigned long start_index, unsigned long finish_index) {  
     void *ida_mem;  
     int t_index;  
     realtype t0, tout, tret, tol;  
     realtype tmin;              /** < The length of a small additional step made after an event is triggered */  
     N_Vector ypret, yret;  
     IntegratorIdaData *enginedata;  
     int i, flag;  
   
     int *rootsfound;            /** < IDA rootfinder reports root index in here */  
     int *rootdir = NULL;                /** < Used to tell IDA to ignore doulve crossings */  
     int *bnd_cond_states;       /** < Record of boundary states so that IDA can tell LRSlv  
                                        how to evaluate a boundary crossing */  
     int *bnd_not_set;  
     int all_bnds_set = 1;  
   
     int need_to_reconfigure;    /** < Flag to indicate system rebuild after crossing */  
     int need_to_reinteg = 0;    /** < Flag for when crossings happen on or very close to timesteps */  
     int skipping_output;        /** < Flag to skip output to reporter */  
     int qrslv_ind, lrslv_ind;  
   
     int after_root = 0;  
   
 #if defined(SOLVE_DEBUG) || defined(IDA_BND_DEBUG)  
     char *relname;  
 #endif  
   
     //CONSOLE_DEBUG("STARTING IDA...");  
     /* Setup boundary list */  
     enginedata = integrator_ida_enginedata(integ);  
     enginedata->bndlist = slv_get_solvers_bnd_list(integ->system);  
     enginedata->nbnds = slv_get_num_solvers_bnds(integ->system);  
     enginedata->safeeval = SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_SAFEEVAL);  
     //CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);  
   
     qrslv_ind = slv_lookup_client("QRSlv");  
     lrslv_ind = slv_lookup_client("LRSlv");  
   
     bnd_not_set = ASC_NEW_ARRAY(int,enginedata->nbnds);  
   
 #ifdef SOLVE_DEBUG  
     integrator_ida_debug(integ, stderr);  
 #endif  
   
     /* store reference to list of relations (in enginedata) */  
     ida_load_rellist(integ);  
   
     /* create IDA object */  
     ida_mem = IDACreate();  
   
     /* Setup parameter inputs and initial conditions for IDA. */  
     tout = samplelist_get(integ->samples, start_index + 1);  
     /* solve the initial conditions, allocate memory, other stuff... */  
     ida_prepare_integrator(integ, ida_mem, tout);  
   
   
   
   
     /* Initialise boundary condition states if appropriate. Reconfigure if necessary */  
     if(enginedata->nbnds){  
         CONSOLE_DEBUG("Initialising boundary states");  
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR<4  
         ERROR_REPORTER_HERE(ASC_PROG_WARNING, "Warning: boundary detection is"  
                 "unreliable with SUNDIALS pre version 2.4.0. Please update if you"  
                 "wish to use IDA for conditional integration");  
 #endif  
         bnd_cond_states = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);  
   
         /* identify if we're exactly *on* any boundaries currently */  
         for(i = 0; i < enginedata->nbnds; i++) {  
 #ifdef IDA_BND_DEBUG  
             relname = bnd_make_name(integ->system,enginedata->bndlist[i]);  
 #endif  
             bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);  
             bnd_set_ida_first_cross(enginedata->bndlist[i],1);  
             if(bndman_real_eval(enginedata->bndlist[i]) == 0) {  
                 /* if the residual for the boundary is zero (ie looks like we are *on* the boundary?) JP */  
 #ifdef IDA_BND_DEBUG  
                 CONSOLE_DEBUG("Boundary '%s': not set",relname);  
 #endif  
                 bnd_not_set[i] = 1;  
                 all_bnds_set = 0;  
             }else{  
                 bnd_not_set[i] = 0;  
             }  
 #ifdef IDA_BND_DEBUG  
             CONSOLE_DEBUG("Boundary '%s' is %d",relname,bnd_cond_states[i]);  
             ASC_FREE(relname);  
 #endif  
         }  
         CONSOLE_DEBUG("Setting up LRSlv...");  
         if(ida_setup_lrslv(integ,qrslv_ind,lrslv_ind)){  
             ERROR_REPORTER_HERE(ASC_USER_ERROR, "Idaanalyse failed.");  
             return 1;  
         }  
   
     }  
   
   
   
   
     tol = 0.0001*(samplelist_get(integ->samples, finish_index)  
                     - samplelist_get(integ->samples, start_index))  
                 /samplelist_length(integ->samples);  
     tmin = tol;  
   
     /* -- set up the IntegratorReporter */  
     integrator_output_init(integ);  
   
     /* -- store the initial values of all the stuff */  
     integrator_output_write(integ);  
     integrator_output_write_obs(integ);  
   
     /* specify where the returned values should be stored */  
     yret    = ida_bnd_new_zero_NV(integ->n_y);  
     ypret   = ida_bnd_new_zero_NV(integ->n_y);  
   
     /* advance solution in time, return values as yret and derivatives as ypret */  
     integ->currentstep = 1;  
   
     rootdir = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);  
   
     for(t_index = start_index + 1; t_index <= finish_index; ++t_index, ++integ->currentstep) {  
         tout = samplelist_get(integ->samples, t_index);  
         t0 = integrator_get_t(integ);  
   
         asc_assert(tout > t0);  
   
 #ifdef SOLVE_DEBUG  
         CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, tout);  
 #endif  
   
         if(!after_root){  
             /* reset the root-crossing-direction flags */  
             for(i = 0; i < enginedata->nbnds; i++){  
                 rootdir[i] = 0;  
 #ifdef IDA_BND_DEBUG  
                 char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);  
                 CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d]=%d, bndman_calc_satisfied=%d; (trigger dirn=%s)"  
                     ,n, i, bnd_cond_states[i], bndman_calc_satisfied(enginedata->bndlist[i])  
                     ,rootdir[i]==1?"UP":(rootdir[i]==-1?"DOWN":"both")  
                 );  
                 ASC_FREE(n);  
 #endif  
             }  
               
             if(enginedata->nbnds) IDASetRootDirection(ida_mem, rootdir);  
         }  
   
         /**  
          * do { solve routine } while(need_to_reintegrate)  
          * If IDA detects a boundary sufficiently far from the desired output  
          * time tout, then the system is reconfigured and IDASolve is recalled  
          * up to the same tindex.  
          *  
          * @TODO "Sufficiently far" is currently tol = 0.0001, i.e complete arbitrary  
          * Should reflect IDACalcIC's requirements?  
          */  
         do{  
             if(need_to_reinteg && !after_root){  
 #ifdef IDA_BND_DEBUG  
                 CONSOLE_DEBUG("Resuming integration from %f to %f", integrator_get_t(integ), tout);  
 #endif  
                 integrator_output_write(integ);  
             }  
   
 #ifdef IDA_BND_DEBUG  
             if(after_root && (tout - tret) > tol){  
                 CONSOLE_DEBUG("Resuming integration from %f to %f", integrator_get_t(integ), tret + tmin);  
             }  
 #endif  
   
             /* Control flags for boundary crossings */  
             need_to_reinteg = 0;  
             skipping_output = 0;  
   
 #ifdef ASC_SIGNAL_TRAPS  
             Asc_SignalHandlerPushDefault(SIGINT);  
             if(setjmp(g_int_env) == 0) {  
 #endif  
                 if(after_root && (tout - tret) > tmin){  
                     CONSOLE_DEBUG("Continuing (slightly after root)...");  
                     /* after a root, we advance time by a tiny amount to ensure  
                     we traverse past the switching poing -- FIXME a bit hacky? */  
                     flag = IDASolve(ida_mem, tret + tmin, &tret, yret, ypret, IDA_NORMAL);  
                 }else{  
                     CONSOLE_DEBUG("Continuing immediately...");  
                     /* if we didn't pass a root, continue immediately from prev end time */  
                     flag = IDASolve(ida_mem, tout, &tret, yret, ypret, IDA_NORMAL);  
                 }  
                 if(after_root){  
                     CONSOLE_DEBUG("Running logic solver after root...");  
                     ida_log_solve(integ, lrslv_ind);  
                     after_root = 0;  
                 }  
 #ifdef ASC_SIGNAL_TRAPS  
             }else{  
                 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Caught interrupt");  
                 flag = -555;  
             }  
             Asc_SignalHandlerPopDefault(SIGINT);  
 #endif  
   
             if(enginedata->nbnds) {  
                 if(flag == IDA_ROOT_RETURN){  
 #ifdef IDA_BND_DEBUG  
                     CONSOLE_DEBUG("IDA reports root found!");  
 #endif  
                     /* Store the root index */  
                     rootsfound = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);  
   
                     if(IDA_SUCCESS != IDAGetRootInfo(ida_mem, rootsfound)) {  
                         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch boundary-crossing info");  
                         return 14;  
                     }  
   
 #ifdef IDA_BND_DEBUG  
                     /* write out the boundaries that were crossed */  
                     for(i = 0; i < enginedata->nbnds; i++) {  
                         if(rootsfound[i]) {  
                             relname = bnd_make_name(integ->system, enginedata->bndlist[i]);  
                             CONSOLE_DEBUG("Boundary '%s' crossed at time x = %f, direction %s"  
                                 ,relname,tret,(rootsfound>0?"UP":"DOWN")  
                             );  
                             ASC_FREE(relname);  
                         }  
                     }  
 #endif  
   
                       
                     if(all_bnds_set == 0){  
 #ifdef IDA_BND_DEBUG  
                         CONSOLE_DEBUG("Unset bounds exist; evaluate them explicitly...");  
 #endif  
                         all_bnds_set = 1;  
                         for(i = 0; i < enginedata->nbnds; i++){  
                             if(bnd_not_set[i]){  
                                 if(!rootsfound[i]){  
                                     bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);  
 #ifdef IDA_BND_DEBUG  
                                     relname = bnd_make_name(integ->system, enginedata->bndlist[i]);  
                                     CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d] = %d"  
                                         ,relname,i,bnd_cond_states[i]  
                                     );  
                                     ASC_FREE(relname);  
 #endif  
                                 }else all_bnds_set = 0;  
                             }  
                         }  
                     }  
   
                     if(!after_root){  
 #ifdef IDA_BND_DEBUG  
                         CONSOLE_DEBUG("Just 'after_root'...");  
                         for(i=0;i<enginedata->nbnds;++i){  
                             relname = bnd_make_name(integ->system, enginedata->bndlist[i]);  
                             CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d] = %d"  
                                 ,relname,i,bnd_cond_states[i]  
                             );  
                             ASC_FREE(relname);  
                         }  
 #endif  
                         need_to_reconfigure = ida_cross_boundary(integ, rootsfound,  
                             bnd_cond_states, qrslv_ind, lrslv_ind);  
                     }  
                     if(need_to_reconfigure == 2) {  
                         ERROR_REPORTER_HERE(ASC_USER_ERROR,"Analysis after the boundary failed.");  
                         return 1;  
                     }  
                     if(need_to_reconfigure){  
                         after_root = 1;  
                         if (ida_bnd_update_relist(integ) != 0) {  
                             /* system not square, failure */  
                             return 1;  
                         }  
 #ifdef IDA_BND_DEBUG  
                         CONSOLE_DEBUG("Boundary(ies) crossed at x=%f: need to reinitialize solver",tret);  
 #endif  
                         /* so, now we need to restart the integration. we will assume that  
                          everything changes: number of variables, etc, etc, etc. */  
   
                         /* Need to destroy and rebuild system */  
                         /* IDAFree(ida_mem); */  
                         /* ida_mem = IDACreate(); */  
                         /* ida_reinit_integrator(integ, ida_mem, tout); */  
                         /* n_y may have changed */  
                         N_VDestroy_Serial(yret);  
                         N_VDestroy_Serial(ypret);  
   
                         yret = N_VNew_Serial(integ->n_y);  
                         ypret = N_VNew_Serial(integ->n_y);  
                     } /* need to reconfigure */  
   
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4  
                     /* set rootdir to -1*rootsfound to set IDA  
                      * to ignore double crossings */  
                     for(i = 0; i < enginedata->nbnds; i++) {  
                         rootdir[i] = -1*rootsfound[i];  
 #ifdef IDA_BND_DEBUG  
                         char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);  
                         CONSOLE_DEBUG("Set direction=%d for boundary '%s'",rootdir[i],n);  
                         ASC_FREE(n);  
 #endif  
                     }  
                     IDASetRootDirection(ida_mem, rootdir);\  
                     /*^^^ FIXME what about setting root direction for 'normal' roots? */  
 #endif  
                     ASC_FREE(rootsfound);  
                     ida_reinit_integrator(integ, ida_mem, tout);  
                 } /* IDA_ROOT_RETURN */  
   
             } /* nbnds */  
   
             /* Are we sufficiently far from tout to continue  
              * integrating on this timestep? */  
             if(fabs(tret - tout) > tol) need_to_reinteg = 1;  
             else if(after_root){  
                 /* Advance timestep, skip writing output once*/  
                 tout = samplelist_get(integ->samples, t_index + 1);  
                 skipping_output = 1;  
             }  
   
         }while (need_to_reinteg); /* end of solve time step */  
   
         if(flag != IDA_ROOT_RETURN && all_bnds_set == 0) {  
             all_bnds_set = 1;  
             for(i = 0; i < enginedata->nbnds; i++) {  
                 if(bnd_not_set[i]) {  
                     bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);  
 #ifdef IDA_BND_DEBUG  
                     char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);  
                     CONSOLE_DEBUG("Boundary '%s' not set; satisfied=%d",n,bnd_cond_states[i]);  
                     ASC_FREE(n);  
 #endif  
                 }  
             }  
         }  
   
         if(!skipping_output){  
             /* pass the values of everything back to the compiler */  
             integrator_set_t(integ, (double) tret);  
             integrator_set_y(integ, NV_DATA_S(yret));  
             integrator_set_ydot(integ, NV_DATA_S(ypret));  
   
             /* -- store the current values of all the stuff */  
             if(integrator_output_write(integ) == 0) {  
                 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Interrupted at t = %f", tout);  
                 break;  
             }  
             integrator_output_write_obs(integ);  
         }  
   
         if(flag < 0) {  
             ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", tout, flag);  
             break;  
         }  
     }/* loop through next sample timestep */  
   
     ASC_FREE(rootdir);  
   
     /* -- close the IntegratorReporter */  
     integrator_output_close(integ);  
   
     /* get optional outputs */  
 #ifdef STATS_DEBUG  
     IntegratorIdaStats stats;  
     if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)) {  
         integrator_ida_write_stats(&stats);  
     }else{  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");  
     }  
 #endif  
   
     /* free solution memory */  
     N_VDestroy_Serial(yret);  
     N_VDestroy_Serial(ypret);  
   
     /* free bnd states if appropriate */  
     if(enginedata->nbnds) {  
         ASC_FREE(bnd_cond_states);  
     }  
   
     /* free solver memory */  
     IDAFree(&ida_mem);  
   
     if(flag < -500) {  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", tout);  
         return -flag;  
     }  
   
     if(flag < 0) {  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", tout);  
         return 14;  
     }  
   
     /* all done, success */  
     return 0;  
 }  
   
 /*----------------------------------------------  
  STATS  
  */  
   
 /**  
  A simple wrapper to the IDAGetIntegratorStats function. Returns all the  
  status in a struct instead of separately.  
   
  @return IDA_SUCCESS on success.  
  */  
 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s) {  
   
 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2  
   
     int res;  
   
     /*  
      There is an error in the documentation for this function in Sundials 2.2.  
      According the the header file, the hinused stat is not provided.  
      */  
     res = IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals,  
             &s->nlinsetups, &s->netfails, &s->qlast, &s->qcur, &s->hlast,  
             &s->hcur, &s->tcur);  
   
     /* get the missing statistic */  
     IDAGetActualInitStep(ida_mem, &s->hinused);  
   
     return res;  
 #else  
   
     return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups  
             ,&s->netfails, &s->qlast, &s->qcur, &s->hinused  
             ,&s->hlast, &s->hcur, &s->tcur  
     );  
   
 #endif  
 }  
   
 /* vim: set ts=4: */  
Legend:
Removed from v.3031  
changed lines
  Added in v.3032
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22