Created
May 18, 2013 16:17
-
-
Save nazarov-yuriy/5604963 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <stdio.h> | |
#include <cuda_runtime.h> | |
#define DIM 2 //Dimension of the dynamical system <--second order as an example | |
#define EGOROVS_CONST 8 //Egolov's constant for integration step control | |
#define PI 3.1415926535897932384626433832795029L //Pi-magic-number is the base of the Universe | |
#define A1 5.0 | |
#define A2 15.0 | |
#define B2 8.0 | |
#define C1 4.0 | |
#define C2 7.0 | |
#define nu 1.0 | |
#define eps 0.60 | |
#define DELTA 3.0 | |
struct in_data{ | |
double x[DIM]; | |
double NumPointsOnMap; | |
double StepIntegr; | |
double TolIntegr; | |
double TolOfSection; | |
double G; | |
}; | |
struct out_data{ | |
double points[4000]; | |
}; | |
// Device code | |
__device__ double F(int number_of_equation, double time, double y[DIM], double G){ | |
if(number_of_equation == 0) | |
return y[1]*(1/C2-cos(y[0])*cos(y[0])/(A1+B2)-sin(y[0])*sin(y[0])/(A1+A2))-DELTA/C2-eps/C2/nu*sin(nu*time); | |
else | |
return -(G*G-y[1]*y[1])*(1/(A1+A2)-1/(A1+B2))*sin(y[0])*cos(y[0]); | |
} | |
__device__ double Section(double time){ | |
return sin((nu*time)/2);// repetition of the phase: nu*time=nu*time+2*Pi | |
} | |
__device__ double RK45(double t, double *T, double *x, double *X, double StepIntegr, double G){ | |
double k1[DIM],k2[DIM],k3[DIM],k4[DIM]; //arrays for RK45 and ControlTerm evaluation | |
double x_k1[DIM],x_k2[DIM],x_k3[DIM]; | |
int i; | |
double value=0; | |
double temp=0; | |
for (i=0; i<DIM; i++) | |
{ | |
k1[i]=StepIntegr*F(i, t, x, G); | |
x_k1[i]=x[i]+0.5*k1[i]; | |
} | |
for (i=0; i<DIM; i++) | |
{ | |
k2[i]=StepIntegr*F(i,t+0.5*StepIntegr,x_k1, G); | |
x_k2[i]=x[i]+0.5*k2[i]; | |
} | |
for (i=0; i<DIM; i++) | |
{ | |
k3[i]=StepIntegr*F(i,t+0.5*StepIntegr,x_k2, G); | |
x_k3[i]=x[i]+k3[i]; | |
} | |
for (i=0; i<DIM; i++) | |
{ | |
k4[i]=StepIntegr*F(i,t+StepIntegr,x_k3, G); | |
} | |
for (i=0; i<DIM; i++) | |
{ | |
X[i]=x[i]+(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6; | |
} | |
*T=t+StepIntegr; | |
for (i=0;i<DIM;i++)//evaluation of Egorov's control term: | |
{ | |
temp=fabs(k1[i]-k2[i]-k3[i]+k4[i])*2/3; | |
if (temp>value) { | |
value=temp; | |
} | |
} | |
return value; | |
} | |
__device__ void process_trajectory(const struct in_data *input, struct out_data *output){ | |
double t = 0; | |
double x[DIM], X[DIM];//arrays for start and final point of one iteration | |
double T; | |
int NumPointsOnMap_NOW = 0; | |
double rk_temp, StepTemp; | |
double StepIntegr = input->StepIntegr; | |
int i; | |
for(i=0;i<DIM;i++) | |
x[i]=input->x[i]; | |
while (NumPointsOnMap_NOW<input->NumPointsOnMap) //finding of each point of the Poincare-map | |
{//1 | |
//---------------ONE ITERATION WITH CONTROL TERM--------------- | |
rk_temp=RK45(t, &T, x, X, StepIntegr, input->G);//execute of one iteration of Runge-Kutta method | |
while (rk_temp>input->TolIntegr) | |
{//autochanging of StepIntegr_SIZE | |
StepIntegr/=2; | |
rk_temp=RK45(t, &T, x, X, StepIntegr, input->G); | |
} | |
while (rk_temp<=input->TolIntegr/EGOROVS_CONST) | |
{ | |
StepIntegr*=2; | |
rk_temp=RK45(t, &T, x, X, StepIntegr, input->G); | |
}//autochanging of StepIntegr_SIZE */ | |
//----------------END OF ITERATION----------------------------- | |
//-----------INTERSECTION FINDING------------------------------ | |
double SectionBefore, SectionAfter;// values of section function befor and after intersection | |
SectionBefore=Section(t); | |
SectionAfter=Section(T); | |
if (SectionBefore*SectionAfter<0)//then we have intersection on this iteration-interval | |
{//2 | |
if (SectionBefore>0)//(X[3]>x[3]) | |
{// direction of trajectory transpotation control | |
StepTemp=StepIntegr;//Let's save in memory step size befor finding | |
while (fabs(SectionAfter)>input->TolOfSection) //correction of tolerance of point | |
{//3 | |
if (SectionBefore*SectionAfter<0) StepIntegr*=0.5; else StepIntegr*=1.5; | |
RK45(t, &T, x, X, StepIntegr, input->G); | |
SectionAfter=Section(T); | |
//This method corespond to "half-secant" ideology | |
}//3 | |
StepIntegr=StepTemp;//restor step size | |
output->points[NumPointsOnMap_NOW*DIM] = fmod(static_cast<double>(fmod(static_cast<double>(x[0]),static_cast<double>(2*PI))+static_cast<double>(2*PI)),static_cast<double>(2*PI));; | |
output->points[NumPointsOnMap_NOW*DIM+1] = x[1]/input->G; | |
NumPointsOnMap_NOW++; | |
}// direction of trajectory transpotation control | |
}//2 | |
//----CHANGING OF START-FINAL POINTS FOR FURTHER ITERATIONS----- | |
for (i=0;i<DIM;i++) | |
x[i]=X[i]; | |
t=T; | |
}//1 | |
} | |
__global__ void VecAdd(const struct in_data *A, struct out_data *C, int N) { | |
int i = blockDim.x * blockIdx.x + threadIdx.x; | |
if (i < N){ | |
process_trajectory(&A[i], &C[i]); | |
} | |
} | |
#define N 160 | |
// Host code | |
void read_x(FILE *f_ini, double *x){ | |
char temp[500]; | |
int i; | |
fscanf(f_ini,"%s",temp); | |
fscanf(f_ini,"%s",temp); | |
int bufint; | |
fscanf(f_ini,"%d",&bufint); | |
for (i=0;i<DIM;i++) | |
{ | |
fscanf(f_ini,"%lf", &x[i]); | |
} | |
fscanf(f_ini,"%s",temp); | |
} | |
void read_some_params(FILE *f_ini, int *NumTrajectories, int *NumPointsOnMap, double *StepIntegr, double *TolIntegr, double *TolOfSection, double *G){ | |
char temp[500]; | |
fscanf(f_ini,"%s",temp); | |
fscanf(f_ini,"%s",temp); | |
fscanf(f_ini,"%s",temp); | |
fscanf(f_ini,"%d %d %lf %lf %lf\n",NumTrajectories, NumPointsOnMap, StepIntegr, TolIntegr, TolOfSection); | |
fscanf(f_ini,"%s",temp); | |
fscanf(f_ini,"%s",temp); | |
fscanf(f_ini,"%lf \n", G); | |
} | |
int main(int argc, char** argv) { | |
struct in_data* host_in; | |
struct out_data* host_out; | |
struct in_data* device_in; | |
struct out_data* device_out; | |
host_in = (struct in_data*) malloc(N*sizeof(*host_in)); | |
host_out = (struct out_data*) malloc(N*sizeof(*host_out)); | |
// Allocate vectors in device memory | |
cudaMalloc((void**) &device_in, N*sizeof(*device_in)); | |
cudaMalloc((void**) &device_out, N*sizeof(*device_out)); | |
if (cudaGetLastError()) | |
printf("allocation failed\n"); | |
FILE *f_ini, *f_lL_res; | |
int NumTrajectories; | |
double StepIntegr; | |
double TolIntegr, TolOfSection; | |
int NumPointsOnMap; //number of points on Puancare map for one trajectory | |
double G; | |
int k; | |
f_lL_res = fopen("./RESULT/lL.txt","w"); | |
f_ini = fopen("ini.txt","r"); | |
read_some_params(f_ini, &NumTrajectories, &NumPointsOnMap, &StepIntegr, &TolIntegr, &TolOfSection, &G); | |
for(k=1;k<=NumTrajectories;k++){ | |
printf("Processing trajectrory %d\n",k ); | |
host_in[k-1].NumPointsOnMap = NumPointsOnMap; | |
host_in[k-1].StepIntegr = StepIntegr; | |
host_in[k-1].G = G; | |
host_in[k-1].TolIntegr = TolIntegr; | |
host_in[k-1].TolOfSection = TolOfSection; | |
read_x(f_ini, host_in[k-1].x); | |
} | |
//========================================================================== | |
printf("Vector Addition Started\n"); | |
cudaMemcpy(device_in, host_in, N*sizeof(*host_in), cudaMemcpyHostToDevice); | |
VecAdd<<<32, 5>>>(device_in, device_out, N); | |
if (cudaGetLastError()) | |
printf("processinf failed\n"); | |
cudaMemcpy(host_out, device_out, N*sizeof(*host_out), cudaMemcpyDeviceToHost); | |
cudaError_t e; | |
if (e=cudaGetLastError()) | |
printf("host <-copy- device failed: %s\n", cudaGetErrorString(e)); | |
printf("Vector Addition Ended\n"); | |
//========================================================================== | |
for(k=1;k<=NumTrajectories;k++){ | |
for(int i=0; i<NumPointsOnMap; i++){ | |
fprintf(f_lL_res,"%15.10f %15.10f\n", host_out[k-1].points[i*DIM], host_out[k-1].points[i*DIM+1]); | |
} | |
} | |
fclose(f_ini); | |
fclose(f_lL_res); | |
free(host_in); | |
free(host_out); | |
cudaFree(device_in); | |
cudaFree(device_out); | |
cudaDeviceReset(); | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment