/[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 944 by johnpye, Sat Nov 25 10:46:13 2006 UTC revision 945 by johnpye, Sat Nov 25 12:41:03 2006 UTC
# Line 56  Line 56 
56  /* SUNDIALS includes */  /* SUNDIALS includes */
57  #ifdef ASC_WITH_IDA  #ifdef ASC_WITH_IDA
58  # include <sundials/sundials_config.h>  # include <sundials/sundials_config.h>
59    # include <sundials/sundials_dense.h>
60  # include <ida/ida.h>  # include <ida/ida.h>
61  # include <nvector/nvector_serial.h>  # include <nvector/nvector_serial.h>
62  # include <ida/ida_spgmr.h>  # include <ida/ida_spgmr.h>
63    # include <ida/ida_dense.h>
64  # ifndef IDA_SUCCESS  # ifndef IDA_SUCCESS
65  #  error "Failed to include SUNDIALS IDA header file"  #  error "Failed to include SUNDIALS IDA header file"
66  # endif  # endif
# Line 111  void integrator_ida_error(int error_code Line 113  void integrator_ida_error(int error_code
113          , char *msg, void *eh_data          , char *msg, void *eh_data
114  );  );
115    
116    int integrator_ida_djex(long int Neq, realtype tt
117            , N_Vector yy, N_Vector yp, N_Vector rr
118            , realtype c_j, void *jac_data, DenseMat Jac
119            , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
120    );
121    
122  /*-------------------------------------------------------------  /*-------------------------------------------------------------
123    SETUP/TEARDOWN ROUTINES    SETUP/TEARDOWN ROUTINES
124  */  */
# Line 147  IntegratorIdaData *integrator_ida_engine Line 155  IntegratorIdaData *integrator_ida_engine
155  */  */
156    
157  enum ida_parameters{  enum ida_parameters{
158      IDA_PARAM_AUTODIFF      IDA_PARAM_LINSOLVER
159        ,IDA_PARAM_AUTODIFF
160      ,IDA_PARAM_RTOL      ,IDA_PARAM_RTOL
161      ,IDA_PARAM_ATOL      ,IDA_PARAM_ATOL
162      ,IDA_PARAM_ATOLVECT      ,IDA_PARAM_ATOLVECT
# Line 209  int integrator_ida_params_default(Integr Line 218  int integrator_ida_params_default(Integr
218              ,"Scalar relative error tolerance",1              ,"Scalar relative error tolerance",1
219              ,"Value of the scalar relative error tolerance."              ,"Value of the scalar relative error tolerance."
220              " See IDA manual, section 5.5.1"              " See IDA manual, section 5.5.1"
221          }, 1e-5, DBL_MIN, DBL_MAX }          }, 1e-4, DBL_MIN, DBL_MAX }
222        );
223    
224        slv_param_char(p,IDA_PARAM_LINSOLVER
225                ,(SlvParameterInitChar){{"linsolver"
226                ,"Linear solver",1
227                ,"See IDA manual, section 5.5.3."
228            }, "DENSE"}, (char *[]){"DENSE","SPGMR",NULL}
229      );      );
230    
231      asc_assert(p->num_parms == IDA_PARAMS_SIZE);      asc_assert(p->num_parms == IDA_PARAMS_SIZE);
# Line 234  int integrator_ida_solve( Line 250  int integrator_ida_solve(
250      realtype t0, reltol, abstol, t, tret, tout1;      realtype t0, reltol, abstol, t, tret, tout1;
251      N_Vector y0, yp0, abstolvect, ypret, yret;      N_Vector y0, yp0, abstolvect, ypret, yret;
252      IntegratorIdaData *enginedata;      IntegratorIdaData *enginedata;
253        char *linsolver;
254    
255      CONSOLE_DEBUG("STARTING IDA...");      CONSOLE_DEBUG("STARTING IDA...");
256    
# Line 291  int integrator_ida_solve( Line 308  int integrator_ida_solve(
308          CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);          CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
309          abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);          abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
310          flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);          flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
   
