- Main Page
- Namespaces
- Data Structures
- Files

- File List
- Globals

# D:/Basile/VisualCppProjects/hydroMEX3/hydroMEX3/f.cpp

Go to the documentation of this file.

    00001
    00012 /*
    00013 Wrap function for solver in dx=f(t,x,parameters)
    00014 */
    00015
    00016
    00017
    00018 #include "math.h"
    00019 #include "mex.h"
    00020
    00021 #include "string.h"
    00022
    00023 #include "hydro_source/hydro_includes.h"
    00024
    00025 #include "CVODEincludes.h"  // ODE solver includes
    00026
    00027
    00028 #include "f.h"
    00029
    00030
    00031 #ifndef PI_DEF
    00032 #define pi 3.141592653589793
    00033 #define PI 3.141592653589793
    00034 #endif
    00035
    00036
    00037
    00038
    00039 #define Ith(v,i)    NV_Ith_S(v,i-1)       /* Ith numbers components 1..NEQ */
    00040 #define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */
    00041
    00042 /* Problem Constants */
    00043
    00044 #define NEQ   12//18                /* number of equations  */
    00045 #define T0    RCONST(0.0)      /* initial time           */
    00046
    00047
    00048
    00049
    00050
    00051
    00052 //  dy/dt = f(t,y,parameters)
    00053
    00055
    00063 static int f(realtype t, N_Vector y, N_Vector ydot, void *f_data)
    00064 {
    00065         datas *d;
    00066         parameters *P;
    00067         ODE_data *f_dat;
    00068         double *state;
    00069
    00070         f_dat = (ODE_data*)f_data;
    00071
    00072         d = f_dat->d;
    00073         P = f_dat->P;
    00074         state = f_dat->state;
    00075
    00076
    00077
    00078         state[0] = Ith(y,1);
    00079         state[1] = Ith(y,2);
    00080         state[2] = Ith(y,3);
    00081         state[3] = Ith(y,4);
    00082         state[4] = Ith(y,5);
    00083         state[5] = Ith(y,6);
    00084         state[6] = Ith(y,7);
    00085         state[7] = Ith(y,8);
    00086         state[8] = Ith(y,9);
    00087         state[9] = Ith(y,10);
    00088         state[10] = Ith(y,11);
    00089         state[11] = Ith(y,12);
    00090         /*state[12] = Ith(y,13); //unused camera stuff
    00091         state[13] = Ith(y,14);
    00092         state[14] = Ith(y,15);
    00093         state[15] = Ith(y,16);
    00094         state[16] = Ith(y,17);
    00095         state[17] = Ith(y,18);*/
    00096
    00097         data_update_state(d, state);
    00098
    00099         compute_state_dot(d,P,state);
    00100
    00101         Ith(ydot,1) = d->dx;
    00102         Ith(ydot,2) = d->dy;
    00103         Ith(ydot,3) = d->dz;
    00104         Ith(ydot,4) = d->dphi;
    00105         Ith(ydot,5) = d->dtheta;
    00106         Ith(ydot,6) = d->dpsi;
    00107         Ith(ydot,7) = d->ddq[0];
    00108         Ith(ydot,8) = d->ddq[1];
    00109         Ith(ydot,9) = d->ddq[2];
    00110         Ith(ydot,10) = d->ddq[3];
    00111         Ith(ydot,11) = d->ddq[4];
    00112         Ith(ydot,12) = d->ddq[5];
    00113         /*Ith(ydot,13) = d->dx_cam;  //unused camera stuff
    00114         Ith(ydot,14) = d->dy_cam;
    00115         Ith(ydot,15) = d->dz_cam;
    00116         Ith(ydot,16) = d->dd_cam[0];
    00117         Ith(ydot,17) = d->dd_cam[1];
    00118         Ith(ydot,18) = d->dd_cam[2];*/
    00119
    00120 /*
    00121         Ith(ydot,1) = f_dat->d->t;  //integrator test:  dy/dt = t
    00122         Ith(ydot,2) = f_dat->d->t;
    00123         Ith(ydot,3) = f_dat->d->t;
    00124         Ith(ydot,4) = f_dat->d->t;
    00125         Ith(ydot,5) = f_dat->d->t;
    00126         Ith(ydot,6) = f_dat->d->t;
    00127         Ith(ydot,7) = f_dat->d->t;
    00128         Ith(ydot,8) = f_dat->d->t;
    00129         Ith(ydot,9) = f_dat->d->t;
    00130         Ith(ydot,10) = f_dat->d->t;
    00131         Ith(ydot,11) = f_dat->d->t;
    00132         Ith(ydot,12) = f_dat->d->t;
    00133         Ith(ydot,13) = f_dat->d->t;
    00134         Ith(ydot,14) = f_dat->d->t;
    00135         Ith(ydot,15) = f_dat->d->t;
    00136         Ith(ydot,16) = f_dat->d->t;
    00137         Ith(ydot,17) = f_dat->d->t;
    00138         Ith(ydot,18) = f_dat->d->t;
    00139 */
    00140
    00141   return(0);
    00142 }
    00143
    00144
    00145
    00146
    00147
    00149
    00154 void solver_init(ODE_data* f_dat)
    00155 {
    00156         int k;
    00157
    00158         f_dat->y = NULL;
    00159         f_dat->abstol = NULL;
    00160         f_dat->cvode_mem = NULL;
    00161
    00162   /* Create serial vector of length NEQ for I.C. and abstol */
    00163         f_dat->y = N_VNew_Serial(NEQ);
    00164         f_dat->abstol = N_VNew_Serial(NEQ);
    00165
    00166   /* Initialize y */
    00167         for (k=1; k<=NEQ; k++)
    00168         {
    00169                 Ith(f_dat->y,k) = f_dat->state[k-1];
    00170                 //printf("Ith(y,%d) = %f \n",k,Ith(y,k));
    00171         }
    00172
    00173
    00174   /* Set the scalar relative tolerance */
    00175         f_dat->reltol = f_dat->P->reltol;
    00176   /* Set the vector absolute tolerance */
    00177         for (k=1; k<=NEQ; k++)
    00178         {
    00179                 Ith(f_dat->abstol,k) = f_dat->P->abstol[k-1];
    00180         }
    00181
    00182
    00183           /*
    00184      Call CVodeCreate to create the solver memory:
    00185
    00186      CV_BDF     specifies the Backward Differentiation Formula
    00187      CV_NEWTON  specifies a Newton iteration
    00188
    00189      A pointer to the integrator problem memory is returned and stored in cvode_mem.
    00190   */
    00191
    00192         f_dat->cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
    00193
    00194
    00195   /*
    00196      Call CVodeMalloc to initialize the integrator memory:
    00197
    00198      cvode_mem is the pointer to the integrator memory returned by CVodeCreate
    00199      f         is the user's right hand side function in y'=f(t,y)
    00200      T0        is the initial time
    00201      y         is the initial dependent variable vector
    00202      CV_SV     specifies scalar relative and vector absolute tolerances
    00203      &reltol   is a pointer to the scalar relative tolerance
    00204      abstol    is the absolute tolerance vector
    00205   */
    00206
    00207         f_dat->flag = CVodeMalloc(f_dat->cvode_mem, f, T0, f_dat->y, CV_SV, f_dat->reltol, f_dat->abstol);
    00208
    00209         f_dat->flag = CVodeSetFdata(f_dat->cvode_mem, (void*)f_dat);
    00210
    00211         /* Call CVodeRootInit to specify the root function g with 2 components */
    00212         //flag = CVodeRootInit(cvode_mem, 2, g, NULL);
    00213
    00214         /* Call CVDense to specify the CVDENSE dense linear solver */
    00215         f_dat->flag = CVDense(f_dat->cvode_mem, NEQ);
    00216         //flag = CVDiag(cvode_mem);
    00217
    00218         /* Set the Jacobian routine to Jac (user-supplied) */
    00219         //flag = CVDenseSetJacFn(cvode_mem, Jac, NULL);
    00220 }
    00221
    00222
    00223
    00224
    00225
    00226
    00228
    00233 void solve_ODE(ODE_data* f_dat, realtype tout)
    00234 {
    00235         int k, tempk;
    00236         long ti, nst;
    00237
    00238         //integrate to next time value:
    00239         f_dat->flag = CVode(f_dat->cvode_mem, tout, f_dat->y, &f_dat->d->t, CV_NORMAL);
    00240
    00241         //store results and time if necessary:
    00242         if (f_dat->time_param != 0.0) store_state(f_dat->d->t, f_dat->y);
    00243
    00244
    00245         f_dat->flag = CVodeGetNumSteps(f_dat->cvode_mem, &nst);
    00246         f_dat->CumulNumOfSteps = (double)nst;
    00247 }
    00248
    00249
    00251
    00255 void solver_free(ODE_data* f_dat)
    00256 {
    00257         /* Free y vector */
    00258         N_VDestroy_Serial(f_dat->y);
    00259
    00260         /* Free integrator memory */
    00261         CVodeFree(&f_dat->cvode_mem);
    00262 }