Smoke: test commit of PCG
authorDaniel Genrich <daniel.genrich@gmx.net>
Fri, 7 Aug 2009 01:07:45 +0000 (01:07 +0000)
committerDaniel Genrich <daniel.genrich@gmx.net>
Fri, 7 Aug 2009 01:07:45 +0000 (01:07 +0000)
intern/smoke/intern/FLUID_3D.cpp
intern/smoke/intern/FLUID_3D.h
intern/smoke/intern/FLUID_3D_SOLVERS.cpp

index 76d94e1b4b7f8b674ea10ff4338ca1d7b231147b..5f40a823238cc9233b85fd93b8bf66e8ae3216d3 100644 (file)
@@ -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();
 
index 7548f48250d60df81bb83bd22f86d4679153ebb0..3e685a8df16189fcdfd8bd938d92bafe4082ee24 100644 (file)
@@ -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
index 5fd8f72d79d7145e627b0f3a65893c03fa75b49a..c6142bbb4bc84a4a98e49c918c5b854f46220bc1 100644 (file)
 //////////////////////////////////////////////////////////////////////
 
 #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;
 }