/[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 990 by johnpye, Thu Dec 21 04:22:07 2006 UTC revision 991 by johnpye, Thu Dec 21 10:44:32 2006 UTC
# Line 179  enum ida_parameters{ Line 179  enum ida_parameters{
179      ,IDA_PARAM_ATOL      ,IDA_PARAM_ATOL
180      ,IDA_PARAM_ATOLVECT      ,IDA_PARAM_ATOLVECT
181      ,IDA_PARAM_GSMODIFIED      ,IDA_PARAM_GSMODIFIED
182        ,IDA_PARAM_MAXNCF
183          ,IDA_PARAMS_SIZE          ,IDA_PARAMS_SIZE
184  };  };
185    
# Line 277  int integrator_ida_params_default(Integr Line 278  int integrator_ida_params_default(Integr
278    
279      slv_param_bool(p,IDA_PARAM_GSMODIFIED      slv_param_bool(p,IDA_PARAM_GSMODIFIED
280              ,(SlvParameterInitBool){{"gsmodified"              ,(SlvParameterInitBool){{"gsmodified"
281              ,"Use modified Gram-Schmidt orthogonalisation in SPGMR?",2              ,"Gram-Schmidt Orthogonalisation Scheme", 2
282              ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"              ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
283              " 5.5.6.6. Only applies when linsolve=SPGMR is selected."              " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
284          }, TRUE}          }, TRUE}
285      );      );
286    
287        slv_param_int(p,IDA_PARAM_MAXNCF
288                ,(SlvParameterInitInt){{"maxncf"
289                ,"Max nonlinear solver convergence failures per step", 2
290                ,"Maximum number of allowable nonlinear solver convergence failures"
291                " on one step. See IDA manual section 5.5.6.1."
292            }, 10,0,1000 }
293        );
294    
295      asc_assert(p->num_parms == IDA_PARAMS_SIZE);      asc_assert(p->num_parms == IDA_PARAMS_SIZE);
296    
297      CONSOLE_DEBUG("Created %d params", p->num_parms);      CONSOLE_DEBUG("Created %d params", p->num_parms);
# Line 295  int integrator_ida_params_default(Integr Line 304  int integrator_ida_params_default(Integr
304  */  */
305    
306  typedef int IdaFlagFn(void *,int *);  typedef int IdaFlagFn(void *,int *);
307  typedef const char *IdaFlagNameFn(int);  typedef char *IdaFlagNameFn(int);
308    
309  /* return 1 on success */  /* return 1 on success */
310  int integrator_ida_solve(  int integrator_ida_solve(
# Line 347  int integrator_ida_solve( Line 356  int integrator_ida_solve(
356      y0 = N_VNew_Serial(size);      y0 = N_VNew_Serial(size);
357      integrator_get_y(blsys,NV_DATA_S(y0));      integrator_get_y(blsys,NV_DATA_S(y0));
358    
359    #ifdef SOLVE_DEBUG
360      CONSOLE_DEBUG("RETRIEVING yp0");      CONSOLE_DEBUG("RETRIEVING yp0");
361    #endif
362    
363      yp0 = N_VNew_Serial(size);      yp0 = N_VNew_Serial(size);
364      integrator_get_ydot(blsys,NV_DATA_S(yp0));      integrator_get_ydot(blsys,NV_DATA_S(yp0));
365    
366    #ifdef SOLVE_DEBUG
367      N_VPrint_Serial(yp0);      N_VPrint_Serial(yp0);
368      CONSOLE_DEBUG("yp0 is at %p",&yp0);      CONSOLE_DEBUG("yp0 is at %p",&yp0);
369    #endif
370    
371      /* create IDA object */      /* create IDA object */
372      ida_mem = IDACreate();      ida_mem = IDACreate();
# Line 399  int integrator_ida_solve( Line 412  int integrator_ida_solve(
412      if(integrator_get_minstep(blsys)>0){      if(integrator_get_minstep(blsys)>0){
413          ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");          ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
414      }      }
415    
416        CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
417        IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
418    
419      /* there's no capability for setting *minimum* step size in IDA */      /* there's no capability for setting *minimum* step size in IDA */
420    
421    
# Line 430  int integrator_ida_solve( Line 447  int integrator_ida_solve(
447          flagfntype = "IDADENSE";          flagfntype = "IDADENSE";
448          flagfn = &IDADenseGetLastFlag;          flagfn = &IDADenseGetLastFlag;
449          flagnamefn = &IDADenseGetReturnFlagName;          flagnamefn = &IDADenseGetReturnFlagName;
   
450      }else{      }else{
451          /* remaining methods are all SPILS */          /* remaining methods are all SPILS */
452          CONSOLE_DEBUG("IDA SPILS");          CONSOLE_DEBUG("IDA SPILS");
# Line 521  int integrator_ida_solve( Line 537  int integrator_ida_solve(
537          flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);          flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
538      # endif      # endif
539    
540          if(flag!=IDA_SUCCESS){          switch(flag){
541              flag = -999;              case IDA_SUCCESS:
542              flag1 = (flagfn)(ida_mem,&flag);                  CONSOLE_DEBUG("Initial conditions solved OK");
543              if(flag1){                  break;
544                  ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error retrieving linear solver error code (%d)",flag1);  
545              }              case IDA_LSETUP_FAIL:
546                case IDA_LINIT_FAIL:
547              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (Last %s error: %s)",flagfntype,(flagnamefn)(flag));              case IDA_LSOLVE_FAIL:
548              return 0;              case IDA_NO_RECOVERY:
549          }/* else success */                  flag1 = -999;
550                    flag = (flagfn)(ida_mem,&flag1);
551                    if(flag){
552                        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to retrieve error code from %s (err %d)",flagfntype,flag);
553                        return 0;
554                    }
555                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);
556                    return 0;
557    
558          CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");              default:
559                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
560                    return 0;
561            }
562      }else{      }else{
563          CONSOLE_DEBUG("Not solving initial conditions (see IDA parameter 'calcic')");          CONSOLE_DEBUG("Not solving initial conditions (see IDA parameter 'calcic')");
564      }      }
# Line 627  int integrator_ida_fex(realtype tt, N_Ve Line 653  int integrator_ida_fex(realtype tt, N_Ve
653      char *relname;      char *relname;
654  #ifdef FEX_DEBUG  #ifdef FEX_DEBUG
655      char *varname;      char *varname;
656        char diffname[30];
657  #endif  #endif
658    
659      blsys = (IntegratorSystem *)res_data;      blsys = (IntegratorSystem *)res_data;
# Line 698  int integrator_ida_fex(realtype tt, N_Ve Line 725  int integrator_ida_fex(realtype tt, N_Ve
725  #ifdef FEX_DEBUG  #ifdef FEX_DEBUG
726      /* output residuals to console */      /* output residuals to console */
727      CONSOLE_DEBUG("RESIDUAL OUTPUT");      CONSOLE_DEBUG("RESIDUAL OUTPUT");
728      fprintf(stderr,"index\t%20s\t%20s\t%s\n","y","ydot","resid");      fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
729      for(i=0; i<blsys->n_y; ++i){      for(i=0; i<blsys->n_y; ++i){
730          varname = var_make_name(blsys->system,blsys->y[i]);          varname = var_make_name(blsys->system,blsys->y[i]);
731          fprintf(stderr,"%d\t%10s=%10f\t",i,varname,NV_Ith_S(yy,i));          fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
732          if(blsys->ydot[i]){          if(blsys->ydot[i]){
733              varname = var_make_name(blsys->system,blsys->ydot[i]);              varname = var_make_name(blsys->system,blsys->ydot[i]);
734              fprintf(stderr,"%10s=%10f\t",varname,NV_Ith_S(yp,i));              fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
735          }else{          }else{
736              fprintf(stderr,"diff(%4s)=%10f\t",varname,NV_Ith_S(yp,i));              snprintf(diffname,99,"diff(%s)",varname);
737                fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
738          }          }
739          ASC_FREE(varname);          ASC_FREE(varname);
740          relname = rel_make_name(blsys->system,enginedata->rellist[i]);          relname = rel_make_name(blsys->system,enginedata->rellist[i]);

Legend:
Removed from v.990  
changed lines
  Added in v.991

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