cosmotool/sample/simple3DFilter.cpp

204 lines
5.6 KiB
C++
Raw Normal View History

#include "omptl/algorithm"
#include <cassert>
#include "yorick.hpp"
#include "sphSmooth.hpp"
#include "mykdtree.hpp"
#include "miniargs.hpp"
#include <H5Cpp.h>
#include "hdf5_array.hpp"
#include <iostream>
using namespace std;
using namespace CosmoTool;
#define N_SPH 16
struct VCoord{
float v[3];
};
template<int i>
ComputePrecision getVelocity(const VCoord& v)
{
return v.v[i];
}
ComputePrecision getUnity(const VCoord& v)
{
return 1.0;
}
typedef SPHSmooth<VCoord> MySmooth;
typedef MySmooth::SPHTree MyTree;
typedef MyTree::Cell MyCell;
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;
int Nres;
MiniArgDesc args[] = {
{ "INPUT DATA1", &fname1, MINIARG_STRING },
{ "RADIUS LIMIT", &rLimit, MINIARG_DOUBLE },
{ "BOXSIZE", &boxsize, MINIARG_DOUBLE },
{ "RESOLUTION", &Nres, MINIARG_INT },
{ "CX", &cx, MINIARG_DOUBLE },
{ "CY", &cy, MINIARG_DOUBLE },
{ "CZ", &cz, MINIARG_DOUBLE },
{ 0, 0, MINIARG_NULL }
};
if (!parseMiniArgs(argc, argv, args))
return 1;
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;
array3_type bins(boost::extents[Nres][Nres][Nres]);
rLimit2 = rLimit*rLimit;
hdf5_read_array(in_f, "particles", v1_data);
assert(v1_data.shape()[1] == 6);
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.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][j];
for (int k = 0; k < 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;
long rx = floor((allCells_1[i].coord[0]+cx)*Nres/boxsize+0.5);
long ry = floor((allCells_1[i].coord[1]+cy)*Nres/boxsize+0.5);
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;
bins[rx][ry][rz]++;
}
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;
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;
#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)
{
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;
#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);
}
}
}
}
//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;
};