guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of bitbucket.org:glavaux/cosmotool

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-06-18 09:12:01 +02:00
commit ff01447ab6
4 changed files with 89 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
sample/gadgetToArray.cpp
View File

@ -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];

118
sample/simple3DFilter.cpp
View File

@ -7,6 +7,8 @@
#include <H5Cpp.h>
#include "hdf5_array.hpp"
#include <iostream>
#include <boost/format.hpp>
#include <boost/bind.hpp>
using namespace std;
using namespace CosmoTool;
@ -17,8 +19,13 @@ struct VCoord{
float v[3];
};
template<int i>
ComputePrecision getVelocity(const VCoord& v)
using boost::format;
using boost::str;
typedef boost::multi_array<float, 2> array_type;
typedef boost::multi_array<float, 3> array3_type;
typedef boost::multi_array<float, 4> array4_type;
ComputePrecision getVelocity(const VCoord& v, int i)
{
return v.v[i];
}
@ -32,11 +39,51 @@ typedef SPHSmooth<VCoord> MySmooth;
typedef MySmooth::SPHTree MyTree;
typedef MyTree::Cell MyCell;
template<typename FuncT>
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<float, 2> array_type;
typedef boost::multi_array<float, 3> array3_type;
typedef boost::multi_array<float, 4> 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;
};

18
src/sphSmooth.hpp
View File

@ -82,10 +82,14 @@ namespace CosmoTool
return internal.currentCenter;
}
template<typename FuncT>
ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
computeParticleValue fun, SPHState *state = 0);
FuncT fun, SPHState *state = 0);
template<typename FuncT>
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<typename FuncT>
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<typename FuncT>
ComputePrecision computeWValue(const typename SPHTree::coords & c,
SPHCell& cell,
CoordType d,
computeParticleValue fun, SPHState *state);
FuncT fun, SPHState *state);
template<typename FuncT>
void runUnrollNode(SPHNode *node,
runParticleValue fun);
FuncT fun);
};
template<class ValType1, class ValType2, int Ndims>

12
src/sphSmooth.tcc
View File

@ -21,10 +21,11 @@ SPHSmooth<ValType,Ndims>::~SPHSmooth()
}
template<typename ValType, int Ndims>
template<typename FuncT>
ComputePrecision SPHSmooth<ValType,Ndims>::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<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords
}
template<typename ValType, int Ndims>
template<typename FuncT>
ComputePrecision
SPHSmooth<ValType,Ndims>::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<typename ValType, int Ndims>
template<typename FuncT>
ComputePrecision SPHSmooth<ValType,Ndims>::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<ValType,Ndims>::computeInterpolatedValue(const typena
}
template<typename ValType, int Ndims>
void SPHSmooth<ValType,Ndims>::runForEachNeighbour(runParticleValue fun, SPHState *state)
template<typename FuncT>
void SPHSmooth<ValType,Ndims>::runForEachNeighbour(FuncT fun, SPHState *state)
{
if (state == 0)
state = &internal;