249 lines
6.3 KiB
C++
249 lines
6.3 KiB
C++
#include <iostream>
|
|
|
|
#include <vector>
|
|
#include <sys/time.h>
|
|
|
|
#include "DKSBase.h"
|
|
|
|
#include <vector_types.h>
|
|
#include "cuda_runtime.h"
|
|
|
|
|
|
using namespace std;
|
|
|
|
typedef struct {
|
|
int label;
|
|
unsigned localID;
|
|
double Rincol[3];
|
|
double Pincol[3];
|
|
} PART_SMALL;
|
|
|
|
typedef struct {
|
|
double x;
|
|
double y;
|
|
double z;
|
|
} Vector;
|
|
|
|
PART_SMALL initPartSmall(int d) {
|
|
|
|
PART_SMALL p;
|
|
p.label = 0;
|
|
p.localID = d;
|
|
|
|
p.Rincol[0] = 0.0;
|
|
p.Rincol[1] = 0.0;
|
|
p.Rincol[2] = 0.02;
|
|
|
|
p.Pincol[0] = 0.0;
|
|
p.Pincol[1] = 0.0;
|
|
p.Pincol[2] = 3.9920183237269791e-01;
|
|
|
|
return p;
|
|
}
|
|
|
|
Vector initVector() {
|
|
Vector tmp;
|
|
tmp.x = 0.5;
|
|
tmp.y = 0.5;
|
|
tmp.z = 0.5;
|
|
|
|
return tmp;
|
|
}
|
|
|
|
void printPart(PART_SMALL p) {
|
|
cout << "label: " << p.label << ", ";
|
|
cout << "localid: " << p.localID << ",";
|
|
cout << "Rincol: " << p.Rincol[0] << ", " << p.Rincol[1] << ", " << p.Rincol[2] << ", ";
|
|
cout << "Pincol: " << p.Pincol[0] << ", " << p.Pincol[1] << ", " << p.Pincol[2];
|
|
cout << endl;
|
|
}
|
|
|
|
void printVector(Vector v) {
|
|
cout << v.x << "\t" << v.y << "\t" << v.z << endl;
|
|
}
|
|
|
|
void initParts(PART_SMALL *p, int N) {
|
|
for (int i = 0; i < N; i++)
|
|
p[i] = initPartSmall(i);
|
|
}
|
|
|
|
void printParts(PART_SMALL *p, int N) {
|
|
for (int i = 0; i < N; i++)
|
|
printPart(p[i]);
|
|
cout << endl;
|
|
}
|
|
|
|
void initVectors(Vector *v, int N) {
|
|
for (int i = 0; i < N; i++)
|
|
v[i] = initVector();
|
|
}
|
|
|
|
void printVectors(Vector *v, int N) {
|
|
for (int i = 0; i < N; i++)
|
|
printVector(v[i]);
|
|
cout << endl;
|
|
}
|
|
|
|
|
|
void initParams(double *data) {
|
|
data[0] = 0.0;//2.0000000000000000e-02;
|
|
data[1] = 1.0;//1.0000000000000000e-02;
|
|
data[2] = 2.2100000000000000e+00;
|
|
data[3] = 6.0000000000000000e+00;
|
|
data[4] = 1.2010700000000000e+01;
|
|
data[5] = 2.6010000000000000e+00;
|
|
data[6] = 1.7010000000000000e+03;
|
|
data[7] = 1.2790000000000000e+03;
|
|
data[8] = 1.6379999999999999e-02;
|
|
data[9] = 1.9321266968325795e-01;
|
|
data[10] = 7.9000000000000000e+01;
|
|
data[11] = 1.0000000000000002e-12;
|
|
|
|
}
|
|
|
|
void printDouble(double *data, int N) {
|
|
for (int i = 0; i < N; i++)
|
|
std::cout << data[i] << "\t";
|
|
std::cout << std::endl;
|
|
}
|
|
|
|
int main(int argc, char *argv[]) {
|
|
|
|
int loop = 10;
|
|
int numpart = 1e5;
|
|
char *api_name = new char[10];
|
|
char *device_name = new char[10];
|
|
strcpy(api_name, "Cuda");
|
|
strcpy(device_name, "-gpu");
|
|
|
|
for (int i = 1; i < argc; i++) {
|
|
|
|
if (argv[i] == string("-mic")) {
|
|
strcpy(api_name, "OpenMP");
|
|
strcpy(device_name, "-mic");
|
|
}
|
|
|
|
if (argv[i] == string("-npart")) {
|
|
numpart = atoi(argv[i+1]);
|
|
i++;
|
|
}
|
|
|
|
if (argv[i] == string("-loop")) {
|
|
loop = atoi(argv[i+1]);
|
|
i++;
|
|
}
|
|
|
|
}
|
|
|
|
cout << "=========================BEGIN TEST=========================" << endl;
|
|
cout << "Use api: " << api_name << "\t" << device_name << endl;
|
|
cout << "Number of particles: " << numpart << endl;
|
|
cout << "Number of loops: " << loop << endl;
|
|
cout << "------------------------------------------------------------" << endl;
|
|
|
|
//init part vector to test mc
|
|
PART_SMALL *parts = new PART_SMALL[numpart];
|
|
initParts(parts, numpart);
|
|
|
|
double *params = new double[12];
|
|
initParams(params);
|
|
|
|
//init dks
|
|
int ierr;
|
|
DKSBase base;
|
|
base.setAPI(api_name, strlen(api_name));
|
|
base.setDevice(device_name, strlen(api_name));
|
|
base.initDevice();
|
|
|
|
//init random
|
|
base.callInitRandoms(numpart);
|
|
|
|
//**test collimator physics and sort***//
|
|
void *part_ptr, *param_ptr;
|
|
|
|
//allocate memory for particles
|
|
part_ptr = base.allocateMemory<PART_SMALL>(numpart, ierr);
|
|
param_ptr = base.allocateMemory<double>(12, ierr);
|
|
|
|
//transfer data to device
|
|
base.writeData<PART_SMALL>(part_ptr, parts, numpart);
|
|
base.writeData<double>(param_ptr, params, 12);
|
|
|
|
int numaddback;
|
|
//test calls to do some first executions
|
|
base.callCollimatorPhysics2(part_ptr, param_ptr, numpart);
|
|
base.callCollimatorPhysicsSort(part_ptr, numpart, numaddback);
|
|
base.syncDevice();
|
|
//std::cout << "particles to add back: " << numaddback << std::endl;
|
|
|
|
struct timeval timeStart, timeEnd;
|
|
std::cout << "Start MC" << std::endl;
|
|
|
|
gettimeofday(&timeStart, NULL);
|
|
for (int i = 0; i < loop; i++) {
|
|
base.callCollimatorPhysics2(part_ptr, param_ptr, numpart);
|
|
base.callCollimatorPhysicsSort(part_ptr, numpart, numaddback);
|
|
base.syncDevice();
|
|
}
|
|
gettimeofday(&timeEnd, NULL);
|
|
|
|
std::cout << "addback: " << numaddback << std::endl;
|
|
|
|
std::cout << "End MC" << std::endl;
|
|
double t = ( (timeEnd.tv_sec - timeStart.tv_sec) * 1000000 +
|
|
(timeEnd.tv_usec - timeStart.tv_usec));
|
|
|
|
std::cout << "Time for " << loop << " MC runs: " << t * 1e-6 << "s" << std::endl;
|
|
std::cout << "Average time for MC run: " << t * 1e-6 / loop << std::endl;
|
|
|
|
//read data from device
|
|
base.readData<PART_SMALL>(part_ptr, parts, numpart);
|
|
|
|
//free memory
|
|
base.freeMemory<PART_SMALL>(part_ptr, numpart);
|
|
base.freeMemory<double>(param_ptr, 12);
|
|
|
|
|
|
std::cout << std::fixed << std::setprecision(4);
|
|
for (int i = 0; i < 10; i++) {
|
|
std::cout << parts[i].label << "\t"
|
|
<< parts[i].Rincol[0] << "\t" << parts[i].Rincol[1] << "\t"
|
|
<< parts[i].Rincol[2] << "\t" << parts[i].Pincol[0] << "\t"
|
|
<< parts[i].Pincol[1] << "\t" << parts[i].Pincol[2] << "\t"
|
|
<< std::endl;
|
|
}
|
|
|
|
std:: cout << "..." << std::endl;
|
|
|
|
for (int i = numpart - 10; i < numpart; i++) {
|
|
std::cout << parts[i].label << "\t"
|
|
<< parts[i].Rincol[0] << "\t" << parts[i].Rincol[1] << "\t"
|
|
<< parts[i].Rincol[2] << "\t" << parts[i].Pincol[0] << "\t"
|
|
<< parts[i].Pincol[1] << "\t" << parts[i].Pincol[2] << "\t"
|
|
<< std::endl;
|
|
}
|
|
|
|
double arx = 0, ary = 0, arz = 0;
|
|
double apx = 0, apy = 0, apz = 0;
|
|
for (int i = 0; i < numpart; i++) {
|
|
|
|
arx += sqrt(parts[i].Rincol[0] * parts[i].Rincol[0]) / numpart;
|
|
ary += sqrt(parts[i].Rincol[1] * parts[i].Rincol[1]) / numpart;
|
|
arz += sqrt(parts[i].Rincol[2] * parts[i].Rincol[2]) / numpart;
|
|
|
|
apx += sqrt(parts[i].Pincol[0] * parts[i].Pincol[0]) / numpart;
|
|
apy += sqrt(parts[i].Pincol[1] * parts[i].Pincol[1]) / numpart;
|
|
apz += sqrt(parts[i].Pincol[2] * parts[i].Pincol[2]) / numpart;
|
|
|
|
}
|
|
|
|
std::cout << std::fixed << std::setprecision(10);
|
|
std::cout << "R (" << arx << ", " << ary << ", " << arz << ") " << std::endl
|
|
<< "P (" << apx << ", " << apy << ", " << apz << ") " << std::endl;
|
|
|
|
|
|
cout << "==========================END TEST==========================" << endl;
|
|
return 0;
|
|
|
|
}
|