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

Diff of /trunk/base/generic/solver/ida.c

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

revision 943 by johnpye, Sat Nov 25 05:26:47 2006 UTC revision 944 by johnpye, Sat Nov 25 10:46:13 2006 UTC
# Line 148  IntegratorIdaData *integrator_ida_engine Line 148  IntegratorIdaData *integrator_ida_engine
148    
149  enum ida_parameters{  enum ida_parameters{
150      IDA_PARAM_AUTODIFF      IDA_PARAM_AUTODIFF
151        ,IDA_PARAM_RTOL
152        ,IDA_PARAM_ATOL
153        ,IDA_PARAM_ATOLVECT
154      ,IDA_PARAMS_SIZE      ,IDA_PARAMS_SIZE
155  };  };
156    
# Line 184  int integrator_ida_params_default(Integr Line 187  int integrator_ida_params_default(Integr
187          }, TRUE}          }, TRUE}
188      );      );
189    
190        slv_param_bool(p,IDA_PARAM_ATOLVECT
191                ,(SlvParameterInitBool){{"atolvect"
192                ,"Use 'ode_atol' values as specified?",1
193                ,"If TRUE, values of 'ode_atol' are taken from your model and used "
194                " in the integration. If FALSE, a scalar absolute tolerance value"
195                " is shared by all variables. See IDA manual, section 5.5.1"
196            }, TRUE }
197        );
198    
199        slv_param_real(p,IDA_PARAM_ATOL
200                ,(SlvParameterInitReal){{"atol"
201                ,"Scalar absolute error tolerance",1
202                ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
203                " See IDA manual, section 5.5.1"
204            }, 1e-5, DBL_MIN, DBL_MAX }
205        );
206    
207        slv_param_real(p,IDA_PARAM_RTOL
208                ,(SlvParameterInitReal){{"rtol"
209                ,"Scalar relative error tolerance",1
210                ,"Value of the scalar relative error tolerance."
211                " See IDA manual, section 5.5.1"
212            }, 1e-5, DBL_MIN, DBL_MAX }
213        );
214    
215      asc_assert(p->num_parms == IDA_PARAMS_SIZE);      asc_assert(p->num_parms == IDA_PARAMS_SIZE);
216    
217      CONSOLE_DEBUG("Created %d params", p->num_parms);      CONSOLE_DEBUG("Created %d params", p->num_parms);
# Line 203  int integrator_ida_solve( Line 231  int integrator_ida_solve(
231  ){  ){
232      void *ida_mem;      void *ida_mem;
233      int size, flag, t_index;      int size, flag, t_index;
234      realtype t0, reltol, t, tret, tout1;      realtype t0, reltol, abstol, t, tret, tout1;
235      N_Vector y0, yp0, abstol, ypret, yret;      N_Vector y0, yp0, abstolvect, ypret, yret;
236      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
237    
238      CONSOLE_DEBUG("STARTING IDA...");      CONSOLE_DEBUG("STARTING IDA...");
# Line 225  int integrator_ida_solve( Line 253  int integrator_ida_solve(
253          return 0; /* failure */          return 0; /* failure */
254      }      }
255    
     CONSOLE_DEBUG("RETRIEVING t0");  
   
256      /* retrieve initial values from the system */      /* retrieve initial values from the system */
257      t0 = samplelist_get(blsys->samples,start_index);      
258        /** @TODO fix this, the starting time != first sample */
259        t0 = integrator_get_t(blsys);
260        CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
261    
262      CONSOLE_DEBUG("RETRIEVING y0");      CONSOLE_DEBUG("RETRIEVING y0");
263    
# Line 246  int integrator_ida_solve( Line 275  int integrator_ida_solve(
275      /* create IDA object */      /* create IDA object */
276      ida_mem = IDACreate();      ida_mem = IDACreate();
277    
278      /* retrieve the absolute tolerance values for each variable */      /* relative error tolerance */  
279      abstol = N_VNew_Serial(size);      reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
     N_VConst(0.1,abstol); /** @TODO fill in the abstol values from the model */  
     reltol = 0.001;  
