add checkGravity()

src/camr.cc

@@ -274,6 +274,7 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
         if (!count_only) {
           var= new double[ngrida*twotondim*nvarh];
           if (is_gravity && (req_bits&POT_BIT || req_bits&ACC_BIT)) {
+            checkGravity(ngrida,ilevel,icpu,&ngridfile[0][0]);
             varg= new double[ngrida*twotondim*nvarg];
           }
         }
@@ -337,7 +338,10 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
                   assert(n==ngrida*8);
 
                 }
-                else grav.skipBlock();
+                else {
+                  int n=grav.skipBlock();
+                  std::cerr << "n="<<n<<"\n";
+                }
               }
             }
           }
@@ -493,4 +497,59 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
   } //for (int icpu=0 .... 
   return nbody;
 }
+//
+// checkGravity
+//
+bool CAmr::checkGravity(const int ngrida, const int ilevel, const int icpu, int *gridfile)
+{
+  bool status=true;
+  // MUST verify this
+  assert(nvarg==ndim);
+  assert(ngrida>0);
+
+  // pointer to gridfile
+  int (*pgf)[nlevelmax][ncpu+nboundary] = (int (*)[nlevelmax][ncpu+nboundary]) gridfile;
+
+  //double * varg = new double[ngrida*twotondim]; // allocate temporary variable
+  int sdouble = sizeof(double); // sizeof double
+
+  // save grav file position
+  std::streampos position = grav.getPos();
+  try {
+    for (int j=0; j<(nboundary+ncpu); j++) {
+      grav.skipBlock(2);
+      if ((*pgf)[ilevel][j]>0) {
+        // Read grav variables
+        for (int ind=0;ind<twotondim;ind++) {
+          for (int ivar=0; ivar<nvarg; ivar++) {
+            if (j==icpu&& ngrida>0) {
+              //int n=grav.readDataBlock((char *) &varg[0]);
+              int n=grav.skipBlock();
+              //std::cerr << "n="<<n<<"\n";
+              if (n!=ngrida*sdouble) {
+                // old RAMSES code detected, nvarg=ndim and should be ndim+1
+                throw -1;
+              }
+              assert(n==ngrida*8);
+            }
+            else {
+              grav.skipBlock();
+            }
+          }
+        }
+      }
+    }
+  } catch (const int e) {
+    std::cerr << "CAmr::checkGravity, old RAMSES code detected, fixing [nvarg=ndim+1]...\n";
+    nvarg++;
+    status=false;
+    assert(nvarg==ndim+1);
+  }
+
+  //delete [] varg;
+  // get back to file position
+  grav.setPos(position);
+  return status;
+
+}
 } // end of namespace ramses

src/camr.h

@@ -80,6 +80,7 @@ private:
     int lmin,lmax;
     CFortIO  amr, hydro, grav;
     int readHeader();
+    bool checkGravity(const int ngrida,const int ilevel,const int icpu, int * gridfile);
     // amr header variables
     static const double XH, mH, kB;
     int ncpu, ndim, nx, ny ,nz, nlevelmax, ngridmax, nboundary, ngrid_current;

src/cfortio.h

@@ -43,6 +43,12 @@ public:
       return len1;
     } else return 1;
   }
+  inline std::streampos getPos() {
+    return in.tellg();
+  }
+  inline void setPos(const std::streampos position) {
+    in.seekg(position);
+  }
   inline int readData(char * ptr,const size_t size_bytes,const  int items) {
     if (!fake_reading) {
       // get data from file
@@ -79,16 +85,19 @@ public:
     else return 1;
   }
   // skip Block
-  inline void skipBlock(int n=1) {
+  inline int skipBlock(int n=1) {
+    int cpt=0;
     if (!fake_reading) {
       for (int i=0;i<n;i++) {
         int len1 = readFRecord();
+        cpt+=len1;
         in.seekg(len1,std::ios::cur);
         int len2 = readFRecord();
         if (len2==len1) ; // remove warning....
         assert(in.good() && len1==len2);
-      }
+      }      
     }
+    return cpt;
 }
 private:
   std::ifstream in;