From 247fc4db6b8cdb55f7335245f4056ad1690be333 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 17 Mar 2011 13:39:20 -0500 Subject: [PATCH] Fixed to correctly support velocities in gadget. Fixed on-disk format on kdtree --- src/loadGadget.cpp | 28 ++++++++++++++++++++-------- src/loadGadget.hpp | 2 +- src/loadRamses.cpp | 1 - src/mykdtree.tcc | 34 +++++++++++++++++++++++++++++----- 4 files changed, 50 insertions(+), 15 deletions(-) diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index a8048e6..c6a2370 100644 --- a/src/loadGadget.cpp +++ b/src/loadGadget.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -63,12 +64,13 @@ PurePositionData *CosmoTool::loadGadgetPosition(const char *fname) -SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags) +SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, int GadgetFormat) { SimuData *data; int p, n; UnformattedRead *f; GadgetHeader h; + float velmul; if (id >= 0) { int k = snprintf(0, 0, "%s.%d", fname, id)+1; @@ -113,9 +115,9 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags) h.flag_cooling = f->readInt32(); h.num_files = f->readInt32(); data->BoxSize = h.BoxSize = f->readReal64(); - h.Omega0 = f->readReal64(); - h.OmegaLambda = f->readReal64(); - h.HubbleParam = f->readReal64(); + data->Omega_M = h.Omega0 = f->readReal64(); + data->Omega_Lambda = h.OmegaLambda = f->readReal64(); + data->Hubble = h.HubbleParam = f->readReal64(); f->endCheckpoint(true); for(int k=0; k<6; k++) @@ -125,6 +127,14 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags) } data->NumPart = NumPart; data->TotalNumPart = NumPartTotal; + if (GadgetFormat == 1) + velmul = sqrt(h.time); + else if (GadgetFormat == 2) + velmul = 1/(h.time); + else { + cerr << "unknown gadget format" << endl; + abort(); + } } catch (const InvalidUnformattedAccess& e) { @@ -197,9 +207,10 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags) f->beginCheckpoint(); for(int k = 0, p = 0; k < 6; k++) { for(int n = 0; n < h.npart[k]; n++) { - data->Vel[0][p] = f->readReal32(); - data->Vel[1][p] = f->readReal32(); - data->Vel[2][p] = f->readReal32(); + // 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; p++; } } @@ -213,7 +224,8 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags) return 0; } - // TODO: FIX THE UNITS OF THESE FUNKY VELOCITIES !!! + // THE VELOCITIES ARE IN PHYSICAL COORDINATES +/// // TODO: FIX THE UNITS OF THESE FUNKY VELOCITIES !!! } else { // Skip velocities f->skip(NumPart*3*sizeof(float)+2*4); diff --git a/src/loadGadget.hpp b/src/loadGadget.hpp index 59a46d4..e92e21f 100644 --- a/src/loadGadget.hpp +++ b/src/loadGadget.hpp @@ -8,7 +8,7 @@ namespace CosmoTool { PurePositionData *loadGadgetPosition(const char *fname); - SimuData *loadGadgetMulti(const char *fname, int id, int flags); + SimuData *loadGadgetMulti(const char *fname, int id, int flags, int GadgetFormat = 1); }; diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp index 74ee575..d243c1a 100644 --- a/src/loadRamses.cpp +++ b/src/loadRamses.cpp @@ -243,7 +243,6 @@ int readInfoFile(const char *basename, int outputId, InfoData& info) iss >> unitLength; infoMap[keyword] = unitLength; - } } diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc index e87b7c2..2a233a4 100644 --- a/src/mykdtree.tcc +++ b/src/mykdtree.tcc @@ -446,6 +446,7 @@ namespace CosmoTool { { char id[KDTREE_DISK_SIGNATURE_LEN]; long nodesUsed, numCells; + long rootId; }; template @@ -456,6 +457,7 @@ namespace CosmoTool { strncpy(h.id, KDTREE_DISK_SIGNATURE, KDTREE_DISK_SIGNATURE_LEN); h.nodesUsed = lastNode; h.numCells = numNodes; + h.rootId = root - nodes; o.write((char*)&h, sizeof(h)); for (long i = 0; i < lastNode; i++) @@ -464,13 +466,17 @@ namespace CosmoTool { node_on_disk.cell_id = nodes[i].value - base_cell; if (nodes[i].children[0] == 0) - node_on_disk.children_node[0] = -1; + node_on_disk.children_node[0] = -1L; else node_on_disk.children_node[0] = nodes[i].children[0] - nodes; + assert((node_on_disk.children_node[0] == -1) || ((node_on_disk.children_node[0] >= 0) && (node_on_disk.children_node[0] < lastNode))); + if (nodes[i].children[1] == 0) - node_on_disk.children_node[1] = -1; + node_on_disk.children_node[1] = -1L; else node_on_disk.children_node[1] = nodes[i].children[1] - nodes; + assert((node_on_disk.children_node[1] == -1) || ((node_on_disk.children_node[1] >= 0) && (node_on_disk.children_node[1] < lastNode))); + memcpy(node_on_disk.minBound, nodes[i].minBound, sizeof(coords)); memcpy(node_on_disk.maxBound, nodes[i].maxBound, sizeof(coords)); @@ -489,10 +495,15 @@ namespace CosmoTool { in.read((char *)&h, sizeof(h)); if (!in || strncmp(h.id, KDTREE_DISK_SIGNATURE, KDTREE_DISK_SIGNATURE_LEN) != 0) - throw InvalidOnDiskKDTree(); + { + std::cerr << "KDTree Signature invalid" << std::endl; + throw InvalidOnDiskKDTree(); + } - if (h.numCells != Ncells || h.nodesUsed < 0) + if (h.numCells != Ncells || h.nodesUsed < 0) { + std::cerr << "The number of cells has changed (" << h.numCells << " != " << Ncells << ") or nodesUsed=" << h.nodesUsed << std::endl; throw InvalidOnDiskKDTree(); + } base_cell = cells; nodes = new Node[h.nodesUsed]; @@ -505,11 +516,18 @@ namespace CosmoTool { in.read((char *)&node_on_disk, sizeof(node_on_disk)); - if (!in || node_on_disk.cell_id > numNodes || node_on_disk.cell_id < 0 || + if (!in) { + std::cerr << "End-of-file reached" << std::endl; + delete[] nodes; + throw InvalidOnDiskKDTree(); + } + if (node_on_disk.cell_id > numNodes || node_on_disk.cell_id < 0 || node_on_disk.children_node[0] > lastNode || node_on_disk.children_node[0] < -1 || node_on_disk.children_node[1] > lastNode || node_on_disk.children_node[1] < -1) { delete[] nodes; + std::cerr << "Invalid cell id or children node id invalid" << std::endl; + std::cerr << node_on_disk.cell_id << std::endl << node_on_disk.children_node[0] << std::endl << node_on_disk.children_node[1] << std::endl; throw InvalidOnDiskKDTree(); } @@ -535,9 +553,15 @@ namespace CosmoTool { if (c != N) { delete[] nodes; + std::cerr << "Coordinates of the cell inconsistent with the boundaries" << std::endl + << " X=" << nodes[i].value->coord[0] << " B=[" << nodes[i].minBound[0] << "," << nodes[i].maxBound[0] << "]" << std::endl + << " Y=" << nodes[i].value->coord[1] << " B=[" << nodes[i].minBound[1] << "," << nodes[i].maxBound[1] << "]" << std::endl + << " Z=" << nodes[i].value->coord[2] << " B=[" << nodes[i].minBound[2] << "," << nodes[i].maxBound[2] << "]" << std::endl; throw InvalidOnDiskKDTree(); } } + + root = &nodes[h.rootId]; sortingHelper = new Cell *[Ncells]; for (uint32_t i = 0; i < Ncells; i++)