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

#include <iostream>
#include <fstream>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "../def.h"
#include "../model.h"
#include "../nr.h"

Go to the source code of this file.

Functions

void tabulateValues (Glob *globs, GlobExp *ex, double t0, double t1, double *t, long n, double *state, double **y, double ***dmds, double ***dmdp, double **dyds, double **dydp)
void simInit (GlobExp *ex, Glob *globs)


Function Documentation

void simInit GlobExp ex,
Glob globs
 

Definition at line 21 of file simInit.cc.

References dmatrix(), dvector(), GlobExp::nPoints, GlobExp::nvar, GlobExp::y0, and GlobExp::yTry.

00022 {
00023   long i,j;
00024   long nvar=ex->nvar;
00025   long nPoints=ex->nPoints;
00026   long idum=-42;
00027   double *t=dvector(1,2);
00028   double **y=dmatrix(1,2,1,nvar);
00029   
00030    if (ex->y0) 
00031     {
00032       for (i=1; i<=nvar; ++i) 
00033         ex->yTry[1][i]=ex->y0[i];
00034     } 
00035    else
00036      {
00037         for (i=1; i<=nvar; ++i) 
00038           {
00039             ex->yTry[1][i]=DefYValues[i-1];
00040           }
00041      }
00042 
00043    //integrating 
00044 
00045    for(i=2;i<=nPoints;i++)
00046      {
00047        for(j=1;j<=nvar;j++)
00048          ex->yTry[i][j]=ex->yTry[i-1][j];
00049        t[1]=ex->mesh[i-1];
00050        t[2]=ex->mesh[i];
00051        tabulateValues(globs,ex,t[1],t[2],t,2,ex->yTry[i],y,NULL,NULL,NULL,NULL);
00052       for(j=1;j<=nvar;j++)
00053         ex->yTry[i][j]+=globs->pert*ex->yTry[i][j]*gasdev(&idum);
00054      }
00055    
00056 
00057 
00058    free_dvector(t,1,2);
00059    free_dmatrix(y,1,2,1,nvar);
00060 }

void tabulateValues Glob globs,
GlobExp ex,
double  t0,
double  t1,
double *  t,
long  n,
double *  state,
double **  y,
double ***  dmds,
double ***  dmdp,
double **  dyds,
double **  dydp
 

Definition at line 108 of file intODE.cc.

References call_odessa(), dfill0(), dmatrix(), Glob::doP, dvector(), Glob::eps, initYt(), integrateRK(), Glob::integrator, max(), Glob::maxstp, Glob::npar, GlobExp::nvar, GlobExp::par, Glob::rkqs_ign, sensDerivs(), Glob::stiff, and TRUE.

