[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [48476] branches/smoke2/intern/smoke/ intern: WIP comit: put preconditioners in a seperate function, use a loop macro

Daniel Genrich daniel.genrich at gmx.net
Mon Jul 2 00:39:27 CEST 2012


Revision: 48476
          http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=48476
Author:   genscher
Date:     2012-07-01 22:39:14 +0000 (Sun, 01 Jul 2012)
Log Message:
-----------
WIP comit: put preconditioners in a seperate function, use a loop macro

Modified Paths:
--------------
    branches/smoke2/intern/smoke/intern/FLUID_3D.cpp
    branches/smoke2/intern/smoke/intern/FLUID_3D.h
    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-07-01 22:19:19 UTC (rev 48475)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D.cpp	2012-07-01 22:39:14 UTC (rev 48476)
@@ -411,6 +411,17 @@
 	for (int i = 0; i < _totalCells; i++)
 	{
 		_xForce[i] = _yForce[i] = _zForce[i] = 0.0f;
+
+		/*
+		if(_density[i] < FLT_EPSILON)
+		{
+			// _density[i] = 0.0f;
+
+			_xVelocity[i] =
+			_yVelocity[i] = 
+			_zVelocity[i] = 0.0f;
+		}
+		*/
 	}
 
 }
@@ -767,7 +778,7 @@
 {
 	int x, y, z;
 	size_t index;
-
+#if USE_NEW_CG == 1
 	VectorXf p(_totalCells);
 
 	ArrayXd A0(_totalCells);
@@ -799,7 +810,7 @@
 	b.setZero(linearIndex);
 	SparseMatrix<float,RowMajor> A(linearIndex, linearIndex);
 	A.reserve(VectorXi::Constant(linearIndex, 7));
-
+#endif
 	// float *_pressure = new float[_totalCells];
 	float *_divergence   = new float[_totalCells];
 
@@ -1566,36 +1577,7 @@
 
 	setZeroBorder(_density, res, zBegin, zEnd);
 	setZeroBorder(_heat, res, zBegin, zEnd);
-#if 0
-	{
-	const size_t index_ = _slabSize + _xRes + 1;
-	int bb=0;
-	int bt=0;
 
-	if (zBegin == 0) {bb = 1;}
-	if (zEnd == _zRes) {bt = 1;}
-	
-	for (int z = zBegin + bb; z < zEnd - bt; z++)
-	{
-		size_t index = index_ +(z-1)*_slabSize;
-
-		for (int y = 1; y < _yRes - 1; y++, index += 2)
-		{
-			for (int x = 1; x < _xRes - 1; x++, index++)
-			{
-				// clean custom velocities from moving obstacles again
-				if (_obstacles[index])
-				{
-					_xVelocity[index] =
-					_yVelocity[index] =
-					_zVelocity[index] = 0.0f;
-				}
-			}
-		}
-	}
-	}
-#endif
-
 	/*int begin=zBegin * _slabSize;
 	int end=begin + (zEnd - zBegin) * _slabSize;
   for (int x = begin; x < end; x++)

Modified: branches/smoke2/intern/smoke/intern/FLUID_3D.h
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D.h	2012-07-01 22:19:19 UTC (rev 48475)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D.h	2012-07-01 22:39:14 UTC (rev 48476)
@@ -176,6 +176,14 @@
 		void solvePressurePre(VectorXf &b, SparseMatrix<float,RowMajor> &A, ArrayXd &gti, VectorXf &result);
 #else
 		void solvePressurePre(float* field, float* b, unsigned char* skip);
+
+		// diagonal preconditioner
+		void precond_init_diag(float *precond, float* field, unsigned char* skip, size_t index);
+		void precond_apply_diag(float *dest, float* source, float *precond, unsigned char* skip, size_t index);
+
+		// modified incomplete cholesky preconditioner
+		void precond_init_mic(float *precond, float* field, unsigned char* skip, size_t index);
+		void precond_apply_mic(float *dest, float* source, float *precond, unsigned char* skip, size_t index);
 #endif
 
 		void solvePressureJacobian(float* p, float* d, unsigned char* ob);

Modified: branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp	2012-07-01 22:19:19 UTC (rev 48475)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp	2012-07-01 22:39:14 UTC (rev 48476)
@@ -31,9 +31,6 @@
 #include <cstring>
 #define SOLVER_ACCURACY 1e-06
 
-#include <iostream>
-#include <fstream>
-
 //////////////////////////////////////////////////////////////////////
 // solve the heat equation with CG
 //////////////////////////////////////////////////////////////////////
@@ -128,37 +125,37 @@
 							alpha += _direction[index] * _q[index];
 						}
 
-				if (fabs(alpha) > 0.0f)
-					alpha = deltaNew / alpha;
+						if (fabs(alpha) > 0.0f)
+							alpha = deltaNew / alpha;
 
-				float deltaOld = deltaNew;
-				deltaNew = 0.0f;
+						float deltaOld = deltaNew;
+						deltaNew = 0.0f;
 
-				maxR = 0.0f;
+						maxR = 0.0f;
 
-				index = _slabSize + _xRes + 1;
-				for (z = 1; z < _zRes - 1; z++, index += twoxr)
-					for (y = 1; y < _yRes - 1; y++, index += 2)
-						for (x = 1; x < _xRes - 1; x++, index++)
-						{
-							field[index] += alpha * _direction[index];
+						index = _slabSize + _xRes + 1;
+						for (z = 1; z < _zRes - 1; z++, index += twoxr)
+							for (y = 1; y < _yRes - 1; y++, index += 2)
+								for (x = 1; x < _xRes - 1; x++, index++)
+								{
+									field[index] += alpha * _direction[index];
 
-							_residual[index] -= alpha * _q[index];
-							maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
+									_residual[index] -= alpha * _q[index];
+									maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
 
-							deltaNew += _residual[index] * _residual[index];
-						}
+									deltaNew += _residual[index] * _residual[index];
+								}
 
-				float beta = deltaNew / deltaOld;
+								float beta = deltaNew / deltaOld;
 
-				index = _slabSize + _xRes + 1;
-				for (z = 1; z < _zRes - 1; z++, index += twoxr)
-					for (y = 1; y < _yRes - 1; y++, index += 2)
-						for (x = 1; x < _xRes - 1; x++, index++)
-							_direction[index] = _residual[index] + beta * _direction[index];
+								index = _slabSize + _xRes + 1;
+								for (z = 1; z < _zRes - 1; z++, index += twoxr)
+									for (y = 1; y < _yRes - 1; y++, index += 2)
+										for (x = 1; x < _xRes - 1; x++, index++)
+											_direction[index] = _residual[index] + beta * _direction[index];
 
 
-				i++;
+								i++;
 			}
 			// cout << i << " iterations converged to " << maxR << endl;
 
@@ -168,253 +165,103 @@
 			if (_Acenter)  delete[] _Acenter;
 }
 
+#define SMOKE_LOOP \
+	index = _slabSize + _xRes + 1; \
+	for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes) \
+	for (y = 1; y < _yRes - 1; y++, index += 2) \
+	for (x = 1; x < _xRes - 1; x++, index++) 
 
-#if 1
-void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
+void FLUID_3D::precond_init_diag(float *precond, float* source, unsigned char* skip, size_t index)
 {
-  int x, y, z, index;
-  float *_q, *_h, *_residual, *_direction;
+	// if the cell is a variable
+	float Acenter = 0.0f;
+	if (!skip[index])
+	{
+		// set the matrix to the Poisson stencil in order
+		if (!skip[index + 1]) Acenter += 1.;
+		if (!skip[index - 1]) Acenter += 1.;
+		if (!skip[index + _xRes]) Acenter += 1.;
+		if (!skip[index - _xRes]) Acenter += 1.;
+		if (!skip[index + _slabSize]) Acenter += 1.;
+		if (!skip[index - _slabSize]) Acenter += 1.;
+	}
 
-  // i = 0
-  int i = 0;
+	// P^-1
+	if(Acenter < 1.0)
+		precond[index] = 0.0;
+	else
+		precond[index] = 1.0 / Acenter;
+}
 
- _residual  = new float[_totalCells]; // set 0
- _h  = new float[_totalCells]; // set 0
- _q  = new float[_totalCells]; // set 0
- _direction  = new float[_totalCells]; // set 0
+void FLUID_3D::precond_apply_diag(float *dest, float* source, float *precond, unsigned char* skip, size_t index)
+{
+	dest[index] = source[index] * precond[index];
+}
 
- memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
- memset(_h, 0, sizeof(float)*_xRes*_yRes*_zRes);
- memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
- memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
+#define SQUARE(a) ((a)*(a))
+void FLUID_3D::precond_init_mic(float *precond, float* source, unsigned char* skip, size_t index)
+{
+	const Real tau = 0.97;
+	const Real sigma = 0.25;
 
-  // r = b - Ax
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-      {
-        // if the cell is a variable
-        float Acenter = 0.0f;
-        if (!skip[index])
-        {
-          // set the matrix to the Poisson stencil in order
-          if (!skip[index + 1]) Acenter += 1.;
-          if (!skip[index - 1]) Acenter += 1.;
-          if (!skip[index + _xRes]) Acenter += 1.;
-          if (!skip[index - _xRes]) Acenter += 1.;
-          if (!skip[index + _slabSize]) Acenter += 1.;
-          if (!skip[index - _slabSize]) Acenter += 1.;
-        }
-        
-        _residual[index] = b[index] - (Acenter * field[index] +  
-          field[index - 1] * (skip[index - 1] ? 0.0 : -1.0f)+ 
-          field[index + 1] * (skip[index + 1] ? 0.0 : -1.0f)+
-          field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f)+ 
-          field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
-          field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f)+ 
-          field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
-        _residual[index] = (skip[index]) ? 0.0f : _residual[index];
-      }
+	// TODO
+}
 
-  // d = r
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-        _direction[index] = _residual[index];
-
-  // deltaNew = transpose(r) * r
-  float deltaNew = 0.0f;
-  index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
-        deltaNew += _residual[index] * _residual[index];
-
-  // delta0 = deltaNew
-  float delta0 = deltaNew;
-
-  // While deltaNew > (eps^2) * delta0
-  const float eps  = SOLVER_ACCURACY;
-  float maxR = 2.0f * eps;
-  while ((i < _iterations) && (maxR > eps))
-  {
-    // q = Ad
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-        {
-          // if the cell is a variable
-          float Acenter = 0.0f;
-          if (!skip[index])
-          {
-            // set the matrix to the Poisson stencil in order
-            if (!skip[index + 1]) Acenter += 1.;
-            if (!skip[index - 1]) Acenter += 1.;
-            if (!skip[index + _xRes]) Acenter += 1.;
-            if (!skip[index - _xRes]) Acenter += 1.;
-            if (!skip[index + _slabSize]) Acenter += 1.;
-            if (!skip[index - _slabSize]) Acenter += 1.;
-          }
-          
-          _q[index] = Acenter * _direction[index] +  
-            _direction[index - 1] * (skip[index - 1] ? 0.0 : -1.0f) + 
-            _direction[index + 1] * (skip[index + 1] ? 0.0 : -1.0f) +
-            _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f) + 
-            _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
-            _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f) + 

@@ Diff output truncated at 10240 characters. @@



More information about the Bf-blender-cvs mailing list