From 34eac2721a29bdb7ad802c3bb0377bfee57e3563 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Fri, 13 Jun 2014 18:28:05 +0200 Subject: [PATCH 1/3] Small Fixes / enhancement to project_tool --- python/project_tool.hpp | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/python/project_tool.hpp b/python/project_tool.hpp index 6f8c1d7..c185c44 100644 --- a/python/project_tool.hpp +++ b/python/project_tool.hpp @@ -5,14 +5,14 @@ template static T project_tool(T *vertex_value, T *u, T *u0) { T ret0 = 0; - for (int i = 0; i < 8; i++) + for (unsigned int i = 0; i < 8; i++) { - int c[3] = { i & 1, (i>>1)&1, (i>>2)&1 }; + unsigned int c[3] = { i & 1, (i>>1)&1, (i>>2)&1 }; int epsilon[3]; T ret = 0; for (int q = 0; q < 3; q++) - epsilon[q] = (2*c[q]-1); + epsilon[q] = 2*c[q]-1; for (int q = 0; q < ProdType::numProducts; q++) ret += ProdType::product(u, u0, epsilon, q); @@ -27,7 +27,8 @@ static T project_tool(T *vertex_value, T *u, T *u0) template static inline T get_u0(const T& u0, int epsilon) { - return (epsilon > 0) ? u0 : (1-u0); + return (1-epsilon)/2 + epsilon*u0; +// return (epsilon > 0) ? u0 : (1-u0); } template @@ -39,7 +40,7 @@ struct ProductTerm0 { T a = 1; - for (int r = 0; r < 3; r++) + for (unsigned int r = 0; r < 3; r++) a *= get_u0(u0[r], epsilon[r]); return a; } @@ -54,14 +55,14 @@ struct ProductTerm1 static T product(T *u, T *u0, int *epsilon, int q) { T a = 1; - double G[3]; + T G[3]; - for (int r = 0; r < 3; r++) + for (unsigned int r = 0; r < 3; r++) { G[r] = get_u0(u0[r], epsilon[r]); } - double F[3] = { G[1]*G[2], G[0]*G[2], G[0]*G[1] }; + T F[3] = { G[1]*G[2], G[0]*G[2], G[0]*G[1] }; return F[q] * u[q] * epsilon[q]; } @@ -75,14 +76,14 @@ struct ProductTerm2 static inline T product(T *u, T *u0, int *epsilon, int q) { T a = 1; - double G[3]; + T G[3]; - for (int r = 0; r < 3; r++) + for (unsigned int r = 0; r < 3; r++) { G[r] = get_u0(u0[r], epsilon[r]); } - double F[3] = { epsilon[1]*epsilon[2]*u[1]*u[2], epsilon[0]*epsilon[2]*u[0]*u[2], epsilon[0]*epsilon[1]*u[0]*u[1] }; + T F[3] = { epsilon[1]*epsilon[2]*u[1]*u[2], epsilon[0]*epsilon[2]*u[0]*u[2], epsilon[0]*epsilon[1]*u[0]*u[1] }; return F[q] * G[q]; } From d23f79cbaf7f3e027f7cabbf11a896b881ec7e82 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 14 Jun 2014 17:18:28 +0200 Subject: [PATCH 2/3] Templatized SPHSmooth. Reduced memory consumption of simple3DFilter --- sample/simple3DFilter.cpp | 118 ++++++++++++++++++++------------------ src/sphSmooth.hpp | 18 ++++-- src/sphSmooth.tcc | 12 ++-- 3 files changed, 84 insertions(+), 64 deletions(-) diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp index b3a74d2..345834c 100644 --- a/sample/simple3DFilter.cpp +++ b/sample/simple3DFilter.cpp @@ -7,6 +7,8 @@ #include #include "hdf5_array.hpp" #include +#include +#include using namespace std; using namespace CosmoTool; @@ -17,8 +19,13 @@ struct VCoord{ float v[3]; }; -template -ComputePrecision getVelocity(const VCoord& v) +using boost::format; +using boost::str; +typedef boost::multi_array array_type; +typedef boost::multi_array array3_type; +typedef boost::multi_array array4_type; + +ComputePrecision getVelocity(const VCoord& v, int i) { return v.v[i]; } @@ -32,11 +39,51 @@ typedef SPHSmooth MySmooth; typedef MySmooth::SPHTree MyTree; typedef MyTree::Cell MyCell; +template +void computeInterpolatedField(MyTree *tree1, double boxsize, int Nres, double cx, double cy, double cz, + array3_type& bins, array3_type& arr, FuncT func, double rLimit2) +{ +#pragma omp parallel + { + MySmooth smooth1(tree1, N_SPH); + +#pragma omp for schedule(dynamic) + for (int rz = 0; rz < Nres; rz++) + { + double pz = (rz)*boxsize/Nres-cz; + + cout << format("[%d] %d / %d") % omp_get_thread_num() % rz % Nres << endl; + for (int ry = 0; ry < Nres; ry++) + { + double py = (ry)*boxsize/Nres-cy; + for (int rx = 0; rx < Nres; rx++) + { + double px = (rx)*boxsize/Nres-cx; + + MyTree::coords c = { px, py, pz }; + + double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2]; + if (r2 > rLimit2) + { + arr[rx][ry][rz] = 0; + continue; + } + + uint32_t numInCell = bins[rx][ry][rz]; + if (numInCell > N_SPH) + smooth1.fetchNeighbours(c, numInCell); + else + smooth1.fetchNeighbours(c); + + arr[rx][ry][rz] = smooth1.computeSmoothedValue(c, func); + } + } + } + } +} + int main(int argc, char **argv) { - typedef boost::multi_array array_type; - typedef boost::multi_array array3_type; - typedef boost::multi_array array4_type; char *fname1, *fname2; double rLimit, boxsize, rLimit2, cx, cy, cz; @@ -105,8 +152,6 @@ int main(int argc, char **argv) cout << "Creating smoothing filter..." << endl; - array3_type out_den_1(boost::extents[Nres][Nres][Nres]); - array4_type out_v3d_1(boost::extents[Nres][Nres][Nres][3]); // array3_type out_rad_1(boost::extents[Nres][Nres][Nres]); cout << "Weighing..." << endl; @@ -149,55 +194,18 @@ int main(int argc, char **argv) } cout << "Interpolating..." << endl; -#pragma omp parallel - { - MySmooth smooth1(&tree1, N_SPH); - -#pragma omp for schedule(dynamic) - for (int rz = 0; rz < Nres; rz++) - { - double pz = (rz)*boxsize/Nres-cz; - cout << rz << " / " << Nres << endl; - for (int ry = 0; ry < Nres; ry++) - { - double py = (ry)*boxsize/Nres-cy; - for (int rx = 0; rx < Nres; rx++) - { - double px = (rx)*boxsize/Nres-cx; - - MyTree::coords c = { px, py, pz }; - - double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2]; - if (r2 > rLimit2) - { - out_v3d_1[rx][ry][rz][0] = 0; - out_v3d_1[rx][ry][rz][1] = 0; - out_v3d_1[rx][ry][rz][2] = 0; - out_den_1[rx][ry][rz] = 0; - //out_rad_1[rx][ry][rz] = -1; - continue; - } - - uint32_t numInCell = bins[rx][ry][rz]; - if (numInCell > N_SPH) - smooth1.fetchNeighbours(c, numInCell); - else - smooth1.fetchNeighbours(c); - - //out_rad_1[rx][ry][rz] = smooth1.getSmoothingLen(); - out_v3d_1[rx][ry][rz][0] = smooth1.computeSmoothedValue(c, getVelocity<0>); - out_v3d_1[rx][ry][rz][1] = smooth1.computeSmoothedValue(c, getVelocity<1>); - out_v3d_1[rx][ry][rz][2] = smooth1.computeSmoothedValue(c, getVelocity<2>); - out_den_1[rx][ry][rz] = smooth1.computeSmoothedValue(c, getUnity); - } - } - } + array3_type interpolated(boost::extents[Nres][Nres][Nres]); + + computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz, + bins, interpolated, getUnity, rLimit2); + hdf5_write_array(out_f, "density", interpolated); + //out_f.flush(); + for (int i = 0; i < 3; i++) { + computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz, + bins, interpolated, boost::bind(getVelocity, _1, i), rLimit2); + hdf5_write_array(out_f, str(format("p%d") % i), interpolated); } - - //hdf5_write_array(out_f, "radii", out_rad_1); - hdf5_write_array(out_f, "velocity", out_v3d_1); - hdf5_write_array(out_f, "density", out_den_1); - + return 0; }; diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp index d116dd4..bcd98a5 100644 --- a/src/sphSmooth.hpp +++ b/src/sphSmooth.hpp @@ -82,10 +82,14 @@ namespace CosmoTool return internal.currentCenter; } + template ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c, - computeParticleValue fun, SPHState *state = 0); + FuncT fun, SPHState *state = 0); + + template ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c, - computeParticleValue fun, SPHState *state = 0); + FuncT fun, SPHState *state = 0); + ComputePrecision getMaxDistance(const typename SPHTree::coords& c, SPHNode *node) const; @@ -101,7 +105,8 @@ namespace CosmoTool } // END - void runForEachNeighbour(runParticleValue fun, SPHState *state = 0); + template + void runForEachNeighbour(FuncT fun, SPHState *state = 0); void addGridSite(const typename SPHTree::coords& c); bool hasNeighbours() const; @@ -127,12 +132,15 @@ namespace CosmoTool uint32_t maxNgb; SPHTree *tree; + template ComputePrecision computeWValue(const typename SPHTree::coords & c, SPHCell& cell, CoordType d, - computeParticleValue fun, SPHState *state); + FuncT fun, SPHState *state); + + template void runUnrollNode(SPHNode *node, - runParticleValue fun); + FuncT fun); }; template diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc index 1ae5486..99094cd 100644 --- a/src/sphSmooth.tcc +++ b/src/sphSmooth.tcc @@ -21,10 +21,11 @@ SPHSmooth::~SPHSmooth() } template +template ComputePrecision SPHSmooth::computeWValue(const typename SPHTree::coords& c, SPHCell& cell, CoordType d, - computeParticleValue fun, SPHState *state) + FuncT fun, SPHState *state) { CoordType weight; @@ -117,9 +118,10 @@ SPHSmooth::fetchNeighboursOnVolume(const typename SPHTree::coords } template +template ComputePrecision SPHSmooth::computeSmoothedValue(const typename SPHTree::coords& c, - computeParticleValue fun, SPHState *state) + FuncT fun, SPHState *state) { if (state == 0) state = &internal; @@ -144,8 +146,9 @@ ComputePrecision interpolateOne(const ValType& t) // WARNING ! Cell's weight must be 1 !!! template +template ComputePrecision SPHSmooth::computeInterpolatedValue(const typename SPHTree::coords& c, - computeParticleValue fun, SPHState *state) + FuncT fun, SPHState *state) { if (state == 0) state = &internal; @@ -164,7 +167,8 @@ ComputePrecision SPHSmooth::computeInterpolatedValue(const typena } template -void SPHSmooth::runForEachNeighbour(runParticleValue fun, SPHState *state) +template +void SPHSmooth::runForEachNeighbour(FuncT fun, SPHState *state) { if (state == 0) state = &internal; From 19b23fbe4fd9bc2d272559d72bd09884809cc943 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 14 Jun 2014 17:26:49 +0200 Subject: [PATCH 3/3] Enforce periodicity --- sample/gadgetToArray.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/sample/gadgetToArray.cpp b/sample/gadgetToArray.cpp index 6f646c6..386fa01 100644 --- a/sample/gadgetToArray.cpp +++ b/sample/gadgetToArray.cpp @@ -75,6 +75,11 @@ int main(int argc, char **argv) parts[q][0] = p->Pos[0][i]/1000; parts[q][1] = p->Pos[1][i]/1000; parts[q][2] = p->Pos[2][i]/1000; + for (int j = 0; j < 3; j++) + { + while (parts[q][j] < 0) parts[q][j] += L0; + while (parts[q][j] >= L0) parts[q][j] -= L0; + } parts[q][3] = p->Vel[0][i]; parts[q][4] = p->Vel[1][i]; parts[q][5] = p->Vel[2][i];