[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [47781] branches/smoke2/intern/smoke/ intern: Fix Eigen solver, can be enabled using define in fluid3d.h
Daniel Genrich
daniel.genrich at gmx.net
Tue Jun 12 12:50:29 CEST 2012
Revision: 47781
http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=47781
Author: genscher
Date: 2012-06-12 10:50:23 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Fix Eigen solver, can be enabled using define in fluid3d.h
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-06-12 09:58:38 UTC (rev 47780)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D.cpp 2012-06-12 10:50:23 UTC (rev 47781)
@@ -1001,11 +1001,15 @@
*/
// copyBorderAll(_pressure, 0, _zRes);
- VectorXf result;
+ VectorXf result(linearIndex);
result.setZero(linearIndex);
// solve Poisson equation
- solvePressurePre(_pressure, _divergence, _obstacles, b, A, gti, result);
+#if USE_NEW_CG == 1
+ solvePressurePre(b, A, gti, result);
+#else
+ solvePressurePre(_pressure, _divergence, _obstacles);
+#endif
for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
{
@@ -1048,24 +1052,24 @@
{
float vMask[3] = {1.0f, 1.0f, 1.0f}, vObst[3] = {0, 0, 0};
float vR = 0.0f, vL = 0.0f, vT = 0.0f, vB = 0.0f, vD = 0.0f, vU = 0.0f;
-/*
+#if USE_NEW_CG == 1
float pC = result[gti(FINDEX(x, y, z))]; // center
float pR = result[gti(FINDEX(x + 1, y, z))]; // right
float pL = result[gti(FINDEX(x - 1, y, z))]; // left
- float pT = result[gti(FINDEX(x, y + 1, z))]; // top
- float pB = result[gti(FINDEX(x, y - 1, z))]; // bottom
- float pU = result[gti(FINDEX(x, y, z + 1))]; // Up
- float pD = result[gti(FINDEX(x, y, z - 1))]; // Down
-*/
+ float pU = result[gti(FINDEX(x, y + 1, z))]; // top
+ float pD = result[gti(FINDEX(x, y - 1, z))]; // bottom
+ 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
- float pT = _pressure[index + _xRes]; // top
- float pB = _pressure[index - _xRes]; // bottom
- float pU = _pressure[index + _slabSize]; // Up
- float pD = _pressure[index - _slabSize]; // Down
-
+ float pU = _pressure[index + _xRes]; // Up
+ float pD = _pressure[index - _xRes]; // Down
+ float pT = _pressure[index + _slabSize]; // top
+ float pB = _pressure[index - _slabSize]; // bottom
+#endif
if(!_obstacles[index])
{
// DG TODO: What if obstacle is left + right and one is moving?
Modified: branches/smoke2/intern/smoke/intern/FLUID_3D.h
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D.h 2012-06-12 09:58:38 UTC (rev 47780)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D.h 2012-06-12 10:50:23 UTC (rev 47781)
@@ -47,6 +47,8 @@
#include <Eigen/Sparse>
using namespace Eigen;
+#define USE_NEW_CG 0
+
// Fluid index
#define FINDEX(indexX, indexY, indexZ) ((indexZ) * _xRes * _yRes + (indexY) * _xRes + (indexX))
@@ -170,7 +172,11 @@
void diffuseHeat();
void solvePressure(float* field, float* b, unsigned char* skip);
- void solvePressurePre(float* field, float* b, unsigned char* skip, VectorXf &myb, SparseMatrix<float,RowMajor> &A, ArrayXd >i, VectorXf &result);
+#if USE_NEW_CG == 1
+ void solvePressurePre(VectorXf &b, SparseMatrix<float,RowMajor> &A, ArrayXd >i, VectorXf &result);
+#else
+ void solvePressurePre(float* field, float* b, unsigned char* skip);
+#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-06-12 09:58:38 UTC (rev 47780)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp 2012-06-12 10:50:23 UTC (rev 47781)
@@ -167,51 +167,126 @@
if (_q) delete[] _q;
if (_Acenter) delete[] _Acenter;
}
+#if USE_NEW_CG == 1
+void FLUID_3D::solvePressurePre(VectorXf &b, SparseMatrix<float,RowMajor> &A, ArrayXd >i, VectorXf &result)
+{
+ unsigned int arraySize = b.size();
-void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip, VectorXf &myb, SparseMatrix<float,RowMajor> &A, ArrayXd >i, VectorXf &result)
+ // 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;
+}
+
+#else
+// Original code
+void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
{
int x, y, z;
size_t index;
float *_q, *_Precond, *_h, *_residual, *_direction;
- unsigned int arraySize = myb.size();
// i = 0
int i = 0;
- result.setZero(arraySize);
-
_residual = new float[_totalCells]; // set 0
- VectorXf residual(arraySize);
-
_direction = new float[_totalCells]; // set 0
- VectorXf direction(arraySize);
-
_q = new float[_totalCells]; // set 0
- VectorXf q(arraySize);
-
_h = new float[_totalCells]; // set 0
- VectorXf h(arraySize);
-
_Precond = new float[_totalCells]; // set 0
- SparseMatrix<float,RowMajor> Pi(arraySize, arraySize);
- Pi.reserve(VectorXi::Constant(arraySize, 1));
memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
- residual.setZero(arraySize);
-
memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
- q.setZero(arraySize);
-
memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
- direction.setZero(arraySize);
-
memset(_h, 0, sizeof(float)*_xRes*_yRes*_zRes);
- h.setZero(arraySize);
-
memset(_Precond, 0, sizeof(float)*_xRes*_yRes*_zRes);
float deltaNew = 0.0f;
- float mydeltaNew = 0.0f;
// r = b - Ax
index = _slabSize + _xRes + 1;
@@ -256,34 +331,6 @@
deltaNew += _residual[index] * _direction[index];
}
- // r = b - A * x
- residual = myb - 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)
- {
- // myfile << "(" << it.row() << ", " << it.col() << ") = " << it.value();
- 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;
-
- // ???
- mydeltaNew = residual.dot(direction);
- printf("DELTA_NEW old: %.12f, new: %.12f\n", deltaNew, mydeltaNew);
-
-
// While deltaNew > (eps^2) * delta0
const float eps = SOLVER_ACCURACY;
//while ((i < _iterations) && (deltaNew > eps*delta0))
@@ -293,7 +340,6 @@
{
float alpha = 0.0f;
- float myalpha = 0.0f;
index = _slabSize + _xRes + 1;
for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
@@ -328,33 +374,18 @@
alpha += _direction[index] * _q[index];
}
- // ???
- q = A.selfadjointView<Upper>() * direction;
-
- // alpha = ...??
- myalpha = direction.dot(q);
-
- printf("ALPHA old: %.12f, new: %.12f\n", alpha, myalpha);
- printf("DELTA_NEW old: %.12f, new: %.12f\n", deltaNew, mydeltaNew);
-
if (fabs(alpha) > 0.0f)
{
alpha = deltaNew / alpha;
- myalpha = mydeltaNew / myalpha;
}
- printf("ALPHA old: %.12f, new: %.12f\n", alpha, myalpha);
-
float deltaOld = deltaNew;
- float mydeltaOld = mydeltaNew;
deltaNew = 0.0f;
- mydeltaNew = 0.0f;
maxR = 0.0;
float tmp;
- float mytmp;
// x = x + alpha * d
index = _slabSize + _xRes + 1;
@@ -374,32 +405,9 @@
}
- // x = x + alpha * d
- result += myalpha * direction;
-
- // ????
- residual -= myalpha * q;
-
- // ????
- h = Pi.selfadjointView<Upper>() * residual;
-
- // ????
- // mytmp = residual.dot(h);
-
- // ????
- mydeltaNew = residual.dot(h);
-
- printf("DELTA_NEW old: %.12f, new: %.12f\n", deltaNew, mydeltaNew);
-
- // ????
- // DG TODO: maxR
-
// beta = deltaNew / deltaOld
- float beta = deltaNew / deltaOld;
- float mybeta = mydeltaNew / mydeltaOld;
+ float beta = deltaNew / deltaOld;
- printf("BETA old: %.12f, new: %.12f\n\n", beta, mybeta);
-
// d = h + beta * d
index = _slabSize + _xRes + 1;
for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
@@ -407,8 +415,6 @@
for (x = 1; x < _xRes - 1; x++, index++)
_direction[index] = _h[index] + beta * _direction[index];
- direction = h + mybeta * direction;
-
// i = i + 1
i++;
}
@@ -420,6 +426,8 @@
if (_direction) delete[] _direction;
if (_q) delete[] _q;
}
+#endif
+
#if 0
void FLUID_3D::solvePressureJacobian(float* p, float* d, unsigned char* ob)
{
More information about the Bf-blender-cvs
mailing list