2009-01-08 16:19:01 +01:00
|
|
|
#include <regex.h>
|
2008-12-02 18:00:28 +01:00
|
|
|
#include <cmath>
|
|
|
|
#include <cstring>
|
|
|
|
#include <cstdlib>
|
|
|
|
#include <cassert>
|
|
|
|
#include <sstream>
|
|
|
|
#include <iostream>
|
|
|
|
#include <iomanip>
|
|
|
|
#include <string>
|
|
|
|
#include <fstream>
|
|
|
|
#include "loadRamses.hpp"
|
|
|
|
#include "load_data.h"
|
2009-01-08 16:19:01 +01:00
|
|
|
#include <map>
|
2008-12-02 18:00:28 +01:00
|
|
|
|
|
|
|
using namespace std;
|
|
|
|
|
|
|
|
#if 0
|
|
|
|
#define read_buf(b, n) \
|
|
|
|
{ \
|
|
|
|
int k; \
|
|
|
|
control_size -= n; \
|
|
|
|
for (k = (n-1); k >= 0; k--) \
|
|
|
|
infile.read(&b[k],1); \
|
|
|
|
}
|
|
|
|
#else
|
|
|
|
#define read_buf(b, n) \
|
|
|
|
{ \
|
|
|
|
int k; \
|
|
|
|
control_size -= n; \
|
|
|
|
for (k = 0; k < n; k++) \
|
|
|
|
infile.read(&b[k], 1); \
|
|
|
|
}
|
|
|
|
#endif
|
|
|
|
|
|
|
|
#define read_int(i) \
|
|
|
|
{ \
|
|
|
|
char *o = (char*)&(i); \
|
|
|
|
read_buf(o, 4); \
|
|
|
|
}
|
|
|
|
|
|
|
|
#define read_real(f) \
|
|
|
|
{ \
|
|
|
|
char *o = (char*)&(f); \
|
|
|
|
read_buf(o, 4); \
|
|
|
|
}
|
|
|
|
|
|
|
|
#define read_characters(c, n) { \
|
|
|
|
int k; \
|
|
|
|
control_size -= n; \
|
|
|
|
infile.read(c, n); \
|
|
|
|
}
|
|
|
|
|
|
|
|
#define push_dummy_control(id) \
|
|
|
|
{ int control_size = 0;
|
|
|
|
|
|
|
|
#define pop_dummy_control() }
|
|
|
|
|
|
|
|
#define push_control(id) \
|
|
|
|
{ \
|
|
|
|
int control_size = 0; \
|
|
|
|
int control_size2 = 0; \
|
|
|
|
char *intbuf = (char*)&control_size; \
|
|
|
|
read_int(control_size);
|
|
|
|
|
|
|
|
#define pop_control(id) \
|
|
|
|
assert(control_size == 0); \
|
|
|
|
read_int(control_size2); \
|
|
|
|
}
|
|
|
|
|
|
|
|
GadgetData *CosmoTool::loadRamses(const char *name)
|
|
|
|
{
|
|
|
|
GadgetData *gd = (GadgetData *)malloc(sizeof(GadgetData));
|
|
|
|
int id = 1;
|
|
|
|
uint32_t totPart = 0;
|
|
|
|
|
|
|
|
cout << "Detecting number of files and particles..." << endl;
|
|
|
|
while (1)
|
|
|
|
{
|
|
|
|
ostringstream ss_fname;
|
|
|
|
ss_fname << name << setfill('0') << setw(5) << id;
|
|
|
|
string fname = ss_fname.str();
|
|
|
|
|
|
|
|
cout << " ... " << fname << endl;
|
|
|
|
|
|
|
|
ifstream infile(fname.c_str());
|
|
|
|
if (!infile)
|
|
|
|
break;
|
|
|
|
|
|
|
|
int nCpu, ndim, nPar;
|
|
|
|
|
|
|
|
push_control(CPU);
|
|
|
|
read_int(nCpu);
|
|
|
|
pop_control(CPU);
|
|
|
|
|
|
|
|
push_control(DIM);
|
|
|
|
read_int(ndim);
|
|
|
|
pop_control(DIM);
|
|
|
|
|
|
|
|
push_control(NPAR);
|
|
|
|
read_int(nPar);
|
|
|
|
pop_control(NPAR);
|
|
|
|
|
|
|
|
cout << " NUMCPU=" << nCpu << " NUMDIM=" << ndim << " NumPart=" << nPar << endl;
|
|
|
|
|
|
|
|
totPart += nPar;
|
|
|
|
id++;
|
|
|
|
}
|
|
|
|
|
|
|
|
assert (totPart <= ((~(size_t)0)/sizeof(ParticleState)));
|
|
|
|
size_t memSize = sizeof(ParticleState)*(size_t)totPart;
|
|
|
|
cout << " Needing " << memSize / 1024 << " kBytes" << endl;
|
|
|
|
gd->particles = (ParticleState *)malloc(memSize);
|
|
|
|
assert(gd->particles != 0);
|
|
|
|
gd->NumPart = totPart;
|
|
|
|
gd->ntot_withmasses = totPart;
|
|
|
|
gd->header.npart[0] = totPart;
|
|
|
|
for (int k = 1; k < 6; k++)
|
|
|
|
{
|
|
|
|
gd->header.npart[k] = 0;
|
|
|
|
gd->header.mass[k] = 0;
|
|
|
|
}
|
|
|
|
gd->header.num_files = id;
|
|
|
|
gd->header.BoxSize = 200.0 * 1000; // kPc
|
|
|
|
gd->header.Omega0 = 0.30; // ????
|
|
|
|
gd->header.OmegaLambda = 0.70; // ????
|
|
|
|
gd->header.HubbleParam = 0.70; // ????
|
|
|
|
|
|
|
|
cout << " Total number part=" << totPart << endl
|
|
|
|
<< "Loading particles ..." << endl;
|
|
|
|
|
|
|
|
uint32_t curPos = 0;
|
|
|
|
id = 1;
|
|
|
|
while (1)
|
|
|
|
{
|
|
|
|
ostringstream ss_fname;
|
|
|
|
ss_fname << name << setfill('0') << setw(5) << id;
|
|
|
|
string fname = ss_fname.str();
|
|
|
|
|
|
|
|
cout << " ... " << id;
|
|
|
|
cout.flush();
|
|
|
|
|
|
|
|
ifstream infile(fname.c_str());
|
|
|
|
if (!infile)
|
|
|
|
break;
|
|
|
|
|
|
|
|
int nCpu, ndim, nPar;
|
|
|
|
|
|
|
|
push_control(CPU);
|
|
|
|
read_int(nCpu);
|
|
|
|
pop_control(CPU);
|
|
|
|
|
|
|
|
push_control(DIM);
|
|
|
|
read_int(ndim);
|
|
|
|
pop_control(DIM);
|
|
|
|
|
|
|
|
push_control(NPAR);
|
|
|
|
read_int(nPar);
|
|
|
|
pop_control(NPAR);
|
|
|
|
|
|
|
|
ParticleState *s = &gd->particles[curPos];
|
|
|
|
|
|
|
|
for (uint32_t i = nPar; i > 0; i--,s++)
|
|
|
|
s->Mass = 1.0;
|
|
|
|
|
|
|
|
s = &gd->particles[curPos];
|
|
|
|
push_control(POSX);
|
|
|
|
for (uint32_t i = nPar; i > 0; i--)
|
|
|
|
{
|
|
|
|
read_real(s->Pos[0]);
|
|
|
|
s->Pos[0] *= gd->header.BoxSize;
|
|
|
|
s++;
|
|
|
|
}
|
|
|
|
pop_control(POSX);
|
|
|
|
|
|
|
|
s = &gd->particles[curPos];
|
|
|
|
push_control(POSY);
|
|
|
|
for (uint32_t i = nPar; i > 0; i--)
|
|
|
|
{
|
|
|
|
read_real(s->Pos[1]);
|
|
|
|
s->Pos[1] *= gd->header.BoxSize;
|
|
|
|
s++;
|
|
|
|
}
|
|
|
|
pop_control(POSY);
|
|
|
|
|
|
|
|
s = &gd->particles[curPos];
|
|
|
|
push_control(POSZ);
|
|
|
|
for (uint32_t i = nPar; i > 0; i--)
|
|
|
|
{
|
|
|
|
read_real(s->Pos[2]);
|
|
|
|
s->Pos[2] *= gd->header.BoxSize;
|
|
|
|
s++;
|
|
|
|
}
|
|
|
|
pop_control(POSZ);
|
|
|
|
|
|
|
|
push_control(VELX);
|
|
|
|
for (uint32_t i = 0; i < nPar; i++)
|
|
|
|
{
|
|
|
|
read_real(gd->particles[curPos+i].Vel[0]);
|
|
|
|
}
|
|
|
|
pop_control(VELX);
|
|
|
|
|
|
|
|
push_control(VELY);
|
|
|
|
for (uint32_t i = 0; i < nPar; i++)
|
|
|
|
{
|
|
|
|
read_real(gd->particles[curPos+i].Vel[1]);
|
|
|
|
}
|
|
|
|
pop_control(VELY);
|
|
|
|
|
|
|
|
push_control(VELZ);
|
|
|
|
for (uint32_t i = 0; i < nPar; i++)
|
|
|
|
{
|
|
|
|
read_real(gd->particles[curPos+i].Vel[2]);
|
|
|
|
}
|
|
|
|
pop_control(VELZ);
|
|
|
|
|
|
|
|
curPos += nPar;
|
|
|
|
|
|
|
|
id++;
|
|
|
|
}
|
|
|
|
|
|
|
|
cout << endl;
|
|
|
|
|
|
|
|
return gd;
|
|
|
|
}
|
|
|
|
|
2009-01-08 16:19:01 +01:00
|
|
|
typedef struct
|
|
|
|
{
|
|
|
|
double unitLength;
|
|
|
|
double aexp;
|
|
|
|
double boxSize;
|
|
|
|
} InfoData;
|
|
|
|
|
|
|
|
int readInfoFile(const char *basename, int outputId, InfoData& info)
|
|
|
|
{
|
|
|
|
ostringstream ss_fname;
|
|
|
|
ss_fname << basename << "/info_" << setfill('0') << setw(5) << outputId << ".txt";
|
|
|
|
|
|
|
|
cout << "Opening info file " << ss_fname.str() << endl;
|
|
|
|
ifstream infile(ss_fname.str().c_str());
|
|
|
|
if (!infile)
|
|
|
|
return 0;
|
|
|
|
|
|
|
|
int err;
|
|
|
|
regex_t unit_l_rx;
|
|
|
|
|
|
|
|
// const char *pattern = "^unit_l[ ]*=[ ]*([0-9\\.E+\\-]+)";
|
|
|
|
const char *pattern = "^([A-Za-z_]+)[ ]*=[ ]*([0-9\\.E+\\-]+)";
|
|
|
|
|
|
|
|
err = regcomp (&unit_l_rx, pattern, REG_EXTENDED);
|
|
|
|
cout << unit_l_rx.re_nsub << endl;
|
|
|
|
if (err)
|
|
|
|
{
|
|
|
|
char errString[255];
|
|
|
|
regerror(err, &unit_l_rx, errString, sizeof(errString));
|
|
|
|
cout << errString << endl;
|
|
|
|
abort();
|
|
|
|
}
|
|
|
|
|
|
|
|
map<string,double> infoMap;
|
|
|
|
string line;
|
|
|
|
while (getline(infile, line))
|
|
|
|
{
|
|
|
|
regmatch_t allMatch[4];
|
|
|
|
if (!regexec(&unit_l_rx, line.c_str(), 4, allMatch, 0))
|
|
|
|
{
|
|
|
|
uint32_t start0 = allMatch[1].rm_so, end0 = allMatch[1].rm_eo;
|
|
|
|
uint32_t start1 = allMatch[2].rm_so, end1 = allMatch[2].rm_eo;
|
|
|
|
|
|
|
|
string keyword = line.substr(start0, end0-start0);
|
|
|
|
istringstream iss(line.substr(start1, end1-start1));
|
|
|
|
double unitLength;
|
|
|
|
|
|
|
|
iss >> unitLength;
|
|
|
|
|
|
|
|
infoMap[keyword] = unitLength;
|
|
|
|
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
regfree(&unit_l_rx);
|
|
|
|
|
|
|
|
info.unitLength = infoMap["unit_l"];
|
|
|
|
info.aexp = infoMap["aexp"];
|
|
|
|
info.boxSize = infoMap["boxlen"];
|
|
|
|
|
|
|
|
return 1;
|
|
|
|
}
|
|
|
|
|
|
|
|
PurePositionData *CosmoTool::loadRamsesPosition(const char *basename, int outputId)
|
2008-12-02 18:00:28 +01:00
|
|
|
{
|
|
|
|
PurePositionData *gd = (PurePositionData *)malloc(sizeof(PurePositionData));
|
|
|
|
int id = 1;
|
|
|
|
uint32_t totPart = 0;
|
|
|
|
int nCpu = 0;
|
2009-01-08 16:19:01 +01:00
|
|
|
InfoData info;
|
|
|
|
|
|
|
|
static const double CM_IN_MPC = 3.08e24;
|
|
|
|
|
|
|
|
if (!readInfoFile(basename, outputId, info))
|
|
|
|
return 0;
|
|
|
|
|
|
|
|
double L0 = info.boxSize*info.unitLength/CM_IN_MPC/info.aexp;
|
|
|
|
cout << "L0=" << L0 << " Mpc" << endl;
|
2008-12-02 18:00:28 +01:00
|
|
|
|
|
|
|
cout << "Detecting number of files and particles..." << endl;
|
|
|
|
while (1)
|
|
|
|
{
|
|
|
|
ostringstream ss_fname;
|
2009-01-08 16:19:01 +01:00
|
|
|
ss_fname << basename << "/part_" << setfill('0') << setw(5) << outputId << ".out" << setfill('0') << setw(5) << id;
|
2008-12-02 18:00:28 +01:00
|
|
|
string fname = ss_fname.str();
|
|
|
|
|
|
|
|
cout << " ... " << fname << endl;
|
|
|
|
|
|
|
|
ifstream infile(fname.c_str());
|
|
|
|
if (!infile)
|
|
|
|
break;
|
|
|
|
|
|
|
|
int ndim, nPar;
|
|
|
|
|
|
|
|
push_control(CPU);
|
|
|
|
read_int(nCpu);
|
|
|
|
pop_control(CPU);
|
|
|
|
|
|
|
|
push_control(DIM);
|
|
|
|
read_int(ndim);
|
|
|
|
pop_control(DIM);
|
|
|
|
|
|
|
|
push_control(NPAR);
|
|
|
|
read_int(nPar);
|
|
|
|
pop_control(NPAR);
|
|
|
|
|
|
|
|
cout << " NUMCPU=" << nCpu << " NUMDIM=" << ndim << " NumPart=" << nPar << endl;
|
|
|
|
|
|
|
|
totPart += nPar;
|
|
|
|
id++;
|
|
|
|
}
|
|
|
|
|
|
|
|
assert (totPart <= ((~(size_t)0)/sizeof(FCoordinates)));
|
|
|
|
size_t memSize = sizeof(FCoordinates)*(size_t)(totPart+totPart/nCpu);
|
|
|
|
cout << " Needing " << memSize / 1024 << " kBytes" << endl;
|
|
|
|
gd->pos = (FCoordinates *)malloc(sizeof(FCoordinates)*totPart);
|
|
|
|
assert(gd->pos != 0);
|
|
|
|
gd->NumPart = totPart;
|
2009-01-08 16:19:01 +01:00
|
|
|
gd->BoxSize = L0*1000;
|
2008-12-02 18:00:28 +01:00
|
|
|
|
|
|
|
cout << " Total number part=" << totPart << endl
|
|
|
|
<< "Loading particles ..." << endl;
|
|
|
|
|
|
|
|
uint32_t curPos = 0;
|
|
|
|
id = 1;
|
|
|
|
while (1)
|
|
|
|
{
|
|
|
|
ostringstream ss_fname;
|
2009-01-08 16:19:01 +01:00
|
|
|
ss_fname << basename << "/part_" << setfill('0') << setw(5) << outputId << ".out" << setfill('0') << setw(5) << id;
|
|
|
|
|
2008-12-02 18:00:28 +01:00
|
|
|
string fname = ss_fname.str();
|
|
|
|
int *idP;
|
|
|
|
|
|
|
|
cout << " ... " << id;
|
|
|
|
cout.flush();
|
|
|
|
|
|
|
|
ifstream infile(fname.c_str());
|
|
|
|
if (!infile)
|
|
|
|
break;
|
|
|
|
|
|
|
|
int nCpu, ndim, nPar;
|
|
|
|
|
|
|
|
push_control(CPU);
|
|
|
|
read_int(nCpu);
|
|
|
|
pop_control(CPU);
|
|
|
|
|
|
|
|
push_control(DIM);
|
|
|
|
read_int(ndim);
|
|
|
|
pop_control(DIM);
|
|
|
|
|
|
|
|
push_control(NPAR);
|
|
|
|
read_int(nPar);
|
|
|
|
pop_control(NPAR);
|
|
|
|
|
|
|
|
FCoordinates *s = new FCoordinates[nPar];
|
|
|
|
|
|
|
|
push_control(POSX);
|
|
|
|
for (uint32_t i = 0; i < nPar; i++)
|
|
|
|
{
|
|
|
|
read_real(s[i][0]);
|
|
|
|
s[i][0] *= gd->BoxSize;
|
|
|
|
}
|
|
|
|
pop_control(POSX);
|
|
|
|
|
|
|
|
push_control(POSY);
|
|
|
|
for (uint32_t i = 0; i < nPar; i++)
|
|
|
|
{
|
|
|
|
read_real(s[i][1]);
|
|
|
|
s[i][1] *= gd->BoxSize;
|
|
|
|
}
|
|
|
|
pop_control(POSY);
|
|
|
|
|
|
|
|
push_control(POSZ);
|
|
|
|
for (uint32_t i = 0; i < nPar; i++)
|
|
|
|
{
|
|
|
|
read_real(s[i][2]);
|
|
|
|
s[i][2] *= gd->BoxSize;
|
|
|
|
}
|
|
|
|
pop_control(POSZ);
|
|
|
|
|
|
|
|
// SKIP VELOCITIES
|
|
|
|
for (int k = 0; k < 3; k++) {
|
|
|
|
push_control(V);
|
|
|
|
for (uint32_t i = nPar; i > 0; i--)
|
|
|
|
{
|
|
|
|
float dummyF;
|
|
|
|
read_real(dummyF);
|
|
|
|
}
|
|
|
|
pop_control(V);
|
|
|
|
}
|
|
|
|
|
|
|
|
float minMass = INFINITY;
|
|
|
|
push_control(MASS);
|
|
|
|
for (uint32_t i = nPar; i > 0; i--)
|
|
|
|
{
|
|
|
|
float dummyF;
|
|
|
|
read_real(dummyF);
|
|
|
|
if (dummyF < minMass) minMass = dummyF;
|
|
|
|
}
|
|
|
|
pop_control(MASS);
|
|
|
|
|
|
|
|
push_control(ID);
|
|
|
|
for (uint32_t i = 0; i < nPar; i++)
|
|
|
|
{
|
|
|
|
int id;
|
|
|
|
read_int(id);
|
|
|
|
id--;
|
|
|
|
assert(id >= 0);
|
|
|
|
assert(id < totPart);
|
|
|
|
memcpy(&gd->pos[id], &s[i], sizeof(FCoordinates));
|
|
|
|
}
|
|
|
|
pop_control(ID);
|
|
|
|
|
|
|
|
delete[] s;
|
|
|
|
|
|
|
|
curPos += nPar;
|
|
|
|
|
|
|
|
id++;
|
|
|
|
}
|
|
|
|
|
|
|
|
cout << endl;
|
|
|
|
|
|
|
|
return gd;
|
|
|
|
}
|
|
|
|
|
|
|
|
#if 0
|
|
|
|
int main(int argc, char **argv)
|
|
|
|
{
|
|
|
|
GadgetData *gd = loadRamses(argv[1]);
|
|
|
|
}
|
|
|
|
|
|
|
|
#endif
|
|
|
|
|
|
|
|
#if 0
|
|
|
|
int main(int argc, char **argv)
|
|
|
|
{
|
|
|
|
GadgetData *gd = loadRamses(argv[1]);
|
|
|
|
}
|
|
|
|
|
|
|
|
#endif
|