diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp index d15566a..250bb3d 100644 --- a/src/loadRamses.cpp +++ b/src/loadRamses.cpp @@ -268,7 +268,7 @@ CosmoTool::PurePositionData *CosmoTool::loadRamsesPosition(const char *basename, return 0; double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC); - double L0 = info.boxSize*info.unitLength/CM_IN_MPC/info.aexp; + double L0 = info.boxSize*info.unitLength*hubble/100/CM_IN_MPC/info.aexp; if (!quiet) { cout << "L0=" << L0 << " Mpc" << endl; @@ -436,7 +436,7 @@ CosmoTool::PhaseSpaceData *CosmoTool::loadRamsesPhase(const char *basename, int if (!readInfoFile(basename, outputId, info)) return 0; - double unit_vel = info.unitLength/info.unit_t; + double unit_vel = info.unitLength/info.unit_t/1e5; double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC); double L0 = info.boxSize*info.unitLength/CM_IN_MPC/info.aexp; if (!quiet) { @@ -596,3 +596,165 @@ CosmoTool::PhaseSpaceData *CosmoTool::loadRamsesPhase(const char *basename, int return gd; } + + +CosmoTool::PhaseSpaceDataID *CosmoTool::loadRamsesPhase1(const char *basename, int outputId, int cpuid, bool quiet) +{ + PhaseSpaceDataID *gd = (PhaseSpaceDataID *)malloc(sizeof(PhaseSpaceDataID)); + int id = 1; + uint32_t totPart = 0; + int nCpu = 0; + InfoData info; + + static const double CM_IN_MPC = 3.08e24; + + if (!readInfoFile(basename, outputId, info)) + return 0; + + double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC); + double L0 = info.boxSize*info.unitLength*hubble/100/CM_IN_MPC/info.aexp; + double unit_vel = 100*L0/info.aexp; + if (!quiet) { + cout << "L0=" << L0 << " Mpc" << endl; + cout << "H=" << hubble << " km/s/Mpc" << endl; + cout << "unitvel=" << unit_vel << " km/s" << endl; + } + + if (!quiet) + cout << "Detecting number of files and particles..." << endl; + id = cpuid; + { + ostringstream ss_fname; + ss_fname << basename << "/part_" << setfill('0') << setw(5) << outputId << ".out" << setfill('0') << setw(5) << id; + string fname = ss_fname.str(); + + if (!quiet) + cout << " ... " << fname << endl; + + + try + { + UnformattedRead infile(fname); + + int ndim, nPar; + + infile.beginCheckpoint(); + nCpu = max(1,infile.readInt32()); + infile.endCheckpoint(); + + infile.beginCheckpoint(); + ndim = infile.readInt32(); + infile.endCheckpoint(); + + infile.beginCheckpoint(); + nPar = infile.readInt32(); + infile.endCheckpoint(); + + if (!quiet) + cout << " NUMCPU=" << nCpu << " NUMDIM=" << ndim << " NumPart=" << nPar << endl; + + totPart += nPar; + id++; + } + catch (const NoSuchFileException& e) + { + return 0; + } + } + + assert (totPart <= ((~(size_t)0)/sizeof(FCoordinates))); + size_t memSize = sizeof(FCoordinates)*(size_t)(totPart+totPart/nCpu); + if (!quiet) + cout << " Needing " << memSize / 1024 << " kBytes" << endl; + gd->pos = (FCoordinates *)malloc(sizeof(FCoordinates)*totPart); + assert(gd->pos != 0); + gd->vel = (FCoordinates *)malloc(sizeof(FCoordinates)*totPart); + assert(gd->vel != 0); + gd->ID = (int *)malloc(sizeof(int)*totPart); + assert(gd->ID != 0); + gd->NumPart = totPart; + gd->BoxSize = L0*1000; + gd->hubble = hubble; + + uint32_t curPos = 0; + id = cpuid; + { + ostringstream ss_fname; + ss_fname << basename << "/part_" << setfill('0') << setw(5) << outputId << ".out" << setfill('0') << setw(5) << id; + + string fname = ss_fname.str(); + int *idP; + + if (!quiet) + (cout << " ... " << id).flush(); + + try + { + UnformattedRead infile(fname); + + int nCpu, ndim, nPar; + + infile.beginCheckpoint(); + nCpu = infile.readInt32(); + infile.endCheckpoint(); + + infile.beginCheckpoint(); + ndim = infile.readInt32(); + infile.endCheckpoint(); + + infile.beginCheckpoint(); + nPar = infile.readInt32(); + infile.endCheckpoint(); + + gd->pos = new FCoordinates[nPar]; + gd->vel = new FCoordinates[nPar]; + + for (int k = 0; k < 3; k++) + { + infile.beginCheckpoint(); + for (uint32_t i = 0; i < nPar; i++) + { + gd->pos[i][k] = infile.readReal32()*gd->BoxSize; + } + infile.endCheckpoint(); + } + + // SKIP VELOCITIES + for (int k = 0; k < 3; k++) { + infile.beginCheckpoint(); + for (uint32_t i = 0; i < nPar; i++) + { + gd->vel[i][k] = infile.readReal32()*unit_vel; + } + infile.endCheckpoint(); + } + + float minMass = INFINITY; + infile.beginCheckpoint(); + for (uint32_t i = nPar; i > 0; i--) + { + float dummyF = infile.readReal32(); + if (dummyF < minMass) minMass = dummyF; + } + infile.endCheckpoint(); + + infile.beginCheckpoint(); + for (uint32_t i = 0; i < nPar; i++) + { + gd->ID[i] = infile.readInt32(); + } + infile.endCheckpoint(); + + curPos += nPar; + } + catch (const NoSuchFileException& e) + { + return 0; + } + } + + if (!quiet) + cout << endl; + + return gd; +} diff --git a/src/loadRamses.hpp b/src/loadRamses.hpp index 55f79e3..d725a70 100644 --- a/src/loadRamses.hpp +++ b/src/loadRamses.hpp @@ -8,6 +8,8 @@ namespace CosmoTool { GadgetData *loadRamses(const char *name, bool quiet = false); PurePositionData *loadRamsesPosition(const char *fname, int id, bool quiet = false); PhaseSpaceData *loadRamsesPhase(const char *fname, int id, bool quiet = false); + + PhaseSpaceDataID *loadRamsesPhase1(const char *fname, int id, int cpuid, bool quiet = false); }; #endif diff --git a/src/load_data.hpp b/src/load_data.hpp index a153815..c4bc611 100644 --- a/src/load_data.hpp +++ b/src/load_data.hpp @@ -72,6 +72,15 @@ namespace CosmoTool { FCoordinates *vel; }; + struct PhaseSpaceDataID { + unsigned int NumPart; + double hubble; + double BoxSize; + FCoordinates *pos; + FCoordinates *vel; + int *ID; + }; + void writeGadget(GadgetData *data, const char *fname); GadgetData *loadGadget(const char *fname); void freeGadget(GadgetData *data);