/home/peifer/diffit/modules/integrateCVODES.cc

Go to the documentation of this file.
00001 #include <math.h>
00002 #include <iostream>
00003 #include "../nr.h"
00004 #include "../def.h"
00005 #include "../model.h"
00006 #include <string.h>
00007 #include "../libCVODES/sundialstypes.h"   
00008 #include "../libCVODES/cvodes.h"          
00009 #include "../libCVODES/cvdense.h"         
00010 #include "../libCVODES/nvector_serial.h"  
00011 #include "../libCVODES/dense.h"           
00012 #include "../libCVODES/sundialsmath.h"
00013 
00014 using namespace std;
00015 
00016 /* Accessor macros */
00017 
00018 #define Ith(v,i)    NV_Ith_S(v,i-1)       // i-th vector component
00019 #define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) // (i,j)-th matrix component
00020 
00021 #define RTOL  RCONST(1e-8) 
00022 #define ATOL RCONST(1e-8) 
00023 
00024 
00025 //user data
00026 typedef struct {
00027   double *p;
00028   long nvar;
00029   int sensi;
00030 } *UserData;
00031 
00032 
00033 //prototypes 
00034 void f(realtype t, N_Vector y, N_Vector ydot, void *f_data);
00035 
00036 static int ewt(N_Vector y, N_Vector w, void *e_data);
00037 
00038 void sensDerivs (double t, double *Yt, double *Ft, double *p);
00039 
00040 void integrateCVODES(Glob *globs,double ystart[], double p[], long nvar, double x1, 
00041                      double x2,double eps,int sensi)
00042 {
00043   
00044   void *cvode_mem;
00045   UserData data;
00046   realtype t,t0,t1;
00047   N_Vector y;
00048   long i;
00049 
00050   cvode_mem = NULL;
00051   data      = NULL;
00052   y         = NULL;
00053 
00054  
00055   data = (UserData) malloc(sizeof *data);
00056   
00057   //user data
00058   data->p=p;
00059   data->nvar=nvar;
00060   data->sensi=sensi;
00061   
00062   //  Initial conditions 
00063   y = N_VNew_Serial(nvar);
00064   for(i=1;i<=nvar;i++)
00065     Ith(y,i)=ystart[i];
00066 
00067   t0= RCONST(x1);
00068   t1= RCONST(x2);
00069   t = RCONST(t1);
00070 
00071   /* Create CVODES object */
00072   cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
00073   /* Allocate space for CVODES */
00074   realtype atol=ATOL;
00075   CVodeMalloc(cvode_mem, f, t0, y, CV_SS, RTOL, &atol);
00076   /* Attach user data */
00077   CVodeSetFdata(cvode_mem, data);
00078   /* Attach linear solver */
00079   CVDense(cvode_mem, nvar);
00080   
00081   //integrate
00082   CVode(cvode_mem, t1, y, &t, CV_NORMAL);
00083 
00084   //copy output back
00085   for(i=1;i<=nvar;i++)
00086     ystart[i]=Ith(y,i);
00087 
00088   //free memory
00089   N_VDestroy_Serial(y);
00090   free(data);
00091   CVodeFree(cvode_mem);
00092 }
00093 
00094 /*
00095  * f routine. Compute f(t,y). 
00096  */
00097 
00098 void f(realtype t, N_Vector y, N_Vector ydot, void *f_data)
00099 {
00100   UserData data;
00101   data = (UserData) f_data;
00102   long i,nvar=data->nvar;
00103   double *ydot_nr=(double*)malloc((nvar+1)*sizeof(double));
00104   double *y_nr=(double*)malloc((nvar+1)*sizeof(double));
00105 
00106   for(i=1;i<=nvar;i++)
00107     y_nr[i]=Ith(y,i);
00108   
00109   if(data->sensi==TRUE)
00110     sensDerivs((double)t,y_nr,ydot_nr,data->p);
00111   else
00112     ode((double)t,y_nr,ydot_nr,data->p);
00113   
00114   for(i=1;i<=nvar;i++)
00115       Ith(ydot,i)=ydot_nr[i];
00116 
00117   free(y_nr);
00118   free(ydot_nr);
00119 
00120 }
00121 
00122  

Generated on Mon Jan 29 17:09:12 2007 for Diffit by  doxygen 1.4.6