From afd706c1dc2d70a66f3cf721df38b20a85cecacd Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@curvaton.astro.uiuc.edu>
Date: Tue, 29 Sep 2009 16:32:59 -0500
Subject: [PATCH] Added phase loading support in Ramses

---
 src/loadRamses.cpp | 179 +++++++++++++++++++++++++++++++++++++++++----
 src/loadRamses.hpp |   2 +-
 src/load_data.hpp  |   7 ++
 3 files changed, 174 insertions(+), 14 deletions(-)

diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp
index a7e8870..84d2b82 100644
--- a/src/loadRamses.cpp
+++ b/src/loadRamses.cpp
@@ -416,18 +416,171 @@ CosmoTool::PurePositionData *CosmoTool::loadRamsesPosition(const char *basename,
   return gd;
 }
 
-#if 0
-int main(int argc, char **argv)
+CosmoTool::PhaseSpaceData *CosmoTool::loadRamsesPhase(const char *basename, int outputId, bool quiet)
 {
-  GadgetData *gd = loadRamses(argv[1]);
+  PhaseSpaceData *gd = (PhaseSpaceData *)malloc(sizeof(PhaseSpaceData));
+  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 L0 = info.boxSize*info.unitLength/CM_IN_MPC/info.aexp;
+  if (!quiet)
+    cout << "L0=" << L0 << " Mpc" << endl;
+
+  if (!quiet)
+    cout << "Detecting number of files and particles..." << endl;
+  while (1)
+    {
+      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)
+	{
+	  break;
+	}
+
+    }
+
+  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->NumPart = totPart;
+  gd->BoxSize = L0*1000; 
+
+  if (!quiet)
+    cout << " Total number part=" << totPart << endl
+	 << "Loading particles ..." << endl;
+  
+  uint32_t curPos = 0;
+  id = 1;
+  while (1)
+    {
+      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();
+	  
+
+	  FCoordinates *s = new FCoordinates[nPar];
+	  FCoordinates *vel = new FCoordinates[nPar];
+	  
+	  for (int k = 0; k < 3; k++)
+	    {
+	      infile.beginCheckpoint();
+	      for (uint32_t i = 0; i < nPar; i++)
+		{
+		  s[i][k] = infile.readReal32();
+		  s[i][k] *= gd->BoxSize;
+		}
+	      infile.endCheckpoint();
+	    }
+
+	  // SKIP VELOCITIES
+	  for (int k = 0; k < 3; k++) {	
+	    infile.beginCheckpoint();
+	    for (uint32_t i = 0; i < nPar; i++)
+	      {
+		vel[i][k] = infile.readReal32()*gd->BoxSize/1000*100;
+	      }
+	    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++)
+	    {
+	      int id = infile.readInt32();
+	      id--;
+	      assert(id >= 0);
+	      assert(id < totPart);
+	      memcpy(&gd->pos[id], &s[i], sizeof(FCoordinates));
+	      memcpy(&gd->vel[id], &vel[i], sizeof(FCoordinates));
+	    }
+	  infile.endCheckpoint();
+      
+	  delete[] vel;
+	  delete[] s;
+	  
+	  curPos += nPar;
+	} catch (const NoSuchFileException& e)
+	{
+	  break;
+	}
+      id++;
+    }
+
+  if (!quiet)
+    cout << endl;
+
+  return gd;
 }
-
-#endif
-
-#if 0
-int main(int argc, char **argv)
-{
-  GadgetData *gd = loadRamses(argv[1]);
-}
-
-#endif
diff --git a/src/loadRamses.hpp b/src/loadRamses.hpp
index 7ad966c..55f79e3 100644
--- a/src/loadRamses.hpp
+++ b/src/loadRamses.hpp
@@ -7,7 +7,7 @@ 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);
 };
 
 #endif
diff --git a/src/load_data.hpp b/src/load_data.hpp
index 7321cf1..be2f088 100644
--- a/src/load_data.hpp
+++ b/src/load_data.hpp
@@ -63,6 +63,13 @@ namespace CosmoTool {
     FCoordinates *pos;
   };
 
+  struct PhaseSpaceData {
+    unsigned int NumPart;
+    double BoxSize;
+    FCoordinates *pos;
+    FCoordinates *vel;
+  };
+
   void writeGadget(GadgetData *data, const char *fname);
   GadgetData *loadGadget(const char *fname);
   void freeGadget(GadgetData *data);