diff --git a/CMakeLists.txt b/CMakeLists.txt index e2e45f6..1830463 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,6 +31,8 @@ find_library(NETCDF_LIBRARY netcdf) find_library(GSL_LIBRARY gsl) find_library(GSLCBLAS_LIBRARY gslcblas) +find_package(Boost 1.53) + if (ENABLE_SHARP) SET(SHARP_SOURCE ${CMAKE_SOURCE_DIR}/external/sharp) SET(DEP_BUILD ${CMAKE_SOURCE_DIR}/external/sharp/auto) diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 2053060..606c8d1 100644 --- a/sample/CMakeLists.txt +++ b/sample/CMakeLists.txt @@ -1,5 +1,7 @@ SET(tolink ${GSL_LIBRARIES} CosmoTool ${CosmoTool_LIBS}) -include_directories(${CMAKE_SOURCE_DIR}/src ${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH}) +include_directories(${CMAKE_SOURCE_DIR}/src) +include_directories(${FFTW3_INCLUDE_DIRS} +${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH}) IF(SHARP_INCLUDE_PATH) include_directories(BEFORE ${SHARP_INCLUDE_PATH}) @@ -11,9 +13,6 @@ target_link_libraries(testBQueue ${tolink}) add_executable(testInterpolate testInterpolate.cpp) target_link_libraries(testInterpolate ${tolink}) -#add_executable(testSmooth testSmooth.cpp) -#target_link_libraries(testSmooth ${tolink}) - add_executable(testkd testkd.cpp) target_link_libraries(testkd ${tolink}) @@ -34,7 +33,7 @@ target_link_libraries(testPool ${tolink}) if (HDF5_FOUND) add_executable(testReadFlash testReadFlash.cpp) - target_link_libraries(testReadFlash ${tolink}) + target_link_libraries(testReadFlash ${tolink}) endif (HDF5_FOUND) @@ -69,3 +68,15 @@ endif (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) add_executable(test_cosmopower test_cosmopower.cpp) target_link_libraries(test_cosmopower ${tolink}) +if (Boost_FOUND) + include_directories(${Boost_INCLUDE_DIRS}) + + add_executable(testSmooth testSmooth.cpp) + target_link_libraries(testSmooth ${tolink}) + + add_executable(simple3DFilter simple3DFilter.cpp) + target_link_libraries(simple3DFilter ${tolink}) +endif (Boost_FOUND) + +add_executable(gadgetToDensity gadgetToDensity.cpp) +target_link_libraries(gadgetToDensity ${tolink}) diff --git a/sample/gadgetToDensity.cpp b/sample/gadgetToDensity.cpp new file mode 100644 index 0000000..d17f263 --- /dev/null +++ b/sample/gadgetToDensity.cpp @@ -0,0 +1,62 @@ +#include +#include +#include +#include "cic.hpp" +#include "loadGadget.hpp" +#include "miniargs.hpp" +#include "yorick.hpp" + +using namespace std; +using namespace CosmoTool; + +int main(int argc, char **argv) +{ + uint32_t res; + char *fname; + int id; + + MiniArgDesc desc[] = { + { "SNAPSHOT", &fname, MINIARG_STRING }, + { "ID", &id, MINIARG_INT }, + { "RESOLUTION", &res, MINIARG_INT }, + { 0, 0, MINIARG_NULL } + }; + + if (!parseMiniArgs(argc, argv, desc)) + return 1; + + SimuData *p = loadGadgetMulti(fname, 0, 0); + double L0 = p->BoxSize/1000; + CICFilter filter(res, L0); + + delete p; + + try { + for (int cpuid=0;;cpuid++) { + p = loadGadgetMulti(fname, cpuid, NEED_POSITION); + for (uint32_t i = 0; i < p->NumPart; i++) + { + CICParticles a; + + a.mass = 1.0; + a.coords[0] = p->Pos[0][i]/1000; + a.coords[1] = p->Pos[1][i]/1000; + a.coords[2] = p->Pos[2][i]/1000; + filter.putParticles(&a, 1); + } + + delete p; + } + } catch (const NoSuchFileException& e) {} + + CICType *denField; + uint32_t Ntot; + filter.getDensityField(denField, Ntot); + + cout << "L0=" << L0 << endl; + cout << "Saving density field" << endl; + uint32_t dimList[] = { res, res, res}; + saveArray("densityField.nc", denField, dimList, 3); + + return 0; +} diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp new file mode 100644 index 0000000..b7d1498 --- /dev/null +++ b/sample/simple3DFilter.cpp @@ -0,0 +1,193 @@ +#include +#include "yorick.hpp" +#include "sphSmooth.hpp" +#include "mykdtree.hpp" +#include "miniargs.hpp" + +using namespace std; +using namespace CosmoTool; + +#define N_SPH 16 + +struct VCoord{ + float v[3]; +}; + +template +ComputePrecision getVelocity(const VCoord& v) +{ + return v.v[i]; +} + +ComputePrecision getUnity(const VCoord& v) +{ + return 1.0; +} + +typedef SPHSmooth MySmooth; +typedef MySmooth::SPHTree MyTree; +typedef MyTree::Cell MyCell; + +int main(int argc, char **argv) +{ + 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; + + float *v1_data; + uint32_t *dimList; + uint32_t rank; + uint32_t N1_points, N2_points; + int *bins = new int[Nres*Nres*Nres]; + + rLimit2 = rLimit*rLimit; + + loadArray(fname1, v1_data, dimList, rank); + assert(rank == 2); + assert(dimList[1] == 6); + + N1_points = dimList[0]; + delete[] dimList; + + 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; + + 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]; + for (int k = 0; k < 3; k++) + allCells_1[i].val.pValue.v[k] = v1_data[i*6 + 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*Nres + rz*Nres*Nres]++; + } + delete[] v1_data; + uint32_t dims[3] = { Nres, Nres, Nres } ; + saveArray("num_in_cell.nc", bins, dims, 3); + + 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); + + 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 }; + + double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2]; + if (r2 > rLimit2) + { + continue; + } + + 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); + } + } + } + + + + + cout << "Interpolating..." << 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 }; + + 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; + } + + uint32_t numInCell = bins[rx + ry*Nres + rz*Nres*Nres]; + if (numInCell > N_SPH) + smooth1.fetchNeighbours(c, numInCell); + else + smooth1.fetchNeighbours(c); + + 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); + } + } + } + + return 0; +}; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8901027..611deaa 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -9,6 +9,7 @@ SET(CosmoTool_SRCS miniargs.cpp growthFactor.cpp cosmopower.cpp + cic.cpp ) IF(FOUND_NETCDF3) diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc index b9c9d2e..1ae5486 100644 --- a/src/sphSmooth.tcc +++ b/src/sphSmooth.tcc @@ -47,12 +47,12 @@ SPHSmooth::fetchNeighbours(const typename SPHTree::coords& c, uin if (requested > maxNgb) { maxNgb = requested; - internal.ngb = new SPHCell *[maxNgb]; - internal.distances = new CoordType[maxNgb]; + internal.ngb = boost::shared_ptr(new SPHCell *[maxNgb]); + internal.distances = boost::shared_ptr(new CoordType[maxNgb]); } memcpy(internal.currentCenter, c, sizeof(c)); - tree->getNearestNeighbours(c, requested, internal.ngb, internal.distances); + tree->getNearestNeighbours(c, requested, internal.ngb.get(), internal.distances.get()); internal.currentNgb = 0; for (uint32_t i = 0; i < requested && internal.ngb[i] != 0; i++,internal.currentNgb++)