[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [47820] branches/smoke2/intern/smoke/ intern: Disable new CG for now

[Daniel Genrich]

Revision: 47820
          http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=47820
Author:   genscher
Date:     2012-06-13 10:12:55 +0000 (Wed, 13 Jun 2012)
Log Message:
-----------
Disable new CG for now

Modified Paths:
--------------
    branches/smoke2/intern/smoke/intern/FLUID_3D.cpp
    branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp

Modified: branches/smoke2/intern/smoke/intern/FLUID_3D.cpp
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D.cpp	2012-06-13 10:03:39 UTC (rev 47819)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D.cpp	2012-06-13 10:12:55 UTC (rev 47820)
@@ -972,12 +972,18 @@
 		if (!_obstacles[FINDEX(x,y,z)])
 			b[gti(FINDEX(x, y, z))] = _divergence[FINDEX(x,y,z)];
 
+	// copyBorderAll(_pressure, 0, _zRes);
+
+	// solve Poisson equation
+#if USE_NEW_CG == 1
 	A.makeCompressed();
-/*
-	ConjugateGradient< SparseMatrix<float, RowMajor>, Lower > solver;
-	solver.setMaxIterations(100);
-	solver.setTolerance(1e-05);
 
+	VectorXf result(linearIndex);
+	result.setZero(linearIndex);
+	ConjugateGradient< SparseMatrix<float, RowMajor>, Upper > solver;
+	solver.setMaxIterations(300);
+	solver.setTolerance(1e-06);
+
 	solver.compute(A);
 	if(solver.info() != Success) 
 	{
@@ -986,7 +992,7 @@
 	}
 	else
 	{
-		VectorXf result = solver.solve(b);
+		result = solver.solve(b);
 		if(solver.info() != Success) 
 		{
 		  // solving failed
@@ -998,20 +1004,21 @@
 		std::cout << "#iterations:     " << solver.iterations() << std::endl;
 		std::cout << "estimated error: " << solver.error()      << std::endl;
 	}
-*/
-	// copyBorderAll(_pressure, 0, _zRes);
-
-	VectorXf result(linearIndex);
-	result.setZero(linearIndex);
-
-	// solve Poisson equation
-#if USE_NEW_CG == 1
-	solvePressurePre(b, A, gti, result);
 #else
 	solvePressurePre(_pressure, _divergence, _obstacles);
 #endif
 
+#if 0
+	for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
 	{
+		float value = (_pressure[i] - result[i]);
+		if(value > 0.01)
+			printf("error: p: %f, b: %f\n", _pressure[i], result[i]);
+	}
+#endif
+
+#if 0
+	{
 		float maxvalue = 0;
 		for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
 		{
@@ -1032,7 +1039,7 @@
 		}
 		// printf("Max pressure: %f, dx: %f\n", maxvalue, _dx);
 	}
-
+#endif
 	// DG TODO: check this function, for now this is done in the next function
 	// setObstaclePressure(_pressure, 0, _zRes);
 
@@ -1054,7 +1061,7 @@
 				float pT = result[gti(FINDEX(x, y, z + 1))]; // Up
 				float pB = result[gti(FINDEX(x, y, z - 1))]; // Down
 #else
-				
+
 				float pC = _pressure[index]; // center
 				float pR = _pressure[index + 1]; // right
 				float pL = _pressure[index - 1]; // left

Modified: branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp	2012-06-13 10:03:39 UTC (rev 47819)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp	2012-06-13 10:12:55 UTC (rev 47820)
@@ -167,104 +167,9 @@
 			if (_q)       delete[] _q;
 			if (_Acenter)  delete[] _Acenter;
 }
-#if USE_NEW_CG == 1
-void FLUID_3D::solvePressurePre(VectorXf &b, SparseMatrix<float,RowMajor> &A, ArrayXd &gti, VectorXf &result)
-{
-	unsigned int arraySize = b.size();
 
-	// i = 0
-	int i = 0;
-
-	result.setZero(arraySize);
-
-	VectorXf residual(arraySize);
-	VectorXf direction(arraySize);
-	VectorXf q(arraySize);
-	VectorXf h(arraySize);
-
-	SparseMatrix<float,RowMajor> Pi(arraySize, arraySize);
-	Pi.reserve(VectorXi::Constant(arraySize, 1));
-
-	residual.setZero(arraySize);
-	q.setZero(arraySize);
-	direction.setZero(arraySize);
-	h.setZero(arraySize);
-
-	float deltaNew = 0.0f;
-
-	// r = b - A * x
-	residual = b -  A.selfadjointView<Upper>() * result;
-
-	// z = M^-1 * r
-	for (int k=0; k<A.outerSize(); ++k)
-		for (SparseMatrix<float,RowMajor>::InnerIterator it(A,k); it; ++it)
-		{
-			if(it.row() == it.col())
-			{
-				float value = 0.0f;
-
-				if(it.value() < 1.0)
-					value = 0.0f;
-				else
-					value = 1.0f / it.value();
-				
-				Pi.insert(it.row(), it.col()) = value;
-			}
-		}
-
-	direction = Pi.selfadjointView<Upper>() * residual;
-	deltaNew = residual.dot(direction);
-
-	// While deltaNew > (eps^2) * delta0
-	const float eps  = SOLVER_ACCURACY;
-	//while ((i < _iterations) && (deltaNew > eps*delta0))
-	float maxR = 2.0f * eps;
-	// while (i < _iterations)
-	while ((i < _iterations) && (maxR > 0.001*eps))
-	{
-
-		float alpha = 0.0f;
-
-		// ???
-		q = A.selfadjointView<Upper>() * direction;
-
-		// alpha = ...??
-		alpha = direction.dot(q);
-
-		if (fabs(alpha) > 0.0f)
-		{
-			alpha = deltaNew / alpha;
-		}
-
-		float deltaOld = deltaNew;
-
-		deltaNew = 0.0f;
-
-		result += alpha * direction;
-		residual -= alpha * q;
-		h = Pi.selfadjointView<Upper>() * residual;
-		deltaNew = residual.dot(h);
-
-		maxR = 0.0;
-		for(unsigned int j = 0; j < arraySize; j++)
-		{
-			if(maxR < residual[j])
-				maxR = residual[j];
-		}
-
-		// beta = deltaNew / deltaOld
-		float beta = deltaNew / deltaOld;
-
-		direction = h + beta * direction;
-
-		// i = i + 1
-		i++;
-	}
-	cout << i << " iterations converged to " << sqrt(maxR) << endl;
-}
-
+#if USE_NEW_CG == 1
 #else
-// Original code
 void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
 {
 	int x, y, z;
@@ -427,7 +332,6 @@
 	if (_q)       delete[] _q;
 }
 #endif
-
 #if 0
 void FLUID_3D::solvePressureJacobian(float* p, float* d, unsigned char* ob)
 {