Page MenuHomeHEPForge
Files F8308620

ScannerSRGEtools.cpp
No OneTemporary
Actions

ScannerSRGEtools.cpp
View Options

#include "ScannerSRGEtools.h"
using namespace std;
////////////////////////////////////////////////////////////////////////
//////////// RGE evolution routines definition //////////
////////////////////////////////////////////////////////////////////////
int func(double x, const double y[], double dydx[],void *params_point){
internal_paramsRGE par =*( (internal_paramsRGE *) params_point);
vector<double> yref(y,y+par.dim);///........??? Vector references figure out a way to assign data in array to reference of vector double without copy. If necessary create reference inside params which is prepared once in the begining and then points to the correct location in memory
vector<double> dydxref(dydx,dydx+par.dim); // Same here!!!!
(par.derivatives)(x,yref,dydxref,par.params);
for(int i=0;i<par.dim;++i)
dydx[i]=dydxref[i];
return GSL_SUCCESS;
}
int dummyjac (double x, const double y[], double *dydx_dy, double dydx_dr[], void *params_point){
cout <<"ScannerSRGEtools.cpp: in dummyjac, this was not supposed to be called at all!!!!"<< endl;
exit(EXIT_FAILURE);
return GSL_SUCCESS;
}
// Basic evolution routine.
void EvolveDirect(const double ti,const double tf,const std::vector<double> & yi,std::vector<double> & yf,void (*derivatives)(const double,const std::vector<double> & y,std::vector<double> & dydx,void *),void * params,const double AbsError,const double RelErr){
//USING integrators which do not require jacobian!!!!!!!!
const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
gsl_odeiv_step * s = gsl_odeiv_step_alloc(T,yi.size());
gsl_odeiv_control * c = gsl_odeiv_control_y_new(AbsError,RelErr);
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc(yi.size());
internal_paramsRGE par={yi.size(),&params,derivatives};
gsl_odeiv_system sys = {func,dummyjac,yi.size(),&par};
double * ytemp;
ytemp= new double [yi.size()];
double x=ti;
double h = 1e-5;
yf=yi;
while (x < tf){
int status = gsl_odeiv_evolve_apply(e,c,s,&sys,&x,tf,&h,&yf[0]);
if(status != GSL_SUCCESS){
std::clog<< "ScannerSRGEtools.cpp: In void EvolveDirect(const double ti, const double tf, const std::vector<double> & yi, std::vector<double> & yf, void (*derivatives)(const double , const std::vector<double> & y,std::vector<double> & dydx , void *), void * params, const double AbsError, const double RelErr): GSL failed."<< endl;
break;
}
}
gsl_odeiv_evolve_free (e);
gsl_odeiv_control_free (c);
gsl_odeiv_step_free (s);
//delete [] ytemp;
}
//Basic evolution that stops if check fails
double EvolvePlusCheckDirect(const double ti,const double tf,const std::vector<double> & yi,std::vector<double> & yf,void (*derivatives)(const double,const std::vector<double> & y,std::vector<double> & dydx,void *),void * params1,bool (*check)(const std::vector<double> & y,void *),void * params2,const double AbsError,const double RelErr){
//USING integrators which do not require jacobian!!!!!!!!
const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
gsl_odeiv_step * s = gsl_odeiv_step_alloc(T,yi.size());
gsl_odeiv_control * c = gsl_odeiv_control_y_new(AbsError,RelErr);
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc(yi.size());
internal_paramsRGE par={yi.size(),&params1,derivatives};
gsl_odeiv_system sys = {func,dummyjac,yi.size(),&par};
double t=ti;
double h = 1e-5;
yf=yi;
while (t < tf && !isNaNorInf(yf) && check(yf,params2)){
int status = gsl_odeiv_evolve_apply(e,c,s,&sys,&t,tf,&h,&yf[0]);
if(status != GSL_SUCCESS){
std::clog<< "ScannerSRGEtools.cpp: In double EvolvePlusCheckDirect(const double ti,const double tf,const std::vector<double> & yi,std::vector<double> & yf,void (*derivatives)(const double,const std::vector<double> & y,std::vector<double> & dydx,void *),void * params1,bool (*check)(const std::vector<double> & y,void *),void * params2,const double AbsError,const double RelErr): GSL failed."<< endl;
break;
}
}
gsl_odeiv_evolve_free (e);
gsl_odeiv_control_free (c);
gsl_odeiv_step_free (s);
return t;
}
//Basic evolution that stops if check fails, keeping the middle points
double EvolvePlusCheckKeepAll(const double ti,const double tf,std::vector<double> & t_evolve,const std::vector<double> & yi,std::vector< std::vector<double> > & y_evolve,void (*derivatives)(const double,const std::vector<double> & y,std::vector<double> & dydx,void *),void * params1,bool (*check)(const std::vector<double> & y,void *),void * params2,const double AbsError,const double RelErr){
//USING integrators which do not require jacobian!!!!!!!!
const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
gsl_odeiv_step * s = gsl_odeiv_step_alloc(T,yi.size());
gsl_odeiv_control * c = gsl_odeiv_control_y_new(AbsError,RelErr);
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc(yi.size());
internal_paramsRGE par={yi.size(),&params1,derivatives};
gsl_odeiv_system sys = {func,dummyjac,yi.size(),&par};
double h=1e-5;
t_evolve.push_back(ti);
y_evolve.push_back(yi);
size_t index=0;
while (t_evolve[index]<tf && check(y_evolve[index],params2)){
int status=gsl_odeiv_evolve_apply(e,c,s,&sys,&t_evolve[index],tf,&h,&y_evolve[index][0]);
if(status!=GSL_SUCCESS){
std::clog<< "ScannerSRGEtools.cpp: In double EvolvePlusCheckKeepAll(const double ti,const double tf,const std::vector<double> & yi,std::vector<double> & yf,void (*derivatives)(const double,const std::vector<double> & y,std::vector<double> & dydx,void *),void * params1,bool (*check)(const std::vector<double> & y,void *),void * params2,const double AbsError,const double RelErr): GSL failed."<< endl;
break;
}
if(isNaNorInf(y_evolve[index]))
break;
t_evolve.push_back(t_evolve[index]);
y_evolve.push_back(y_evolve[index]);
++index;
}
gsl_odeiv_evolve_free (e);
gsl_odeiv_control_free (c);
gsl_odeiv_step_free (s);
return t_evolve.back();
}
//Evolution storing couplings at npoints points; if yevolution is not an empty vector, it starts writing after the last element.
bool EvolveBySteps(const double ti,const double tf,std::vector<double> & tevolution,const std::vector<double> & yi,std::vector< std::vector<double> > & yevolution,void (*derivatives)(const double ,const std::vector<double> & y,std::vector<double> & dydx , void *),void * params, const double AbsError,const double RelErr,const int npoints){
double step=(tf-ti)/npoints; //time step
double t=ti+step; //time after each evolution
int index=yevolution.size(); //so it starts writing on the 1st empty slot
std::clog<<"inside evolvebysteps"<<endl;
//only the 1st step uses yi
yevolution.push_back(vector<double>(yi.size())); //creates a vector of the apropriate size
EvolveDirect(ti,t,yi,yevolution[index],derivatives,params,AbsError,RelErr);
tevolution.push_back(t); //stores the time evolution
if(isNaNorInf(yevolution[index])){
std::clog<<"ScannerS Warning: 1 point could not be evolved with the RGE's until the final time. Integration produced a NaN or inf before that."<<endl;
return false; //returns false if it wasn't able to fully evolve the couplings
}
++index;
//after that every step used the previous one as initial state
for(int i=1;i!=npoints;++i){
yevolution.push_back(vector<double>(yi.size())); //creates a vector of the apropriate size
EvolveDirect(t,t+step,yevolution[index-1],yevolution[index],derivatives,params,AbsError,RelErr);
tevolution.push_back(t+step); //stores the time evolution
t+=step;
if(isNaNorInf(yevolution[index])){
std::clog<<"ScannerS Warning: 1 point could not be evolved with the RGE's until the final time. Integration produced a NaN or inf before that."<<endl;
return false; //returns false if it wasn't able to fully evolve the couplings
}
++index;
}
return true; //if the point can be fully evolved return true
}
// Evolution with nchecks points, stopping when checks fail. Returns t correpsonding to failure.
double EvolvePlusCheckDirect(const double ti,const double tf,const std::vector<double> & yi,std::vector<double> & yf,void (*derivatives)(const double,const std::vector<double> & y,std::vector<double> & dydx,void *),void * params1,bool (*check)(const std::vector<double> & y,void *),void * params2,const double AbsError,const double RelErr,const int npoints){
double step=(tf-ti)/npoints; //time step
double t=ti+step; //time after each evolution
vector<double> ytemp(yf.size());
//only the 1st step uses yi
EvolveDirect(ti,t,yi,yf,derivatives,params1,AbsError,RelErr);
if(isNaNorInf(yf)){
std::clog<<"ScannerS Warning: 1 point could not be evolved with the RGE's until the final time. Integration produced a NaN or inf before that."<<endl;
return t; //returns false if it wasn't able to fully evolve the couplings
}
if(!check(yf,params2)){
return t;
}
//after that every step used the previous one as initial state
for(int i=1;i!=npoints;++i){
ytemp=yf;
EvolveDirect(t,t+step,ytemp,yf,derivatives,params1,AbsError,RelErr);
t+=step;
if(isNaNorInf(yf)){
std::clog<<"ScannerS Warning: 1 point could not be evolved with the RGE's until the final time. Integration produced a NaN or inf before that."<<endl;
return t; //returns false if it wasn't able to fully evolve the couplings
}
if(!check(yf,params2)){
return t;
}
}
return t;
}
// Evolution storing couplings at nchecks points, stopping when checks fail. Returns t correpsonding to failure.
double EvolvePlusCheckBySteps(const double ti,const double tf,std::vector<double> & tevolution,const std::vector<double> & yi,std::vector< std::vector<double> > & yevolution,void (*derivatives)(const double ,const std::vector<double> & y,std::vector<double> & dydx,void *),void * params1,bool (*check)(const std::vector<double> & y , void *),void * params2,const double AbsError,const double RelErr,const int npoints){
double step=(tf-ti)/npoints; //time step
double t=ti+step; //time after each evolution
int index=yevolution.size(); //so it starts writing on the 1st empty slot
//only the 1st step uses yi
yevolution.push_back(vector<double>(yi.size())); //creates a vector of the apropriate size
EvolveDirect(ti,t,yi,yevolution[index],derivatives,params1,AbsError,RelErr);
tevolution.push_back(t); //stores the time evolution
if(isNaNorInf(yevolution[index])){
std::clog<<"ScannerS Warning: 1 point could not be evolved with the RGE's until the final time. Integration produced a NaN or inf before that."<<endl;
return t; //returns false if it wasn't able to fully evolve the couplings
}
if(!check(yevolution[index],params2)){
return t;
}
++index;
//after that every step used the previous one as initial state
for(int i=1;i!=npoints;++i){
yevolution.push_back(vector<double>(yi.size())); //creates a vector of the apropriate size
EvolveDirect(t,t+step,yevolution[index-1],yevolution[index],derivatives,params1,AbsError,RelErr);
tevolution.push_back(t+step); //stores the time evolution
t+=step;
if(isNaNorInf(yevolution[index])){
std::clog<<"ScannerS Warning: 1 point could not be evolved with the RGE's until the final time. Integration produced a NaN or inf before that."<<endl;
return t; //returns false if it wasn't able to fully evolve the couplings
}
if(!check(yevolution[index],params2)){
return t;
}
++index;
}
return t;
}
//-------------------------------------------------//
const double NaN_doubles[]={sqrt(-1),-sqrt(-1),1./0.,-1./0.};
const unsigned char *point2nan=(unsigned char*)&NaN_doubles;
bool isNaNorInf(const double var){
//This function basicly servers the same purpose as (!isfinite()), but since the c++ standards do NOT demand compilers to handle "not numbers" when optimized (thus rendering functions like isfinite(),isnan(),etc... useless) we decided to build our own
const unsigned char* myvar=(unsigned char*)&var;
for(int i=0;i!=4;++i){
bool nancheck=true;
for(int j=0;j!=sizeof(double);++j)
if(*(myvar+j) != *(point2nan+j+i*sizeof(double)))
nancheck=false;
if(nancheck)
return true;
}
return false;
}
bool isNaNorInf(const std::vector<double> & y){
//checks if the vector has any nan or inf members
for(int i=0;i<y.size();++i)
if(isNaNorInf(y[i]))
return true;
return false;
}

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 12:38 PM (1 d, 18 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4020147
Default Alt Text
ScannerSRGEtools.cpp (12 KB)

Event Timeline