From b668f9947dbde3280902406a5eac9bce6c23c9f1 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 9 Jun 2015 20:01:56 +0200 Subject: [PATCH] fixed include ordering in h5_readFlash. Expanded limits of kdtree. Better support for NETCDF4/HDF5. --- external/external_build.cmake | 6 ++++- python/cosmotool/grafic.py | 7 ++++-- src/h5_readFlash.cpp | 3 --- src/h5_readFlash.hpp | 11 ++++----- src/kdtree_splitters.hpp | 12 +++++----- src/mykdtree.hpp | 26 ++++++++++++---------- src/mykdtree.tcc | 42 ++++++++++++++++++++--------------- 7 files changed, 60 insertions(+), 47 deletions(-) diff --git a/external/external_build.cmake b/external/external_build.cmake index 8725fbc..2fd6380 100644 --- a/external/external_build.cmake +++ b/external/external_build.cmake @@ -12,6 +12,7 @@ SET(GSL_URL "ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz" CACHE STRING "URL to dow mark_as_advanced(FFTW_URL EIGEN_URL HDF5_URL NETCDF_URL BOOST_URL GSL_URL) + IF(ENABLE_OPENMP) IF (NOT OPENMP_FOUND) MESSAGE(ERROR "No known compiler option for enabling OpenMP") @@ -25,6 +26,8 @@ ENDIF(ENABLE_OPENMP) SET(BUILD_PREFIX ${CMAKE_BINARY_DIR}/external_build) SET(EXT_INSTALL ${CMAKE_BINARY_DIR}/ext_install) +SET(CONFIGURE_LIBS) +SET(CONFIGURE_CPP_FLAGS) if (ENABLE_SHARP) SET(SHARP_SOURCE ${CMAKE_SOURCE_DIR}/external/sharp) @@ -72,6 +75,7 @@ if (INTERNAL_HDF5) SET(ENV{HDF5_ROOT} ${HDF5_BIN_DIR}) SET(HDF5_ROOTDIR ${HDF5_BIN_DIR}) SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${HDF5_BIN_DIR}/lib") + SET(CONFIGURE_LIBS "${CONFIGURE_LIBS} -ldl") else(INTERNAL_HDF5) find_path(HDF5_INCLUDE_PATH hdf5.h) find_library(HDF5_LIBRARY hdf5) @@ -92,7 +96,7 @@ if (INTERNAL_NETCDF) SET(NETCDF_BIN_DIR ${EXT_INSTALL}) SET(CONFIGURE_CPP_FLAGS "${CONFIGURE_CPP_FLAGS} -I${NETCDF_BIN_DIR}/include") SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${NETCDF_BIN_DIR}/lib") - SET(EXTRA_NC_FLAGS CPPFLAGS=${CONFIGURE_CPP_FLAGS} LDFLAGS=${CONFIGURE_LDFLAGS}) + SET(EXTRA_NC_FLAGS CPPFLAGS=${CONFIGURE_CPP_FLAGS} LIBS=${CONFIGURE_LIBS} LDFLAGS=${CONFIGURE_LDFLAGS}) ExternalProject_Add(netcdf DEPENDS ${hdf5_built} PREFIX ${BUILD_PREFIX}/netcdf-prefix diff --git a/python/cosmotool/grafic.py b/python/cosmotool/grafic.py index 9310172..c23073b 100644 --- a/python/cosmotool/grafic.py +++ b/python/cosmotool/grafic.py @@ -72,7 +72,7 @@ def writeWhitePhase(filename, field): f.write(struct.pack("IIIIII", checkPoint, Nx, Ny, Nz, 0, checkPoint)) field = field.reshape(field.shape, order='F') - checkPoint = struct.pack("I", 4*Nx*Ny) + checkPoint = struct.pack("I", 4*Ny*Nz) for i in range(Nx): f.write(checkPoint) f.write(field[i].astype(np.float32).tostring()) @@ -91,7 +91,10 @@ def readWhitePhase(filename): if struct.unpack("I", f.read(4))[0] != checkPoint_ref: raise ValueError("Invalid unformatted access") - a[i, :, :] = np.fromfile(f, dtype=np.float32, count=Ny*Nz).reshape((Ny, Nz), order='F') + b = np.fromfile(f, dtype=np.float32, count=Ny*Nz).reshape((Ny, Nz)) + if i==0: + print(b) + a[i, : ,:] = b if struct.unpack("I", f.read(4))[0] != checkPoint_ref: raise ValueError("Invalid unformatted access") diff --git a/src/h5_readFlash.cpp b/src/h5_readFlash.cpp index 2ace735..3d6d205 100644 --- a/src/h5_readFlash.cpp +++ b/src/h5_readFlash.cpp @@ -212,7 +212,6 @@ void h5_read_runtime_parameters *numPart = *numPart * *numPart * *numPart; } } - } @@ -253,7 +252,6 @@ void h5_read_flash3_particles (H5File* file, return; } - /* first determine how many particle properties are present in the input data file, and determine which of these are the properties we are interested in */ @@ -424,7 +422,6 @@ void h5_read_flash3_particles (H5File* file, //status = H5Tclose(datatype); //status = H5Sclose(dataspace); //status = H5Dclose(dataset); - } diff --git a/src/h5_readFlash.hpp b/src/h5_readFlash.hpp index 3c6982d..56f290f 100644 --- a/src/h5_readFlash.hpp +++ b/src/h5_readFlash.hpp @@ -45,12 +45,13 @@ This file has been developped by P. M. Sutter. * The functions accept the PARAMESH data through arguments. */ -#include -#include -#include -#include -#include "hdf5_flash.h" +#include +#include +//#include +#include +#include #include "H5Cpp.h" +#include "hdf5_flash.h" using namespace H5; diff --git a/src/kdtree_splitters.hpp b/src/kdtree_splitters.hpp index 8f0e22e..d1f00a9 100644 --- a/src/kdtree_splitters.hpp +++ b/src/kdtree_splitters.hpp @@ -48,17 +48,17 @@ namespace CosmoTool typedef typename KDDef::CoordType ctype; - void check_splitting(KDCell **cells, uint32_t Ncells, int axis, uint32_t split_index, ctype midCoord) + void check_splitting(KDCell **cells, NodeIntType Ncells, int axis, NodeIntType split_index, ctype midCoord) { ctype delta = std::numeric_limits::max(); assert(split_index < Ncells); assert(axis < N); - for (uint32_t i = 0; i < split_index; i++) + for (NodeIntType i = 0; i < split_index; i++) { assert(cells[i]->coord[axis] <= midCoord); delta = min(midCoord-cells[i]->coord[axis], delta); } - for (uint32_t i = split_index+1; i < Ncells; i++) + for (NodeIntType i = split_index+1; i < Ncells; i++) { assert(cells[i]->coord[axis] > midCoord); delta = min(cells[i]->coord[axis]-midCoord, delta); @@ -67,7 +67,7 @@ namespace CosmoTool assert (std::abs(cells[split_index]->coord[axis]-midCoord) <= delta); } - void operator()(KDCell **cells, uint32_t Ncells, uint32_t& split_index, int axis, coords minBound, coords maxBound) + void operator()(KDCell **cells, NodeIntType Ncells, NodeIntType& split_index, int axis, coords minBound, coords maxBound) { if (Ncells == 1) { @@ -76,9 +76,9 @@ namespace CosmoTool } ctype midCoord = 0.5*(maxBound[axis]+minBound[axis]); - uint32_t below = 0, above = Ncells-1; + NodeIntType below = 0, above = Ncells-1; ctype delta_min = std::numeric_limits::max(); - uint32_t idx_min = std::numeric_limits::max(); + NodeIntType idx_min = std::numeric_limits::max(); while (below < above) { diff --git a/src/mykdtree.hpp b/src/mykdtree.hpp index 10fac70..e84dea2 100644 --- a/src/mykdtree.hpp +++ b/src/mykdtree.hpp @@ -46,6 +46,8 @@ knowledge of the CeCILL license and that you accept its terms. namespace CosmoTool { + typedef uint64_t NodeIntType; + template struct KDDef { @@ -84,7 +86,7 @@ namespace CosmoTool { KDTreeNode *children[2]; typename KDDef::KDCoordinates minBound, maxBound; #ifdef __KD_TREE_NUMNODES - uint32_t numNodes; + NodeIntType numNodes; #endif }; @@ -122,7 +124,7 @@ namespace CosmoTool { template struct KD_default_cell_splitter { - void operator()(KDCell **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef::KDCoordinates minBound, typename KDDef::KDCoordinates maxBound); + void operator()(KDCell **cells, NodeIntType Ncells, NodeIntType& split_index, int axis, typename KDDef::KDCoordinates minBound, typename KDDef::KDCoordinates maxBound); }; template > @@ -136,7 +138,7 @@ namespace CosmoTool { CellSplitter splitter; - KDTree(Cell *cells, uint32_t Ncells); + KDTree(Cell *cells, NodeIntType Ncells); ~KDTree(); void setPeriodic(bool on, CoordType replicate) @@ -175,14 +177,14 @@ namespace CosmoTool { void optimize(); Node *getAllNodes() { return nodes; } - uint32_t getNumNodes() const { return lastNode; } + NodeIntType getNumNodes() const { return lastNode; } - uint32_t countActives() const; + NodeIntType countActives() const; #ifdef __KD_TREE_NUMNODES - uint32_t getNumberInNode(const Node *n) const { return n->numNodes; } + NodeIntType getNumberInNode(const Node *n) const { return n->numNodes; } #else - uint32_t getNumberInNode(const Node *n) const { + NodeIntType getNumberInNode(const Node *n) const { if (n == 0) return 0; return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]); @@ -190,15 +192,15 @@ namespace CosmoTool { #endif #ifdef __KD_TREE_SAVE_ON_DISK - KDTree(std::istream& i, Cell *cells, uint32_t Ncells) + KDTree(std::istream& i, Cell *cells, NodeIntType Ncells) throw (InvalidOnDiskKDTree); void saveTree(std::ostream& o) const; #endif protected: Node *nodes; - uint32_t numNodes; - uint32_t lastNode; + NodeIntType numNodes; + NodeIntType lastNode; Node *root; Cell **sortingHelper; @@ -208,7 +210,7 @@ namespace CosmoTool { coords replicate; Node *buildTree(Cell **cell0, - uint32_t NumCells, + NodeIntType NumCells, uint32_t depth, coords minBound, coords maxBound); @@ -231,7 +233,7 @@ namespace CosmoTool { }; template - uint32_t gatherActiveCells(KDCell **cells, uint32_t numCells); + NodeIntType gatherActiveCells(KDCell **cells, NodeIntType numCells); }; diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc index 62c74da..e9ce02c 100644 --- a/src/mykdtree.tcc +++ b/src/mykdtree.tcc @@ -31,7 +31,7 @@ namespace CosmoTool { } template - KDTree::KDTree(Cell *cells, uint32_t Ncells) + KDTree::KDTree(Cell *cells, NodeIntType Ncells) { periodic = false; @@ -40,7 +40,7 @@ namespace CosmoTool { nodes = new Node[numNodes]; sortingHelper = new Cell *[Ncells]; - for (uint32_t i = 0; i < Ncells; i++) + for (NodeIntType i = 0; i < Ncells; i++) sortingHelper[i] = &cells[i]; optimize(); @@ -52,7 +52,7 @@ namespace CosmoTool { coords absoluteMin, absoluteMax; std::cout << "Optimizing the tree..." << std::endl; - uint32_t activeCells = gatherActiveCells(sortingHelper, numNodes); + NodeIntType activeCells = gatherActiveCells(sortingHelper, numNodes); std::cout << " number of active cells = " << activeCells << std::endl; lastNode = 0; @@ -62,7 +62,7 @@ namespace CosmoTool { absoluteMax[i] = -std::numeric_limits::max(); } // Find min and max corner - for (uint32_t i = 0; i < activeCells; i++) + for (NodeIntType i = 0; i < activeCells; i++) { KDCell *cell = sortingHelper[i]; @@ -226,11 +226,11 @@ namespace CosmoTool { } template - uint32_t gatherActiveCells(KDCell **cells, - uint32_t Ncells) + NodeIntType gatherActiveCells(KDCell **cells, + NodeIntType Ncells) { - uint32_t swapId = Ncells-1; - uint32_t i = 0; + NodeIntType swapId = Ncells-1; + NodeIntType i = 0; #if __KD_TREE_ACTIVE_CELLS == 1 while (!cells[swapId]->active && swapId > 0) @@ -253,7 +253,10 @@ namespace CosmoTool { } template - void KD_default_cell_splitter::operator()(KDCell **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef::KDCoordinates minBound, typename KDDef::KDCoordinates maxBound) + void KD_default_cell_splitter::operator()(KDCell **cells, NodeIntType Ncells, + NodeIntType& split_index, int axis, + typename KDDef::KDCoordinates minBound, + typename KDDef::KDCoordinates maxBound) { CellCompare compare(axis); omptl::sort(cells,cells+Ncells,compare); // std::sort(cells, cells+Ncells, compare); @@ -263,7 +266,7 @@ namespace CosmoTool { template KDTreeNode *KDTree::buildTree(Cell **cell0, - uint32_t Ncells, + NodeIntType Ncells, uint32_t depth, coords minBound, coords maxBound) @@ -273,11 +276,14 @@ namespace CosmoTool { Node *node; int axis = depth % N; - uint32_t mid; + NodeIntType mid; coords tmpBound; + NodeIntType nodeId; -#pragma omp critical - node = &nodes[lastNode++]; +#pragma omp atomic + nodeId = (this->lastNode)++; + + node = &nodes[nodeId]; // Isolate the environment splitter(cell0, Ncells, mid, axis, minBound, maxBound); @@ -315,10 +321,10 @@ namespace CosmoTool { } template - uint32_t KDTree::countActives() const + NodeIntType KDTree::countActives() const { - uint32_t numActive = 0; - for (uint32_t i = 0; i < lastNode; i++) + NodeIntType numActive = 0; + for (NodeIntType i = 0; i < lastNode; i++) { #if __KD_TREE_ACTIVE_CELLS == 1 if (nodes[i].value->active) @@ -604,7 +610,7 @@ namespace CosmoTool { } template - KDTree::KDTree(std::istream& in, Cell *cells, uint32_t Ncells) + KDTree::KDTree(std::istream& in, Cell *cells, NodeIntType Ncells) throw (InvalidOnDiskKDTree) { KDTreeHeader h; @@ -683,7 +689,7 @@ namespace CosmoTool { root = &nodes[h.rootId]; sortingHelper = new Cell *[Ncells]; - for (uint32_t i = 0; i < Ncells; i++) + for (NodeIntType i = 0; i < Ncells; i++) sortingHelper[i] = &cells[i]; } #endif