guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 15:21:11 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Imported newer version of some cosmotool files

Browse source
This commit is contained in:
Guilhem Lavaux 2016-01-26 23:07:47 +01:00
parent 90fd74d907
commit 555e0a4d5a
6 changed files with 260 additions and 80 deletions
Download patch file Download diff file Expand all files Collapse all files

30
external/cosmotool/src/fortran.cpp vendored
View file

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

13
external/cosmotool/src/fortran.hpp vendored
View file

@ -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

211
external/cosmotool/src/loadGadget.cpp vendored
View file

@ -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 <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <boost/function.hpp>
#include <boost/bind.hpp>
#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<typename T>
T myRead64(UnformattedRead *f) { return f->readReal64(); }
template<typename T>
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<double> readToDouble;
boost::function0<float> 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<std::string, int64_t> 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<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);
}
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<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++) {
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++) {

2
external/cosmotool/src/loadGadget.hpp vendored
View file

@ -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

76
external/cosmotool/src/loadSimu.hpp vendored
View file

@ -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 <sys/types.h>
#include <map>
#include <string>
@ -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<std::string, std::pair<void *, FreeFunction> > 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<typename T>
T *as(const std::string& n)
{
AttributeMap::iterator i = attributes.find(n);
if (i == attributes.end())
return 0;
return reinterpret_cast<T *>(i->second.first);
AttributeMap::iterator i = attributes.find(n);
if (i == attributes.end())
return 0;
return reinterpret_cast<T *>(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);
}
};

8
external/cosmotool/src/load_data.hpp vendored
View file

@ -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