diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp index b7d1498..1932e97 100644 --- a/sample/simple3DFilter.cpp +++ b/sample/simple3DFilter.cpp @@ -3,6 +3,9 @@ #include "sphSmooth.hpp" #include "mykdtree.hpp" #include "miniargs.hpp" +#include +#include "hdf5_array.hpp" +#include using namespace std; using namespace CosmoTool; @@ -30,6 +33,10 @@ typedef MyTree::Cell MyCell; 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; int Nres; @@ -48,35 +55,34 @@ int main(int argc, char **argv) if (!parseMiniArgs(argc, argv, args)) return 1; - float *v1_data; - uint32_t *dimList; - uint32_t rank; + H5::H5File in_f(fname1, 0); + H5::H5File out_f("fields.h5", H5F_ACC_TRUNC); + array_type v1_data; uint32_t N1_points, N2_points; - int *bins = new int[Nres*Nres*Nres]; + + array3_type bins(boost::extents[Nres][Nres][Nres]); rLimit2 = rLimit*rLimit; - loadArray(fname1, v1_data, dimList, rank); - assert(rank == 2); - assert(dimList[1] == 6); + hdf5_read_array(in_f, "particles", v1_data); + assert(v1_data.shape()[1] == 6); - N1_points = dimList[0]; - delete[] dimList; + N1_points = v1_data.shape()[0]; cout << "Got " << N1_points << " in the first file." << endl; MyCell *allCells_1 = new MyCell[N1_points]; for (long i = 0; i < Nres*Nres*Nres; i++) - bins[i] = 0; + bins.data()[i] = 0; cout << "Shuffling data in cells..." << endl; for (int i = 0 ; i < N1_points; i++) { 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++) - 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].val.weight = 0.0; @@ -85,109 +91,114 @@ int main(int argc, char **argv) long rz = floor((allCells_1[i].coord[2]+cz)*Nres/boxsize+0.5); 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; - uint32_t dims[3] = { Nres, Nres, Nres } ; - saveArray("num_in_cell.nc", bins, dims, 3); + v1_data.resize(boost::extents[1][1]); + + hdf5_write_array(out_f, "num_in_cell", bins); cout << "Building trees..." << endl; MyTree tree1(allCells_1, N1_points); cout << "Creating smoothing filter..." << endl; - MySmooth smooth1(&tree1, N_SPH); - uint32_t outDimList[3] = { Nres, Nres, Nres }; - uint32_t outDimList2[4] = { 3, Nres, Nres, Nres }; - ProgressiveOutput out_den_1 = - ProgressiveOutput::saveArrayProgressive("density.nc", outDimList, 3); - ProgressiveOutput out_vel_1 = - ProgressiveOutput::saveArrayProgressive("v3d.nc", outDimList2, 4); - ProgressiveOutput out_rad_1 = - ProgressiveOutput::saveArrayProgressive("rad.nc", outDimList, 3); + 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; - 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 }; +#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; - double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2]; - if (r2 > rLimit2) - { - continue; - } + 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 }; - uint32_t numInCell = bins[rx + Nres*ry + Nres*Nres*rz]; - if (numInCell > N_SPH) - smooth1.fetchNeighbours(c, numInCell); - else - smooth1.fetchNeighbours(c); - smooth1.addGridSite(c); - } - } - } + double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2]; + if (r2 > rLimit2) + { + continue; + } + uint32_t numInCell = bins[rx][ry][rz]; + if (numInCell > N_SPH) + smooth1.fetchNeighbours(c, numInCell); + else + smooth1.fetchNeighbours(c); +#pragma omp critical + smooth1.addGridSite(c); + } + } + } + } - - cout << "Interpolating..." << endl; - for (int rz = 0; rz < Nres; rz++) - { - double pz = (rz)*boxsize/Nres-cz; +#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 }; + 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_vel_1.put(0); - out_vel_1.put(0); - out_vel_1.put(0); - out_den_1.put(0); - out_rad_1.put(0); - continue; - } + 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*Nres + rz*Nres*Nres]; - if (numInCell > N_SPH) - smooth1.fetchNeighbours(c, numInCell); - else - smooth1.fetchNeighbours(c); + uint32_t numInCell = bins[rx][ry][rz]; + if (numInCell > N_SPH) + smooth1.fetchNeighbours(c, numInCell); + else + smooth1.fetchNeighbours(c); - float val; + float val; - out_rad_1.put(smooth1.getSmoothingLen()); - val = smooth1.computeSmoothedValue(c, getVelocity<0>); - out_vel_1.put(val); - val = smooth1.computeSmoothedValue(c, getVelocity<1>); - out_vel_1.put(val); - val = smooth1.computeSmoothedValue(c, getVelocity<2>); - out_vel_1.put(val); - val = smooth1.computeSmoothedValue(c, getUnity); - out_den_1.put(val); - } - } - } + 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); + } + } + } + } + + 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; };