aaron/ispc
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

added type definition for real

Browse Source
This commit is contained in:
Evghenii
2014-02-02 12:07:53 +01:00
parent 522257343b
commit 796942e6fa
2 changed files with 273 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

138
examples_ptx/nbody/plummer.h Normal file
View File

@@ -0,0 +1,138 @@
#pragma once
#include <vector>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include "vector3.h"
struct Plummer{
std::vector<double> mass;
std::vector<dvec3> pos, vel;
Plummer(
unsigned long n,
unsigned int seed = 19810614,
const char *filename = "plummer.dat")
: mass(n), pos(n), vel(n)
{
{
std::ifstream ifs(filename);
if(!ifs.fail()){
unsigned long ntmp, stmp;
ifs.read((char *)&ntmp, sizeof(unsigned long));
ifs.read((char *)&stmp, sizeof(unsigned long));
if(n == ntmp && seed == stmp){
ifs.read((char *)&mass[0], n*sizeof(double));
ifs.read((char *)& pos[0], n*sizeof(dvec3));
ifs.read((char *)& vel[0], n*sizeof(dvec3));
if(!ifs.fail()){
fprintf(stdout, "plummer : read from %s\n", filename);
}
return;
}
}
}
srand(seed);
unsigned long i = 0;
while(i < n){
double X1 = my_rand();
double X2 = my_rand();
double X3 = my_rand();
double R = 1.0/sqrt( (pow(X1,-2.0/3.0) - 1.0) );
if(R < 100.0) {
double Z = (1.0 - 2.0*X2)*R;
double X = sqrt(R*R - Z*Z) * cos(2.0*M_PI*X3);
double Y = sqrt(R*R - Z*Z) * sin(2.0*M_PI*X3);
double Ve = sqrt(2.0)*pow( (1.0 + R*R), -0.25 );
double X4 = 0.0;
double X5 = 0.0;
while( 0.1*X5 >= X4*X4*pow( (1.0-X4*X4), 3.5) ) {
X4 = my_rand(); X5 = my_rand();
}
double V = Ve*X4;
double X6 = my_rand();
double X7 = my_rand();
double Vz = (1.0 - 2.0*X6)*V;
double Vx = sqrt(V*V - Vz*Vz) * cos(2.0*M_PI*X7);
double Vy = sqrt(V*V - Vz*Vz) * sin(2.0*M_PI*X7);
double conv = 3.0*M_PI/16.0;
X *= conv; Y *= conv; Z *= conv;
Vx /= sqrt(conv); Vy /= sqrt(conv); Vz /= sqrt(conv);
double M = 1.0;
mass[i] = M/n;
pos[i][0] = X;
pos[i][1] = Y;
pos[i][2] = Z;
vel[i][0] = Vx;
vel[i][1] = Vy;
vel[i][2] = Vz;
/*
tmp_i = ldiv(i, 256);
if(tmp_i.rem == 0) printf("i = %d \n", i);
*/
ldiv_t tmp_i = ldiv(i, n/64);
if(tmp_i.rem == 0) {
printf(".");
fflush(stdout);
}
i++;
}
} // while (i<n)
double mcm = 0.0;
double xcm[3], vcm[3];
for(int k=0;k<3;k++) {
xcm[k] = 0.0; vcm[k] = 0.0;
} /* k */
for(i=0; i<n; i++) {
mcm += mass[i];
for(int k=0;k<3;k++) {
xcm[k] += mass[i] * pos[i][k];
vcm[k] += mass[i] * vel[i][k];
} /* k */
} /* i */
for(int k=0;k<3;k++) {
xcm[k] /= mcm; vcm[k] /= mcm;
} /* k */
for(i=0; i<n; i++) {
for(int k=0;k<3;k++) {
pos[i][k] -= xcm[k];
vel[i][k] -= vcm[k];
} /* k */
} /* i */
printf("\n");
{
std::ofstream ofs(filename);
if(!ofs.fail()){
unsigned long ntmp = n;
unsigned long stmp = seed;
ofs.write((char *)&ntmp, sizeof(unsigned long));
ofs.write((char *)&stmp, sizeof(unsigned long));
ofs.write((char *)&mass[0], n*sizeof(double));
ofs.write((char *)& pos[0], n*sizeof(dvec3));
ofs.write((char *)& vel[0], n*sizeof(dvec3));
if(!ofs.fail()){
fprintf(stdout, "plummer : wrote to %s\n", filename);
}
}
}
}
static double my_rand(void) {
return rand()/(1. + RAND_MAX);
}
};

