various updates for double precision compiles

- additional USE64 define to identify enabled 64b support
- inverted real grid scale , to be in line with other displays
- new uni file format version, includes build info string and time stamp
- added support for loading/saving uni files when compiling with double precision (uni files are still always float!)
- cleanup of wlt setup , smaller cg fixes
- added uv coord test case, smaller fixes for grids in test case handling
- added build info to particle uni file format

--HG--
branch : particleSys

source/conjugategrad.cpp

@@ -132,9 +132,9 @@ void ApplyPreconditionModifiedIncompCholesky2(Grid<Real>& dst, Grid<Real>& Var1,
         if (!flags.isFluid(i,j,k)) continue;
         const Real p = Aprecond(i,j,k);
         dst(i,j,k) = p * (Var1(i,j,k)
-                 - dst(i-1,j,k) * Ai(i-1,j,k) * p
-                 - dst(i,j-1,k) * Aj(i,j-1,k) * p
-                 - dst(i,j,k-1) * Ak(i,j,k-1) * p);
+                 - dst(i-1,j,k) * Ai(i-1,j,k) * Aprecond(i-1,j,k)
+                 - dst(i,j-1,k) * Aj(i,j-1,k) * Aprecond(i,j-1,k)
+                 - dst(i,j,k-1) * Ak(i,j,k-1) * Aprecond(i,j,k-1) );
     }
 
     // backward substitution
