Undo revision 23130 which was a merge with 2.5, a messy one because I did something...
[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
23 #include "FLUID_3D.h"
24 #include <cstring>
25 #define SOLVER_ACCURACY 1e-06
26
27 void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
28 {
29         int x, y, z;
30         size_t index;
31
32         // i = 0
33         int i = 0;
34
35         memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
36         memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
37         memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
38         memset(_h, 0, sizeof(float)*_xRes*_yRes*_zRes);
39         memset(_Precond, 0, sizeof(float)*_xRes*_yRes*_zRes);
40
41         // r = b - Ax
42         index = _slabSize + _xRes + 1;
43         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
44                 for (y = 1; y < _yRes - 1; y++, index += 2)
45                   for (x = 1; x < _xRes - 1; x++, index++)
46                   {
47                         // if the cell is a variable
48                         float Acenter = 0.0f;
49                         if (!skip[index])
50                         {
51                           // set the matrix to the Poisson stencil in order
52                           if (!skip[index + 1]) Acenter += 1.;
53                           if (!skip[index - 1]) Acenter += 1.;
54                           if (!skip[index + _xRes]) Acenter += 1.;
55                           if (!skip[index - _xRes]) Acenter += 1.;
56                           if (!skip[index + _slabSize]) Acenter += 1.;
57                           if (!skip[index - _slabSize]) Acenter += 1.;
58                         }
59                     
60                         _residual[index] = b[index] - (Acenter * field[index] +  
61                           field[index - 1] * (skip[index - 1] ? 0.0 : -1.0f)+ 
62                           field[index + 1] * (skip[index + 1] ? 0.0 : -1.0f)+
63                           field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f)+ 
64                           field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
65                           field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f)+ 
66                           field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
67                         _residual[index] = (skip[index]) ? 0.0f : _residual[index];
68
69                         // P^-1
70                         if(Acenter < 1.0)
71                                 _Precond[index] = 0.0;
72                         else
73                                 _Precond[index] = 1.0 / Acenter;
74
75                         // p = P^-1 * r
76                         _direction[index] = _residual[index] * _Precond[index];
77                   }
78
79         // deltaNew = transpose(r) * p
80         float deltaNew = 0.0f;
81         index = _slabSize + _xRes + 1;
82         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
83                 for (y = 1; y < _yRes - 1; y++, index += 2)
84                   for (x = 1; x < _xRes - 1; x++, index++)
85                         deltaNew += _residual[index] * _direction[index];
86
87         // delta0 = deltaNew
88         // float delta0 = deltaNew;
89
90   // While deltaNew > (eps^2) * delta0
91   const float eps  = SOLVER_ACCURACY;
92   //while ((i < _iterations) && (deltaNew > eps*delta0))
93   float maxR = 2.0f * eps;
94   // while (i < _iterations)
95   while ((i < _iterations) && (maxR > 0.001*eps))
96   {
97     // (s) q = Ad (p)
98     index = _slabSize + _xRes + 1;
99     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
100       for (y = 1; y < _yRes - 1; y++, index += 2)
101         for (x = 1; x < _xRes - 1; x++, index++)
102         {
103           // if the cell is a variable
104           float Acenter = 0.0f;
105           if (!skip[index])
106           {
107             // set the matrix to the Poisson stencil in order
108             if (!skip[index + 1]) Acenter += 1.;
109             if (!skip[index - 1]) Acenter += 1.;
110             if (!skip[index + _xRes]) Acenter += 1.;
111             if (!skip[index - _xRes]) Acenter += 1.;
112             if (!skip[index + _slabSize]) Acenter += 1.;
113             if (!skip[index - _slabSize]) Acenter += 1.;
114           }
115           
116           _q[index] = Acenter * _direction[index] +  
117             _direction[index - 1] * (skip[index - 1] ? 0.0 : -1.0f) + 
118             _direction[index + 1] * (skip[index + 1] ? 0.0 : -1.0f) +
119             _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f) + 
120             _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
121             _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f) + 
122             _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f);
123           _q[index] = (skip[index]) ? 0.0f : _q[index];
124         }
125
126     // alpha = deltaNew / (transpose(d) * q)
127     float alpha = 0.0f;
128     index = _slabSize + _xRes + 1;
129     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
130       for (y = 1; y < _yRes - 1; y++, index += 2)
131         for (x = 1; x < _xRes - 1; x++, index++)
132           alpha += _direction[index] * _q[index];
133     if (fabs(alpha) > 0.0f)
134       alpha = deltaNew / alpha;
135
136     // x = x + alpha * d
137     index = _slabSize + _xRes + 1;
138     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
139       for (y = 1; y < _yRes - 1; y++, index += 2)
140         for (x = 1; x < _xRes - 1; x++, index++)
141           field[index] += alpha * _direction[index];
142
143     // r = r - alpha * q
144         maxR = 0.0;
145     index = _slabSize + _xRes + 1;
146     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
147       for (y = 1; y < _yRes - 1; y++, index += 2)
148         for (x = 1; x < _xRes - 1; x++, index++)
149         {
150           _residual[index] -= alpha * _q[index];
151                   // maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
152         }
153
154         // if(maxR <= eps)
155         //      break;
156
157         // h = P^-1 * r
158          index = _slabSize + _xRes + 1;
159     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
160       for (y = 1; y < _yRes - 1; y++, index += 2)
161         for (x = 1; x < _xRes - 1; x++, index++)
162                 {
163                         _h[index] = _Precond[index] * _residual[index];
164                 }
165
166     // deltaOld = deltaNew
167     float deltaOld = deltaNew;
168
169     // deltaNew = transpose(r) * h
170     deltaNew = 0.0f;
171     index = _slabSize + _xRes + 1;
172     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
173       for (y = 1; y < _yRes - 1; y++, index += 2)
174         for (x = 1; x < _xRes - 1; x++, index++)
175                 {
176           deltaNew += _residual[index] * _h[index];
177                   maxR = (_residual[index]* _h[index] > maxR) ? _residual[index]* _h[index] : maxR;
178                 }
179
180     // beta = deltaNew / deltaOld
181     float beta = deltaNew / deltaOld;
182
183     // d = h + beta * d
184     index = _slabSize + _xRes + 1;
185     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
186       for (y = 1; y < _yRes - 1; y++, index += 2)
187         for (x = 1; x < _xRes - 1; x++, index++)
188           _direction[index] = _h[index] + beta * _direction[index];
189
190     // i = i + 1
191     i++;
192   }
193   // cout << i << " iterations converged to " << sqrt(maxR) << endl;
194 }
195
196 //////////////////////////////////////////////////////////////////////
197 // solve the poisson equation with CG
198 //////////////////////////////////////////////////////////////////////
199 void FLUID_3D::solvePressure(float* field, float* b, unsigned char* skip)
200 {
201         int x, y, z;
202         size_t index;
203
204         // i = 0
205         int i = 0;
206
207         memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
208         memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
209         memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
210
211   // r = b - Ax
212   index = _slabSize + _xRes + 1;
213   for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
214     for (y = 1; y < _yRes - 1; y++, index += 2)
215       for (x = 1; x < _xRes - 1; x++, index++)
216       {
217         // if the cell is a variable
218         float Acenter = 0.0f;
219         if (!skip[index])
220         {
221           // set the matrix to the Poisson stencil in order
222           if (!skip[index + 1]) Acenter += 1.;
223           if (!skip[index - 1]) Acenter += 1.;
224           if (!skip[index + _xRes]) Acenter += 1.;
225           if (!skip[index - _xRes]) Acenter += 1.;
226           if (!skip[index + _slabSize]) Acenter += 1.;
227           if (!skip[index - _slabSize]) Acenter += 1.;
228         }
229         
230         _residual[index] = b[index] - (Acenter * field[index] +  
231           field[index - 1] * (skip[index - 1] ? 0.0 : -1.0f)+ 
232           field[index + 1] * (skip[index + 1] ? 0.0 : -1.0f)+
233           field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f)+ 
234           field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
235           field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f)+ 
236           field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
237         _residual[index] = (skip[index]) ? 0.0f : _residual[index];
238       }
239         
240
241   // d = r
242   index = _slabSize + _xRes + 1;
243   for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
244     for (y = 1; y < _yRes - 1; y++, index += 2)
245       for (x = 1; x < _xRes - 1; x++, index++)
246         _direction[index] = _residual[index];
247
248   // deltaNew = transpose(r) * r
249   float deltaNew = 0.0f;
250   index = _slabSize + _xRes + 1;
251   for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
252     for (y = 1; y < _yRes - 1; y++, index += 2)
253       for (x = 1; x < _xRes - 1; x++, index++)
254         deltaNew += _residual[index] * _residual[index];
255
256   // delta0 = deltaNew
257   // float delta0 = deltaNew;
258
259   // While deltaNew > (eps^2) * delta0
260   const float eps  = SOLVER_ACCURACY;
261   float maxR = 2.0f * eps;
262   while ((i < _iterations) && (maxR > eps))
263   {
264     // q = Ad
265     index = _slabSize + _xRes + 1;
266     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
267       for (y = 1; y < _yRes - 1; y++, index += 2)
268         for (x = 1; x < _xRes - 1; x++, index++)
269         {
270           // if the cell is a variable
271           float Acenter = 0.0f;
272           if (!skip[index])
273           {
274             // set the matrix to the Poisson stencil in order
275             if (!skip[index + 1]) Acenter += 1.;
276             if (!skip[index - 1]) Acenter += 1.;
277             if (!skip[index + _xRes]) Acenter += 1.;
278             if (!skip[index - _xRes]) Acenter += 1.;
279             if (!skip[index + _slabSize]) Acenter += 1.;
280             if (!skip[index - _slabSize]) Acenter += 1.;
281           }
282           
283           _q[index] = Acenter * _direction[index] +  
284             _direction[index - 1] * (skip[index - 1] ? 0.0 : -1.0f) + 
285             _direction[index + 1] * (skip[index + 1] ? 0.0 : -1.0f) +
286             _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f) + 
287             _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
288             _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f) + 
289             _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f);
290           _q[index] = (skip[index]) ? 0.0f : _q[index];
291         }
292
293     // alpha = deltaNew / (transpose(d) * q)
294     float alpha = 0.0f;
295     index = _slabSize + _xRes + 1;
296     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
297       for (y = 1; y < _yRes - 1; y++, index += 2)
298         for (x = 1; x < _xRes - 1; x++, index++)
299           alpha += _direction[index] * _q[index];
300     if (fabs(alpha) > 0.0f)
301       alpha = deltaNew / alpha;
302
303     // x = x + alpha * d
304     index = _slabSize + _xRes + 1;
305     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
306       for (y = 1; y < _yRes - 1; y++, index += 2)
307         for (x = 1; x < _xRes - 1; x++, index++)
308           field[index] += alpha * _direction[index];
309
310     // r = r - alpha * q
311     maxR = 0.0f;
312     index = _slabSize + _xRes + 1;
313     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
314       for (y = 1; y < _yRes - 1; y++, index += 2)
315         for (x = 1; x < _xRes - 1; x++, index++)
316         {
317           _residual[index] -= alpha * _q[index];
318           maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
319         }
320
321     // deltaOld = deltaNew
322     float deltaOld = deltaNew;
323
324     // deltaNew = transpose(r) * r
325     deltaNew = 0.0f;
326     index = _slabSize + _xRes + 1;
327     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
328       for (y = 1; y < _yRes - 1; y++, index += 2)
329         for (x = 1; x < _xRes - 1; x++, index++)
330           deltaNew += _residual[index] * _residual[index];
331
332     // beta = deltaNew / deltaOld
333     float beta = deltaNew / deltaOld;
334
335     // d = r + beta * d
336     index = _slabSize + _xRes + 1;
337     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
338       for (y = 1; y < _yRes - 1; y++, index += 2)
339         for (x = 1; x < _xRes - 1; x++, index++)
340           _direction[index] = _residual[index] + beta * _direction[index];
341
342     // i = i + 1
343     i++;
344   }
345   // cout << i << " iterations converged to " << maxR << endl;
346 }
347
348 //////////////////////////////////////////////////////////////////////
349 // solve the heat equation with CG
350 //////////////////////////////////////////////////////////////////////
351 void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip)
352 {
353         int x, y, z;
354         size_t index;
355         const float heatConst = _dt * _heatDiffusion / (_dx * _dx);
356
357         // i = 0
358         int i = 0;
359
360         memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
361         memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
362         memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
363
364   // r = b - Ax
365   index = _slabSize + _xRes + 1;
366   for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
367     for (y = 1; y < _yRes - 1; y++, index += 2)
368       for (x = 1; x < _xRes - 1; x++, index++)
369       {
370         // if the cell is a variable
371         float Acenter = 1.0f;
372         if (!skip[index])
373         {
374           // set the matrix to the Poisson stencil in order
375           if (!skip[index + 1]) Acenter += heatConst;
376           if (!skip[index - 1]) Acenter += heatConst;
377           if (!skip[index + _xRes]) Acenter += heatConst;
378           if (!skip[index - _xRes]) Acenter += heatConst;
379           if (!skip[index + _slabSize]) Acenter += heatConst;
380           if (!skip[index - _slabSize]) Acenter += heatConst;
381         }
382         
383         _residual[index] = b[index] - (Acenter * field[index] + 
384           field[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
385           field[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
386           field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
387           field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
388           field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
389           field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
390         _residual[index] = (skip[index]) ? 0.0f : _residual[index];
391       }
392
393   // d = r
394   index = _slabSize + _xRes + 1;
395   for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
396     for (y = 1; y < _yRes - 1; y++, index += 2)
397       for (x = 1; x < _xRes - 1; x++, index++)
398         _direction[index] = _residual[index];
399
400   // deltaNew = transpose(r) * r
401   float deltaNew = 0.0f;
402   index = _slabSize + _xRes + 1;
403   for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
404     for (y = 1; y < _yRes - 1; y++, index += 2)
405       for (x = 1; x < _xRes - 1; x++, index++)
406         deltaNew += _residual[index] * _residual[index];
407
408   // delta0 = deltaNew
409   // float delta0 = deltaNew;
410
411   // While deltaNew > (eps^2) * delta0
412   const float eps  = SOLVER_ACCURACY;
413   float maxR = 2.0f * eps;
414   while ((i < _iterations) && (maxR > eps))
415   {
416     // q = Ad
417     index = _slabSize + _xRes + 1;
418     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
419       for (y = 1; y < _yRes - 1; y++, index += 2)
420         for (x = 1; x < _xRes - 1; x++, index++)
421         {
422           // if the cell is a variable
423           float Acenter = 1.0f;
424           if (!skip[index])
425           {
426             // set the matrix to the Poisson stencil in order
427             if (!skip[index + 1]) Acenter += heatConst;
428             if (!skip[index - 1]) Acenter += heatConst;
429             if (!skip[index + _xRes]) Acenter += heatConst;
430             if (!skip[index - _xRes]) Acenter += heatConst;
431             if (!skip[index + _slabSize]) Acenter += heatConst;
432             if (!skip[index - _slabSize]) Acenter += heatConst;
433           }
434
435           _q[index] = (Acenter * _direction[index] + 
436             _direction[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
437             _direction[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
438             _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
439             _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
440             _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
441             _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
442          
443           _q[index] = (skip[index]) ? 0.0f : _q[index];
444         }
445
446     // alpha = deltaNew / (transpose(d) * q)
447     float alpha = 0.0f;
448     index = _slabSize + _xRes + 1;
449     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
450       for (y = 1; y < _yRes - 1; y++, index += 2)
451         for (x = 1; x < _xRes - 1; x++, index++)
452           alpha += _direction[index] * _q[index];
453     if (fabs(alpha) > 0.0f)
454       alpha = deltaNew / alpha;
455
456     // x = x + alpha * d
457     index = _slabSize + _xRes + 1;
458     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
459       for (y = 1; y < _yRes - 1; y++, index += 2)
460         for (x = 1; x < _xRes - 1; x++, index++)
461           field[index] += alpha * _direction[index];
462
463     // r = r - alpha * q
464     maxR = 0.0f;
465     index = _slabSize + _xRes + 1;
466     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
467       for (y = 1; y < _yRes - 1; y++, index += 2)
468         for (x = 1; x < _xRes - 1; x++, index++)
469         {
470           _residual[index] -= alpha * _q[index];
471           maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
472         }
473
474     // deltaOld = deltaNew
475     float deltaOld = deltaNew;
476
477     // deltaNew = transpose(r) * r
478     deltaNew = 0.0f;
479     index = _slabSize + _xRes + 1;
480     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
481       for (y = 1; y < _yRes - 1; y++, index += 2)
482         for (x = 1; x < _xRes - 1; x++, index++)
483           deltaNew += _residual[index] * _residual[index];
484
485     // beta = deltaNew / deltaOld
486     float beta = deltaNew / deltaOld;
487
488     // d = r + beta * d
489     index = _slabSize + _xRes + 1;
490     for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
491       for (y = 1; y < _yRes - 1; y++, index += 2)
492         for (x = 1; x < _xRes - 1; x++, index++)
493           _direction[index] = _residual[index] + beta * _direction[index];
494
495     // i = i + 1
496     i++;
497   }
498   // cout << i << " iterations converged to " << maxR << endl;
499 }
500