fixed include ordering in h5_readFlash. Expanded limits of kdtree. Better support for NETCDF4/HDF5.

This commit is contained in:
Guilhem Lavaux 2015-06-09 20:01:56 +02:00
parent 59b1b12be8
commit b668f9947d
7 changed files with 60 additions and 47 deletions

View File

@ -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) mark_as_advanced(FFTW_URL EIGEN_URL HDF5_URL NETCDF_URL BOOST_URL GSL_URL)
IF(ENABLE_OPENMP) IF(ENABLE_OPENMP)
IF (NOT OPENMP_FOUND) IF (NOT OPENMP_FOUND)
MESSAGE(ERROR "No known compiler option for enabling OpenMP") 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(BUILD_PREFIX ${CMAKE_BINARY_DIR}/external_build)
SET(EXT_INSTALL ${CMAKE_BINARY_DIR}/ext_install) SET(EXT_INSTALL ${CMAKE_BINARY_DIR}/ext_install)
SET(CONFIGURE_LIBS)
SET(CONFIGURE_CPP_FLAGS)
if (ENABLE_SHARP) if (ENABLE_SHARP)
SET(SHARP_SOURCE ${CMAKE_SOURCE_DIR}/external/sharp) SET(SHARP_SOURCE ${CMAKE_SOURCE_DIR}/external/sharp)
@ -72,6 +75,7 @@ if (INTERNAL_HDF5)
SET(ENV{HDF5_ROOT} ${HDF5_BIN_DIR}) SET(ENV{HDF5_ROOT} ${HDF5_BIN_DIR})
SET(HDF5_ROOTDIR ${HDF5_BIN_DIR}) SET(HDF5_ROOTDIR ${HDF5_BIN_DIR})
SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${HDF5_BIN_DIR}/lib") SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${HDF5_BIN_DIR}/lib")
SET(CONFIGURE_LIBS "${CONFIGURE_LIBS} -ldl")
else(INTERNAL_HDF5) else(INTERNAL_HDF5)
find_path(HDF5_INCLUDE_PATH hdf5.h) find_path(HDF5_INCLUDE_PATH hdf5.h)
find_library(HDF5_LIBRARY hdf5) find_library(HDF5_LIBRARY hdf5)
@ -92,7 +96,7 @@ if (INTERNAL_NETCDF)
SET(NETCDF_BIN_DIR ${EXT_INSTALL}) SET(NETCDF_BIN_DIR ${EXT_INSTALL})
SET(CONFIGURE_CPP_FLAGS "${CONFIGURE_CPP_FLAGS} -I${NETCDF_BIN_DIR}/include") SET(CONFIGURE_CPP_FLAGS "${CONFIGURE_CPP_FLAGS} -I${NETCDF_BIN_DIR}/include")
SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${NETCDF_BIN_DIR}/lib") 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 ExternalProject_Add(netcdf
DEPENDS ${hdf5_built} DEPENDS ${hdf5_built}
PREFIX ${BUILD_PREFIX}/netcdf-prefix PREFIX ${BUILD_PREFIX}/netcdf-prefix

View File

@ -72,7 +72,7 @@ def writeWhitePhase(filename, field):
f.write(struct.pack("IIIIII", checkPoint, Nx, Ny, Nz, 0, checkPoint)) f.write(struct.pack("IIIIII", checkPoint, Nx, Ny, Nz, 0, checkPoint))
field = field.reshape(field.shape, order='F') 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): for i in range(Nx):
f.write(checkPoint) f.write(checkPoint)
f.write(field[i].astype(np.float32).tostring()) 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: if struct.unpack("I", f.read(4))[0] != checkPoint_ref:
raise ValueError("Invalid unformatted access") 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: if struct.unpack("I", f.read(4))[0] != checkPoint_ref:
raise ValueError("Invalid unformatted access") raise ValueError("Invalid unformatted access")

View File