311      }      }
312    
313      if(flag==IDA_MEM_NULL){      if(flag==IDA_MEM_NULL){
# Line 316  int integrator_ida_solve( Line 332  int integrator_ida_solve(
332      CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");      CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
333    
334      /* attach linear solver module, using the default value of maxl */      /* attach linear solver module, using the default value of maxl */
335      flag = IDASpgmr(ida_mem, 0);      linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
336      if(flag==IDASPILS_MEM_NULL){      if(strcmp(linsolver,"SPGMR")==0){
337          ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");          flag = IDASpgmr(ida_mem, 0);
         return 0;  
     }else if(flag==IDASPILS_MEM_FAIL){  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");  
         return 0;  
     }/* else success */  
   
     /* assign the J*v function */  
     if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){  
         CONSOLE_DEBUG("USING AUTODIFF");  
         flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);  
338          if(flag==IDASPILS_MEM_NULL){          if(flag==IDASPILS_MEM_NULL){
339              ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");              ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
340              return 0;              return 0;
341          }else if(flag==IDASPILS_LMEM_NULL){          }else if(flag==IDASPILS_MEM_FAIL){
342              ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");              ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
343              return 0;              return 0;
344          }/* else success */          }/* else success */
345    
346            /* assign the J*v function */
347            if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
348                CONSOLE_DEBUG("USING AUTODIFF");
349                flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
350                if(flag==IDASPILS_MEM_NULL){
351                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
352                    return 0;
353                }else if(flag==IDASPILS_LMEM_NULL){
354                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
355                    return 0;
356                }/* else success */
357            }else{
358                CONSOLE_DEBUG("USING NUMERICAL DIFF");
359            }
360        }else if(strcmp(linsolver,"DENSE")==0){
361            CONSOLE_DEBUG("USING IDADEENSE SOLVER, size = %d",size);
362            flag = IDADense(ida_mem, size);
363            switch(flag){
364                case IDADENSE_SUCCESS: break;
365                case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 0;
366                case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 0;
367                case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;
368                default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;
369            }
370            
371            if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
372                CONSOLE_DEBUG("USING AUTODIFF");
373                flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
374                switch(flag){
375                    case IDADENSE_SUCCESS: break;
376                    default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 0;
377                }
378            }else{
379                CONSOLE_DEBUG("USING NUMERICAL DIFF");
380            }
381      }else{      }else{
382          CONSOLE_DEBUG("USING NUMERICAL DIFF");          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
383      }                return 0;
384        }
385    
386      /* set linear solver optional inputs...      /* set linear solver optional inputs...
387    
# Line 505  int integrator_ida_fex(realtype tt, N_Ve Line 548  int integrator_ida_fex(realtype tt, N_Ve
548  }  }
549    
550  /**  /**
551        Dense Jacobian evaluation. Only suitable for small problems!
552    */
553    int integrator_ida_djex(long int Neq, realtype tt
554            , N_Vector yy, N_Vector yp, N_Vector rr
555            , realtype c_j, void *jac_data, DenseMat Jac
556            , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
557    ){
558        IntegratorSystem *blsys;
559        IntegratorIdaData *enginedata;
560        realtype *yval;
561    
562        blsys = (IntegratorSystem *)jac_data;
563        enginedata = integrator_ida_enginedata(blsys);
564    
565        /* pass the values of everything back to the compiler */
566        integrator_set_t(blsys, (double)tt);
567        integrator_set_y(blsys, NV_DATA_S(yy));
568        integrator_set_ydot(blsys, NV_DATA_S(yp));
569    
570        yval = NV_DATA_S(yy);
571    
572        /* AWFUL HACK! In an attempt to prove that 'IDADense' works as expected,
573        I've pulled the jacobians straight from idadenx.c (an IDA example) */
574    # define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1)
575        IJth(Jac,1,1) = RCONST(-0.04) - c_j;
576        IJth(Jac,2,1) = RCONST(0.04);
577        IJth(Jac,3,1) = 1.0;
578        IJth(Jac,1,2) = RCONST(1.0e4)*yval[2];
579        IJth(Jac,2,2) = RCONST(-1.0e4)*yval[2] - RCONST(6.0e7)*yval[1] - c_j;
580        IJth(Jac,3,2) = 1.0;
581        IJth(Jac,1,3) = RCONST(1.0e4)*yval[1];
582        IJth(Jac,2,3) = RCONST(-1.0e4)*yval[1];
583        IJth(Jac,3,3) = 1.0;
584    #undef IJth
585        return(0);
586    }
587    
588    /**
589      Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)      Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
590    
591      Given tt, yy, yp, rr and v, we need to evaluate and return Jv.      Given tt, yy, yp, rr and v, we need to evaluate and return Jv.

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

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