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 >i, 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. @@
_______________________________________________
Bf-blender-cvs mailing list
[email protected]
http://lists.blender.org/mailman/listinfo/bf-blender-cvs