@ -212,7 +212,6 @@ void h5_read_runtime_parameters
*numPart = *numPart * *numPart * *numPart; *numPart = *numPart * *numPart * *numPart;
} }
} }
} }
@ -253,7 +252,6 @@ void h5_read_flash3_particles (H5File* file,
return; return;
} }
/* first determine how many particle properties are /* first determine how many particle properties are
present in the input data file, and determine which of these present in the input data file, and determine which of these
are the properties we are interested in */ are the properties we are interested in */
@ -424,7 +422,6 @@ void h5_read_flash3_particles (H5File* file,
//status = H5Tclose(datatype); //status = H5Tclose(datatype);
//status = H5Sclose(dataspace); //status = H5Sclose(dataspace);
//status = H5Dclose(dataset); //status = H5Dclose(dataset);
} }

View File

@ -45,12 +45,13 @@ This file has been developped by P. M. Sutter.
* The functions accept the PARAMESH data through arguments. * The functions accept the PARAMESH data through arguments.
*/ */
#include <stdio.h> #include <cstdio>
#include <stdlib.h> #include <cstdlib>
#include <stddef.h> //#include <stddef.h>
#include <string.h> #include <string>
#include "hdf5_flash.h" #include <cstring>
#include "H5Cpp.h" #include "H5Cpp.h"
#include "hdf5_flash.h"
using namespace H5; using namespace H5;

View File