00110 {
00111   // integrate from t0 to t1, storing values at t[1..n]
00112   // t[n] may be ==t1 or <t1
00113   // in the latter case we will integrate from t[n] to t1 in the
00114   // final loop k=n+1 without storing any data points
00115   // If dmds!=NULL, compute derivatives as well
00116   // and store them in dmds,dmdp,dyds,dydp
00117   // with respect to observation, see odin.doc: Observation function
00118   
00119   //BEGIN von Felix veraendert
00120   //   double hmin=max(fabs(t0),fabs(t1))*1e-8;
00121   double hmin = max (fabs (t0), fabs (t1)) * 1e-15;
00122   // Das legt die minimale Schrittweite, die ODESSA erlaubt wird fest.
00123   // hmin sollte so gross sein, dass es noch nicht zu Abschneidefehlern
00124   // beim Aufaddieren zum Zeitpunkt kommt, deshalb die max(...)-Konstruktion.
00125   // Der Faktor 1e-15 gewaehrleistet, dass noch richtig kleine Schritte gemacht
00126   // werden duerfen, aber man trotzdem von der rel. Maschinengenauigkeit (im Bereich
00127   // von 1e-18 noch einigermassen weit entfernt bleibt. 
00128   // Man koennte das in Zukunft noch in einen externen Parameter stecken.
00129   //END von Felix veraendert
00130 
00131   long nPExp=globs->npar;
00132   long nVar=ex->nvar;
00133   long nobs=NOBS;
00134   int generic;
00135   double *p=ex->par;
00136   double eps=globs->eps;  
00137   double h1 = (t1 - t0) * .01;
00138   long nok,nbad;
00139 
00140   initInt(globs,ex);
00141   if (n > 0 && (t0 > t[1] || t1 < t[n]))
00142     {
00143       cerr << "tabVal: t[1..n] exceeds [t0;t1]" << endl;
00144       throw 1;
00145     }
00146 
00147   if (globs->integrator==3)
00148     {                           // said -odessa
00149 #ifdef ODESSA
00150       double &rtol = eps;
00151 //BEGIN von Felix veraendert
00152 //       double atol=eps*1e-6;  
00153       double atol = eps;
00154 // Fehlerkontrolle fuer ODESSA: Das Spiel mit relativem und absolutem LOKALEN(!)
00155 // Fehler sollte man unbedingt noch in externe Parameter stecken. Doku dazu
00156 // siehe odessa.f: ATOL und RTOL
00157 //END von Felix veraendert
00158       // odeint uses TINY=1e-30 in a similar way as odessa's atol, 
00159       // but atol<1e-16 or so results in step size underflow
00160       
00161       int odeint_MAXSTP=globs->maxstp;
00162       int d_flag=1;
00163       int stif=globs->stiff;
00164       char doPStr[globs->npar+1];
00165       
00166       for(long i=1;i<=globs->npar;i++)
00167         doPStr[i]=(char)globs->doP[i];
00168       
00169       call_odessa(globs,ex,nVar, nPExp, doPStr,
00170          d_flag, odeint_MAXSTP, stif,
00171          TRUE, n, t, t0, t1, state, y,
00172          dmds, dmdp, dyds, dydp, p, rtol, atol, hmin, t1 - t0, h1);
00173 
00174       return;
00175 #else
00176       cerr << "-odessa used but not linked" << endl;
00177       throw 1;
00178 #endif // ODESSA
00179     }                           // if odessa
00180       double *gy=dvector(1,nobs);
00181       long i, j, k;
00182       double xx1 = t0, xx2;
00183       
00184       if (dmds)
00185         {                           // compute sensitivities
00186           // create extended state vector Y
00187           double **Yt = dmatrix (0, nPExp + nVar, 1, nVar);
00188           initYt (Yt[0], state);
00189           
00190           // chain rule for extra observables
00191           double **dgdy, **dgdp;
00192           
00193           dgdy = dmatrix (1, nobs, 1, nVar);
00194           dgdp = dmatrix (1, nobs, 1, nPExp);
00195           dfill0 (dgdy[1], nobs * nVar);
00196           dfill0 (dgdp[1], nobs * nPExp);
00197           
00198           
00199           for (k = 1; k <= n + 1; k++)
00200             {
00201               xx2 = (k > n) ? t1 :  // end point
00202                 t[k];               // data point
00203               
00204               globs->rkqs_ign = (nPExp + nVar) * nVar;
00205               if (xx2 > xx1 && nVar > 0)
00206                 {                 
00207                   if(globs->integrator==1)
00208                     {
00209                       //Runge-Kutta integation
00210                       integrateRK(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,
00211                                   eps, h1, hmin, &nok, &nbad, sensDerivs);
00212                     }
00213                   else if(globs->integrator==2)
00214                     {
00215                       cerr << "CVODES not optimized skip integrator" << endl;
00216                       throw 1;
00217             
00218                       //integrateCVODES(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,eps,TRUE);
00219                     }
00220                 }
00221               // write state to y and sensitivities to dm/dp and dm/ds
00222               if (k <= n)
00223                 {
00224                   dcopy (y[k], Yt[0], nVar);
00225                   generic=observation(globs,ex,xx2, y[k],gy, p, dgdy, dgdp);
00226                   for(i=1;i<=nobs;i++)
00227                     y[k][i]=gy[i];
00228                   
00229                   for (i = 1; i <= nobs; i++)
00230                     {
00231                       for (j = 1; j <= nPExp; j++)
00232                         {
00233                           double &dest = dmdp[k][i][j];
00234                           if (!generic)
00235                             {
00236                               dest = dgdp[i][j];
00237                               for (long l = 1; l <= nVar; l++)
00238                                 dest += dgdy[i][l] * Yt[j][l];
00239                             }
00240                           else
00241                             dest = Yt[j][i];
00242                         }           // for j
00243                       for (j = 1; j <= nVar; j++)
00244                         {
00245                           double &dest = dmds[k][i][j];
00246                           if (!generic)
00247                             {
00248                               dest = 0;
00249                               for (long l = 1; l <= nVar; l++)
00250                                 dest += dgdy[i][l] * Yt[nPExp + j][l];
00251                             }
00252                           else
00253                             dest = Yt[nPExp + j][i];
00254                         }           // for j
00255                     }               // for i
00256                 }                   // if k<=n
00257               xx1 = xx2;
00258             }                       // for k
00259           
00260           // copy state and its derivatives to state, dydp and dyds
00261           dcopy (state, Yt[0], nVar);       // original state vector
00262           // write sensitivities to dy/dp and dy/ds
00263           for (i = 1; i <= nVar; i++)
00264             {
00265               for (j = 1; j <= nPExp; j++)
00266                 dydp[i][j] = Yt[j][i];
00267               for (j = 1; j <= nVar; j++)
00268                 dyds[i][j] = Yt[nPExp + j][i];
00269             }
00270           
00271           free_dmatrix (dgdp,1 ,nobs , 1, nPExp);
00272           free_dmatrix (dgdy,1 ,nobs, 1, nVar);
00273           free_dmatrix (Yt, 0, nPExp + nVar, 1, nVar);
00274         }
00275       else
00276         {                           // !dmds
00277           for (k = 1; k <= n + 1; k++)
00278             {
00279               xx2 = (k <= n) ? t[k] : t1;
00280               
00281               globs->rkqs_ign = 0;
00282               if (xx2 > xx1 && nVar > 0)
00283                 { 
00284                   if(globs->integrator==1)
00285                     {
00286                       //Runge-Kutta integration
00287                       integrateRK(globs,state, p, nVar, xx1, xx2,
00288                                   eps, h1, hmin, &nok, &nbad, ode);
00289                     }
00290                   else if(globs->integrator==2)
00291                     {
00292                       cerr << "CVODES not optimized skip integrator" << endl;
00293                       throw 1;
00294 
00295                       //integrateCVODES(globs,state, p, nVar, xx1, xx2,eps,FALSE);
00296                     }
00297                 }
00298               if (k <= n)
00299                 {
00300                   dcopy (y[k], state, nVar);
00301                   generic=observation (globs,ex,xx2, y[k],gy,p, NULL, NULL);
00302                   for(i=1;i<=nobs;i++)
00303                     y[k][i]=gy[i];
00304                 }
00305               xx1 = xx2;
00306             }                       // for k
00307         }           
00308       // if dmds else
00309       free_dvector(gy,1,nobs); 
00310 }


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