From 0b40de23d09610bf98e576c14ffddfd133974726 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 26 Feb 2015 10:36:52 +0100 Subject: [PATCH] Added some support for double precision gadget snapshot --- python/_cosmo_power.pyx | 6 ++-- src/cosmopower.cpp | 17 ++++++---- src/cosmopower.hpp | 2 +- src/loadGadget.cpp | 72 ++++++++++++++++++++++++++++------------- src/loadSimu.hpp | 66 ++++++++++++++++++------------------- src/load_data.hpp | 6 +++- 6 files changed, 102 insertions(+), 67 deletions(-) diff --git a/python/_cosmo_power.pyx b/python/_cosmo_power.pyx index 2e87c5c..701788f 100644 --- a/python/_cosmo_power.pyx +++ b/python/_cosmo_power.pyx @@ -44,7 +44,7 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool": void setFunction(CosmoFunction) void updateCosmology() void updatePhysicalCosmology() - void normalize(double) + void normalize(double,double) void setNormalization(double) double power(double) @@ -75,7 +75,7 @@ cdef class CosmologyPower: self.power.updateCosmology() - def normalize(self,s8,k_max=-1): + def normalize(self,s8,k_min=-1,k_max=-1): """normalize(self, sigma8) Compute the normalization of the power spectrum using sigma8. @@ -84,7 +84,7 @@ cdef class CosmologyPower: sigma8 (float): standard deviation of density field smoothed at 8 Mpc/h """ self.power.SIGMA8 = s8 - self.power.normalize(k_max) + self.power.normalize(k_min, k_max) def setFunction(self,funcname): diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp index 5104148..d81cb85 100644 --- a/src/cosmopower.cpp +++ b/src/cosmopower.cpp @@ -123,6 +123,9 @@ double CosmoPower::powerEfstathiou(double k) void CosmoPower::updateHuWigglesConsts() { + double f_b = OMEGA_B / OMEGA_0; + double f_c = OMEGA_C / OMEGA_0; + double k_silk = 1.6 * pow(OMEGA_B * h * h, 0.52) * pow(OmegaEff, 0.73) * (1 + pow(10.4 * OmegaEff, -0.95)); double z_eq = 2.50e4 * OmegaEff * pow(Theta_27, -4); //double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75))); @@ -139,14 +142,14 @@ void CosmoPower::updateHuWigglesConsts() double a1 = pow(46.9 * OmegaEff, 0.670) * (1 + pow(32.1 * OmegaEff, -0.532)); double a2 = pow(12.0 * OmegaEff, 0.424) * (1 + pow(45.0 * OmegaEff, -0.582)); - double alpha_c = pow(a1, -OMEGA_B/ OMEGA_0) * pow(a2, -pow(OMEGA_B / OMEGA_0, 3)); + double alpha_c = pow(a1, -f_b) * pow(a2, -pow(f_b, 3)); double b1_betac = 0.944 * 1/(1 + pow(458 * OmegaEff, -0.708)); double b2_betac = pow(0.395 * OmegaEff, -0.0266); - double beta_c = 1/ ( 1 + b1_betac * (pow(OMEGA_C / OMEGA_0, b2_betac) - 1) ); + double beta_c = 1/ ( 1 + b1_betac * (pow(f_c, b2_betac) - 1) ); double alpha_b = 2.07 * k_eq * s * pow(1 + R_d, -0.75) * powG((1 + z_eq)/(1 + z_d)); - double beta_b = 0.5 + OMEGA_B / OMEGA_0 + (3 - 2 * OMEGA_B / OMEGA_0) * sqrt(pow(17.2 * OmegaEff, 2) + 1); + double beta_b = 0.5 + f_b + (3 - 2 * f_b) * sqrt(pow(17.2 * OmegaEff, 2) + 1); double beta_node = 8.41 * pow(OmegaEff, 0.435); ehu.k_silk = k_silk; @@ -262,16 +265,18 @@ double CosmoPower::integrandNormalize(double x) return power(k)*k*k*f*f/(x*x); } -void CosmoPower::normalize(double k_max) +void CosmoPower::normalize(double k_min, double k_max) { double normVal = 0; double abserr; gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); gsl_function f; - double x_min = 0; + double x_min = 0, x_max = 1; if (k_max > 0) x_min = 1/(1+k_max); + if (k_min > 0) + x_max = 1/(1+k_min); f.function = gslPowSpecNorm; f.params = this; @@ -286,7 +291,7 @@ void CosmoPower::normalize(double k_max) } // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr); - gsl_integration_qag(&f, x_min, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr); + gsl_integration_qag(&f, x_min, x_max, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr); gsl_integration_workspace_free(w); normVal /= (2*M_PI*M_PI); diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp index f1ba45d..b32c6a7 100644 --- a/src/cosmopower.hpp +++ b/src/cosmopower.hpp @@ -92,7 +92,7 @@ namespace CosmoTool { void updateCosmology(); void updatePhysicalCosmology(); - void normalize(double k_max = -1); + void normalize(double k_min = -1, double k_max = -1); void setNormalization(double A_K); void updateHuWigglesConsts(); diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index cbd458d..61bb922 100644 --- a/src/loadGadget.cpp +++ b/src/loadGadget.cpp @@ -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,6 +67,13 @@ 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); ssize_t NumPart = 0, NumPartTotal = 0; @@ -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,9 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, UnformattedRead *f; GadgetHeader h; float velmul; + boost::function0 readToDouble; + boost::function0 readToSingle; + long float_size; if (id >= 0) { int k = snprintf(0, 0, "%s.%d", fname, id)+1; @@ -124,9 +142,18 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, velmul = 1/(h.time); else { cerr << "unknown gadget format" << endl; - abort(); + abort(); + } + + if (h.flag_doubleprecision) { + readToDouble = boost::bind(myRead64, f); + readToSingle = boost::bind(myRead64, f); + float_size = sizeof(double); + } else { + readToDouble = boost::bind(myRead32, f); + readToSingle = boost::bind(myRead32, f); + float_size = sizeof(float); } - NumPart = data->NumPart; NumPartTotal = data->TotalNumPart; } @@ -144,28 +171,28 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, 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 { f->beginCheckpoint(); 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(); + data->Pos[0][p] = readToSingle(); + data->Pos[1][p] = readToSingle(); + data->Pos[2][p] = readToSingle(); p++; } } @@ -173,15 +200,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) { @@ -202,9 +230,9 @@ 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 - data->Vel[0][p] = f->readReal32()*velmul; - data->Vel[1][p] = f->readReal32()*velmul; - data->Vel[2][p] = f->readReal32()*velmul; + data->Vel[0][p] = readToSingle()*velmul; + data->Vel[1][p] = readToSingle()*velmul; + data->Vel[2][p] = readToSingle()*velmul; p++; } } @@ -222,7 +250,7 @@ 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 @@ -280,7 +308,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, if (h.mass[k] == 0) { for(int n = 0; n < h.npart[k]; n++) { - data->Mass[l++] = f->readReal32(); + data->Mass[l++] = readToSingle(); } } else { for(int n = 0; n < h.npart[k]; n++) diff --git a/src/loadSimu.hpp b/src/loadSimu.hpp index b098351..6f848ce 100644 --- a/src/loadSimu.hpp +++ b/src/loadSimu.hpp @@ -46,6 +46,7 @@ namespace CosmoTool 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 { @@ -88,49 +89,46 @@ namespace CosmoTool 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() { - 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); - } + 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/src/load_data.hpp b/src/load_data.hpp index 6643eaa..36a2422 100644 --- a/src/load_data.hpp +++ b/src/load_data.hpp @@ -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