@ -48,17 +48,17 @@ namespace CosmoTool
typedef typename KDDef<N,CType>::CoordType ctype; typedef typename KDDef<N,CType>::CoordType ctype;
void check_splitting(KDCell<N,ValType,CType> **cells, uint32_t Ncells, int axis, uint32_t split_index, ctype midCoord) void check_splitting(KDCell<N,ValType,CType> **cells, NodeIntType Ncells, int axis, NodeIntType split_index, ctype midCoord)
{ {
ctype delta = std::numeric_limits<ctype>::max(); ctype delta = std::numeric_limits<ctype>::max();
assert(split_index < Ncells); assert(split_index < Ncells);
assert(axis < N); 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); assert(cells[i]->coord[axis] <= midCoord);
delta = min(midCoord-cells[i]->coord[axis], delta); 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); assert(cells[i]->coord[axis] > midCoord);
delta = min(cells[i]->coord[axis]-midCoord, delta); delta = min(cells[i]->coord[axis]-midCoord, delta);
@ -67,7 +67,7 @@ namespace CosmoTool
assert (std::abs(cells[split_index]->coord[axis]-midCoord) <= delta); assert (std::abs(cells[split_index]->coord[axis]-midCoord) <= delta);
} }
void operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, coords minBound, coords maxBound) void operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells, NodeIntType& split_index, int axis, coords minBound, coords maxBound)
{ {
if (Ncells == 1) if (Ncells == 1)
{ {
@ -76,9 +76,9 @@ namespace CosmoTool
} }
ctype midCoord = 0.5*(maxBound[axis]+minBound[axis]); 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<ctype>::max(); ctype delta_min = std::numeric_limits<ctype>::max();
uint32_t idx_min = std::numeric_limits<uint32_t>::max(); NodeIntType idx_min = std::numeric_limits<NodeIntType>::max();
while (below < above) while (below < above)
{ {

View File

@ -46,6 +46,8 @@ knowledge of the CeCILL license and that you accept its terms.
namespace CosmoTool { namespace CosmoTool {
typedef uint64_t NodeIntType;
template<int N, typename CType = ComputePrecision> template<int N, typename CType = ComputePrecision>
struct KDDef struct KDDef
{ {
@ -84,7 +86,7 @@ namespace CosmoTool {
KDTreeNode<N,ValType,CType> *children[2]; KDTreeNode<N,ValType,CType> *children[2];
typename KDDef<N,CType>::KDCoordinates minBound, maxBound; typename KDDef<N,CType>::KDCoordinates minBound, maxBound;
#ifdef __KD_TREE_NUMNODES #ifdef __KD_TREE_NUMNODES
uint32_t numNodes; NodeIntType numNodes;
#endif #endif
}; };
@ -122,7 +124,7 @@ namespace CosmoTool {
template<int N, typename ValType, typename CType = ComputePrecision> template<int N, typename ValType, typename CType = ComputePrecision>
struct KD_default_cell_splitter struct KD_default_cell_splitter
{ {
void operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound); void operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells, NodeIntType& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound);
}; };
template<int N, typename ValType, typename CType = ComputePrecision, typename CellSplitter = KD_default_cell_splitter<N,ValType,CType> > template<int N, typename ValType, typename CType = ComputePrecision, typename CellSplitter = KD_default_cell_splitter<N,ValType,CType> >
@ -136,7 +138,7 @@ namespace CosmoTool {
CellSplitter splitter; CellSplitter splitter;
KDTree(Cell *cells, uint32_t Ncells); KDTree(Cell *cells, NodeIntType Ncells);
~KDTree(); ~KDTree();
void setPeriodic(bool on, CoordType replicate) void setPeriodic(bool on, CoordType replicate)
@ -175,14 +177,14 @@ namespace CosmoTool {
void optimize(); void optimize();
Node *getAllNodes() { return nodes; } 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 #ifdef __KD_TREE_NUMNODES
uint32_t getNumberInNode(const Node *n) const { return n->numNodes; } NodeIntType getNumberInNode(const Node *n) const { return n->numNodes; }
#else #else
uint32_t getNumberInNode(const Node *n) const { NodeIntType getNumberInNode(const Node *n) const {
if (n == 0) if (n == 0)
return 0; return 0;
return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]); return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]);
@ -190,15 +192,15 @@ namespace CosmoTool {
#endif #endif
#ifdef __KD_TREE_SAVE_ON_DISK #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); throw (InvalidOnDiskKDTree);
void saveTree(std::ostream& o) const; void saveTree(std::ostream& o) const;
#endif #endif
protected: protected:
Node *nodes; Node *nodes;
uint32_t numNodes; NodeIntType numNodes;
uint32_t lastNode; NodeIntType lastNode;
Node *root; Node *root;
Cell **sortingHelper; Cell **sortingHelper;
@ -208,7 +210,7 @@ namespace CosmoTool {
coords replicate; coords replicate;
Node *buildTree(Cell **cell0, Node *buildTree(Cell **cell0,
uint32_t NumCells, NodeIntType NumCells,
uint32_t depth, uint32_t depth,
coords minBound, coords minBound,
coords maxBound); coords maxBound);
@ -231,7 +233,7 @@ namespace CosmoTool {
}; };
template<int N, typename T, typename CType> template<int N, typename T, typename CType>
uint32_t gatherActiveCells(KDCell<N,T,CType> **cells, uint32_t numCells); NodeIntType gatherActiveCells(KDCell<N,T,CType> **cells, NodeIntType numCells);
}; };

View File

@ -31,7 +31,7 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType, typename CellSplitter> template<int N, typename ValType, typename CType, typename CellSplitter>
KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, uint32_t Ncells) KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, NodeIntType Ncells)
{ {
periodic = false; periodic = false;
@ -40,7 +40,7 @@ namespace CosmoTool {
nodes = new Node[numNodes]; nodes = new Node[numNodes];
sortingHelper = new Cell *[Ncells]; sortingHelper = new Cell *[Ncells];
for (uint32_t i = 0; i < Ncells; i++) for (NodeIntType i = 0; i < Ncells; i++)
sortingHelper[i] = &cells[i]; sortingHelper[i] = &cells[i];
optimize(); optimize();
@ -52,7 +52,7 @@ namespace CosmoTool {
coords absoluteMin, absoluteMax; coords absoluteMin, absoluteMax;
std::cout << "Optimizing the tree..." << std::endl; 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; std::cout << " number of active cells = " << activeCells << std::endl;
lastNode = 0; lastNode = 0;
@ -62,7 +62,7 @@ namespace CosmoTool {
absoluteMax[i] = -std::numeric_limits<typeof (absoluteMax[0])>::max(); absoluteMax[i] = -std::numeric_limits<typeof (absoluteMax[0])>::max();
} }
// Find min and max corner // Find min and max corner
for (uint32_t i = 0; i < activeCells; i++) for (NodeIntType i = 0; i < activeCells; i++)
{ {
KDCell<N,ValType,CType> *cell = sortingHelper[i]; KDCell<N,ValType,CType> *cell = sortingHelper[i];
@ -226,11 +226,11 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType> template<int N, typename ValType, typename CType>
uint32_t gatherActiveCells(KDCell<N,ValType,CType> **cells, NodeIntType gatherActiveCells(KDCell<N,ValType,CType> **cells,
uint32_t Ncells) NodeIntType Ncells)
{ {
uint32_t swapId = Ncells-1; NodeIntType swapId = Ncells-1;
uint32_t i = 0; NodeIntType i = 0;
#if __KD_TREE_ACTIVE_CELLS == 1 #if __KD_TREE_ACTIVE_CELLS == 1
while (!cells[swapId]->active && swapId > 0) while (!cells[swapId]->active && swapId > 0)
@ -253,7 +253,10 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType> template<int N, typename ValType, typename CType>
void KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound) void KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells,
NodeIntType& split_index, int axis,
typename KDDef<N,CType>::KDCoordinates minBound,
typename KDDef<N,CType>::KDCoordinates maxBound)
{ {
CellCompare<N,ValType,CType> compare(axis); CellCompare<N,ValType,CType> compare(axis);
omptl::sort(cells,cells+Ncells,compare); // std::sort(cells, cells+Ncells, compare); omptl::sort(cells,cells+Ncells,compare); // std::sort(cells, cells+Ncells, compare);
@ -263,7 +266,7 @@ namespace CosmoTool {
template<int N, typename ValType, typename CType, typename CellSplitter> template<int N, typename ValType, typename CType, typename CellSplitter>
KDTreeNode<N,ValType,CType> *KDTree<N,ValType,CType,CellSplitter>::buildTree(Cell **cell0, KDTreeNode<N,ValType,CType> *KDTree<N,ValType,CType,CellSplitter>::buildTree(Cell **cell0,
uint32_t Ncells, NodeIntType Ncells,
uint32_t depth, uint32_t depth,
coords minBound, coords minBound,
coords maxBound) coords maxBound)
@ -273,11 +276,14 @@ namespace CosmoTool {
Node *node; Node *node;
int axis = depth % N; int axis = depth % N;
uint32_t mid; NodeIntType mid;
coords tmpBound; coords tmpBound;
NodeIntType nodeId;
#pragma omp critical #pragma omp atomic
node = &nodes[lastNode++]; nodeId = (this->lastNode)++;
node = &nodes[nodeId];
// Isolate the environment // Isolate the environment
splitter(cell0, Ncells, mid, axis, minBound, maxBound); splitter(cell0, Ncells, mid, axis, minBound, maxBound);
@ -315,10 +321,10 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType, typename CellSplitter> template<int N, typename ValType, typename CType, typename CellSplitter>
uint32_t KDTree<N,ValType,CType,CellSplitter>::countActives() const NodeIntType KDTree<N,ValType,CType,CellSplitter>::countActives() const
{ {
uint32_t numActive = 0; NodeIntType numActive = 0;
for (uint32_t i = 0; i < lastNode; i++) for (NodeIntType i = 0; i < lastNode; i++)
{ {
#if __KD_TREE_ACTIVE_CELLS == 1 #if __KD_TREE_ACTIVE_CELLS == 1
if (nodes[i].value->active) if (nodes[i].value->active)
@ -604,7 +610,7 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType, typename CellSplitter> template<int N, typename ValType, typename CType, typename CellSplitter>
KDTree<N,ValType,CType,CellSplitter>::KDTree(std::istream& in, Cell *cells, uint32_t Ncells) KDTree<N,ValType,CType,CellSplitter>::KDTree(std::istream& in, Cell *cells, NodeIntType Ncells)
throw (InvalidOnDiskKDTree) throw (InvalidOnDiskKDTree)
{ {
KDTreeHeader h; KDTreeHeader h;
@ -683,7 +689,7 @@ namespace CosmoTool {
root = &nodes[h.rootId]; root = &nodes[h.rootId];
sortingHelper = new Cell *[Ncells]; sortingHelper = new Cell *[Ncells];
for (uint32_t i = 0; i < Ncells; i++) for (NodeIntType i = 0; i < Ncells; i++)
sortingHelper[i] = &cells[i]; sortingHelper[i] = &cells[i];
} }
#endif #endif