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