Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#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
You can’t perform that action at this time.