280    
281      /* allocate internal memory */      /* allocate internal memory */
282      flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstol);      if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
283            /* vector of absolute tolerances */
284            CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
285            abstolvect = N_VNew_Serial(size);
286            integrator_get_atol(blsys,NV_DATA_S(abstolvect));
287    
288            flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
289        }else{
290            /* scalar absolute tolerance (one value for all) */
291            CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
292            abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
293            flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
294    
295        }
296    
297      if(flag==IDA_MEM_NULL){      if(flag==IDA_MEM_NULL){
298          ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
299          return 0;          return 0;
# Line 307  int integrator_ida_solve( Line 348  int integrator_ida_solve(
348    
349      /* correct initial values, given derivatives */      /* correct initial values, given derivatives */
350      blsys->currentstep=0;      blsys->currentstep=0;
351      t_index=start_index+1;      t_index=start_index;
352      tout1 = samplelist_get(blsys->samples, t_index);      tout1 = samplelist_get(blsys->samples, t_index);
353    
354        /* CONSOLE_DEBUG("Giving t value %f to IDACalcIC", tout1);*/
355    
356  #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3  #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
357        /* note the new API from version 2.3 and onwards */
358      flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);      flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
359  #else  #else
360      flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);      flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
# Line 495  int integrator_ida_jvex(realtype tt, N_V Line 539  int integrator_ida_jvex(realtype tt, N_V
539      var_filter_t filter;      var_filter_t filter;
540      int count;      int count;
541    
542      fprintf(stderr,"\n--------------\n");      /* fprintf(stderr,"\n--------------\n"); */
543      CONSOLE_DEBUG("EVALUTING JACOBIAN...");      /* CONSOLE_DEBUG("EVALUTING JACOBIAN..."); */
544    
545      blsys = (IntegratorSystem *)jac_data;      blsys = (IntegratorSystem *)jac_data;
546      enginedata = integrator_ida_enginedata(blsys);      enginedata = integrator_ida_enginedata(blsys);
# Line 517  int integrator_ida_jvex(realtype tt, N_V Line 561  int integrator_ida_jvex(realtype tt, N_V
561      filter.matchbits = VAR_SVAR;      filter.matchbits = VAR_SVAR;
562      filter.matchvalue = VAR_SVAR;      filter.matchvalue = VAR_SVAR;
563    
564      CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v));      /* CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v)); */
565      for(i=0; i<NV_LENGTH_S(v); ++i){      /* for(i=0; i<NV_LENGTH_S(v); ++i){
566          CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));          CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));
567      }      }*/
568    
569      Asc_SignalHandlerPush(SIGFPE,SIG_IGN);      Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
570      if (setjmp(g_fpe_env)==0) {      if (setjmp(g_fpe_env)==0) {
# Line 528  int integrator_ida_jvex(realtype tt, N_V Line 572  int integrator_ida_jvex(realtype tt, N_V
572                  i< enginedata->nrels && relptr != NULL;                  i< enginedata->nrels && relptr != NULL;
573                  ++i, ++relptr                  ++i, ++relptr
574          ){          ){
575              fprintf(stderr,"\n");              /* fprintf(stderr,"\n"); */
576              relname = rel_make_name(blsys->system, *relptr);              relname = rel_make_name(blsys->system, *relptr);
577              CONSOLE_DEBUG("RELATION %d '%s'",i,relname);              /* CONSOLE_DEBUG("RELATION %d '%s'",i,relname); */
578              ASC_FREE(relname);              ASC_FREE(relname);
579    
580              /* get derivatives for this particular relation */              /* get derivatives for this particular relation */
581              status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);              status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
582              CONSOLE_DEBUG("Got derivatives against %d matching variables", count);              /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
583    
584              for(j=0;j<count;++j){              for(j=0;j<count;++j){
585                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
586                  CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname);                  /* CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname); */
587                  ASC_FREE(varname);                  ASC_FREE(varname);
588              }              }
589    
590              if(!status){              if(!status){
591                  CONSOLE_DEBUG("Derivatives for relation %d OK",i);                  /* CONSOLE_DEBUG("Derivatives for relation %d OK",i); */
592              }else{              }else{
593                  CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);                  CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
594                  break;                  break;
# Line 561  int integrator_ida_jvex(realtype tt, N_V Line 605  int integrator_ida_jvex(realtype tt, N_V
605                  /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */                  /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
606                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);                  varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
607                  if(varname){                  if(varname){
608                      CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);                      /* CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]); */
609                      ASC_FREE(varname);                      ASC_FREE(varname);
610                  }else{                  }else{
611                      CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);                      CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
612                  }                  }
613                                    
614                  var_yindex = blsys->y_id[variables[j]];                  var_yindex = blsys->y_id[variables[j]];
615                  CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex);                  /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */
616    
617                  if(var_yindex >= 0){                  if(var_yindex >= 0){
618                      CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex));                      /* CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
619                      Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);                      Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
620                  }else{                  }else{
621                      var_yindex = -var_yindex-1;                      var_yindex = -var_yindex-1;
622                      CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex));                      /* CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
623                      Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;                      Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
624                  }                  }
625              }              }
626    
627              NV_Ith_S(Jv,i) = Jv_i;              NV_Ith_S(Jv,i) = Jv_i;
628              CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i);              /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */
629    
630              if(status){              if(status){
631                  /* presumably some error_reporter will already have been made*/                  /* presumably some error_reporter will already have been made*/

Legend:
Removed from v.943  
changed lines
  Added in v.944

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