diff --git a/CMakeLists.txt b/CMakeLists.txt index b9d2edf..1b3c80e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -39,13 +39,13 @@ MESSAGE(STATUS "Using the target ${CosmoTool_local} to build python module") include(${CMAKE_SOURCE_DIR}/external/external_build.cmake) IF(YORICK_SUPPORT) - IF(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf AND ${NETCDFCPP_LIBRARY} MATCHES "netcdf_c\\+\\+4") + IF((EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf AND ${NETCDFCPP_LIBRARY} MATCHES "netcdf_c\\+\\+4") OR (INTERNAL_NETCDF)) SET(FOUND_NETCDF4 1) FILE(WRITE ${CMAKE_BINARY_DIR}/src/ctool_netcdf_ver.hpp "#define NETCDFCPP4 1") - ELSE(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf AND ${NETCDFCPP_LIBRARY} MATCHES "netcdf_c\\+\\+4") + ELSE((EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf AND ${NETCDFCPP_LIBRARY} MATCHES "netcdf_c\\+\\+4") OR (INTERNAL_NETCDF)) SET(FOUND_NETCDF3 1) FILE(WRITE ${CMAKE_BINARY_DIR}/src/ctool_netcdf_ver.hpp "#undef NETCDFCPP4") - ENDIF(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf AND ${NETCDFCPP_LIBRARY} MATCHES "netcdf_c\\+\\+4") + ENDIF((EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf AND ${NETCDFCPP_LIBRARY} MATCHES "netcdf_c\\+\\+4") OR (INTERNAL_NETCDF)) ENDIF(YORICK_SUPPORT) find_program(CYTHON cython) diff --git a/external/external_build.cmake b/external/external_build.cmake index 170cf8f..8725fbc 100644 --- a/external/external_build.cmake +++ b/external/external_build.cmake @@ -5,7 +5,7 @@ OPTION(ENABLE_OPENMP "Set to Yes if Healpix and/or you need openMP" OFF) SET(FFTW_URL "http://www.fftw.org/fftw-3.3.3.tar.gz" CACHE URL "URL to download FFTW from") SET(EIGEN_URL "http://bitbucket.org/eigen/eigen/get/3.1.4.tar.gz" CACHE URL "URL to download Eigen from") SET(GENGETOPT_URL "ftp://ftp.gnu.org/gnu/gengetopt/gengetopt-2.22.5.tar.gz" CACHE STRING "URL to download gengetopt from") -SET(HDF5_URL "http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.9/src/hdf5-1.8.9.tar.gz" CACHE STRING "URL to download HDF5 from") +SET(HDF5_URL "http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.15/src/hdf5-1.8.15.tar.gz" CACHE STRING "URL to download HDF5 from") SET(NETCDF_URL "http://www.unidata.ucar.edu/downloads/netcdf/ftp/netcdf-4.1.3.tar.gz" CACHE STRING "URL to download NetCDF from") SET(BOOST_URL "http://sourceforge.net/projects/boost/files/boost/1.49.0/boost_1_49_0.tar.gz/download" CACHE STRING "URL to download Boost from") SET(GSL_URL "ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz" CACHE STRING "URL to download GSL from ") @@ -142,6 +142,7 @@ if (INTERNAL_BOOST) ) set(Boost_INCLUDE_DIRS ${BOOST_SOURCE_DIR} CACHE STRING "Boost path" FORCE) set(Boost_LIBRARIES ${BOOST_SOURCE_DIR}/stage/lib/libboost_python.a CACHE STRING "Boost libraries" FORCE) + set(Boost_FOUND YES) ELSE (INTERNAL_BOOST) find_package(Boost 1.53) diff --git a/python/cosmotool/simu.py b/python/cosmotool/simu.py index 8aa7722..3dd4cf6 100644 --- a/python/cosmotool/simu.py +++ b/python/cosmotool/simu.py @@ -1,3 +1,4 @@ +import warnings from _cosmotool import * class SimulationBare(PySimulationBase): @@ -18,6 +19,10 @@ class SimulationBare(PySimulationBase): self.Hubble = s.getHubble() self.Omega_M = s.getOmega_M() self.Omega_Lambda = s.getOmega_Lambda() + try: + self.masses = s.getMasses().copy() if s.getMasses() is not None else None + except Exception as e: + warnings.warn("Unexpected exception: " + repr(e)) def merge(self, other): @@ -48,6 +53,11 @@ class SimulationBare(PySimulationBase): self.positions = _safe_merge(self.positions, other.getPositions()) self.velocities = _safe_merge(self.velocities, other.getVelocities()) self.identifiers = _safe_merge0(self.identifiers, other.getIdentifiers()) + try: + self.masses = _safe_merge0(self.masses, other.getMasses()) + except Exception as e: + warnings.warn("Unexpected exception: " + repr(e)); + self.masses = None def getPositions(self): return self.positions @@ -58,6 +68,9 @@ class SimulationBare(PySimulationBase): def getIdentifiers(self): return self.identifiers + def getMasses(self): + return self.masses + def getTime(self): return self.time diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py index f5653b3..86bc55d 100644 --- a/python_sample/icgen/borgicgen.py +++ b/python_sample/icgen/borgicgen.py @@ -93,7 +93,7 @@ def do_supergenerate(density, density_out=None, mulfac=None,zero_fill=False,Pk=N def build_Pk(): ik = np.fft.fftfreq(Ns, d=L/Ns)*2*np.pi k = ne.evaluate('sqrt(kx**2 + ky**2 + kz**2)', {'kx':ik[:,None,None], 'ky':ik[None,:,None], 'kz':ik[None,None,:(Ns/2+1)]}) - return Pk.compute(k)*(h*L)**3 + return Pk.compute(k)*L**3 print np.where(np.isnan(density_out))[0].size Nz = np.count_nonzero(cond) @@ -185,7 +185,7 @@ def whitify(density, L, cosmo, supergenerate=1, zero_fill=False, func='HU_WIGGLE def build_Pk(): ik = np.fft.fftfreq(N, d=L/N)*2*np.pi k = np.sqrt(ik[:,None,None]**2 + ik[None,:,None]**2 + ik[None,None,:(N/2+1)]**2) - return p.compute(k)*cosmo['h']**3*L**3 + return p.compute(k)*L**3 Pk = build_Pk() Pk[0,0,0]=1 diff --git a/sample/gadgetToArray.cpp b/sample/gadgetToArray.cpp index d91b95e..4995277 100644 --- a/sample/gadgetToArray.cpp +++ b/sample/gadgetToArray.cpp @@ -66,12 +66,15 @@ int main(int argc, char **argv) SimuData *p = loadGadgetMulti(fname, 0, 0); double L0 = p->BoxSize/MPC; + cout << "Will read " << p->TotalNumPart << " particles" << endl; array_type parts(boost::extents[p->TotalNumPart][7]); uint64_t q = 0; try { for (int cpuid=0;;cpuid++) { + cout << " = CPU " << cpuid << " = " << endl; p = loadGadgetMulti(fname, cpuid, NEED_POSITION|NEED_VELOCITY|NEED_MASS); + cout << " = DONE LOAD, COPYING IN PLACE" << endl; for (uint32_t i = 0; i < p->NumPart; i++) { for (int j = 0; j < 3; j++) @@ -87,11 +90,13 @@ int main(int argc, char **argv) parts[q][6] = p->Mass[i]; q++; } + cout << " = DONE (q=" << q << ")" << endl; delete p; } } catch (const NoSuchFileException& e) {} + cout << " ++ WRITING ++" << endl; hdf5_write_array(f, "particles", parts); return 0; diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp index 9073546..114ea09 100644 --- a/sample/simple3DFilter.cpp +++ b/sample/simple3DFilter.cpp @@ -123,10 +123,12 @@ int main(int argc, char **argv) MyCell *allCells_1 = new MyCell[N1_points]; +#pragma omp parallel for schedule(static) for (long i = 0; i < Nres*Nres*Nres; i++) bins.data()[i] = 0; cout << "Shuffling data in cells..." << endl; +#pragma omp parallel for schedule(static) for (int i = 0 ; i < N1_points; i++) { for (int j = 0; j < 3; j++) @@ -144,6 +146,7 @@ int main(int argc, char **argv) if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres) continue; +#pragma omp atomic update bins[rx][ry][rz]++; } v1_data.resize(boost::extents[1][1]); @@ -168,7 +171,7 @@ int main(int argc, char **argv) { double pz = (rz)*boxsize/Nres-cz; - cout << rz << " / " << Nres << endl; + (cout << rz << " / " << Nres << endl).flush(); for (int ry = 0; ry < Nres; ry++) { double py = (ry)*boxsize/Nres-cy; @@ -193,6 +196,7 @@ int main(int argc, char **argv) smooth1.addGridSite(c); } } + (cout << " Done " << rz << endl).flush(); } } @@ -204,7 +208,7 @@ int main(int argc, char **argv) bins, interpolated, getMass, rLimit2); hdf5_write_array(out_f, "density", interpolated); //out_f.flush(); - for (int i = 0; i < 0; i++) { + for (int i = 0; i < 3; i++) { computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz, bins, interpolated, boost::bind(getVelocity, _1, i), rLimit2); hdf5_write_array(out_f, str(format("p%d") % i), interpolated); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 659d4cc..a6f241b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -3,7 +3,6 @@ SET(CosmoTool_SRCS fortran.cpp interpolate.cpp load_data.cpp - loadGadget.cpp loadRamses.cpp powerSpectrum.cpp miniargs.cpp @@ -12,6 +11,14 @@ SET(CosmoTool_SRCS cic.cpp ) +IF (Boost_FOUND) + include_directories(${Boost_INCLUDE_DIRS}) + SET(CosmoTool_SRCS + ${CosmoTool_SRCS} + loadGadget.cpp + ) +ENDIF (Boost_FOUND) + IF (ENABLE_OPENMP) ENDIF (ENABLE_OPENMP) diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp index 3aa2854..9848597 100644 --- a/src/cosmopower.cpp +++ b/src/cosmopower.cpp @@ -99,7 +99,7 @@ static double powC(double q, double alpha_c) static double T_tilde_0(double q, double alpha_c, double beta_c) { - static const double c_E = 2.718282; //M_E; + static const double c_E = M_E; double a = log(c_E + 1.8 * beta_c * q); return a / ( a + powC(q, alpha_c) * q * q); } diff --git a/src/fortran.cpp b/src/fortran.cpp index f0b9c86..367bf40 100644 --- a/src/fortran.cpp +++ b/src/fortran.cpp @@ -112,7 +112,7 @@ void UnformattedRead::beginCheckpoint() checkPointRef = (cSize == Check_32bits) ? 4 : 8; checkPointAccum = 0; - checkPointRef = (cSize == Check_32bits) ? readInt32() : readInt64(); + checkPointRef = (cSize == Check_32bits) ? readUint32() : readInt64(); checkPointAccum = 0; if (f->eof()) @@ -144,7 +144,7 @@ void UnformattedRead::endCheckpoint(bool autodrop) void UnformattedRead::readOrderedBuffer(void *buffer, int size) throw (InvalidUnformattedAccess) { - if ((checkPointAccum+size) > checkPointRef) + if ((checkPointAccum+(uint64_t)size) > checkPointRef) throw InvalidUnformattedAccess(); f->read((char *)buffer, size); @@ -186,6 +186,20 @@ float UnformattedRead::readReal32() return a.f; } +uint32_t UnformattedRead::readUint32() + throw (InvalidUnformattedAccess) +{ + union + { + char b[4]; + uint32_t i; + } a; + + readOrderedBuffer(&a, 4); + + return a.i; +} + int32_t UnformattedRead::readInt32() throw (InvalidUnformattedAccess) { diff --git a/src/fortran.hpp b/src/fortran.hpp index 85e7618..81414ad 100644 --- a/src/fortran.hpp +++ b/src/fortran.hpp @@ -73,6 +73,8 @@ namespace CosmoTool // Todo implement primitive description void setOrdering(Ordering o); void setCheckpointSize(CheckpointSize cs); + + uint64_t getBlockSize() const { return checkPointRef; } void beginCheckpoint() throw (InvalidUnformattedAccess,EndOfFileException); @@ -83,6 +85,8 @@ namespace CosmoTool throw (InvalidUnformattedAccess); float readReal32() throw (InvalidUnformattedAccess); + uint32_t readUint32() + throw (InvalidUnformattedAccess); int32_t readInt32() throw (InvalidUnformattedAccess); int64_t readInt64() diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index 61bb922..e307dbf 100644 --- a/src/loadGadget.cpp +++ b/src/loadGadget.cpp @@ -146,10 +146,12 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, } if (h.flag_doubleprecision) { + cout << "Gadget format with double precision" << endl; readToDouble = boost::bind(myRead64, f); readToSingle = boost::bind(myRead64, f); float_size = sizeof(double); } else { + cout << "Gadget format with single precision" << endl; readToDouble = boost::bind(myRead32, f); readToSingle = boost::bind(myRead32, f); float_size = sizeof(float); @@ -188,8 +190,19 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, try { f->beginCheckpoint(); + if (f->getBlockSize() != NumPart*float_size*3) { + // Check that single would work + if (f->getBlockSize() == NumPart*sizeof(float)*3) { + // Change to single + cout << "Change back to single. Buggy header." << endl; + readToDouble = boost::bind(myRead32, f); + readToSingle = boost::bind(myRead32, f); + float_size = sizeof(float); + } + } for(int k = 0, p = 0; k < 6; k++) { for(int n = 0; n < h.npart[k]; n++) { +// if ((n%1000000)==0) cout << n << endl; data->Pos[0][p] = readToSingle(); data->Pos[1][p] = readToSingle(); data->Pos[2][p] = readToSingle(); @@ -230,6 +243,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, for(int k = 0, p = 0; k < 6; k++) { for(int n = 0; n < h.npart[k]; n++) { // THIS IS GADGET 1 +// if ((n%1000000)==0) cout << n << endl; data->Vel[0][p] = readToSingle()*velmul; data->Vel[1][p] = readToSingle()*velmul; data->Vel[2][p] = readToSingle()*velmul; @@ -292,6 +306,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, if (loadflags & NEED_MASS) { bool do_load = false; + cout << "=> Mass" << endl; for (int k = 0; k < 6; k++) { do_load = do_load || ((h.mass[k] == 0)&&(h.npart[k]>0)); @@ -308,11 +323,13 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, if (h.mass[k] == 0) { for(int n = 0; n < h.npart[k]; n++) { +// if ((n%1000000)==0) cout << n << endl; data->Mass[l++] = readToSingle(); } } else { for(int n = 0; n < h.npart[k]; n++) { + if ((n%1000000)==0) cout << n << endl; data->Mass[l++] = h.mass[k]; } } diff --git a/src/load_data.cpp b/src/load_data.cpp index 9fd2552..da7ca28 100644 --- a/src/load_data.cpp +++ b/src/load_data.cpp @@ -41,7 +41,7 @@ knowledge of the CeCILL license and that you accept its terms. using namespace CosmoTool; //#define LARGE_CONTROL -#define LITTLE_ENDIAN +//#define LITTLE_ENDIAN #define NEW(t,n) ((t *)malloc(sizeof(t)*n)) #define SKIP(f) fread(&dummy,sizeof(dummy),1,f); diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp index bcd98a5..c364621 100644 --- a/src/sphSmooth.hpp +++ b/src/sphSmooth.hpp @@ -60,9 +60,10 @@ namespace CosmoTool typedef void (*runParticleValue)(ValType& t); public: + typedef SPHCell *P_SPHCell; struct SPHState { - boost::shared_ptr ngb; + boost::shared_ptr ngb; boost::shared_ptr distances; typename SPHTree::coords currentCenter; int currentNgb; diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc index 99094cd..45a3c29 100644 --- a/src/sphSmooth.tcc +++ b/src/sphSmooth.tcc @@ -48,15 +48,15 @@ SPHSmooth::fetchNeighbours(const typename SPHTree::coords& c, uin if (requested > maxNgb) { maxNgb = requested; - internal.ngb = boost::shared_ptr(new SPHCell *[maxNgb]); + internal.ngb = boost::shared_ptr(new P_SPHCell[maxNgb]); internal.distances = boost::shared_ptr(new CoordType[maxNgb]); } memcpy(internal.currentCenter, c, sizeof(c)); - tree->getNearestNeighbours(c, requested, internal.ngb.get(), internal.distances.get()); + tree->getNearestNeighbours(c, requested, (SPHCell **)internal.ngb.get(), (CoordType*)internal.distances.get()); internal.currentNgb = 0; - for (uint32_t i = 0; i < requested && internal.ngb[i] != 0; i++,internal.currentNgb++) + for (uint32_t i = 0; i < requested && (internal.ngb)[i] != 0; i++,internal.currentNgb++) { internal.distances[i] = sqrt(internal.distances[i]); d2 = internal.distances[i]; @@ -191,9 +191,8 @@ void SPHSmooth::addGridSite(const typename SPHTree::coords& c) for (uint32_t i = 0; i < internal.currentNgb; i++) { ComputePrecision d = internal.distances[i]; - SPHCell& cell = *internal.ngb[i]; + SPHCell& cell = *(internal.ngb[i]); cell.val.weight += getKernel(d/internal.smoothRadius) / r3; - } }