/[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 908 by johnpye, Wed Oct 25 13:07:12 2006 UTC revision 909 by johnpye, Thu Oct 26 12:44:41 2006 UTC
# Line 74  Line 74 
74  */  */
75  typedef struct{  typedef struct{
76      struct rel_relation **rellist;   /**< NULL terminated list of rels */      struct rel_relation **rellist;   /**< NULL terminated list of rels */
77        struct var_variable **varlist;   /**< NULL terminated list of rels */
78      int nrels;      int nrels;
79      int safeeval;                    /**< whether to pass the 'safe' flag to relman_eval */      int safeeval;                    /**< whether to pass the 'safe' flag to relman_eval */
80  } IntegratorIdaData;  } IntegratorIdaData;
# Line 84  typedef struct{ Line 85  typedef struct{
85  /* residual function forward declaration */  /* residual function forward declaration */
86  int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);  int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
87    
88    int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
89            , N_Vector v, N_Vector Jv, realtype c_j
90            , void *jac_data, N_Vector tmp1, N_Vector tmp2
91    );
92    
93  /* error handler forward declaration */  /* error handler forward declaration */
94  void integrator_ida_error(int error_code  void integrator_ida_error(int error_code
95          , const char *module, const char *function          , const char *module, const char *function
# Line 98  void integrator_ida_create(IntegratorSys Line 104  void integrator_ida_create(IntegratorSys
104      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
105      enginedata = ASC_NEW(IntegratorIdaData);      enginedata = ASC_NEW(IntegratorIdaData);
106      enginedata->rellist = NULL;      enginedata->rellist = NULL;
107        enginedata->varlist = NULL;
108      enginedata->safeeval = 1;      enginedata->safeeval = 1;
109      blsys->enginedata = (void *)enginedata;      blsys->enginedata = (void *)enginedata;
110  }  }
# Line 143  int integrator_ida_solve( Line 150  int integrator_ida_solve(
150      /* store reference to list of relations (in enginedata) */      /* store reference to list of relations (in enginedata) */
151      enginedata->nrels = slv_get_num_solvers_rels(blsys->system);      enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
152      enginedata->rellist = slv_get_solvers_rel_list(blsys->system);      enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
153        enginedata->varlist = slv_get_solvers_var_list(blsys->system);
154      CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);      CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
155      CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);      CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
156      size = blsys->n_y;      size = blsys->n_y;
# Line 203  int integrator_ida_solve( Line 211  int integrator_ida_solve(
211          return 0;          return 0;
212      }/* else success */      }/* else success */
213    
214      /* set linear solver optional inputs... */      /* assign the J*v function */
215        flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
216        if(flag==IDASPILS_MEM_NULL){
217            ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
218            return 0;
219        }else if(flag==IDASPILS_LMEM_NULL){
220            ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
221            return 0;
222        }/* else success */
223    
224        /* set linear solver optional inputs...
225    
226            ...nothing here at the moment...
227    
228        */
229    
230      /* correct initial values, given derivatives */      /* correct initial values, given derivatives */
231      blsys->currentstep=0;      blsys->currentstep=0;
# Line 356  int integrator_ida_fex(realtype tt, N_Ve Line 377  int integrator_ida_fex(realtype tt, N_Ve
377  }  }
378    
379  /**  /**
380      @TODO implement calculation of the jacobian, right?      Function to evaluate the product J*v, in the form required for IDA.
381    
382        Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
383    
384        @param tt current value of the independent variable (time, t)
385        @param yy current value of the dependent variable vector, y(t).
386        @param yp current value of y'(t).
387        @param rr current value of the residual vector F(t, y, y').
388        @param v  the vector by which the Jacobian must be multiplied to the right.
389        @param Jv the output vector computed
390        @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
391        @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
392        @param tmp1 @see tmp2
393        @param tmp2 (as well as tmp1) pointers to memory allocated for variables of type N_Vector for use here as temporary storage or work space.
394        @return 0 on success
395  */  */
396    int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
397            , N_Vector v, N_Vector Jv, realtype c_j
398            , void *jac_data, N_Vector tmp1, N_Vector tmp2
399    ){
400        IntegratorSystem *blsys;
401        IntegratorIdaData *enginedata;
402        int i, j, is_error=0;
403        struct rel_relation** relptr;
404        char *relname, *varname;
405        int status;
406        double Jv_i;
407    
408        int *variables;
409        double *derivatives;
410        var_filter_t filter;
411        int count;
412    
413        CONSOLE_DEBUG("EVALUTING JACOBIAN...");
414    
415        blsys = (IntegratorSystem *)jac_data;
416        enginedata = integrator_ida_enginedata(blsys);
417    
418        /* pass the values of everything back to the compiler */
419        integrator_set_t(blsys, (double)tt);
420        integrator_set_y(blsys, NV_DATA_S(yy));
421        integrator_set_ydot(blsys, NV_DATA_S(yp));
422        /* no real use for residuals (rr) here, I don't think? */
423    
424        /* allocate space for returns from relman_diff2 */
425        variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
426        derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
427    
428        /* evaluate the derivatives... */
429        /* J = dG_dy = dF_dy + alpha * dF_dyp */
430    
431        filter.matchbits = VAR_SVAR;
432        filter.matchvalue = VAR_SVAR;
433    
434        Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
435        if (setjmp(g_fpe_env)==0) {
436            for(i=0, relptr = enginedata->rellist;
437                    i< enginedata->nrels && relptr != NULL;
438                    ++i, ++relptr
439            ){
440                /* get derivatives for this particular relation */
441                status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
442    
443                CONSOLE_DEBUG("Got derivatives against %d matching variables", count);
444    
445                relname = rel_make_name(blsys->system, *relptr);
446                if(!status){
447                    fprintf(stderr,"\n\n");
448                    CONSOLE_DEBUG("Derivatives for relation %d '%s' OK",i,relname);
449                }else{
450                    CONSOLE_DEBUG("ERROR calculating derivatives for relation %d '%s'",i,relname);
451                    break;
452                }
453                ASC_FREE(relname);
454    
455                Jv_i = 0;
456                for(j=0; j < count; ++j){
457                    CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
458                    varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
459                    if(varname){
460                        CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
461                        ASC_FREE(varname);
462                    }else{
463                        CONSOLE_DEBUG("Variable %d: derivative = %f",variables[j],derivatives[j]);
464                    }
465    
466                    Jv_i += derivatives[j] * NV_Ith_S(v,variables[j]);
467                }
468    
469                NV_Ith_S(Jv,i) = Jv_i;
470                if(status){
471                    /* presumably some error_reporter will already have been made*/
472                    is_error = 1;
473                }
474            }
475        }else{
476            CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
477        }
478        Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
479    
480        if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
481        return is_error;
482    }
483    
484  /*----------------------------------------------  /*----------------------------------------------
485    ERROR REPORTING    ERROR REPORTING

Legend:
Removed from v.908  
changed lines
  Added in v.909

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