diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py index 28f597e..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) 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/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/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; - } }