/[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 991 by johnpye, Thu Dec 21 10:44:32 2006 UTC revision 992 by johnpye, Thu Dec 21 13:37:45 2006 UTC
# Line 110  typedef struct{ Line 110  typedef struct{
110      struct bnd_boundary **bndlist;   /**< NULL-terminated list of boundaries, for use in the root-finding  code */      struct bnd_boundary **bndlist;   /**< NULL-terminated list of boundaries, for use in the root-finding  code */
111      int nrels;      int nrels;
112      int safeeval;                    /**< whether to pass the 'safe' flag to relman_eval */      int safeeval;                    /**< whether to pass the 'safe' flag to relman_eval */
113        N_Vector pp;                     /**< Preconditioner values */
114  } IntegratorIdaData;  } IntegratorIdaData;
115    
116  /*-------------------------------------------------------------  /*-------------------------------------------------------------
# Line 135  int integrator_ida_djex(long int Neq, re Line 136  int integrator_ida_djex(long int Neq, re
136          , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3          , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
137  );  );
138    
139    int integrator_ida_psetup(realtype tt,
140             N_Vector yy, N_Vector yp, N_Vector rr,
141             realtype c_j, void *prec_data,
142             N_Vector tmp1, N_Vector tmp2,
143             N_Vector tmp3
144    );
145    
146    int integrator_ida_psolve(realtype tt,
147             N_Vector yy, N_Vector yp, N_Vector rr,
148             N_Vector rvec, N_Vector zvec,
149             realtype c_j, realtype delta, void *prec_data,
150             N_Vector tmp
151    );
152    
153    static const IDASpilsPrecSetupFn psetup1 = &integrator_ida_psetup;
154    static const IDASpilsPrecSolveFn psolve1 = &integrator_ida_psolve;
155    
156  /*-------------------------------------------------------------  /*-------------------------------------------------------------
157    SETUP/TEARDOWN ROUTINES    SETUP/TEARDOWN ROUTINES
158  */  */
# Line 145  void integrator_ida_create(IntegratorSys Line 163  void integrator_ida_create(IntegratorSys
163      enginedata->rellist = NULL;      enginedata->rellist = NULL;
164      enginedata->varlist = NULL;      enginedata->varlist = NULL;
165      enginedata->safeeval = 0;      enginedata->safeeval = 0;
166        enginedata->pp = NULL;
167      blsys->enginedata = (void *)enginedata;      blsys->enginedata = (void *)enginedata;
168      integrator_ida_params_default(blsys);      integrator_ida_params_default(blsys);
169  }  }
# Line 152  void integrator_ida_create(IntegratorSys Line 171  void integrator_ida_create(IntegratorSys
171  void integrator_ida_free(void *enginedata){  void integrator_ida_free(void *enginedata){
172      CONSOLE_DEBUG("DELETING IDA ENGINE DATA");      CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
173      IntegratorIdaData *d = (IntegratorIdaData *)enginedata;      IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
174        if(d->pp)N_VDestroy_Serial(d->pp);
175      /* note, we don't own the rellist, so don't need to free it */      /* note, we don't own the rellist, so don't need to free it */
176      ASC_FREE(d);      ASC_FREE(d);
177  }  }
# Line 180  enum ida_parameters{ Line 200  enum ida_parameters{
200      ,IDA_PARAM_ATOLVECT      ,IDA_PARAM_ATOLVECT
201      ,IDA_PARAM_GSMODIFIED      ,IDA_PARAM_GSMODIFIED
202      ,IDA_PARAM_MAXNCF      ,IDA_PARAM_MAXNCF
203        ,IDA_PARAM_PREC
204          ,IDA_PARAMS_SIZE          ,IDA_PARAMS_SIZE
205  };  };
206    
# Line 187  enum ida_parameters{ Line 208  enum ida_parameters{
208      Here the full set of parameters is defined, along with upper/lower bounds,      Here the full set of parameters is defined, along with upper/lower bounds,
209      etc. The values are stuck into the blsys->params structure.      etc. The values are stuck into the blsys->params structure.
210    
211        To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
212        list. Then add a slv_param_*(...) statement below to define the type, description and range
213        for the new parameter.
214    
215      @return 0 on success      @return 0 on success
216  */  */
217  int integrator_ida_params_default(IntegratorSystem *blsys){  int integrator_ida_params_default(IntegratorSystem *blsys){
# Line 292  int integrator_ida_params_default(Integr Line 317  int integrator_ida_params_default(Integr
317          }, 10,0,1000 }          }, 10,0,1000 }
318      );      );
319    
320        slv_param_char(p,IDA_PARAM_PREC
321                ,(SlvParameterInitChar){{"prec"
322                ,"Preconditioner",1
323                ,"See IDA manual, section section 5.6.8."
324            },"NONE"}, (char *[]){"NONE","DIAG",NULL}
325        );
326    
327      asc_assert(p->num_parms == IDA_PARAMS_SIZE);      asc_assert(p->num_parms == IDA_PARAMS_SIZE);
328    
329      CONSOLE_DEBUG("Created %d params", p->num_parms);      CONSOLE_DEBUG("Created %d params", p->num_parms);
# Line 322  int integrator_ida_solve( Line 354  int integrator_ida_solve(
354      IdaFlagFn *flagfn;      IdaFlagFn *flagfn;
355      IdaFlagNameFn *flagnamefn;      IdaFlagNameFn *flagnamefn;
356      const char *flagfntype;      const char *flagfntype;
357        char *pname;
358        IDASpilsPrecSetupFn psetup;
359        IDASpilsPrecSolveFn psolve;
360    
361      CONSOLE_DEBUG("STARTING IDA...");      CONSOLE_DEBUG("STARTING IDA...");
362    
# Line 454  int integrator_ida_solve( Line 489  int integrator_ida_solve(
489          maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);          maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);
490          CONSOLE_DEBUG("maxl = %d",maxl);          CONSOLE_DEBUG("maxl = %d",maxl);
491    
492            pname = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_PREC);
493            if(strcmp(pname,"NONE")==0){
494                psetup=NULL;
495                psolve=NULL;
496            }else if(strcmp(pname,"DIAG")==0){
497                psetup=&integrator_ida_psetup;
498                psolve=&integrator_ida_psolve;
499            }
500    
501          if(strcmp(linsolver,"SPGMR")==0){          if(strcmp(linsolver,"SPGMR")==0){
502              CONSOLE_DEBUG("IDA SPGMR");              CONSOLE_DEBUG("IDA SPGMR");
503              flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */              flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
# Line 468  int integrator_ida_solve( Line 512  int integrator_ida_solve(
512              return 0;              return 0;
513          }          }
514    
515            enginedata->pp = N_VNew_Serial(blsys->n_y);
516            IDASpilsSetPreconditioner(ida_mem,psetup,psolve,(void *)blsys);
517            CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
518    
519          flagfntype = "IDASPILS";          flagfntype = "IDASPILS";
520          flagfn = &IDASpilsGetLastFlag;          flagfn = &IDASpilsGetLastFlag;
521          flagnamefn = &IDASpilsGetReturnFlagName;          flagnamefn = &IDASpilsGetReturnFlagName;
# Line 1056  int integrator_ida_jvex(realtype tt, N_V Line 1104  int integrator_ida_jvex(realtype tt, N_V
1104  }  }
1105    
1106  /*----------------------------------------------  /*----------------------------------------------
1107      PRECONDITIONER
1108    */
1109    
1110    
1111    /**
1112        EXPERIMENTAL. A diagonal preconditioner for use with IDA Krylov solvers
1113    
1114        'setup' function.
1115    */
1116    int integrator_ida_psetup(realtype tt,
1117             N_Vector yy, N_Vector yp, N_Vector rr,
1118             realtype c_j, void *p_data,
1119             N_Vector tmp1, N_Vector tmp2,
1120             N_Vector tmp3
1121    ){
1122        int i;
1123        IntegratorSystem *blsys;
1124        IntegratorIdaData *data;
1125        blsys = (IntegratorSystem *)p_data;
1126        data = blsys->enginedata;
1127    
1128        double y = 1; /* derivative of y[i] in relation i */
1129        double ydot = 1; /* derivative of ydot[i] in relation i */
1130    
1131        N_VConst(1, data->pp);
1132        for(i=0; i<blsys->n_y; ++i){
1133            /* @TODO calculate y, ydot here */
1134            NV_Ith_S(data->pp, i) = 1 / (y + c_j * ydot);
1135        }
1136        return 0;
1137    };
1138    
1139    /**
1140        EXPERIMENTAL. A diagonal preconditioner for use with IDA Krylov solvers
1141    
1142        'solve' function.
1143    */
1144    int integrator_ida_psolve(realtype tt,
1145             N_Vector yy, N_Vector yp, N_Vector rr,
1146             N_Vector rvec, N_Vector zvec,
1147             realtype c_j, realtype delta, void *p_data,
1148             N_Vector tmp
1149    ){
1150        IntegratorSystem *blsys;
1151        IntegratorIdaData *data;
1152        blsys = (IntegratorSystem *)p_data;
1153        data = blsys->enginedata;
1154    
1155        N_VProd(data->pp, rvec, zvec);
1156        return 0;
1157    };
1158    
1159    /*----------------------------------------------
1160    ERROR REPORTING    ERROR REPORTING
1161  */  */
1162  /**  /**

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

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