From 555e0a4d5a0e63b672f4f270993b127b62e06bc5 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 26 Jan 2016 23:07:47 +0100 Subject: [PATCH] Imported newer version of some cosmotool files --- external/cosmotool/src/fortran.cpp | 30 +++- external/cosmotool/src/fortran.hpp | 13 +- external/cosmotool/src/loadGadget.cpp | 211 +++++++++++++++++++++----- external/cosmotool/src/loadGadget.hpp | 2 +- external/cosmotool/src/loadSimu.hpp | 76 +++++----- external/cosmotool/src/load_data.hpp | 8 +- 6 files changed, 260 insertions(+), 80 deletions(-) diff --git a/external/cosmotool/src/fortran.cpp b/external/cosmotool/src/fortran.cpp index 79641ac..e497ff1 100644 --- a/external/cosmotool/src/fortran.cpp +++ b/external/cosmotool/src/fortran.cpp @@ -1,5 +1,5 @@ /*+ -This is CosmoTool (./src/fortran.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) +This is CosmoTool (./src/fortran.cpp) -- Copyright (C) Guilhem Lavaux (2007-2014) guilhem.lavaux@gmail.com @@ -71,6 +71,16 @@ UnformattedRead::~UnformattedRead() delete f; } +int64_t UnformattedRead::position() const +{ + return f->tellg(); +} + +void UnformattedRead::seek(int64_t pos) +{ + f->seekg(pos, istream::beg); +} + // Todo implement primitive description void UnformattedRead::setOrdering(Ordering o) { @@ -112,7 +122,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 +154,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 +196,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/external/cosmotool/src/fortran.hpp b/external/cosmotool/src/fortran.hpp index d63e8d3..48b4a1d 100644 --- a/external/cosmotool/src/fortran.hpp +++ b/external/cosmotool/src/fortran.hpp @@ -1,5 +1,5 @@ /*+ -This is CosmoTool (./src/fortran.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) +This is CosmoTool (./src/fortran.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014) guilhem.lavaux@gmail.com @@ -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() @@ -91,6 +95,11 @@ namespace CosmoTool void skip(int64_t off) throw (InvalidUnformattedAccess); + int64_t position() const; + void seek(int64_t pos); + + void readOrderedBuffer(void *buffer, int size) + throw (InvalidUnformattedAccess); protected: bool swapOrdering; CheckpointSize cSize; @@ -98,8 +107,6 @@ namespace CosmoTool uint64_t checkPointAccum; std::ifstream *f; - void readOrderedBuffer(void *buffer, int size) - throw (InvalidUnformattedAccess); }; class UnformattedWrite: public FortranTypes diff --git a/external/cosmotool/src/loadGadget.cpp b/external/cosmotool/src/loadGadget.cpp index f9f1072..9f77bb4 100644 --- a/external/cosmotool/src/loadGadget.cpp +++ b/external/cosmotool/src/loadGadget.cpp @@ -1,5 +1,5 @@ /*+ -This is CosmoTool (./src/loadGadget.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) +This is CosmoTool (./src/loadGadget.cpp) -- Copyright (C) Guilhem Lavaux (2007-2014) guilhem.lavaux@gmail.com @@ -38,6 +38,8 @@ knowledge of the CeCILL license and that you accept its terms. #include #include #include +#include +#include #include "load_data.hpp" #include "loadGadget.hpp" #include "fortran.hpp" @@ -65,9 +67,16 @@ void loadGadgetHeader(UnformattedRead *f, GadgetHeader& h, SimuData *data, int i data->Omega_M = h.Omega0 = f->readReal64(); data->Omega_Lambda = h.OmegaLambda = f->readReal64(); data->Hubble = h.HubbleParam = f->readReal64(); + (int)f->readInt32(); // stellarage + (int)f->readInt32(); // metals + for (int i = 0; i < 6; i++) + h.npartTotal[i] |= ((unsigned long)f->readInt32()) << 32; + (int)f->readInt32(); // entropy instead of u + h.flag_doubleprecision = f->readInt32(); + f->endCheckpoint(true); - long NumPart = 0, NumPartTotal = 0; + ssize_t NumPart = 0, NumPartTotal = 0; for(int k=0; k<6; k++) { NumPart += h.npart[k]; @@ -77,6 +86,12 @@ void loadGadgetHeader(UnformattedRead *f, GadgetHeader& h, SimuData *data, int i data->TotalNumPart = NumPartTotal; } +template +T myRead64(UnformattedRead *f) { return f->readReal64(); } + +template +T myRead32(UnformattedRead *f) { return f->readReal32(); } + SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, int GadgetFormat, SimuFilter filter) @@ -86,6 +101,14 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, UnformattedRead *f; GadgetHeader h; float velmul; + boost::function0 readToDouble; + boost::function0 readToSingle; + long float_size; + + if (GadgetFormat > 2) { + cerr << "Unknown gadget format" << endl; + return 0; + } if (id >= 0) { int k = snprintf(0, 0, "%s.%d", fname, id)+1; @@ -96,7 +119,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, if (f == 0) return 0; - delete out_fname; + delete[] out_fname; } else { @@ -106,27 +129,68 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, } + typedef std::map BlockMap; + BlockMap blockTable; + + if (GadgetFormat == 2) { + int64_t startBlock = 0; + char block[5]; + int32_t blocksize; + + try { + while (true) { + f->beginCheckpoint(); + f->readOrderedBuffer(block, 4); + block[4] = 0; + blocksize = f->readInt32(); + f->endCheckpoint(); + blockTable[block] = f->position(); + f->skip(blocksize); + } + } catch (EndOfFileException&) {} + + f->seek(startBlock); + } + data = new SimuData; if (data == 0) { delete f; return 0; } - long NumPart = 0, NumPartTotal = 0; + ssize_t NumPart = 0, NumPartTotal = 0; +#define ENSURE(name) { \ + if (GadgetFormat == 2) { \ + BlockMap::iterator iter = blockTable.find(name); \ + if (iter == blockTable.end()) { \ + std::cerr << "GADGET2: Cannot find block named '" << name << "'" << endl; \ + if (data) delete data; \ + delete f; \ + return 0; \ + } \ + f->seek(iter->second); \ + } \ + } + try { + ENSURE("HEAD"); loadGadgetHeader(f, h, data, id); - if (GadgetFormat == 1) - velmul = sqrt(h.time); - else if (GadgetFormat == 2) - velmul = 1/(h.time); - else { - cerr << "unknown gadget format" << endl; - abort(); + velmul = sqrt(h.time); + + 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); } - NumPart = data->NumPart; NumPartTotal = data->TotalNumPart; } @@ -138,34 +202,47 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, return 0; } + if (loadflags & NEED_TYPE) { int p = 0; - + data->type = new int[data->NumPart]; for (int k = 0; k < 6; k++) - for (int n = 0; n < h.npart[k]; n++,p++) - data->type[p] = k; + for (int n = 0; n < h.npart[k]; n++,p++) + data->type[p] = k; } if (loadflags & NEED_POSITION) { for (int i = 0; i < 3; i++) { - data->Pos[i] = new float[data->NumPart]; - if (data->Pos[i] == 0) { - delete data; - return 0; - } + data->Pos[i] = new float[data->NumPart]; + if (data->Pos[i] == 0) { + delete data; + return 0; } + } try { + ENSURE("POS "); 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++) { - data->Pos[0][p] = f->readReal32(); - data->Pos[1][p] = f->readReal32(); - data->Pos[2][p] = f->readReal32(); +// if ((n%1000000)==0) cout << n << endl; + data->Pos[0][p] = readToSingle(); + data->Pos[1][p] = readToSingle(); + data->Pos[2][p] = readToSingle(); p++; } } @@ -173,15 +250,16 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, } catch (const InvalidUnformattedAccess& e) { - cerr << "Invalid format while reading positions" << endl; - delete f; - delete data; - return 0; + cerr << "Invalid format while reading positions" << endl; + delete f; + delete data; + return 0; } } else { + long float_size = (h.flag_doubleprecision) ? sizeof(double) : sizeof(float); // Skip positions - f->skip(NumPart * 3 * sizeof(float) + 2*4); + f->skip(NumPart * 3 * float_size + 2*4); } if (loadflags & NEED_VELOCITY) { @@ -198,13 +276,15 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, try { + ENSURE("VEL "); f->beginCheckpoint(); for(int k = 0, p = 0; k < 6; k++) { for(int n = 0; n < h.npart[k]; n++) { // THIS IS GADGET 1 - data->Vel[0][p] = f->readReal32()*velmul; - data->Vel[1][p] = f->readReal32()*velmul; - data->Vel[2][p] = f->readReal32()*velmul; +// 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; p++; } } @@ -222,13 +302,14 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, /// // TODO: FIX THE UNITS OF THESE FUNKY VELOCITIES !!! } else { // Skip velocities - f->skip(NumPart*3*sizeof(float)+2*4); + f->skip(NumPart*3*float_size+2*4); } // Skip ids if (loadflags & NEED_GADGET_ID) { try { + ENSURE("ID "); f->beginCheckpoint(); data->Id = new long[data->NumPart]; if (data->Id == 0) @@ -261,11 +342,66 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, f->skip(h.npart[k]*4); } + if (loadflags & NEED_MASS) { + bool do_load = false; + + for (int k = 0; k < 6; k++) + { + do_load = do_load || ((h.mass[k] == 0)&&(h.npart[k]>0)); + } + + try + { + long l = 0; + if (do_load) { + ENSURE("MASS"); + f->beginCheckpoint(); + } + data->Mass = new float[NumPart]; + for (int k = 0; k < 6; k++) + { + 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]; + } + } + } + if (do_load) + f->endCheckpoint(); + } + catch (const InvalidUnformattedAccess& e) + { + cerr << "Invalid unformatted access while reading ID" << endl; + delete f; + delete data; + return 0; + } + catch (const EndOfFileException& e) + { + for (int k = 0; k < 6; k++) + cerr << "mass[" << k << "] = " << h.mass[k] << endl; + } + } else { + f->skip(2*4); + for (int k = 0; k < 6; k++) + if (h.mass[k] == 0) + f->skip(h.npart[k]*4); + } + delete f; return data; } +#undef ENSURE void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int GadgetFormat) @@ -274,12 +410,16 @@ void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int Gadget int npart[6], npartTotal[6]; float mass[6]; - if (data->Pos[0] == 0 || data->Vel[0] == 0 || data->Id == 0) + if (data->Pos[0] == 0 || data->Vel[0] == 0 || data->Id == 0) { + cerr << "Invalid simulation data array" << endl; return; + } f = new UnformattedWrite(fname); - if (f == 0) + if (f == 0) { + cerr << "Cannot create file" << endl; return; + } for (int i = 0; i < 6; i++) { @@ -321,8 +461,7 @@ void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int Gadget f->endCheckpoint(); float velmul = 1.0; - if (GadgetFormat == 1) - velmul = sqrt(data->time); + velmul = sqrt(data->time); f->beginCheckpoint(); for(int n = 0; n < data->NumPart; n++) { diff --git a/external/cosmotool/src/loadGadget.hpp b/external/cosmotool/src/loadGadget.hpp index e05a187..adbc818 100644 --- a/external/cosmotool/src/loadGadget.hpp +++ b/external/cosmotool/src/loadGadget.hpp @@ -1,5 +1,5 @@ /*+ -This is CosmoTool (./src/loadGadget.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) +This is CosmoTool (./src/loadGadget.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014) guilhem.lavaux@gmail.com diff --git a/external/cosmotool/src/loadSimu.hpp b/external/cosmotool/src/loadSimu.hpp index 1c360a4..2eb5ac8 100644 --- a/external/cosmotool/src/loadSimu.hpp +++ b/external/cosmotool/src/loadSimu.hpp @@ -1,5 +1,5 @@ /*+ -This is CosmoTool (./src/loadSimu.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) +This is CosmoTool (./src/loadSimu.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014) guilhem.lavaux@gmail.com @@ -36,6 +36,7 @@ knowledge of the CeCILL license and that you accept its terms. #ifndef __COSMOTOOLBOX_HPP #define __COSMOTOOLBOX_HPP +#include #include #include @@ -45,6 +46,8 @@ namespace CosmoTool static const int NEED_POSITION = 2; static const int NEED_VELOCITY = 4; static const int NEED_TYPE = 8; + static const int NEED_MASS = 16; + static const int NEED_DOUBLE_PRECISION = 32; struct SimuParticle { @@ -64,6 +67,8 @@ namespace CosmoTool typedef void (*FreeFunction)(void *); typedef std::map > AttributeMap; + bool noAuto; + float BoxSize; float time; float Hubble; @@ -71,59 +76,60 @@ namespace CosmoTool float Omega_M; float Omega_Lambda; - long NumPart; - long TotalNumPart; + ssize_t NumPart; + ssize_t TotalNumPart; long *Id; float *Pos[3]; float *Vel[3]; + float *Mass; int *type; AttributeMap attributes; public: - SimuData() : Id(0),NumPart(0),type(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; } + SimuData() : Mass(0), Id(0),NumPart(0),type(0),noAuto(false) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; } ~SimuData() { - for (int j = 0; j < 3; j++) - { - if (Pos[j]) - delete[] Pos[j]; - if (Vel[j]) - delete[] Vel[j]; - } - if (type) - delete[] type; - if (Id) - delete[] Id; - - for (AttributeMap::iterator i = attributes.begin(); - i != attributes.end(); - ++i) - { - if (i->second.second) - i->second.second(i->second.first); - } + if (!noAuto) { + for (int j = 0; j < 3; j++) { + if (Pos[j]) + delete[] Pos[j]; + if (Vel[j]) + delete[] Vel[j]; + } + if (type) + delete[] type; + if (Id) + delete[] Id; + if (Mass) + delete[] Mass; + } + for (AttributeMap::iterator i = attributes.begin(); + i != attributes.end(); + ++i) { + if (i->second.second) + i->second.second(i->second.first); + } } template T *as(const std::string& n) { - AttributeMap::iterator i = attributes.find(n); - if (i == attributes.end()) - return 0; - - return reinterpret_cast(i->second.first); + AttributeMap::iterator i = attributes.find(n); + if (i == attributes.end()) + return 0; + + return reinterpret_cast(i->second.first); } void new_attribute(const std::string& n, void *p, FreeFunction free_func) { - AttributeMap::iterator i = attributes.find(n); - if (i != attributes.end()) - { - if (i->second.second) - i->second.second(i->second.first); - } - attributes[n] = std::make_pair(p, free_func); + AttributeMap::iterator i = attributes.find(n); + if (i != attributes.end()) { + if (i->second.second) + i->second.second(i->second.first); + } + attributes[n] = std::make_pair(p, free_func); } }; diff --git a/external/cosmotool/src/load_data.hpp b/external/cosmotool/src/load_data.hpp index 6be3bbb..36a2422 100644 --- a/external/cosmotool/src/load_data.hpp +++ b/external/cosmotool/src/load_data.hpp @@ -1,5 +1,5 @@ /*+ -This is CosmoTool (./src/load_data.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) +This is CosmoTool (./src/load_data.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014) guilhem.lavaux@gmail.com @@ -55,7 +55,11 @@ namespace CosmoTool { double Omega0; double OmegaLambda; double HubbleParam; - char fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8]; /* fills to 256 Bytes */ + int flag_doubleprecision; + int flag_ic_info; + float lpt_scalingfactor; + char fill[18]; /*!< fills to 256 Bytes */ + char names[15][2]; }; struct ParticleState