- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /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