135
examples_ptx/nbody/vector3.h Normal file
View File

@@ -0,0 +1,135 @@
#pragma once
#include <iostream>
#include <fstream>
#include <cmath>
template <class REAL> struct vector3{
public:
REAL x, y, z;
vector3(){
x = y = z = REAL(0);
}
vector3(const REAL &r){
x = y = z = r;
}
vector3(const REAL &_x, const REAL &_y, const REAL &_z){
x = _x; y = _y; z = _z;
}
vector3(const REAL *p){
x = p[0]; y = p[1]; z = p[2];
}
~vector3(){}
REAL &operator [](int i){
return (&x)[i];
}
const REAL &operator [](int i) const{
return (&x)[i];
}
template <class real>
operator vector3<real> () const {
return vector3<real> (real(x), real(y), real(z));
}
operator REAL *(){
return &x;
}
REAL (*toPointer())[3]{
return (REAL (*)[3])&x;
}
typedef REAL (*pArrayOfReal3)[3];
operator pArrayOfReal3(){
return toPointer();
}
void outv(std::ostream &ofs = std::cout) const{
ofs << "(" << x << ", " << y << ", " << z << ")" << std::endl;
}
bool are_numbers () const{
// returns false if *this has (a) NaN member(s)
return (norm2() >= REAL(0));
}
REAL norm2() const{
return (*this)*(*this);
}
REAL abs() const{
return std::sqrt(norm2());
}
friend std::ostream &operator << (std::ostream &ofs, const vector3<REAL> &v){
ofs << v.x << " " << v.y << " " << v.z;
return ofs;
}
friend std::istream &operator >> (std::istream &ifs, vector3<REAL> &v){
ifs >> v.x >> v.y >> v.z;
return ifs;
}
const vector3<REAL> operator + (const vector3<REAL> &v) const{
return vector3<REAL> (x+v.x, y+v.y, z+v.z);
}
const inline vector3<REAL> operator - (const vector3<REAL> &v) const{
return vector3<REAL> (x-v.x, y-v.y, z-v.z);
}
const vector3<REAL> operator * (const REAL &s) const{
return vector3<REAL> (x*s, y*s, z*s);
}
friend const vector3<REAL> operator * (const REAL &s, const vector3<REAL> &v){
return v*s;
}
// dot product
const inline REAL operator * (const vector3<REAL> &v) const{
return (x*v.x + y*v.y + z*v.z);
}
// vector product
const vector3<REAL> operator % (const vector3<REAL> &v) const{
return vector3<REAL> (
y*v.z - z*v.y,
z*v.x - x*v.z,
x*v.y - y*v.x);
}
const vector3<REAL> operator / (const REAL &s) const{
REAL r = REAL(1)/s;
return (*this)*r;
}
const vector3<REAL> operator = (const vector3<REAL> &v){
x = v.x; y=v.y; z=v.z;
return *this;
}
const vector3<REAL> operator - (){
return vector3<REAL> (-x, -y, -z);
}
const vector3<REAL> &operator += (const vector3<REAL> &v){
*this = *this + v;
return *this;
}
const vector3<REAL> &operator -= (const vector3<REAL> &v){
*this = *this - v;
return *this;
}
const vector3<REAL> &operator *= (const REAL &s){
*this = *this * s;
return *this;
}
const vector3<REAL> &operator /= (const REAL &s){
*this = *this / s;
return *this;
}
friend const vector3<REAL> maxeach (const vector3<REAL> &a, const vector3<REAL> &b){
return vector3<REAL> (std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z));
}
friend const vector3<REAL> mineach (const vector3<REAL> &a, const vector3<REAL> &b){
return vector3<REAL> (std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z));
}
const vector3<REAL> abseach(){
return vector3<REAL> (std::fabs(x), std::fabs(y), std::fabs(z));
}
};
typedef vector3<double> dvec3;
typedef vector3<float> fvec3;