Added some support for double precision gadget snapshot
This commit is contained in:
Guilhem Lavaux
parent
b7392204ad
commit
0b40de23d0
@ -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):
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -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();
|
||||
|
||||
|
||||
@ -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,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<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,9 @@ 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 (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<double>, f);
|
||||
readToSingle = boost::bind(myRead64<float>, f);
|
||||
float_size = sizeof(double);
|
||||
} else {
|
||||
readToDouble = boost::bind(myRead32<double>, f);
|
||||
readToSingle = boost::bind(myRead32<float>, 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++)
|
||||
|
||||
@ -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<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);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user