@@ -246,23 +246,25 @@ bool GridCg<APPLYMAT>::iterate() {
 
     // compute norm of the residual?
     if(this->mUseResNorm) { 
-        mResNorm = GridSumSqr(mResidual).sum;
+        mResNorm = GridSumSqr(mResidual).sum; 
     } else {
-        // TODO test...
         mResNorm = mResidual.getMaxAbsValue();        
     }
     //if(mIterations % 10 == 9) debMsg("GridCg::Iteration i="<<mIterations<<", resNorm="<<mResNorm<<" accuracy="<<mAccuracy, 1);
 
     // abort here to safe some work...
-    if(mResNorm<mAccuracy) return false;
+    if(mResNorm<mAccuracy) {
+        mSigma = mResNorm; // this will be returned later on to the caller...
+        return false;
+    }
 
     Real sigmaNew = GridDotProduct(mTmp, mResidual);
     Real beta = sigmaNew / mSigma;
 
     // search =  tmp + beta * search
     UpdateSearchVec (mSearch, mTmp, beta);
 
-    debMsg("PB-Cg::iter i="<<mIterations<<" sigma="<<mSigma<<" alpha="<<alpha<<" beta="<<beta<<" ", CG_DEBUGLEVEL);
+    debMsg("PB-Cg::iter i="<<mIterations<<" sigmaNew="<<sigmaNew<<" sigmaLast="<<mSigma<<" alpha="<<alpha<<" beta="<<beta<<" ", CG_DEBUGLEVEL);
     mSigma = sigmaNew;
 
     //debMsg("PB-CG-Norms::p"<<sqrt( GridOpNormNosqrt(mpDst, mpFlags).getValue() ) <<" search"<<sqrt( GridOpNormNosqrt(mpSearch, mpFlags).getValue(), CG_DEBUGLEVEL ) <<" res"<<sqrt( GridOpNormNosqrt(mpResidual, mpFlags).getValue() ) <<" tmp"<<sqrt( GridOpNormNosqrt(mpTmp, mpFlags).getValue() ), CG_DEBUGLEVEL ); // debug

source/fileio.cpp

@@ -230,7 +230,7 @@ void readGridRaw(const string& name, Grid<T>* grid) {
 #	endif
 }
 
-
+//! legacy headers for reading old files
 typedef struct {
     int dimX, dimY, dimZ;
     int frames, elements, elementType, bytesPerElement, bytesPerFrame;
@@ -239,20 +239,84 @@ typedef struct {
 typedef struct {
     int dimX, dimY, dimZ;
     int gridType, elementType, bytesPerElement;
+} UniLegacyHeader2;
+
+//! uni file header 
+typedef struct {
+    int dimX, dimY, dimZ; // grid size
+    int gridType, elementType, bytesPerElement; // data type info
+	char info[256]; // mantaflow build information
+    unsigned long timestamp; // creation time
 } UniHeader;
 
+//! for test run debugging
+PYTHON void printUniFileInfoString(const string& name) {
+#	if NO_ZLIB!=1
+    gzFile gzf = gzopen(name.c_str(), "rb");
+    if (gzf) { 
+		char ID[5]={0,0,0,0,0};
+		gzread(gzf, ID, 4); 
+		if (!strcmp(ID, "MNT2")) {
+			UniHeader head;
+        	assertMsg (gzread(gzf, &head, sizeof(UniHeader)) == sizeof(UniHeader), "can't read file, no header present"); 
+			gzclose(gzf);
+			debMsg("File '"<<name<<"' info: "<< head.info ,1);
+			return; // all good!
+		}
+		gzclose(gzf);
+	}
+#	endif
+	debMsg("File '"<<name<<"', no valid info string found",1);
+}
+
+template <class T>
+void convertDoubleAndWrite(Grid<T>& grid, void* ptr, gzFile& gzf, UniHeader& head) {
+	errMsg("unknown type, not yet supported");
+}
+
+template <>
+void convertDoubleAndWrite(Grid<int>& grid, void* ptr, gzFile& gzf, UniHeader& head) {
+    gzwrite(gzf, &head, sizeof(UniHeader));
+    gzwrite(gzf, ptr, sizeof(int)*head.dimX*head.dimY*head.dimZ);
+}
+
+template <>
+void convertDoubleAndWrite(Grid<double>& grid, void* ptr, gzFile& gzf, UniHeader& head) {
+	head.bytesPerElement = sizeof(float);
+    gzwrite(gzf, &head, sizeof(UniHeader));
+	float* ptrf = (float*)ptr;
+	for(int i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i,++ptrf) {
+		*ptrf = (float)grid[i];
+	} 
+    gzwrite(gzf, ptr, sizeof(float)* head.dimX*head.dimY*head.dimZ);
+}
+
+template <>
+void convertDoubleAndWrite(Grid<Vector3D<double> >& grid, void* ptr, gzFile& gzf, UniHeader& head) {
+	head.bytesPerElement = sizeof(Vector3D<float>);
+    gzwrite(gzf, &head, sizeof(UniHeader));
+	float* ptrf = (float*)ptr;
+	for(int i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i) {
+		for(int c=0; c<3; ++c) { *ptrf = (float)grid[i][c]; ptrf++; }
+	} 
+    gzwrite(gzf, ptr, sizeof(float)*3 *head.dimX*head.dimY*head.dimZ);
+}
+
 template <class T>
 void writeGridUni(const string& name, Grid<T>* grid) {
     cout << "writing grid " << grid->getName() << " to uni file " << name << endl;
 
 #	if NO_ZLIB!=1
-    char ID[5] = "MNT1";
+    char ID[5] = "MNT2";
     UniHeader head;
 	head.dimX = grid->getSizeX();
     head.dimY = grid->getSizeY();
     head.dimZ = grid->getSizeZ();
     head.gridType = grid->getType();
     head.bytesPerElement = sizeof(T);
+	snprintf( head.info, 256, "%s", buildInfoString().c_str() );	
+	MuTime stamp; stamp.get();
+	head.timestamp = stamp.time;
 
     if (grid->getType() & GridBase::TypeInt)
         head.elementType = 0;
@@ -267,18 +331,60 @@ void writeGridUni(const string& name, Grid<T>* grid) {
     if (!gzf) errMsg("can't open file");
 
     gzwrite(gzf, ID, 4);
+	void* ptr = &((*grid)[0]);
+#	if FLOATINGPOINT_PRECISION!=1
+	// always write float values, even if compiled with double precision...
+	Grid<T> temp(grid->getParent());
+	// "misuse" temp grid as storage for floating point values (we have double, so it will always fit)
+	//ptr = &(temp[0]);
+	//float* ptrf = (float*)ptr;
+	convertDoubleAndWrite( *grid, &(temp[0]), gzf, head);
+#	endif
     gzwrite(gzf, &head, sizeof(UniHeader));
-    gzwrite(gzf, &((*grid)[0]), sizeof(T)*head.dimX*head.dimY*head.dimZ);
+    gzwrite(gzf, ptr, sizeof(T)*head.dimX*head.dimY*head.dimZ);
     gzclose(gzf);
 #	else
     cout << "file format not supported without zlib" << endl;
 #	endif
 };
 
+// grid conversion functions for double precision
+template <class T>
+void convertFloatGridToDouble(Grid<T>& grid, void* ptr, int bytesPerElement) {
+	errMsg("unknown type, not yet supported");
+}
+
+template <>
+void convertFloatGridToDouble<int>(Grid<int>& grid, void* ptr, int bytesPerElement) {
+	assertMsg (bytesPerElement == sizeof(int), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(int) );
+	// easy, nothing to do for ints
+    memcpy(ptr, &(grid[0]), sizeof(int) * grid.getSizeX()*grid.getSizeY()*grid.getSizeZ() );
+}
+
+template <>
+void convertFloatGridToDouble<double>(Grid<double>& grid, void* ptr, int bytesPerElement) {
+	assertMsg (bytesPerElement == sizeof(float), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(float) );
+	float* ptrf = (float*)ptr;
+	for(int i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i,++ptrf) {
+		grid[i] = (double)(*ptrf);
+	} 
+}
+
+template <>
+void convertFloatGridToDouble<Vec3>(Grid<Vec3>& grid, void* ptr, int bytesPerElement) {
+	assertMsg (bytesPerElement == sizeof(Vector3D<float>), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(Vector3D<float>) );
+	float* ptrf = (float*)ptr;
+	for(int i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i) {
+		Vec3 v;
+		for(int c=0; c<3; ++c) { v[c] = double(*ptrf); ptrf++; }
+		grid[i] = v;
+	} 
+}
+
 template <class T>
 void readGridUni(const string& name, Grid<T>* grid) {
     cout << "reading grid " << grid->getName() << " from uni file " << name << endl;
-    
+
 #	if NO_ZLIB!=1
     gzFile gzf = gzopen(name.c_str(), "rb");
     if (!gzf) errMsg("can't open file");
@@ -299,21 +405,37 @@ void readGridUni(const string& name, Grid<T>* grid) {
         gzread(gzf, &((*grid)[0]), sizeof(T)*numEl);
     } 
     else if (!strcmp(ID, "MNT1")) {
+        // legacy file format 2
+        UniLegacyHeader2 head;
+        assertMsg (gzread(gzf, &head, sizeof(UniLegacyHeader2)) == sizeof(UniLegacyHeader2), "can't read file, no header present");
+		assertMsg (head.dimX == grid->getSizeX() && head.dimY == grid->getSizeY() && head.dimZ == grid->getSizeZ(), "grid dim doesn't match, "<< Vec3(head.dimX,head.dimY,head.dimZ)<<" vs "<< grid->getSize() );
+		assertMsg (head.gridType == grid->getType(), "grid type doesn't match "<< head.gridType<<" vs "<< grid->getType() );
+		assertMsg (head.bytesPerElement == sizeof(T), "grid element size doesn't match "<< head.bytesPerElement <<" vs "<< sizeof(T) );
+        gzread(gzf, &((*grid)[0]), sizeof(T)*head.dimX*head.dimY*head.dimZ);
+    }
+    else if (!strcmp(ID, "MNT2")) {
         // current file format
         UniHeader head;
         assertMsg (gzread(gzf, &head, sizeof(UniHeader)) == sizeof(UniHeader), "can't read file, no header present");
 		assertMsg (head.dimX == grid->getSizeX() && head.dimY == grid->getSizeY() && head.dimZ == grid->getSizeZ(), "grid dim doesn't match, "<< Vec3(head.dimX,head.dimY,head.dimZ)<<" vs "<< grid->getSize() );
 		assertMsg (head.gridType == grid->getType(), "grid type doesn't match "<< head.gridType<<" vs "<< grid->getType() );
+#		if FLOATINGPOINT_PRECISION!=1
+		// convert float to double
+		Grid<T> temp(grid->getParent());
+		void*  ptr  = &(temp[0]);
+        gzread(gzf, ptr, sizeof(T)*head.dimX*head.dimY*head.dimZ);
+		convertFloatGridToDouble<T>(*grid, ptr, head.bytesPerElement);
+#		else
 		assertMsg (head.bytesPerElement == sizeof(T), "grid element size doesn't match "<< head.bytesPerElement <<" vs "<< sizeof(T) );
         gzread(gzf, &((*grid)[0]), sizeof(T)*head.dimX*head.dimY*head.dimZ);
+#		endif
     }
     gzclose(gzf);
 #	else
     cout << "file format not supported without zlib" << endl;
 #	endif
 };
 
-
 template <class T>
 void writeGridVol(const string& name, Grid<T>* grid) {
 	cout << "writing grid " << grid->getName() << " to vol file " << name << endl;

source/fluidsolver.cpp

@@ -147,6 +147,10 @@ void FluidSolver::printMemInfo() {
 	printf("%s\n", msg.str().c_str() );
 }
 
+PYTHON void printBuildInfo() {
+	debMsg( "Build info: "<<buildInfoString().c_str()<<" ",1);
+}
+
 void FluidSolver::saveMeanTimings(string filename) {
     ofstream ofs(filename.c_str());
     if (!ofs.good())

source/general.cpp

@@ -77,4 +77,46 @@ ostream& operator<<(ostream& os, const MuTime& t) {
     return os;
 }
 
-} // namespace
+std::string buildInfoString() {
+	std::ostringstream infoStr;
+	infoStr << "mantaflow";
+	// TODO , include hg branch info
+
+	// os
+#	ifdef WIN32
+		infoStr << " win";
+#	endif
+#	ifdef __APPLE__
+		infoStr << " mac";
+#	endif
+#	ifdef LINUX
+		infoStr << " linux";
+#	endif
+
+	// 32/64
+#	ifdef USE64
+		infoStr << " 64bit";
+#	else
+		infoStr << " 32bit";
+#	endif
+
+	// fp precision
+#	if FLOATINGPOINT_PRECISION==2
+		infoStr << " fp2";
+#	else
+		infoStr << " fp1";
+#	endif
+
+	// other compile switches
+#	ifdef DEBUG
+		infoStr << " debug";
+#	endif 
+#	ifdef OPENMP
+		infoStr << " omp";
+#	endif
+
+	infoStr << " from "<< __DATE__<<", "<<__TIME__;
+	return infoStr.str();
+}
+
+} // namespace

source/general.h

@@ -86,6 +86,8 @@ struct MuTime {
 };
 std::ostream& operator<< (std::ostream& os, const MuTime& t);
 
+//! generate a string with infos about the current mantaflow build
+std::string buildInfoString();
 
 // Some commonly used math helpers
 template<class T> inline T square(T a) {

source/grid.h

@@ -147,7 +147,7 @@ class Grid : public GridBase {
     PYTHON void copyFrom(const Grid<T>& a) { *this = a; }
 
 	//! debugging helper, print grid from python
-	PYTHON void print(int zSlice=-1,  bool printIndex=false); 
+	PYTHON void printGrid(int zSlice=-1,  bool printIndex=false); 
 
     // common compound operators
     //! get absolute max value in grid (only Real grids)

source/gui/particlepainter.cpp

@@ -111,7 +111,7 @@ void ParticlePainter::paint() {
 	mHavePdata = false;
 	mMaxVal = 0.;
 
-    Real scale = getScale(); // 0.4;
+    Real scale = getScale(); 
 
     glDisable(GL_BLEND);
     glEnable(GL_DEPTH_TEST);

source/particle.cpp

@@ -182,9 +182,12 @@ void BasicParticleSystem::writeParticlesRawVelocityGz(string name) {
 #	endif
 }
 
+//! in line with grid uni header
 typedef struct {
-    int dim;
-    int elementType, bytesPerElement;
+    int dim; // number of partilces
+    int elementType, bytesPerElement; // type id and byte size
+	char info[256]; // mantaflow build information
+    unsigned long timestamp; // creation time
 } UniPartHeader;
 
 template <class T>
@@ -197,6 +200,9 @@ void writeParticlesUni(const string& name, BasicParticleSystem* parts ) {
 	head.dim      = parts->size();
     head.bytesPerElement = sizeof(T);
     head.elementType = 0; // 0 for base data
+	snprintf( head.info, 256, "%s", buildInfoString().c_str() );	
+	MuTime stamp; stamp.get();
+	head.timestamp = stamp.time;
 
     gzFile gzf = gzopen(name.c_str(), "wb1"); // do some compression
     if (!gzf) errMsg("can't open file");
@@ -225,7 +231,7 @@ void readParticlesUni(const string& name, BasicParticleSystem* parts ) {
         // current file format
         UniPartHeader head;
         assertMsg (gzread(gzf, &head, sizeof(UniPartHeader)) == sizeof(UniPartHeader), "can't read file, no header present");
-        assertMsg (head.bytesPerElement == sizeof(T), "particle type doesn't match");
+        assertMsg ( ((head.bytesPerElement == sizeof(T)) && (head.elementType==0) ), "particle type doesn't match");
 
 		// re-allocate all data
 		parts->resizeAll( head.dim );
@@ -366,6 +372,9 @@ void writePdataUni(const string& name, ParticleDataImpl<T>* pdata ) {
 	head.dim      = pdata->size();
     head.bytesPerElement = sizeof(T);
     head.elementType = 1; // 1 for particle data, todo - add sub types?
+	snprintf( head.info, 256, "%s", buildInfoString().c_str() );	
+	MuTime stamp; stamp.get();
+	head.timestamp = stamp.time;
 
     gzFile gzf = gzopen(name.c_str(), "wb1"); // do some compression
     if (!gzf) errMsg("can't open file");
@@ -393,7 +402,7 @@ void readPdataUni(const string& name, ParticleDataImpl<T>* pdata ) {
     if (!strcmp(ID, "PD01")) {
         UniPartHeader head;
         assertMsg (gzread(gzf, &head, sizeof(UniPartHeader)) == sizeof(UniPartHeader), "can't read file, no header present");
-        assertMsg (head.bytesPerElement == sizeof(T), "pdata type doesn't match");
+        assertMsg ( ((head.bytesPerElement == sizeof(T)) && (head.elementType==1) ), "pdata type doesn't match");
         assertMsg (head.dim == pdata->size() , "pdata size doesn't match");
     	int bytes = sizeof(T)*head.dim;
         int readBytes = gzread(gzf, &(pdata->get(0)), sizeof(T)*head.dim);

tools/tests/helperInclude.py

@@ -60,7 +60,7 @@ def getGenRefFileSetting( ):
 def getStrictSetting( ):
 	# check env var whether strict mode enabled
 	ret = int(os.getenv('MANTA_TEST_STRICT', 0))
-	print("Strict-test-setting: " + str(ret))
+	#print("Strict-test-setting: " + str(ret))
 	if(ret>0):
 		return 1
 	return 0