diff --git a/CMakeLists.txt b/CMakeLists.txt index 97575e1..880f616 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,13 +31,15 @@ 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) ExternalProject_Add(sharp SOURCE_DIR ${SHARP_SOURCE} BUILD_IN_SOURCE 1 - CONFIGURE_COMMAND ${SHARP_SOURCE}/configure --prefix=${DEP_BUILD} + CONFIGURE_COMMAND ${SHARP_SOURCE}/configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" --prefix=${DEP_BUILD} BUILD_COMMAND ${CMAKE_MAKE_PROGRAM} INSTALL_COMMAND echo "No install" ) @@ -67,12 +69,28 @@ pkg_check_modules(FFTW3 fftw3>=3.3) pkg_check_modules(FFTW3F fftw3f>=3.3) pkg_check_modules(EIGEN3 eigen3) + +find_library(FFTW3F_LIBRARY_FULL fftw3f PATHS ${FFTW3F_LIBDIR} NO_DEFAULT_PATH) +find_library(FFTW3_LIBRARY_FULL fftw3 PATHS ${FFTW3_LIBDIR} NO_DEFAULT_PATH) + include(FindPackageHandleStandardArgs) set(NETCDF_FIND_REQUIRED TRUE) set(GSL_FIND_REQUIRED TRUE) FIND_PACKAGE_HANDLE_STANDARD_ARGS(NetCDF DEFAULT_MSG NETCDF_LIBRARY NETCDFCPP_LIBRARY NETCDF_INCLUDE_PATH) FIND_PACKAGE_HANDLE_STANDARD_ARGS(GSL DEFAULT_MSG GSL_LIBRARY GSLCBLAS_LIBRARY GSL_INCLUDE_PATH) +IF(ENABLE_OPENMP) + + IF (NOT OPENMP_FOUND) + MESSAGE(ERROR "No known compiler option for enabling OpenMP") + ENDIF(NOT OPENMP_FOUND) + + SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKED_FLAGS} ${OpenMP_C_FLAGS}") +ENDIF(ENABLE_OPENMP) + + set(GSL_LIBRARIES ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY}) diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 969480c..47a74b5 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,14 @@ target_link_libraries(testPool ${tolink}) if (HDF5_FOUND) add_executable(testReadFlash testReadFlash.cpp) - target_link_libraries(testReadFlash ${tolink}) + target_link_libraries(testReadFlash ${tolink}) + + add_executable(testHDF5 testHDF5.cpp) + target_link_libraries(testHDF5 ${tolink}) + + add_executable(gadgetToArray gadgetToArray.cpp) + target_link_libraries(gadgetToArray ${tolink}) + endif (HDF5_FOUND) @@ -48,18 +54,36 @@ add_executable(testBSP testBSP.cpp) target_link_libraries(testBSP ${tolink}) if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) + IF (FFTW3F_LIBRARY_FULL) + SET(FFTW3_LIB ${FFTW3F_LIBRARY_FULL}) + ENDIF (FFTW3F_LIBRARY_FULL) + IF (FFTW3_LIBRARY_FULL) + SET(FFTW3_LIB ${FFTW3_LIB} ${FFTW3_LIBRARY_FULL}) + ENDIF (FFTW3_LIBRARY_FULL) add_executable(test_fft_calls test_fft_calls.cpp) - target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES}) + target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIB}) endif (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) if (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) include_directories(${SHARP_INCLUDE_PATH}) add_executable(test_healpix_calls test_healpix_calls.cpp) target_link_libraries(test_healpix_calls ${tolink} ${SHARP_LIBRARIES}) - set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) + set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS "${OpenMP_CXX_FLAGS}" LINK_FLAGS "${OpenMP_CXX_FLAGS}") add_dependencies(test_healpix_calls sharp) 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/gadgetToArray.cpp b/sample/gadgetToArray.cpp new file mode 100644 index 0000000..7d4669c --- /dev/null +++ b/sample/gadgetToArray.cpp @@ -0,0 +1,57 @@ +#include +#include +#include +#include "cic.hpp" +#include "loadGadget.hpp" +#include "miniargs.hpp" +#include +#include "hdf5_array.hpp" + +using namespace CosmoTool; +using namespace std; + +int main(int argc, char **argv) +{ + typedef boost::multi_array array_type; + uint32_t res; + char *fname; + int id; + + MiniArgDesc desc[] = { + { "SNAPSHOT", &fname, MINIARG_STRING }, + { 0, 0, MINIARG_NULL } + }; + + if (!parseMiniArgs(argc, argv, desc)) + return 1; + + H5::H5File f("density.h5", H5F_ACC_TRUNC); + + + SimuData *p = loadGadgetMulti(fname, 0, 0); + double L0 = p->BoxSize/1000; + array_type parts(boost::extents[p->TotalNumPart][6]); + uint64_t q = 0; + + try { + for (int cpuid=0;;cpuid++) { + p = loadGadgetMulti(fname, cpuid, NEED_POSITION|NEED_VELOCITY); + for (uint32_t i = 0; i < p->NumPart; i++) + { + parts[q][0] = p->Pos[0][i]/1000; + parts[q][1] = p->Pos[1][i]/1000; + parts[q][2] = p->Pos[2][i]/1000; + parts[q][3] = p->Vel[0][i]; + parts[q][4] = p->Vel[1][i]; + parts[q][5] = p->Vel[2][i]; + q++; + } + + delete p; + } + } catch (const NoSuchFileException& e) {} + + hdf5_write_array(f, "particles", parts); + + return 0; +} 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..1932e97 --- /dev/null +++ b/sample/simple3DFilter.cpp @@ -0,0 +1,204 @@ +#include +#include "yorick.hpp" +#include "sphSmooth.hpp" +#include "mykdtree.hpp" +#include "miniargs.hpp" +#include +#include "hdf5_array.hpp" +#include + +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) +{ + 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; + + 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); + + float 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; +}; diff --git a/sample/testHDF5.cpp b/sample/testHDF5.cpp new file mode 100644 index 0000000..82c3230 --- /dev/null +++ b/sample/testHDF5.cpp @@ -0,0 +1,65 @@ +#include +#include +#include "hdf5_array.hpp" + +using namespace std; + +int main() +{ + typedef boost::multi_array array_type; + typedef boost::multi_array array3_type; + typedef boost::multi_array, 2> arrayc_type; + typedef array_type::index index; + + H5::H5File f("test.h5", H5F_ACC_TRUNC); + + H5::Group g = f.createGroup("test_group"); + + array_type A(boost::extents[2][3]); + array_type B; + array3_type C(boost::extents[2][3][4]); + arrayc_type D, E; + + int values = 0; + for (index i = 0; i != 2; i++) + for (index j = 0; j != 3; j++) + A[i][j] = values++; + + CosmoTool::hdf5_write_array(g, "test_data", A); + + CosmoTool::hdf5_read_array(g, "test_data", B); + + int verify = 0; + for (index i = 0; i != 2; i++) + for (index j = 0; j != 3; j++) + if (B[i][j] != verify++) { + std::cout << "Invalid array content" << endl; + abort(); + } + + try + { + CosmoTool::hdf5_read_array(g, "test_data", C); + std::cout << "Did not throw InvalidDimensions" << endl; + abort(); + } + catch (const CosmoTool::InvalidDimensions&) + {} + + D.resize(boost::extents[2][3]); + D = A; + + CosmoTool::hdf5_write_array(g, "test_data_c", D); + + CosmoTool::hdf5_read_array(g, "test_data_c", E); + + verify = 0; + for (index i = 0; i != 2; i++) + for (index j = 0; j != 3; j++) + if (E[i][j].real() != verify++) { + std::cout << "Invalid array content" << endl; + abort(); + } + + return 0; +} diff --git a/sample/testSmooth.cpp b/sample/testSmooth.cpp index 4f3011d..ddb8c3e 100644 --- a/sample/testSmooth.cpp +++ b/sample/testSmooth.cpp @@ -85,14 +85,16 @@ int main() uint32_t dims[] = { NX, NX }; ProgressiveOutput out = ProgressiveOutput::saveArrayProgressive("out.nc", dims, 2); +//#pragma omp parallel for schedule(static) for (uint32_t iy = 0; iy < NX; iy++) { + MySmooth::SPHState state; cout << "iy=" << iy << endl; for (uint32_t ix = 0; ix < NX; ix++) { MyTree::coords c = { 1.0*ix/NX, 1.0*iy/NX }; - smooth.fetchNeighbours(c); - out.put(smooth.computeSmoothedValue(c, unit_fun)); + smooth.fetchNeighbours(c, &state); + out.put(smooth.computeSmoothedValue(c, unit_fun, &state)); } diff --git a/sample/test_fft_calls.cpp b/sample/test_fft_calls.cpp index b98fcea..e49ad82 100644 --- a/sample/test_fft_calls.cpp +++ b/sample/test_fft_calls.cpp @@ -89,5 +89,6 @@ int main(int argc, char **argv) test_2d(128,128); test_2d(128,131); test_2d(128,130); + test_2d(131,128); 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/cosmopower.cpp b/src/cosmopower.cpp index a96a6bb..20026a2 100644 --- a/src/cosmopower.cpp +++ b/src/cosmopower.cpp @@ -212,7 +212,7 @@ double CosmoPower::powerBDM(double k) double CosmoPower::powerTest(double k) { - return 1/(1+k*k); + return normPower;//1/(1+k*k); } /* diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp index ff092ef..42b579a 100644 --- a/src/fourier/base_types.hpp +++ b/src/fourier/base_types.hpp @@ -55,6 +55,7 @@ namespace CosmoTool protected: SpectrumFunction() {} public: + typedef T type; typedef Eigen::Array VecType; typedef Eigen::Map MapType; typedef Eigen::Map ConstMapType; @@ -88,6 +89,7 @@ namespace CosmoTool FourierMap() {} public: + typedef T type; typedef Eigen::Array VecType; typedef Eigen::Map MapType; typedef Eigen::Map ConstMapType; @@ -190,6 +192,7 @@ namespace CosmoTool class MapUtilityFunction { public: + typedef T type; typedef SpectrumFunction Spectrum; typedef boost::shared_ptr Spectrum_ptr; typedef FourierMap > FMap; diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp index 53b0278..4cbf6a3 100644 --- a/src/fourier/details/euclidian_maps.hpp +++ b/src/fourier/details/euclidian_maps.hpp @@ -36,6 +36,7 @@ knowledge of the CeCILL license and that you accept its terms. #ifndef __DETAILS_EUCLIDIAN_MAPS #define __DETAILS_EUCLIDIAN_MAPS +#include namespace CosmoTool { @@ -172,7 +173,7 @@ namespace CosmoTool } template - void get_Kvec(long p, Array2& kvec) const + void get_Kvec_p(long p, Array2& kvec) const { const DimArray& dims = this->getDims(); DimArray d(delta_k.size()); @@ -212,7 +213,7 @@ namespace CosmoTool return std::sqrt(k2); } - double get_K(long p) const + double get_K_p(long p) const { const DimArray& dims = this->getDims(); DimArray d(delta_k.size()); @@ -246,7 +247,7 @@ namespace CosmoTool { long idx = q0+dims[0]*p; assert(idx < this->size()); - result += 2*(conj(d1[idx]) * d2[idx]).real(); + result += T(2)*(std::conj(d1[idx]) * d2[idx]).real(); } } if (even0) @@ -254,8 +255,8 @@ namespace CosmoTool for (long p = 0; p < plane_size; p++) { long q0 = N0*p, q1 = (p+1)*N0-1; - result += conj(d1[q0]) * d2[q0]; - result += conj(d1[q1]) * d2[q1]; + result += T(2)*std::conj(d1[q0]) * d2[q0]; + result += T(2)*std::conj(d1[q1]) * d2[q1]; } } return result; diff --git a/src/fourier/details/euclidian_spectrum_1d.hpp b/src/fourier/details/euclidian_spectrum_1d.hpp index 36b524a..01a7dac 100644 --- a/src/fourier/details/euclidian_spectrum_1d.hpp +++ b/src/fourier/details/euclidian_spectrum_1d.hpp @@ -115,7 +115,7 @@ namespace CosmoTool for (long p = 1; p < rand_map.size(); p++) { - double A_k = std::sqrt(0.5*V*f(rand_map.get_K(p))); + double A_k = std::sqrt(0.5*V*f(rand_map.get_K_p(p))); d[p] = std::complex(gsl_ran_gaussian(rng, A_k), gsl_ran_gaussian(rng, A_k)); } @@ -138,7 +138,7 @@ namespace CosmoTool plane_size *= dims[q]; } - for (long p = 1; p < plane_size/2; p++) + for (long p = 1; p < plane_size/2+1; p++) { long q = (p+1)*dims[0]-1; long q2 = (plane_size-p+1)*dims[0]-1; @@ -147,6 +147,13 @@ namespace CosmoTool d[q] = conj(d[q2]); } + for (long p = 1; p < plane_size/2+1; p++) + { + long q = (p)*dims[0]; + long q2 = (plane_size-p)*dims[0]; + d[q] = conj(d[q2]); + } + if (alleven) { long q = 0; @@ -164,7 +171,7 @@ namespace CosmoTool std::complex *d = m.data(); for (long p = 0; p < m_c.size(); p++) - d[p] *= f(m_c.get_K(p)); + d[p] *= f(m_c.get_K_p(p)); } template @@ -174,7 +181,7 @@ namespace CosmoTool std::complex *d = m.data(); for (long p = 0; p < m_c.size(); p++) - d[p] *= std::sqrt(f(m_c.get_K(p))); + d[p] *= std::sqrt(f(m_c.get_K_p(p))); } template @@ -185,7 +192,7 @@ namespace CosmoTool for (long p = 0; p < m_c.size(); p++) { - T A = f(m_c.get_K(p)); + T A = f(m_c.get_K_p(p)); if (A==0) d[p] = 0; else @@ -201,7 +208,7 @@ namespace CosmoTool for (long p = 0; p < m_c.size(); p++) { - T A = std::sqrt(f(m_c.get_K(p))); + T A = std::sqrt(f(m_c.get_K_p(p))); if (A == 0) d[p] = 0; else diff --git a/src/fourier/details/euclidian_transform.hpp b/src/fourier/details/euclidian_transform.hpp index eaab4fd..e5ef34a 100644 --- a/src/fourier/details/euclidian_transform.hpp +++ b/src/fourier/details/euclidian_transform.hpp @@ -122,12 +122,22 @@ namespace CosmoTool calls::execute(m_synthesis); realMap->scale(1/volume); } + + void synthesis_unnormed() + { + calls::execute(m_synthesis); + } void analysis() { calls::execute(m_analysis); fourierMap->scale(volume/N); } + + void analysis_unnormed() + { + calls::execute(m_analysis); + } void synthesis_conjugate() { diff --git a/src/fourier/details/healpix_transform.hpp b/src/fourier/details/healpix_transform.hpp index b269471..37e79e8 100644 --- a/src/fourier/details/healpix_transform.hpp +++ b/src/fourier/details/healpix_transform.hpp @@ -67,7 +67,7 @@ namespace CosmoTool sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo); } - HealpixFourierTransform(long nSide, long Lmax, long Mmax, const std::valarray& weights, int iterate = 0 ) + HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate, const std::valarray& weights ) : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate) { sharp_make_weighted_healpix_geom_info (nSide, 1, &weights[0], &ginfo); diff --git a/src/hdf5_array.hpp b/src/hdf5_array.hpp new file mode 100644 index 0000000..46673ba --- /dev/null +++ b/src/hdf5_array.hpp @@ -0,0 +1,175 @@ +#ifndef __COSMO_HDF5_ARRAY_HPP +#define __COSMO_HDF5_ARRAY_HPP + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace CosmoTool { + + //!_______________________________________________________________________________________ + //! + //! map types to HDF5 types + //! + //! + //! Leo Goodstadt (04 March 2013), improved with enable_if by Guilhem Lavaux (May 2014) + //!_______________________________________________________________________________________ + + template struct get_hdf5_data_type + { + static H5::PredType type() + { + BOOST_MPL_ASSERT_MSG(0, Unknown_HDF5_data_type, ()); + return H5::PredType::NATIVE_DOUBLE; + } + }; + + #define HDF5_TYPE(tl, thdf5) \ + template struct get_hdf5_data_type >::type> \ + { static H5::PredType type() { return H5::PredType::thdf5; }; } + + HDF5_TYPE(char , NATIVE_CHAR); + HDF5_TYPE(long long , NATIVE_LLONG); + HDF5_TYPE(unsigned long long, NATIVE_ULLONG); + HDF5_TYPE(int8_t , NATIVE_INT8); + HDF5_TYPE(uint8_t , NATIVE_UINT8); + HDF5_TYPE(int16_t , NATIVE_INT16); + HDF5_TYPE(uint16_t , NATIVE_UINT16); + HDF5_TYPE(int32_t , NATIVE_INT32); + HDF5_TYPE(uint32_t , NATIVE_UINT32); + HDF5_TYPE(int64_t , NATIVE_INT64); + HDF5_TYPE(uint64_t , NATIVE_UINT64); + HDF5_TYPE(float , NATIVE_FLOAT); + HDF5_TYPE(double , NATIVE_DOUBLE); + HDF5_TYPE(long double , NATIVE_LDOUBLE); + + #undef HDF5_TYPE + +//!_______________________________________________________________________________________ +//! +//! write_hdf5 multi_array +//! +//! \author leo Goodstadt (04 March 2013) +//! +//!_______________________________________________________________________________________ + template + void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, + const boost::multi_array& data, + const hdf5_data_type& datatype) + { + std::vector dimensions(data.shape(), data.shape() + DIMENSIONS); + H5::DataSpace dataspace(DIMENSIONS, dimensions.data()); + + H5::DataSet dataset = fg.createDataSet(data_set_name, datatype, dataspace); + + dataset.write(data.data(), datatype); + } + + /* HDF5 complex type */ + template + class hdf5_ComplexType + { + public: + H5::CompType type; + + hdf5_ComplexType() + : type(sizeof(std::complex)) + { + get_hdf5_data_type hdf_data_type; + type.insertMember("r", 0, hdf_data_type.type()); + type.insertMember("i", sizeof(T), hdf_data_type.type()); + type.pack(); + } + + static const hdf5_ComplexType *ctype() + { + static hdf5_ComplexType singleton; + + return &singleton; + } + }; + + template + void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array& data ) + { + get_hdf5_data_type hdf_data_type; + hdf5_write_array(fg, data_set_name, data, hdf_data_type.type()); + } + + template + void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array, DIMENSIONS>& data ) + { + hdf5_write_array(fg, data_set_name, data, hdf5_ComplexType::ctype()->type); + } + + + // HDF5 array reader + // + // Author Guilhem Lavaux (May 2014) + + class InvalidDimensions: virtual std::exception { + }; + + template + struct hdf5_extent_gen { + typedef typename boost::detail::multi_array::extent_gen type; + + static inline type build(hsize_t *d) + { + return (hdf5_extent_gen::build(d))[d[r-1]]; + } + }; + + template<> + struct hdf5_extent_gen<0> { + static inline boost::multi_array_types::extent_gen build(hsize_t *d) + { + return boost::extents; + } + }; + + template + void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, + boost::multi_array& data, + const hdf5_data_type& datatype) + { + H5::DataSet dataset = fg.openDataSet(data_set_name); + H5::DataSpace dataspace = dataset.getSpace(); + std::vector dimensions(DIMENSIONS); + + if (dataspace.getSimpleExtentNdims() != DIMENSIONS) + { + throw InvalidDimensions(); + } + + dataspace.getSimpleExtentDims(dimensions.data()); + data.resize(hdf5_extent_gen::build(dimensions.data())); + dataset.read(data.data(), datatype); + } + + template + void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array& data ) + { + get_hdf5_data_type hdf_data_type; + hdf5_read_array(fg, data_set_name, data, hdf_data_type.type()); + } + + template + void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array, DIMENSIONS>& data ) + { + hdf5_read_array(fg, data_set_name, data, hdf5_ComplexType::ctype()->type); + } + +}; + +#endif + + diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp index 3eb909c..5154921 100644 --- a/src/loadRamses.cpp +++ b/src/loadRamses.cpp @@ -309,7 +309,7 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC); double L0 = info.boxSize*info.unitLength*hubble/100/CM_IN_MPC/info.aexp; - double unit_vel = L0*hubble/info.aexp; + double unit_vel = L0*100/info.aexp; while (1) { diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp index f6fe582..b347178 100644 --- a/src/sphSmooth.hpp +++ b/src/sphSmooth.hpp @@ -35,7 +35,7 @@ knowledge of the CeCILL license and that you accept its terms. #ifndef __COSMOTOOL_SPH_SMOOTH_HPP #define __COSMOTOOL_SPH_SMOOTH_HPP - +#include #include "config.hpp" #include "mykdtree.hpp" @@ -60,37 +60,48 @@ namespace CosmoTool typedef void (*runParticleValue)(ValType& t); public: + struct SPHState + { + boost::shared_ptr ngb; + boost::shared_ptr distances; + typename SPHTree::coords currentCenter; + int currentNgb; + ComputePrecision smoothRadius; + }; + + SPHSmooth(SPHTree *tree, uint32_t Nsph); virtual ~SPHSmooth(); - void fetchNeighbours(const typename SPHTree::coords& c); + void fetchNeighbours(const typename SPHTree::coords& c, SPHState *state = 0); + void fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNsph); void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius); const typename SPHTree::coords& getCurrentCenter() const { - return currentCenter; + return internal.currentCenter; } ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c, - computeParticleValue fun); + computeParticleValue fun, SPHState *state = 0); ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c, - computeParticleValue fun); + computeParticleValue fun, SPHState *state = 0); ComputePrecision getMaxDistance(const typename SPHTree::coords& c, SPHNode *node) const; ComputePrecision getSmoothingLen() const { - return smoothRadius; + return internal.smoothRadius; } // TO USE WITH EXTREME CARE ! void setSmoothingLen(ComputePrecision len) { - smoothRadius = len; + internal.smoothRadius = len; } // END - void runForEachNeighbour(runParticleValue fun); + void runForEachNeighbour(runParticleValue fun, SPHState *state = 0); void addGridSite(const typename SPHTree::coords& c); bool hasNeighbours() const; @@ -100,32 +111,26 @@ namespace CosmoTool SPHTree *getTree() { return tree; } void changeNgb(uint32_t newMax) { - delete[] ngb; - delete[] distances; - ngb = new SPHCell *[newMax]; - distances = new CoordType[newMax]; + internal.ngb = boost::shared_ptr(new SPHCell *[newMax]); + internal.distances = boost::shared_ptr(new CoordType[newMax]); maxNgb = newMax; } - uint32_t getCurrent() const { return currentNgb; } + uint32_t getCurrent() const { return internal.currentNgb; } uint32_t getNgb() const { return maxNgb; } protected: - SPHCell **ngb; - CoordType *distances; + SPHState internal; uint32_t Nsph; uint32_t deltaNsph; uint32_t maxNgb; - uint32_t currentNgb; SPHTree *tree; - ComputePrecision smoothRadius; - typename SPHTree::coords currentCenter; ComputePrecision computeWValue(const typename SPHTree::coords & c, SPHCell& cell, CoordType d, - computeParticleValue fun); + computeParticleValue fun, SPHState *state); void runUnrollNode(SPHNode *node, runParticleValue fun); }; diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc index 0cf7de8..1ae5486 100644 --- a/src/sphSmooth.tcc +++ b/src/sphSmooth.tcc @@ -1,4 +1,5 @@ #include +#include "algo.hpp" namespace CosmoTool { @@ -7,29 +8,27 @@ SPHSmooth::SPHSmooth(SPHTree *tree, uint32_t Nsph) { this->Nsph = Nsph; this->tree = tree; - this->currentNgb = 0; + internal.currentNgb = 0; this->maxNgb = Nsph; - ngb = new SPHCell *[maxNgb]; - distances = new CoordType[maxNgb]; + internal.ngb = boost::shared_ptr(new SPHCell *[maxNgb]); + internal.distances = boost::shared_ptr(new CoordType[maxNgb]); } template SPHSmooth::~SPHSmooth() { - delete[] ngb; - delete[] distances; } template ComputePrecision SPHSmooth::computeWValue(const typename SPHTree::coords& c, SPHCell& cell, CoordType d, - computeParticleValue fun) + computeParticleValue fun, SPHState *state) { CoordType weight; - d /= smoothRadius; + d /= state->smoothRadius; weight = getKernel(d); if (cell.val.weight != 0) @@ -47,86 +46,91 @@ SPHSmooth::fetchNeighbours(const typename SPHTree::coords& c, uin if (requested > maxNgb) { - delete[] ngb; - delete[] distances; - maxNgb = requested; - ngb = new SPHCell *[maxNgb]; - distances = new CoordType[maxNgb]; + internal.ngb = boost::shared_ptr(new SPHCell *[maxNgb]); + internal.distances = boost::shared_ptr(new CoordType[maxNgb]); } - memcpy(currentCenter, c, sizeof(c)); - tree->getNearestNeighbours(c, requested, ngb, distances); + memcpy(internal.currentCenter, c, sizeof(c)); + tree->getNearestNeighbours(c, requested, internal.ngb.get(), internal.distances.get()); - currentNgb = 0; - for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++) + internal.currentNgb = 0; + for (uint32_t i = 0; i < requested && internal.ngb[i] != 0; i++,internal.currentNgb++) { - distances[i] = sqrt(distances[i]); - d2 = distances[i]; + internal.distances[i] = sqrt(internal.distances[i]); + d2 = internal.distances[i]; if (d2 > max_dist) max_dist = d2; } - smoothRadius = max_dist / 2; + internal.smoothRadius = max_dist / 2; } template -void -SPHSmooth::fetchNeighbours(const typename SPHTree::coords& c) +void SPHSmooth::fetchNeighbours(const typename SPHTree::coords& c, SPHState *state) { ComputePrecision d2, max_dist = 0; uint32_t requested = Nsph; - memcpy(currentCenter, c, sizeof(c)); - tree->getNearestNeighbours(c, requested, ngb, distances); + if (state != 0) { + state->distances = boost::shared_ptr(new CoordType[Nsph]); + state->ngb = boost::shared_ptr(new SPHCell *[Nsph]); + } else + state = &internal; + + memcpy(state->currentCenter, c, sizeof(c)); + + tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get()); - currentNgb = 0; - for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++) + state->currentNgb = 0; + for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++) { - distances[i] = sqrt(distances[i]); - d2 = distances[i]; + d2 = internal.distances[i] = sqrt(internal.distances[i]); if (d2 > max_dist) - max_dist = d2; + max_dist = d2; } - smoothRadius = max_dist / 2; -} + state->smoothRadius = max_dist / 2; +} + template void SPHSmooth::fetchNeighboursOnVolume(const typename SPHTree::coords& c, - ComputePrecision radius) + ComputePrecision radius) { uint32_t numPart; ComputePrecision d2, max_dist = 0; - memcpy(currentCenter, c, sizeof(c)); + memcpy(internal.currentCenter, c, sizeof(c)); - currentNgb = tree->getIntersection(c, radius, ngb, distances, + internal.currentNgb = tree->getIntersection(c, radius, internal.ngb, internal.distances, maxNgb); - for (uint32_t i = 0; i < currentNgb; i++) + for (uint32_t i = 0; i < internal.currentNgb; i++) { - distances[i] = sqrt(distances[i]); - d2 = distances[i]; + d2 = internal.distances[i] = sqrt(internal.distances[i]); if (d2 > max_dist) - max_dist = d2; + max_dist = d2; } - smoothRadius = max_dist / 2; + internal.smoothRadius = max_dist / 2; } template ComputePrecision SPHSmooth::computeSmoothedValue(const typename SPHTree::coords& c, - computeParticleValue fun) + computeParticleValue fun, SPHState *state) { + if (state == 0) + state = &internal; + ComputePrecision outputValue = 0; ComputePrecision max_dist = 0; - ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius; + ComputePrecision r3 = cube(state->smoothRadius); - for (uint32_t i = 0; i < currentNgb; i++) + for (uint32_t i = 0; i < state->currentNgb; i++) { - outputValue += computeWValue(c, *ngb[i], distances[i], fun); + outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun, state); } return outputValue / r3; @@ -141,27 +145,33 @@ ComputePrecision interpolateOne(const ValType& t) // WARNING ! Cell's weight must be 1 !!! template ComputePrecision SPHSmooth::computeInterpolatedValue(const typename SPHTree::coords& c, - computeParticleValue fun) + computeParticleValue fun, SPHState *state) { + if (state == 0) + state = &internal; + ComputePrecision outputValue = 0; ComputePrecision max_dist = 0; ComputePrecision weight = 0; - for (uint32_t i = 0; i < currentNgb; i++) + for (uint32_t i = 0; i < state->currentNgb; i++) { - outputValue += computeWValue(c, *ngb[i], distances[i], fun); - weight += computeWValue(c, *ngb[i], distances[i], interpolateOne); + outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun); + weight += computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne); } return (outputValue == 0) ? 0 : (outputValue / weight); } template -void SPHSmooth::runForEachNeighbour(runParticleValue fun) +void SPHSmooth::runForEachNeighbour(runParticleValue fun, SPHState *state) { - for (uint32_t i = 0; i < currentNgb; i++) + if (state == 0) + state = &internal; + + for (uint32_t i = 0; i < state->currentNgb; i++) { - fun(ngb[i]); + fun(state->ngb[i]); } } @@ -172,13 +182,13 @@ void SPHSmooth::addGridSite(const typename SPHTree::coords& c) ComputePrecision outputValue = 0; ComputePrecision max_dist = 0; - ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius; + ComputePrecision r3 = cube(internal.smoothRadius); - for (uint32_t i = 0; i < currentNgb; i++) + for (uint32_t i = 0; i < internal.currentNgb; i++) { - ComputePrecision d = distances[i]; - SPHCell& cell = *ngb[i]; - cell.val.weight += getKernel(d/smoothRadius) / r3; + ComputePrecision d = internal.distances[i]; + SPHCell& cell = *internal.ngb[i]; + cell.val.weight += getKernel(d/internal.smoothRadius) / r3; } } @@ -202,7 +212,7 @@ SPHSmooth::getKernel(ComputePrecision x) const template bool SPHSmooth::hasNeighbours() const { - return (currentNgb != 0); + return (internal.currentNgb != 0); } template diff --git a/src/yorick.hpp b/src/yorick.hpp index 39fd20d..59d4656 100644 --- a/src/yorick.hpp +++ b/src/yorick.hpp @@ -218,7 +218,7 @@ namespace CosmoTool template void saveArray(const std::string& fname, - T *array, uint32_t *dimList, uint32_t rank); + const T *array, uint32_t *dimList, uint32_t rank); template void loadArray(const std::string& fname, diff --git a/src/yorick_nc3.cpp b/src/yorick_nc3.cpp index d80781e..8a69d9b 100644 --- a/src/yorick_nc3.cpp +++ b/src/yorick_nc3.cpp @@ -235,7 +235,7 @@ namespace CosmoTool { template void saveArray(const std::string& fname, - T *array, uint32_t *dimList, uint32_t rank) + const T *array, uint32_t *dimList, uint32_t rank) { NcFile f(fname.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4); @@ -300,10 +300,10 @@ namespace CosmoTool { double*& array, uint32_t *&dimList, uint32_t& rank); template void saveArray(const std::string& fname, - int *array, uint32_t *dimList, uint32_t rank); + const int *array, uint32_t *dimList, uint32_t rank); template void saveArray(const std::string& fname, - float *array, uint32_t *dimList, uint32_t rank); + const float *array, uint32_t *dimList, uint32_t rank); template void saveArray(const std::string& fname, - double *array, uint32_t *dimList, uint32_t rank); + const double *array, uint32_t *dimList, uint32_t rank); } diff --git a/src/yorick_nc4.cpp b/src/yorick_nc4.cpp index ed4cd9e..49ae5b5 100644 --- a/src/yorick_nc4.cpp +++ b/src/yorick_nc4.cpp @@ -203,7 +203,7 @@ namespace CosmoTool { template void saveArray(const std::string& fname, - T *array, uint32_t *dimList, uint32_t rank) + const T *array, uint32_t *dimList, uint32_t rank) { NcFile f(fname.c_str(), NcFile::replace); @@ -263,10 +263,10 @@ namespace CosmoTool { double*& array, uint32_t *&dimList, uint32_t& rank); template void saveArray(const std::string& fname, - int *array, uint32_t *dimList, uint32_t rank); + const int *array, uint32_t *dimList, uint32_t rank); template void saveArray(const std::string& fname, - float *array, uint32_t *dimList, uint32_t rank); + const float *array, uint32_t *dimList, uint32_t rank); template void saveArray(const std::string& fname, - double *array, uint32_t *dimList, uint32_t rank); + const double *array, uint32_t *dimList, uint32_t rank); }