poincare.cu

#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();
}