Cleanup: remove redundant doxygen \file argument
[blender.git] / intern / smoke / intern / FLUID_3D_SOLVERS.cpp
1 /** \file \ingroup smoke
2  */
3 //////////////////////////////////////////////////////////////////////
4 // This file is part of Wavelet Turbulence.
5 // 
6 // Wavelet Turbulence is free software: you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation, either version 3 of the License, or
9 // (at your option) any later version.
10 // 
11 // Wavelet Turbulence is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 // GNU General Public License for more details.
15 // 
16 // You should have received a copy of the GNU General Public License
17 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
18 // 
19 // Copyright 2008 Theodore Kim and Nils Thuerey
20 // 
21 // FLUID_3D.cpp: implementation of the FLUID_3D class.
22 //
23 //////////////////////////////////////////////////////////////////////
24 // Both solvers optimized by merging loops and precalculating
25 // stuff used in iteration loop.
26 //              - MiikaH
27 //////////////////////////////////////////////////////////////////////
28
29 #include "FLUID_3D.h"
30 #include <cstring>
31 #define SOLVER_ACCURACY 1e-06
32
33 //////////////////////////////////////////////////////////////////////
34 // solve the heat equation with CG
35 //////////////////////////////////////////////////////////////////////
36 void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip)
37 {
38         int x, y, z;
39         const int twoxr = 2 * _xRes;
40         size_t index;
41         const float heatConst = _dt * _heatDiffusion / (_dx * _dx);
42         float *_q, *_residual, *_direction, *_Acenter;
43
44         // i = 0
45         int i = 0;
46
47         _residual     = new float[_totalCells]; // set 0
48         _direction    = new float[_totalCells]; // set 0
49         _q            = new float[_totalCells]; // set 0
50         _Acenter       = new float[_totalCells]; // set 0
51
52         memset(_residual, 0, sizeof(float)*_totalCells);
53         memset(_q, 0, sizeof(float)*_totalCells);
54         memset(_direction, 0, sizeof(float)*_totalCells);
55         memset(_Acenter, 0, sizeof(float)*_totalCells);
56
57         float deltaNew = 0.0f;
58
59   // r = b - Ax
60   index = _slabSize + _xRes + 1;
61   for (z = 1; z < _zRes - 1; z++, index += twoxr)
62     for (y = 1; y < _yRes - 1; y++, index += 2)
63       for (x = 1; x < _xRes - 1; x++, index++)
64       {
65         // if the cell is a variable
66         _Acenter[index] = 1.0f;
67         if (!skip[index])
68         {
69           // set the matrix to the Poisson stencil in order
70           if (!skip[index + 1]) _Acenter[index] += heatConst;
71           if (!skip[index - 1]) _Acenter[index] += heatConst;
72           if (!skip[index + _xRes]) _Acenter[index] += heatConst;
73           if (!skip[index - _xRes]) _Acenter[index] += heatConst;
74           if (!skip[index + _slabSize]) _Acenter[index] += heatConst;
75           if (!skip[index - _slabSize]) _Acenter[index] += heatConst;
76
77                   _residual[index] = b[index] - (_Acenter[index] * field[index] + 
78           field[index - 1] * (skip[index - 1] ? 0.0f : -heatConst) +
79           field[index + 1] * (skip[index + 1] ? 0.0f : -heatConst) +
80           field[index - _xRes] * (skip[index - _xRes] ? 0.0f : -heatConst) +
81           field[index + _xRes] * (skip[index + _xRes] ? 0.0f : -heatConst) +
82           field[index - _slabSize] * (skip[index - _slabSize] ? 0.0f : -heatConst) +
83           field[index + _slabSize] * (skip[index + _slabSize] ? 0.0f : -heatConst));
84         }
85                 else
86                 {
87         _residual[index] = 0.0f;
88                 }
89
90                 _direction[index] = _residual[index];
91                 deltaNew += _residual[index] * _residual[index];
92       }
93
94
95   // While deltaNew > (eps^2) * delta0
96   const float eps  = SOLVER_ACCURACY;
97   float maxR = 2.0f * eps;
98   while ((i < _iterations) && (maxR > eps))
99   {
100     // q = Ad
101         float alpha = 0.0f;
102
103     index = _slabSize + _xRes + 1;
104     for (z = 1; z < _zRes - 1; z++, index += twoxr)
105       for (y = 1; y < _yRes - 1; y++, index += 2)
106         for (x = 1; x < _xRes - 1; x++, index++)
107         {
108           // if the cell is a variable
109           if (!skip[index])
110           {
111
112                         _q[index] = (_Acenter[index] * _direction[index] + 
113             _direction[index - 1] * (skip[index - 1] ? 0.0f : -heatConst) +
114             _direction[index + 1] * (skip[index + 1] ? 0.0f : -heatConst) +
115             _direction[index - _xRes] * (skip[index - _xRes] ? 0.0f : -heatConst) +
116             _direction[index + _xRes] * (skip[index + _xRes] ? 0.0f : -heatConst) +
117             _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0f : -heatConst) +
118             _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0f : -heatConst));
119           }
120                   else
121                   {
122           _q[index] = 0.0f;
123                   }
124                   alpha += _direction[index] * _q[index];
125         }
126
127     if (fabs(alpha) > 0.0f)
128       alpha = deltaNew / alpha;
129         
130         float deltaOld = deltaNew;
131         deltaNew = 0.0f;
132
133         maxR = 0.0f;
134
135     index = _slabSize + _xRes + 1;
136     for (z = 1; z < _zRes - 1; z++, index += twoxr)
137       for (y = 1; y < _yRes - 1; y++, index += 2)
138         for (x = 1; x < _xRes - 1; x++, index++)
139                 {
140           field[index] += alpha * _direction[index];
141
142                   _residual[index] -= alpha * _q[index];
143           maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
144
145                   deltaNew += _residual[index] * _residual[index];
146                 }
147
148     float beta = deltaNew / deltaOld;
149
150     index = _slabSize + _xRes + 1;
151     for (z = 1; z < _zRes - 1; z++, index += twoxr)
152       for (y = 1; y < _yRes - 1; y++, index += 2)
153         for (x = 1; x < _xRes - 1; x++, index++)
154          _direction[index] = _residual[index] + beta * _direction[index];
155
156         
157     i++;
158   }
159   // cout << i << " iterations converged to " << maxR << endl;
160
161         if (_residual) delete[] _residual;
162         if (_direction) delete[] _direction;
163         if (_q)       delete[] _q;
164         if (_Acenter)  delete[] _Acenter;
165 }
166
167 void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
168 {
169         int x, y, z;
170         size_t index;
171         float *_q, *_Precond, *_h, *_residual, *_direction;
172
173         // i = 0
174         int i = 0;
175
176         _residual     = new float[_totalCells]; // set 0
177         _direction    = new float[_totalCells]; // set 0
178         _q            = new float[_totalCells]; // set 0
179         _h                        = new float[_totalCells]; // set 0
180         _Precond          = new float[_totalCells]; // set 0
181
182         memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
183         memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
184         memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
185         memset(_h, 0, sizeof(float)*_xRes*_yRes*_zRes);
186         memset(_Precond, 0, sizeof(float)*_xRes*_yRes*_zRes);
187
188         float deltaNew = 0.0f;
189
190         // r = b - Ax
191         index = _slabSize + _xRes + 1;
192         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
193                 for (y = 1; y < _yRes - 1; y++, index += 2)
194                   for (x = 1; x < _xRes - 1; x++, index++)
195                   {
196                         // if the cell is a variable
197                         float Acenter = 0.0f;
198                         if (!skip[index])
199                         {
200                           // set the matrix to the Poisson stencil in order
201                           if (!skip[index + 1]) Acenter += 1.0f;
202                           if (!skip[index - 1]) Acenter += 1.0f;
203                           if (!skip[index + _xRes]) Acenter += 1.0f;
204                           if (!skip[index - _xRes]) Acenter += 1.0f;
205                           if (!skip[index + _slabSize]) Acenter += 1.0f;
206                           if (!skip[index - _slabSize]) Acenter += 1.0f;
207
208                           _residual[index] = b[index] - (Acenter * field[index] +  
209                           field[index - 1] * (skip[index - 1] ? 0.0f : -1.0f) +
210                           field[index + 1] * (skip[index + 1] ? 0.0f : -1.0f) +
211                           field[index - _xRes] * (skip[index - _xRes] ? 0.0f : -1.0f)+
212                           field[index + _xRes] * (skip[index + _xRes] ? 0.0f : -1.0f)+
213                           field[index - _slabSize] * (skip[index - _slabSize] ? 0.0f : -1.0f)+
214                           field[index + _slabSize] * (skip[index + _slabSize] ? 0.0f : -1.0f) );
215                         }
216                         else
217                         {
218                         _residual[index] = 0.0f;
219                         }
220
221                         // P^-1
222                         if(Acenter < 1.0f)
223                                 _Precond[index] = 0.0;
224                         else
225                                 _Precond[index] = 1.0f / Acenter;
226
227                         // p = P^-1 * r
228                         _direction[index] = _residual[index] * _Precond[index];
229
230                         deltaNew += _residual[index] * _direction[index];
231                   }
232
233
234   // While deltaNew > (eps^2) * delta0
235   const float eps  = SOLVER_ACCURACY;
236   //while ((i < _iterations) && (deltaNew > eps*delta0))
237   float maxR = 2.0f * eps;
238   // while (i < _iterations)
239   while ((i < _iterations) && (maxR > 0.001f * eps))
240   {
241
242         float alpha = 0.0f;
243
244     index = _slabSize + _xRes + 1;
245     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
246       for (y = 1; y < _yRes - 1; y++, index += 2)
247         for (x = 1; x < _xRes - 1; x++, index++)
248         {
249           // if the cell is a variable
250           float Acenter = 0.0f;
251           if (!skip[index])
252           {
253             // set the matrix to the Poisson stencil in order
254             if (!skip[index + 1]) Acenter += 1.0f;
255             if (!skip[index - 1]) Acenter += 1.0f;
256             if (!skip[index + _xRes]) Acenter += 1.0f;
257             if (!skip[index - _xRes]) Acenter += 1.0f;
258             if (!skip[index + _slabSize]) Acenter += 1.0f;
259             if (!skip[index - _slabSize]) Acenter += 1.0f;
260
261                         _q[index] = Acenter * _direction[index] +  
262             _direction[index - 1] * (skip[index - 1] ? 0.0f : -1.0f) +
263             _direction[index + 1] * (skip[index + 1] ? 0.0f : -1.0f) +
264             _direction[index - _xRes] * (skip[index - _xRes] ? 0.0f : -1.0f) +
265             _direction[index + _xRes] * (skip[index + _xRes] ? 0.0f : -1.0f)+
266             _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0f : -1.0f) +
267             _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0f : -1.0f);
268           }
269                   else
270                   {
271           _q[index] = 0.0f;
272                   }
273
274                   alpha += _direction[index] * _q[index];
275         }
276
277
278     if (fabs(alpha) > 0.0f)
279       alpha = deltaNew / alpha;
280
281         float deltaOld = deltaNew;
282         deltaNew = 0.0f;
283
284         maxR = 0.0;
285
286         float tmp;
287
288     // x = x + alpha * d
289     index = _slabSize + _xRes + 1;
290     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
291       for (y = 1; y < _yRes - 1; y++, index += 2)
292         for (x = 1; x < _xRes - 1; x++, index++)
293                 {
294           field[index] += alpha * _direction[index];
295
296                   _residual[index] -= alpha * _q[index];
297
298                   _h[index] = _Precond[index] * _residual[index];
299
300                   tmp = _residual[index] * _h[index];
301                   deltaNew += tmp;
302                   maxR = (tmp > maxR) ? tmp : maxR;
303
304                 }
305
306
307     // beta = deltaNew / deltaOld
308     float beta = deltaNew / deltaOld;
309
310     // d = h + beta * d
311     index = _slabSize + _xRes + 1;
312     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
313       for (y = 1; y < _yRes - 1; y++, index += 2)
314         for (x = 1; x < _xRes - 1; x++, index++)
315           _direction[index] = _h[index] + beta * _direction[index];
316
317     // i = i + 1
318     i++;
319   }
320   // cout << i << " iterations converged to " << sqrt(maxR) << endl;
321
322         if (_h) delete[] _h;
323         if (_Precond) delete[] _Precond;
324         if (_residual) delete[] _residual;
325         if (_direction) delete[] _direction;
326         if (_q)       delete[] _q;
327 }