More ramses load modules

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;
+}

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

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);