- Main Page
- Namespaces
- Data Structures
- Files

- Search for |

- File List
- Globals

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