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

Parallelized and HDF5ized simple3DFilter

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-05-20 17:16:54 +02:00
parent c88599aa59
commit 113a14e00c
1 changed files with 105 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
sample/simple3DFilter.cpp
View File

@ -3,6 +3,9 @@
#include "sphSmooth.hpp" #include "sphSmooth.hpp"
#include "mykdtree.hpp" #include "mykdtree.hpp"
#include "miniargs.hpp" #include "miniargs.hpp"
#include <H5Cpp.h>
#include "hdf5_array.hpp"
#include <iostream>
using namespace std; using namespace std;
using namespace CosmoTool; using namespace CosmoTool;
@ -30,6 +33,10 @@ typedef MyTree::Cell MyCell;
int main(int argc, char **argv) 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; char *fname1, *fname2;
double rLimit, boxsize, rLimit2, cx, cy, cz; double rLimit, boxsize, rLimit2, cx, cy, cz;
int Nres; int Nres;
@ -48,35 +55,34 @@ int main(int argc, char **argv)
if (!parseMiniArgs(argc, argv, args)) if (!parseMiniArgs(argc, argv, args))
return 1; return 1;
float *v1_data; H5::H5File in_f(fname1, 0);
uint32_t *dimList; H5::H5File out_f("fields.h5", H5F_ACC_TRUNC);
uint32_t rank; array_type v1_data;
uint32_t N1_points, N2_points; uint32_t N1_points, N2_points;
int *bins = new int[Nres*Nres*Nres];
array3_type bins(boost::extents[Nres][Nres][Nres]);
rLimit2 = rLimit*rLimit; rLimit2 = rLimit*rLimit;
loadArray(fname1, v1_data, dimList, rank); hdf5_read_array(in_f, "particles", v1_data);
assert(rank == 2); assert(v1_data.shape()[1] == 6);
assert(dimList[1] == 6);
N1_points = dimList[0]; N1_points = v1_data.shape()[0];
delete[] dimList;
cout << "Got " << N1_points << " in the first file." << endl; cout << "Got " << N1_points << " in the first file." << endl;
MyCell *allCells_1 = new MyCell[N1_points]; MyCell *allCells_1 = new MyCell[N1_points];
for (long i = 0; i < Nres*Nres*Nres; i++) for (long i = 0; i < Nres*Nres*Nres; i++)
bins[i] = 0; bins.data()[i] = 0;
cout << "Shuffling data in cells..." << endl; cout << "Shuffling data in cells..." << endl;
for (int i = 0 ; i < N1_points; i++) for (int i = 0 ; i < N1_points; i++)
{ {
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
allCells_1[i].coord[j] = v1_data[i*6 + j]; allCells_1[i].coord[j] = v1_data[i][j];
for (int k = 0; k < 3; k++) for (int k = 0; k < 3; k++)
allCells_1[i].val.pValue.v[k] = v1_data[i*6 + 3 + k]; allCells_1[i].val.pValue.v[k] = v1_data[i][3+k];
allCells_1[i].active = true; allCells_1[i].active = true;
allCells_1[i].val.weight = 0.0; allCells_1[i].val.weight = 0.0;
@ -87,28 +93,28 @@ int main(int argc, char **argv)
if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres) if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
continue; continue;
bins[rx + ry*Nres + rz*Nres*Nres]++; bins[rx][ry][rz]++;
} }
delete[] v1_data; v1_data.resize(boost::extents[1][1]);
uint32_t dims[3] = { Nres, Nres, Nres } ;
saveArray("num_in_cell.nc", bins, dims, 3); hdf5_write_array(out_f, "num_in_cell", bins);
cout << "Building trees..." << endl; cout << "Building trees..." << endl;
MyTree tree1(allCells_1, N1_points); MyTree tree1(allCells_1, N1_points);
cout << "Creating smoothing filter..." << endl; cout << "Creating smoothing filter..." << endl;
MySmooth smooth1(&tree1, N_SPH);
uint32_t outDimList[3] = { Nres, Nres, Nres }; array3_type out_den_1(boost::extents[Nres][Nres][Nres]);
uint32_t outDimList2[4] = { 3, Nres, Nres, Nres }; array4_type out_v3d_1(boost::extents[Nres][Nres][Nres][3]);
ProgressiveOutput<float> out_den_1 = array3_type out_rad_1(boost::extents[Nres][Nres][Nres]);
ProgressiveOutput<float>::saveArrayProgressive("density.nc", outDimList, 3);
ProgressiveOutput<float> out_vel_1 =
ProgressiveOutput<float>::saveArrayProgressive("v3d.nc", outDimList2, 4);
ProgressiveOutput<float> out_rad_1 =
ProgressiveOutput<float>::saveArrayProgressive("rad.nc", outDimList, 3);
cout << "Weighing..." << endl; cout << "Weighing..." << endl;
#pragma omp parallel
{
MySmooth smooth1(&tree1, N_SPH);
#pragma omp for schedule(dynamic)
for (int rz = 0; rz < Nres; rz++) for (int rz = 0; rz < Nres; rz++)
{ {
double pz = (rz)*boxsize/Nres-cz; double pz = (rz)*boxsize/Nres-cz;
@ -129,20 +135,24 @@ int main(int argc, char **argv)
continue; continue;
} }
uint32_t numInCell = bins[rx + Nres*ry + Nres*Nres*rz]; uint32_t numInCell = bins[rx][ry][rz];
if (numInCell > N_SPH) if (numInCell > N_SPH)
smooth1.fetchNeighbours(c, numInCell); smooth1.fetchNeighbours(c, numInCell);
else else
smooth1.fetchNeighbours(c); smooth1.fetchNeighbours(c);
#pragma omp critical
smooth1.addGridSite(c); smooth1.addGridSite(c);
} }
} }
} }
}
cout << "Interpolating..." << endl; cout << "Interpolating..." << endl;
#pragma omp parallel
{
MySmooth smooth1(&tree1, N_SPH);
#pragma omp for schedule(dynamic)
for (int rz = 0; rz < Nres; rz++) for (int rz = 0; rz < Nres; rz++)
{ {
double pz = (rz)*boxsize/Nres-cz; double pz = (rz)*boxsize/Nres-cz;
@ -160,15 +170,15 @@ int main(int argc, char **argv)
double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2]; double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
if (r2 > rLimit2) if (r2 > rLimit2)
{ {
out_vel_1.put(0); out_v3d_1[rx][ry][rz][0] = 0;
out_vel_1.put(0); out_v3d_1[rx][ry][rz][1] = 0;
out_vel_1.put(0); out_v3d_1[rx][ry][rz][2] = 0;
out_den_1.put(0); out_den_1[rx][ry][rz] = 0;
out_rad_1.put(0); out_rad_1[rx][ry][rz] = -1;
continue; continue;
} }
uint32_t numInCell = bins[rx + ry*Nres + rz*Nres*Nres]; uint32_t numInCell = bins[rx][ry][rz];
if (numInCell > N_SPH) if (numInCell > N_SPH)
smooth1.fetchNeighbours(c, numInCell); smooth1.fetchNeighbours(c, numInCell);
else else
@ -176,18 +186,19 @@ int main(int argc, char **argv)
float val; float val;
out_rad_1.put(smooth1.getSmoothingLen()); out_rad_1[rx][ry][rz] = smooth1.getSmoothingLen();
val = smooth1.computeSmoothedValue(c, getVelocity<0>); out_v3d_1[rx][ry][rz][0] = smooth1.computeSmoothedValue(c, getVelocity<0>);
out_vel_1.put(val); out_v3d_1[rx][ry][rz][1] = smooth1.computeSmoothedValue(c, getVelocity<1>);
val = smooth1.computeSmoothedValue(c, getVelocity<1>); out_v3d_1[rx][ry][rz][2] = smooth1.computeSmoothedValue(c, getVelocity<2>);
out_vel_1.put(val); out_den_1[rx][ry][rz] = smooth1.computeSmoothedValue(c, getUnity);
val = smooth1.computeSmoothedValue(c, getVelocity<2>); }
out_vel_1.put(val);
val = smooth1.computeSmoothedValue(c, getUnity);
out_den_1.put(val);
} }
} }
} }
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; return 0;
}; };