diff options
author | Daniel Genrich <daniel.genrich@gmx.net> | 2009-08-07 05:07:45 +0400 |
---|---|---|
committer | Daniel Genrich <daniel.genrich@gmx.net> | 2009-08-07 05:07:45 +0400 |
commit | 7c8c83f30dd0442a666d1f68c8bfd0cf1b68bd89 (patch) | |
tree | e14cdfd257296a7db6a65fb4b429762c280aaf97 /intern/smoke | |
parent | 917fce65a631f500359146854018522602d5496d (diff) |
Smoke: test commit of PCG
Diffstat (limited to 'intern/smoke')
-rw-r--r-- | intern/smoke/intern/FLUID_3D.cpp | 8 | ||||
-rw-r--r-- | intern/smoke/intern/FLUID_3D.h | 3 | ||||
-rw-r--r-- | intern/smoke/intern/FLUID_3D_SOLVERS.cpp | 166 |
3 files changed, 175 insertions, 2 deletions
diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp index 76d94e1b4b7..5f40a823238 100644 --- a/intern/smoke/intern/FLUID_3D.cpp +++ b/intern/smoke/intern/FLUID_3D.cpp @@ -96,6 +96,8 @@ FLUID_3D::FLUID_3D(int *res, int amplify, float *p0, float dt) : _xVorticity = new float[_totalCells]; _yVorticity = new float[_totalCells]; _zVorticity = new float[_totalCells]; + _h = new float[_totalCells]; + _Precond = new float[_totalCells]; for (int x = 0; x < _totalCells; x++) { @@ -118,6 +120,8 @@ FLUID_3D::FLUID_3D(int *res, int amplify, float *p0, float dt) : _yVorticity[x] = 0.0f; _zVorticity[x] = 0.0f; _residual[x] = 0.0f; + _h[x] = 0.0f; + _Precond[x] = 0.0f; _obstacles[x] = false; } @@ -189,6 +193,8 @@ FLUID_3D::~FLUID_3D() if (_yVorticity) delete[] _yVorticity; if (_zVorticity) delete[] _zVorticity; if (_vorticity) delete[] _vorticity; + if (_h) delete[] _h; + if (_Precond) delete[] _Precond; if (_obstacles) delete[] _obstacles; if (_wTurbulence) delete _wTurbulence; @@ -414,7 +420,7 @@ void FLUID_3D::project() copyBorderAll(_pressure); // solve Poisson equation - solvePressure(_pressure, _divergence, _obstacles); + solvePressurePre(_pressure, _divergence, _obstacles); setObstaclePressure(); diff --git a/intern/smoke/intern/FLUID_3D.h b/intern/smoke/intern/FLUID_3D.h index 7548f48250d..3e685a8df16 100644 --- a/intern/smoke/intern/FLUID_3D.h +++ b/intern/smoke/intern/FLUID_3D.h @@ -96,6 +96,8 @@ class FLUID_3D float* _yVorticity; float* _zVorticity; float* _vorticity; + float* _h; + float* _Precond; unsigned char* _obstacles; // CG fields @@ -128,6 +130,7 @@ class FLUID_3D void project(); void diffuseHeat(); void solvePressure(float* field, float* b, unsigned char* skip); + void solvePressurePre(float* field, float* b, unsigned char* skip); void solveHeat(float* field, float* b, unsigned char* skip); // handle obstacle boundaries diff --git a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp index 5fd8f72d79d..c6142bbb4bc 100644 --- a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp +++ b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp @@ -21,8 +21,171 @@ ////////////////////////////////////////////////////////////////////// #include "FLUID_3D.h" + #define SOLVER_ACCURACY 1e-06 +void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip) +{ + int x, y, z, index; + + // i = 0 + int i = 0; + + // 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]; + + // P^-1 + if(Acenter < 1.0) + _Precond[index] = 0.0; + else + _Precond[index] = 1.0 / Acenter; + + // p = P^-1 * r + _direction[index] = _residual[index] * _Precond[index]; + } + + // deltaNew = transpose(r) * p + 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] * _direction[index]; + + // delta0 = deltaNew + float delta0 = deltaNew; + + // 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 > eps)) + { + // (s) q = Ad (p) + 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) + + _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f); + _q[index] = (skip[index]) ? 0.0f : _q[index]; + } + + // alpha = deltaNew / (transpose(d) * q) + float alpha = 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++) + alpha += _direction[index] * _q[index]; + if (fabs(alpha) > 0.0f) + alpha = deltaNew / alpha; + + // x = x + alpha * d + 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++) + field[index] += alpha * _direction[index]; + + // r = r - alpha * q + maxR = 0.0; + 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++) + { + _residual[index] -= alpha * _q[index]; + // maxR = (_residual[index] > maxR) ? _residual[index] : maxR; + } + + // if(maxR <= eps) + // break; + + // h = P^-1 * 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++) + { + _h[index] = _Precond[index] * _residual[index]; + } + + // deltaOld = deltaNew + float deltaOld = deltaNew; + + // deltaNew = transpose(r) * h + 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] * _h[index]; + maxR = (_residual[index]* _h[index] > maxR) ? _residual[index]* _h[index] : maxR; + } + + // beta = deltaNew / deltaOld + float beta = deltaNew / deltaOld; + + // d = h + beta * d + 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] = _h[index] + beta * _direction[index]; + + // i = i + 1 + i++; + } + cout << i << " iterations converged to " << maxR << endl; +} + ////////////////////////////////////////////////////////////////////// // solve the poisson equation with CG ////////////////////////////////////////////////////////////////////// @@ -61,6 +224,7 @@ void FLUID_3D::solvePressure(float* field, float* b, unsigned char* skip) field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) ); _residual[index] = (skip[index]) ? 0.0f : _residual[index]; } + // d = r index = _slabSize + _xRes + 1; @@ -314,6 +478,6 @@ void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip) // i = i + 1 i++; } - cout << i << " iterations converged to " << maxR << endl; + // cout << i << " iterations converged to " << maxR << endl; } |