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