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

#include <math.h>
#include <iostream>
#include "../nr.h"
#include "../def.h"
#include "../model.h"
#include <string.h>
#include "../libCVODES/sundialstypes.h"
#include "../libCVODES/cvodes.h"
#include "../libCVODES/cvdense.h"
#include "../libCVODES/nvector_serial.h"
#include "../libCVODES/dense.h"
#include "../libCVODES/sundialsmath.h"

Go to the source code of this file.

Defines

#define Ith(v, i)   NV_Ith_S(v,i-1)
#define IJth(A, i, j)   DENSE_ELEM(A,i-1,j-1)
#define RTOL   RCONST(1e-8)
#define ATOL   RCONST(1e-8)

Functions

void f (realtype t, N_Vector y, N_Vector ydot, void *f_data)
void sensDerivs (double t, double *Yt, double *Ft, double *p)
void integrateCVODES (Glob *globs, double ystart[], double p[], long nvar, double x1, double x2, double eps, int sensi)


Define Documentation

#define ATOL   RCONST(1e-8)
 

Definition at line 22 of file integrateCVODES.cc.

Referenced by integrateCVODES().

#define IJth A,
i,
 )     DENSE_ELEM(A,i-1,j-1)
 

Definition at line 19 of file integrateCVODES.cc.

#define Ith v,
 )     NV_Ith_S(v,i-1)
 

Definition at line 18 of file integrateCVODES.cc.

Referenced by f(), and integrateCVODES().

#define RTOL   RCONST(1e-8)
 

Definition at line 21 of file integrateCVODES.cc.

Referenced by integrateCVODES().


Function Documentation

void f realtype  t,
N_Vector  y,
N_Vector  ydot,
void *  f_data
 

Definition at line 98 of file integrateCVODES.cc.

References Ith, sensDerivs(), and TRUE.

Referenced by integrateCVODES(), and objectiveF().

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 }

void integrateCVODES Glob globs,
double  ystart[],
double  p[],
long  nvar,
double  x1,
double  x2,
double  eps,
int  sensi
 

Definition at line 40 of file integrateCVODES.cc.

References ATOL, f(), Ith, and RTOL.

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 }

void sensDerivs double  t,
double *  Yt,
double *  Ft,
double *  p
 

Definition at line 72 of file intODE.cc.

References dfill0(), and dmatrix().

Referenced by f(), and tabulateValues().

00073 {
00074   //  sensitivity equations
00075   // P=(p,y0)
00076   // S=dy/dP
00077   // Y=(y,S)
00078   // F=dY/dt=(f,df/dy*S+df/dP)
00079   // Y and F are stored in column major format and rolled out as Yt and Ft
00080 
00081   long i, j, k;
00082   dfill0 (Ft, (1 + NPARAMS + NEQNS) * NEQNS);       // initialize to zero
00083   // ode, inhomo and jacobi set only non-zero elements
00084   // inhome expects a dmatrix; we create one who's data lies in Ft
00085   double **dfdpt = new double *[NPARAMS] - 1;
00086   for (i = 1; i <= NPARAMS; i++)
00087     dfdpt[i] = Ft + i * NEQNS;
00088   ode (t, Yt, Ft, p);           // first part of Y,F are y,f
00089   inhomo (t, Yt, dfdpt, p);     // fill df/dp into F
00090   delete[]++ dfdpt;
00091 
00092   // now F=df/dP
00093   // add M*S, M=df/dy
00094   double **M = dmatrix (1, NEQNS, 1, NEQNS);
00095   dfill0 (M[1], NEQNS * NEQNS);
00096   jacobi (t, Yt, M, p);         // get M
00097   for (i = 1; i <= NEQNS; i++)
00098     for (j = 1; j <= NPARAMS + NEQNS; j++)
00099       {
00100         double &sum = Ft[j * NEQNS + i];
00101         for (k = 1; k <= NEQNS; k++)
00102           sum += M[k][i] * Yt[j * NEQNS + k];
00103       }                         // for i,j
00104   free_dmatrix (M, 1, NEQNS, 1, NEQNS);
00105  }                               // sensDerivs


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