Fixes for SPH. Workaround for buggy gadget headers. Fixed borg generation for newer API.

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)

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