Merge branch 'master' of bitbucket.org:glavaux/cosmotool

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)

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)

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
     

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

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;

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);

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)
 

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);
 }

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)
 {

src/fortran.hpp

@@ -74,6 +74,8 @@ namespace CosmoTool
     void setOrdering(Ordering o);
     void setCheckpointSize(CheckpointSize cs);
 
+    uint64_t getBlockSize() const { return checkPointRef; }
+    
     void beginCheckpoint()
       throw (InvalidUnformattedAccess,EndOfFileException);
     void endCheckpoint(bool autodrop = false)
@@ -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()

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<double>, f);
         readToSingle = boost::bind(myRead64<float>, f);
         float_size = sizeof(double);
       } else {
+        cout << "Gadget format with single precision" << endl;
         readToDouble = boost::bind(myRead32<double>, f);
         readToSingle = boost::bind(myRead32<float>, 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<double>, f);
+            readToSingle = boost::bind(myRead32<float>, 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];
                 }
             }

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);

src/sphSmooth.hpp

@@ -60,9 +60,10 @@ namespace CosmoTool
     typedef void (*runParticleValue)(ValType& t);
 
   public:
+    typedef SPHCell *P_SPHCell;
     struct SPHState
     {
-      boost::shared_ptr<SPHCell *[]> ngb;
+      boost::shared_ptr<P_SPHCell[]> ngb;
       boost::shared_ptr<CoordType[]> distances;
       typename SPHTree::coords currentCenter;
       int currentNgb;

src/sphSmooth.tcc

@@ -48,15 +48,15 @@ SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uin
   if (requested > maxNgb)
     {
       maxNgb = requested;
-      internal.ngb = boost::shared_ptr<SPHCell*[]>(new SPHCell *[maxNgb]);
+      internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]);
       internal.distances = boost::shared_ptr<CoordType[]>(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<ValType,Ndims>::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;
-
     }
 }