initial commit of the fluid simulator.
[blender.git] / intern / elbeem / intern / lbmfsgrsolver.h
1 /******************************************************************************
2  *
3  * El'Beem - the visual lattice boltzmann freesurface simulator
4  * All code distributed as part of El'Beem is covered by the version 2 of the 
5  * GNU General Public License. See the file COPYING for details.
6  * Copyright 2003-2005 Nils Thuerey
7  *
8  * Combined 2D/3D Lattice Boltzmann standard solver classes
9  *
10  *****************************************************************************/
11
12
13 #ifndef LBMFSGRSOLVER_H
14 #include "utilities.h"
15 #include "arrays.h"
16 #include "lbmdimensions.h"
17 #include "lbmfunctions.h"
18 #include "ntl_scene.h"
19 #include <stdio.h>
20
21 #if PARALLEL==1
22 #include <omp.h>
23 #endif // PARALLEL=1
24 #ifndef PARALLEL
25 #define PARALLEL 0
26 #endif // PARALLEL
27
28 #if ELBEEM_BLENDER==1
29 #include "SDL.h"
30 #include "SDL_thread.h"
31 #include "SDL_mutex.h"
32 extern "C" {
33         void simulateThreadIncreaseFrame(void);
34         extern SDL_mutex *globalBakeLock;
35         extern int        globalBakeState;
36         extern int        globalBakeFrame;
37 }
38 #endif // ELBEEM_BLENDER==1
39
40 #ifndef LBMMODEL_DEFINED
41 // force compiler error!
42 ERROR - define model first!
43 #endif // LBMMODEL_DEFINED
44
45
46 //! debug coordinate accesses and the like? (much slower)
47 #define FSGR_STRICT_DEBUG 0
48
49 //! debug coordinate accesses and the like? (much slower)
50 #define FSGR_OMEGA_DEBUG 0
51
52 //! OPT3D quick LES on/off, only debug/benachmarking
53 #define USE_LES 1
54
55 //! order of interpolation (1/2)
56 #define INTORDER 1
57
58 //! order of interpolation (0=always current/1=interpolate/2=always other)
59 #define TIMEINTORDER 0
60
61 //! refinement border method (1 = small border / 2 = larger)
62 #define REFINEMENTBORDER 1
63
64 // interpolateCellFromCoarse test
65 #define INTCFCOARSETEST 1
66
67 // use optimized 3D code?
68 #if LBMDIM==2
69 #define OPT3D false
70 #else
71 // determine with debugging...
72 #       if FSGR_STRICT_DEBUG==1
73 #               define OPT3D false
74 #       else // FSGR_STRICT_DEBUG==1
75 // usually switch optimizations for 3d on, when not debugging
76 #               define OPT3D true
77 // COMPRT
78 //#             define OPT3D false
79 #       endif // FSGR_STRICT_DEBUG==1
80 #endif
81
82 // enable/disable fine grid compression for finest level
83 #if LBMDIM==3
84 #define COMPRESSGRIDS 1
85 #else 
86 #define COMPRESSGRIDS 0
87 #endif 
88
89 // cell mark debugging
90 #if FSGR_STRICT_DEBUG==10
91 #define debugMarkCell(lev,x,y,z) \
92         errMsg("debugMarkCell",D::mName<<" step: "<<D::mStepCnt<<" lev:"<<(lev)<<" marking "<<PRINT_VEC((x),(y),(z))<<" line "<< __LINE__ ); \
93         debugMarkCellCall((lev),(x),(y),(z));
94 #else // FSGR_STRICT_DEBUG==1
95 #define debugMarkCell(lev,x,y,z) \
96         debugMarkCellCall((lev),(x),(y),(z));
97 #endif // FSGR_STRICT_DEBUG==1
98
99
100 //! threshold for level set fluid generation/isosurface
101 #define LS_FLUIDTHRESHOLD 0.5
102
103 //! invalid mass value for unused mass data
104 #define MASS_INVALID -1000.0
105
106 // empty/fill cells without fluid/empty NB's by inserting them into the full/empty lists?
107 #define FSGR_LISTTRICK          true
108 #define FSGR_LISTTTHRESHEMPTY   0.10
109 #define FSGR_LISTTTHRESHFULL    0.90
110 #define FSGR_MAGICNR            0.025
111 //0.04
112
113 // flag array defines -----------------------------------------------------------------------------------------------
114
115 // lbm testsolver get index define
116 #define _LBMGI(level, ii,ij,ik, is) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
117
118 //! flag array acces macro
119 #define _RFLAG(level,xx,yy,zz,set) mLevel[level].mprsFlags[set][ LBMGI((level),(xx),(yy),(zz),(set)) ]
120 #define _RFLAG_NB(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+D::dfVecX[dir],(yy)+D::dfVecY[dir],(zz)+D::dfVecZ[dir],set) ]
121 #define _RFLAG_NBINV(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+D::dfVecX[D::dfInv[dir]],(yy)+D::dfVecY[D::dfInv[dir]],(zz)+D::dfVecZ[D::dfInv[dir]],set) ]
122
123 // array data layouts
124 // standard array layout  -----------------------------------------------------------------------------------------------
125 #define ALSTRING "Standard Array Layout"
126 //#define _LBMQI(level, ii,ij,ik, is, lunused) ( ((is)*mLevel[level].lOffsz) + (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
127 #define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
128 #define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)])
129 #define _QCELL_NB(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+D::dfVecX[dir],(yy)+D::dfVecY[dir],(zz)+D::dfVecZ[dir],set, l)*dTotalNum +(l)])
130 #define _QCELL_NBINV(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+D::dfVecX[D::dfInv[dir]],(yy)+D::dfVecY[D::dfInv[dir]],(zz)+D::dfVecZ[D::dfInv[dir]],set, l)*dTotalNum +(l)])
131
132 #define QCELLSTEP dTotalNum
133 #define _RACPNT(level, ii,ij,ik, is )  &QCELL(level,ii,ij,ik,is,0)
134 #define _RAC(s,l) (s)[(l)]
135
136 // standard arrays
137 #define CSRC_C    RAC(ccel                                , dC )
138 #define CSRC_E    RAC(ccel + (-1)             *(dTotalNum), dE )
139 #define CSRC_W    RAC(ccel + (+1)             *(dTotalNum), dW )
140 #define CSRC_N    RAC(ccel + (-mLevel[lev].lOffsx)        *(dTotalNum), dN )
141 #define CSRC_S    RAC(ccel + (+mLevel[lev].lOffsx)        *(dTotalNum), dS )
142 #define CSRC_NE   RAC(ccel + (-mLevel[lev].lOffsx-1)      *(dTotalNum), dNE)
143 #define CSRC_NW   RAC(ccel + (-mLevel[lev].lOffsx+1)      *(dTotalNum), dNW)
144 #define CSRC_SE   RAC(ccel + (+mLevel[lev].lOffsx-1)      *(dTotalNum), dSE)
145 #define CSRC_SW   RAC(ccel + (+mLevel[lev].lOffsx+1)      *(dTotalNum), dSW)
146 #define CSRC_T    RAC(ccel + (-mLevel[lev].lOffsy)        *(dTotalNum), dT )
147 #define CSRC_B    RAC(ccel + (+mLevel[lev].lOffsy)        *(dTotalNum), dB )
148 #define CSRC_ET   RAC(ccel + (-mLevel[lev].lOffsy-1)      *(dTotalNum), dET)
149 #define CSRC_EB   RAC(ccel + (+mLevel[lev].lOffsy-1)      *(dTotalNum), dEB)
150 #define CSRC_WT   RAC(ccel + (-mLevel[lev].lOffsy+1)      *(dTotalNum), dWT)
151 #define CSRC_WB   RAC(ccel + (+mLevel[lev].lOffsy+1)      *(dTotalNum), dWB)
152 #define CSRC_NT   RAC(ccel + (-mLevel[lev].lOffsy-mLevel[lev].lOffsx) *(dTotalNum), dNT)
153 #define CSRC_NB   RAC(ccel + (+mLevel[lev].lOffsy-mLevel[lev].lOffsx) *(dTotalNum), dNB)
154 #define CSRC_ST   RAC(ccel + (-mLevel[lev].lOffsy+mLevel[lev].lOffsx) *(dTotalNum), dST)
155 #define CSRC_SB   RAC(ccel + (+mLevel[lev].lOffsy+mLevel[lev].lOffsx) *(dTotalNum), dSB)
156
157 #define XSRC_C(x)    RAC(ccel + (x)                 *dTotalNum, dC )
158 #define XSRC_E(x)    RAC(ccel + ((x)-1)             *dTotalNum, dE )
159 #define XSRC_W(x)    RAC(ccel + ((x)+1)             *dTotalNum, dW )
160 #define XSRC_N(x)    RAC(ccel + ((x)-mLevel[lev].lOffsx)        *dTotalNum, dN )
161 #define XSRC_S(x)    RAC(ccel + ((x)+mLevel[lev].lOffsx)        *dTotalNum, dS )
162 #define XSRC_NE(x)   RAC(ccel + ((x)-mLevel[lev].lOffsx-1)      *dTotalNum, dNE)
163 #define XSRC_NW(x)   RAC(ccel + ((x)-mLevel[lev].lOffsx+1)      *dTotalNum, dNW)
164 #define XSRC_SE(x)   RAC(ccel + ((x)+mLevel[lev].lOffsx-1)      *dTotalNum, dSE)
165 #define XSRC_SW(x)   RAC(ccel + ((x)+mLevel[lev].lOffsx+1)      *dTotalNum, dSW)
166 #define XSRC_T(x)    RAC(ccel + ((x)-mLevel[lev].lOffsy)        *dTotalNum, dT )
167 #define XSRC_B(x)    RAC(ccel + ((x)+mLevel[lev].lOffsy)        *dTotalNum, dB )
168 #define XSRC_ET(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy-1)      *dTotalNum, dET)
169 #define XSRC_EB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy-1)      *dTotalNum, dEB)
170 #define XSRC_WT(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy+1)      *dTotalNum, dWT)
171 #define XSRC_WB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy+1)      *dTotalNum, dWB)
172 #define XSRC_NT(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy-mLevel[lev].lOffsx) *dTotalNum, dNT)
173 #define XSRC_NB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy-mLevel[lev].lOffsx) *dTotalNum, dNB)
174 #define XSRC_ST(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy+mLevel[lev].lOffsx) *dTotalNum, dST)
175 #define XSRC_SB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy+mLevel[lev].lOffsx) *dTotalNum, dSB)
176
177
178
179 #if FSGR_STRICT_DEBUG==1
180
181 #define LBMGI(level,ii,ij,ik, is)                 debLBMGI(level,ii,ij,ik, is)         
182 #define RFLAG(level,xx,yy,zz,set)                 debRFLAG(level,xx,yy,zz,set)            
183 #define RFLAG_NB(level,xx,yy,zz,set, dir)         debRFLAG_NB(level,xx,yy,zz,set, dir)    
184 #define RFLAG_NBINV(level,xx,yy,zz,set, dir)      debRFLAG_NBINV(level,xx,yy,zz,set, dir) 
185
186 #define LBMQI(level,ii,ij,ik, is, l)              debLBMQI(level,ii,ij,ik, is, l)         
187 #define QCELL(level,xx,yy,zz,set,l)               debQCELL(level,xx,yy,zz,set,l)         
188 #define QCELL_NB(level,xx,yy,zz,set, dir,l)       debQCELL_NB(level,xx,yy,zz,set, dir,l)
189 #define QCELL_NBINV(level,xx,yy,zz,set, dir,l)    debQCELL_NBINV(level,xx,yy,zz,set, dir,l)
190 #define RACPNT(level, ii,ij,ik, is )              debRACPNT(level, ii,ij,ik, is )          
191 #define RAC(s,l)                                                debRAC(s,l)                  
192
193 #else // FSGR_STRICT_DEBUG==1
194
195 #define LBMGI(level,ii,ij,ik, is)                 _LBMGI(level,ii,ij,ik, is)         
196 #define RFLAG(level,xx,yy,zz,set)                 _RFLAG(level,xx,yy,zz,set)            
197 #define RFLAG_NB(level,xx,yy,zz,set, dir)         _RFLAG_NB(level,xx,yy,zz,set, dir)    
198 #define RFLAG_NBINV(level,xx,yy,zz,set, dir)      _RFLAG_NBINV(level,xx,yy,zz,set, dir) 
199
200 #define LBMQI(level,ii,ij,ik, is, l)              _LBMQI(level,ii,ij,ik, is, l)         
201 #define QCELL(level,xx,yy,zz,set,l)               _QCELL(level,xx,yy,zz,set,l)         
202 #define QCELL_NB(level,xx,yy,zz,set, dir,l)       _QCELL_NB(level,xx,yy,zz,set, dir, l)
203 #define QCELL_NBINV(level,xx,yy,zz,set, dir,l)    _QCELL_NBINV(level,xx,yy,zz,set, dir,l)
204 #define RACPNT(level, ii,ij,ik, is )              _RACPNT(level, ii,ij,ik, is )          
205 #define RAC(s,l)                                  _RAC(s,l)                  
206
207 #endif // FSGR_STRICT_DEBUG==1
208
209 // general defines -----------------------------------------------------------------------------------------------
210
211 #define TESTFLAG(flag, compflag) ((flag & compflag)==compflag)
212
213 #if LBMDIM==2
214 #define dC 0
215 #define dN 1
216 #define dS 2
217 #define dE 3
218 #define dW 4
219 #define dNE 5
220 #define dNW 6
221 #define dSE 7
222 #define dSW 8
223 #define LBM_DFNUM 9
224 #else
225 // direction indices
226 #define dC 0
227 #define dN 1
228 #define dS 2
229 #define dE 3
230 #define dW 4
231 #define dT 5
232 #define dB 6
233 #define dNE 7
234 #define dNW 8
235 #define dSE 9
236 #define dSW 10
237 #define dNT 11
238 #define dNB 12
239 #define dST 13
240 #define dSB 14
241 #define dET 15
242 #define dEB 16
243 #define dWT 17
244 #define dWB 18
245 #define LBM_DFNUM 19
246 #endif
247 #define dFfrac 19
248 #define dMass 20
249 #define dFlux 21
250 #define dTotalNum 22
251 #define dWB 18
252
253 // default init for dFlux values
254 #define FLUX_INIT 0.5f * (float)(D::cDfNum)
255
256 // only for non DF dir handling!
257 #define dNET 19
258 #define dNWT 20
259 #define dSET 21
260 #define dSWT 22
261 #define dNEB 23
262 #define dNWB 24
263 #define dSEB 25
264 #define dSWB 26
265
266 //! fill value for boundary cells
267 #define BND_FILL 0.0
268
269 #define DFL1 (1.0/ 3.0)
270 #define DFL2 (1.0/18.0)
271 #define DFL3 (1.0/36.0)
272
273 #define OMEGA(l) mLevel[(l)].omega
274
275 #define EQC (  DFL1*(rho - usqr))
276 #define EQN (  DFL2*(rho + uy*(4.5*uy + 3.0) - usqr))
277 #define EQS (  DFL2*(rho + uy*(4.5*uy - 3.0) - usqr))
278 #define EQE (  DFL2*(rho + ux*(4.5*ux + 3.0) - usqr))
279 #define EQW (  DFL2*(rho + ux*(4.5*ux - 3.0) - usqr))
280 #define EQT (  DFL2*(rho + uz*(4.5*uz + 3.0) - usqr))
281 #define EQB (  DFL2*(rho + uz*(4.5*uz - 3.0) - usqr))
282                     
283 #define EQNE ( DFL3*(rho + (+ux+uy)*(4.5*(+ux+uy) + 3.0) - usqr))
284 #define EQNW ( DFL3*(rho + (-ux+uy)*(4.5*(-ux+uy) + 3.0) - usqr))
285 #define EQSE ( DFL3*(rho + (+ux-uy)*(4.5*(+ux-uy) + 3.0) - usqr))
286 #define EQSW ( DFL3*(rho + (-ux-uy)*(4.5*(-ux-uy) + 3.0) - usqr))
287 #define EQNT ( DFL3*(rho + (+uy+uz)*(4.5*(+uy+uz) + 3.0) - usqr))
288 #define EQNB ( DFL3*(rho + (+uy-uz)*(4.5*(+uy-uz) + 3.0) - usqr))
289 #define EQST ( DFL3*(rho + (-uy+uz)*(4.5*(-uy+uz) + 3.0) - usqr))
290 #define EQSB ( DFL3*(rho + (-uy-uz)*(4.5*(-uy-uz) + 3.0) - usqr))
291 #define EQET ( DFL3*(rho + (+ux+uz)*(4.5*(+ux+uz) + 3.0) - usqr))
292 #define EQEB ( DFL3*(rho + (+ux-uz)*(4.5*(+ux-uz) + 3.0) - usqr))
293 #define EQWT ( DFL3*(rho + (-ux+uz)*(4.5*(-ux+uz) + 3.0) - usqr))
294 #define EQWB ( DFL3*(rho + (-ux-uz)*(4.5*(-ux-uz) + 3.0) - usqr))
295
296
297 // this is a bit ugly, but necessary for the CSRC_ access...
298 #define MSRC_C    m[dC ]
299 #define MSRC_N    m[dN ]
300 #define MSRC_S    m[dS ]
301 #define MSRC_E    m[dE ]
302 #define MSRC_W    m[dW ]
303 #define MSRC_T    m[dT ]
304 #define MSRC_B    m[dB ]
305 #define MSRC_NE   m[dNE]
306 #define MSRC_NW   m[dNW]
307 #define MSRC_SE   m[dSE]
308 #define MSRC_SW   m[dSW]
309 #define MSRC_NT   m[dNT]
310 #define MSRC_NB   m[dNB]
311 #define MSRC_ST   m[dST]
312 #define MSRC_SB   m[dSB]
313 #define MSRC_ET   m[dET]
314 #define MSRC_EB   m[dEB]
315 #define MSRC_WT   m[dWT]
316 #define MSRC_WB   m[dWB]
317
318 // this is a bit ugly, but necessary for the ccel local access...
319 #define CCEL_C    RAC(ccel, dC )
320 #define CCEL_N    RAC(ccel, dN )
321 #define CCEL_S    RAC(ccel, dS )
322 #define CCEL_E    RAC(ccel, dE )
323 #define CCEL_W    RAC(ccel, dW )
324 #define CCEL_T    RAC(ccel, dT )
325 #define CCEL_B    RAC(ccel, dB )
326 #define CCEL_NE   RAC(ccel, dNE)
327 #define CCEL_NW   RAC(ccel, dNW)
328 #define CCEL_SE   RAC(ccel, dSE)
329 #define CCEL_SW   RAC(ccel, dSW)
330 #define CCEL_NT   RAC(ccel, dNT)
331 #define CCEL_NB   RAC(ccel, dNB)
332 #define CCEL_ST   RAC(ccel, dST)
333 #define CCEL_SB   RAC(ccel, dSB)
334 #define CCEL_ET   RAC(ccel, dET)
335 #define CCEL_EB   RAC(ccel, dEB)
336 #define CCEL_WT   RAC(ccel, dWT)
337 #define CCEL_WB   RAC(ccel, dWB)
338
339
340 #if PARALLEL==1
341 #define CSMOMEGA_STATS(dlev, domega) 
342 #else // PARALLEL==1
343 #if FSGR_OMEGA_DEBUG==1
344 #define CSMOMEGA_STATS(dlev, domega) \
345         mLevel[dlev].avgOmega += domega; mLevel[dlev].avgOmegaCnt+=1.0; 
346 #else // FSGR_OMEGA_DEBUG==1
347 #define CSMOMEGA_STATS(dlev, domega) 
348 #endif // FSGR_OMEGA_DEBUG==1
349 #endif // PARALLEL==1
350
351
352 // used for main loops and grav init
353 // source set
354 #define SRCS(l) mLevel[(l)].setCurr
355 // target set
356 #define TSET(l) mLevel[(l)].setOther
357
358
359 // complete default stream&collide, 2d/3d
360 /* read distribution funtions of adjacent cells = sweep step */ 
361 #if OPT3D==false 
362
363 #if FSGR_STRICT_DEBUG==1
364 #define MARKCELLCHECK \
365         debugMarkCell(lev,i,j,k); D::mPanic=1;
366 #define STREAMCHECK(ni,nj,nk,nl) \
367         if((m[l] < -1.0) || (m[l]>1.0)) {\
368                 errMsg("STREAMCHECK","Invalid streamed DF l"<<l<<" value:"<<m[l]<<" at "<<PRINT_IJK<<" from "<<PRINT_VEC(ni,nj,nk)<<" nl"<<(nl)<<\
369                                 " nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther)  ); \
370                 MARKCELLCHECK; \
371         }
372 #define COLLCHECK \
373         if( (rho>2.0) || (rho<-1.0) || (ABS(ux)>1.0) || (ABS(uy)>1.0) |(ABS(uz)>1.0) ) {\
374                 errMsg("COLLCHECK","Invalid collision values r:"<<rho<<" u:"PRINT_VEC(ux,uy,uz)<<" at? "<<PRINT_IJK ); \
375                 MARKCELLCHECK; \
376         }
377 #else
378 #define STREAMCHECK(ni,nj,nk,nl) 
379 #define COLLCHECK
380 #endif
381
382 // careful ux,uy,uz need to be inited before!
383
384 #define DEFAULT_STREAM \
385                 m[dC] = RAC(ccel,dC); \
386                 FORDF1 { \
387                         if(NBFLAG( D::dfInv[l] )&CFBnd) { \
388                                 m[l] = RAC(ccel, D::dfInv[l] ); \
389                                 STREAMCHECK(i,j,k, D::dfInv[l]); \
390                         } else { \
391                                 m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
392                                 STREAMCHECK(i+D::dfVecX[D::dfInv[l]], j+D::dfVecY[D::dfInv[l]],k+D::dfVecZ[D::dfInv[l]], l); \
393                         } \
394                 }   
395
396 // careful ux,uy,uz need to be inited before!
397 #define DEFAULT_COLLIDE \
398                         D::collideArrays( m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago, &mDebugOmegaRet ); \
399                         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
400                         FORDF0 { RAC(tcel,l) = m[l]; }   \
401                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
402                         COLLCHECK;
403 #define OPTIMIZED_STREAMCOLLIDE \
404                         m[0] = RAC(ccel,0); \
405                         FORDF1 { /* df0 is set later on... */ \
406                                 /* FIXME CHECK INV ? */\
407                                 if(RFLAG_NBINV(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); D::mPanic=1;  \
408                                 } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
409                                 STREAMCHECK(i+D::dfVecX[D::dfInv[l]], j+D::dfVecY[D::dfInv[l]],k+D::dfVecZ[D::dfInv[l]], l); \
410                         }   \
411                         rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
412                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
413                         D::collideArrays( m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago , &mDebugOmegaRet  ); \
414                         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
415                         FORDF0 { RAC(tcel,l) = m[l]; } \
416                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
417                         COLLCHECK;
418
419 #else  // 3D, opt OPT3D==true
420
421 #define DEFAULT_STREAM \
422                 m[dC] = RAC(ccel,dC); \
423                 /* explicit streaming */ \
424                 if((!nbored & CFBnd)) { \
425                         /* no boundary near?, no real speed diff.? */ \
426                         m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
427                         m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
428                         m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
429                         m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
430                         m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
431                         m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
432                 } else { \
433                         /* explicit streaming */ \
434                         if(NBFLAG(dS )&CFBnd) { m[dN ] = RAC(ccel,dS ); } else { m[dN ] = CSRC_N ; } \
435                         if(NBFLAG(dN )&CFBnd) { m[dS ] = RAC(ccel,dN ); } else { m[dS ] = CSRC_S ; } \
436                         if(NBFLAG(dW )&CFBnd) { m[dE ] = RAC(ccel,dW ); } else { m[dE ] = CSRC_E ; } \
437                         if(NBFLAG(dE )&CFBnd) { m[dW ] = RAC(ccel,dE ); } else { m[dW ] = CSRC_W ; } \
438                         if(NBFLAG(dB )&CFBnd) { m[dT ] = RAC(ccel,dB ); } else { m[dT ] = CSRC_T ; } \
439                         if(NBFLAG(dT )&CFBnd) { m[dB ] = RAC(ccel,dT ); } else { m[dB ] = CSRC_B ; } \
440                         \
441                         if(NBFLAG(dSW)&CFBnd) { m[dNE] = RAC(ccel,dSW); } else { m[dNE] = CSRC_NE; } \
442                         if(NBFLAG(dSE)&CFBnd) { m[dNW] = RAC(ccel,dSE); } else { m[dNW] = CSRC_NW; } \
443                         if(NBFLAG(dNW)&CFBnd) { m[dSE] = RAC(ccel,dNW); } else { m[dSE] = CSRC_SE; } \
444                         if(NBFLAG(dNE)&CFBnd) { m[dSW] = RAC(ccel,dNE); } else { m[dSW] = CSRC_SW; } \
445                         if(NBFLAG(dSB)&CFBnd) { m[dNT] = RAC(ccel,dSB); } else { m[dNT] = CSRC_NT; } \
446                         if(NBFLAG(dST)&CFBnd) { m[dNB] = RAC(ccel,dST); } else { m[dNB] = CSRC_NB; } \
447                         if(NBFLAG(dNB)&CFBnd) { m[dST] = RAC(ccel,dNB); } else { m[dST] = CSRC_ST; } \
448                         if(NBFLAG(dNT)&CFBnd) { m[dSB] = RAC(ccel,dNT); } else { m[dSB] = CSRC_SB; } \
449                         if(NBFLAG(dWB)&CFBnd) { m[dET] = RAC(ccel,dWB); } else { m[dET] = CSRC_ET; } \
450                         if(NBFLAG(dWT)&CFBnd) { m[dEB] = RAC(ccel,dWT); } else { m[dEB] = CSRC_EB; } \
451                         if(NBFLAG(dEB)&CFBnd) { m[dWT] = RAC(ccel,dEB); } else { m[dWT] = CSRC_WT; } \
452                         if(NBFLAG(dET)&CFBnd) { m[dWB] = RAC(ccel,dET); } else { m[dWB] = CSRC_WB; } \
453                 } 
454
455
456
457 #define COLL_CALCULATE_DFEQ(dstarray) \
458                         dstarray[dN ] = EQN ; dstarray[dS ] = EQS ; \
459                         dstarray[dE ] = EQE ; dstarray[dW ] = EQW ; \
460                         dstarray[dT ] = EQT ; dstarray[dB ] = EQB ; \
461                         dstarray[dNE] = EQNE; dstarray[dNW] = EQNW; dstarray[dSE] = EQSE; dstarray[dSW] = EQSW; \
462                         dstarray[dNT] = EQNT; dstarray[dNB] = EQNB; dstarray[dST] = EQST; dstarray[dSB] = EQSB; \
463                         dstarray[dET] = EQET; dstarray[dEB] = EQEB; dstarray[dWT] = EQWT; dstarray[dWB] = EQWB; 
464 #define COLL_CALCULATE_NONEQTENSOR(csolev, srcArray ) \
465                         lcsmqadd  = (srcArray##NE - lcsmeq[ dNE ]); \
466                         lcsmqadd -= (srcArray##NW - lcsmeq[ dNW ]); \
467                         lcsmqadd -= (srcArray##SE - lcsmeq[ dSE ]); \
468                         lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
469                         lcsmqo = (lcsmqadd*    lcsmqadd); \
470                         lcsmqadd  = (srcArray##ET - lcsmeq[  dET ]); \
471                         lcsmqadd -= (srcArray##EB - lcsmeq[  dEB ]); \
472                         lcsmqadd -= (srcArray##WT - lcsmeq[  dWT ]); \
473                         lcsmqadd += (srcArray##WB - lcsmeq[  dWB ]); \
474                         lcsmqo += (lcsmqadd*    lcsmqadd); \
475                         lcsmqadd  = (srcArray##NT - lcsmeq[  dNT ]); \
476                         lcsmqadd -= (srcArray##NB - lcsmeq[  dNB ]); \
477                         lcsmqadd -= (srcArray##ST - lcsmeq[  dST ]); \
478                         lcsmqadd += (srcArray##SB - lcsmeq[  dSB ]); \
479                         lcsmqo += (lcsmqadd*    lcsmqadd); \
480                         lcsmqo *= 2.0; \
481                         lcsmqadd  = (srcArray##E  -  lcsmeq[ dE  ]); \
482                         lcsmqadd += (srcArray##W  -  lcsmeq[ dW  ]); \
483                         lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
484                         lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
485                         lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
486                         lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
487                         lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
488                         lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
489                         lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
490                         lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
491                         lcsmqo += (lcsmqadd*    lcsmqadd); \
492                         lcsmqadd  = (srcArray##N  -  lcsmeq[ dN  ]); \
493                         lcsmqadd += (srcArray##S  -  lcsmeq[ dS  ]); \
494                         lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
495                         lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
496                         lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
497                         lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
498                         lcsmqadd += (srcArray##NT  - lcsmeq[ dNT ]); \
499                         lcsmqadd += (srcArray##NB  - lcsmeq[ dNB ]); \
500                         lcsmqadd += (srcArray##ST  - lcsmeq[ dST ]); \
501                         lcsmqadd += (srcArray##SB  - lcsmeq[ dSB ]); \
502                         lcsmqo += (lcsmqadd*    lcsmqadd); \
503                         lcsmqadd  = (srcArray##T  -  lcsmeq[ dT  ]); \
504                         lcsmqadd += (srcArray##B  -  lcsmeq[ dB  ]); \
505                         lcsmqadd += (srcArray##NT -  lcsmeq[ dNT ]); \
506                         lcsmqadd += (srcArray##NB -  lcsmeq[ dNB ]); \
507                         lcsmqadd += (srcArray##ST -  lcsmeq[ dST ]); \
508                         lcsmqadd += (srcArray##SB -  lcsmeq[ dSB ]); \
509                         lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
510                         lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
511                         lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
512                         lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
513                         lcsmqo += (lcsmqadd*    lcsmqadd); \
514                         lcsmqo = sqrt(lcsmqo); /* FIXME check effect of sqrt*/ \
515
516 //                      COLL_CALCULATE_CSMOMEGAVAL(csolev, lcsmomega); 
517
518 // careful - need lcsmqo 
519 #define COLL_CALCULATE_CSMOMEGAVAL(csolev, dstomega ) \
520                         dstomega =  1.0/\
521                                         ( 3.0*( mLevel[(csolev)].lcnu+mLevel[(csolev)].lcsmago_sqr*(\
522                                                         -mLevel[(csolev)].lcnu + sqrt( mLevel[(csolev)].lcnu*mLevel[(csolev)].lcnu + 18.0*mLevel[(csolev)].lcsmago_sqr* lcsmqo ) \
523                                                         / (6.0*mLevel[(csolev)].lcsmago_sqr)) \
524                                                 ) +0.5 ); 
525
526 #define DEFAULT_COLLIDE_LES \
527                         rho = + MSRC_C  + MSRC_N  \
528                                 + MSRC_S  + MSRC_E  \
529                                 + MSRC_W  + MSRC_T  \
530                                 + MSRC_B  + MSRC_NE \
531                                 + MSRC_NW + MSRC_SE \
532                                 + MSRC_SW + MSRC_NT \
533                                 + MSRC_NB + MSRC_ST \
534                                 + MSRC_SB + MSRC_ET \
535                                 + MSRC_EB + MSRC_WT \
536                                 + MSRC_WB; \
537                         \
538                         ux += MSRC_E - MSRC_W \
539                                 + MSRC_NE - MSRC_NW \
540                                 + MSRC_SE - MSRC_SW \
541                                 + MSRC_ET + MSRC_EB \
542                                 - MSRC_WT - MSRC_WB ;  \
543                         \
544                         uy += MSRC_N - MSRC_S \
545                                 + MSRC_NE + MSRC_NW \
546                                 - MSRC_SE - MSRC_SW \
547                                 + MSRC_NT + MSRC_NB \
548                                 - MSRC_ST - MSRC_SB ;  \
549                         \
550                         uz += MSRC_T - MSRC_B \
551                                 + MSRC_NT - MSRC_NB \
552                                 + MSRC_ST - MSRC_SB \
553                                 + MSRC_ET - MSRC_EB \
554                                 + MSRC_WT - MSRC_WB ;  \
555                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
556                         COLL_CALCULATE_DFEQ(lcsmeq); \
557                         COLL_CALCULATE_NONEQTENSOR(lev, MSRC_)\
558                         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
559                         CSMOMEGA_STATS(lev,lcsmomega); \
560                         \
561                         RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
562                         \
563                         RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ]; \
564                         RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ]; \
565                         RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
566                         RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ]; \
567                         RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ]; \
568                         RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
569                         \
570                         RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
571                         RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
572                         RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
573                         RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
574                         RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
575                         RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
576                         RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
577                         RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
578                         RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
579                         RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
580                         RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
581                         RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; 
582
583 #define DEFAULT_COLLIDE_NOLES \
584                         rho = + MSRC_C  + MSRC_N  \
585                                 + MSRC_S  + MSRC_E  \
586                                 + MSRC_W  + MSRC_T  \
587                                 + MSRC_B  + MSRC_NE \
588                                 + MSRC_NW + MSRC_SE \
589                                 + MSRC_SW + MSRC_NT \
590                                 + MSRC_NB + MSRC_ST \
591                                 + MSRC_SB + MSRC_ET \
592                                 + MSRC_EB + MSRC_WT \
593                                 + MSRC_WB; \
594                         \
595                         ux += MSRC_E - MSRC_W \
596                                 + MSRC_NE - MSRC_NW \
597                                 + MSRC_SE - MSRC_SW \
598                                 + MSRC_ET + MSRC_EB \
599                                 - MSRC_WT - MSRC_WB ;  \
600                         \
601                         uy += MSRC_N - MSRC_S \
602                                 + MSRC_NE + MSRC_NW \
603                                 - MSRC_SE - MSRC_SW \
604                                 + MSRC_NT + MSRC_NB \
605                                 - MSRC_ST - MSRC_SB ;  \
606                         \
607                         uz += MSRC_T - MSRC_B \
608                                 + MSRC_NT - MSRC_NB \
609                                 + MSRC_ST - MSRC_SB \
610                                 + MSRC_ET - MSRC_EB \
611                                 + MSRC_WT - MSRC_WB ;  \
612                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
613                         \
614                         RAC(tcel,dC ) = (1.0-OMEGA(lev))*MSRC_C  + OMEGA(lev)*EQC ; \
615                         \
616                         RAC(tcel,dN ) = (1.0-OMEGA(lev))*MSRC_N  + OMEGA(lev)*EQN ; \
617                         RAC(tcel,dS ) = (1.0-OMEGA(lev))*MSRC_S  + OMEGA(lev)*EQS ; \
618                         RAC(tcel,dE ) = (1.0-OMEGA(lev))*MSRC_E  + OMEGA(lev)*EQE ; \
619                         RAC(tcel,dW ) = (1.0-OMEGA(lev))*MSRC_W  + OMEGA(lev)*EQW ; \
620                         RAC(tcel,dT ) = (1.0-OMEGA(lev))*MSRC_T  + OMEGA(lev)*EQT ; \
621                         RAC(tcel,dB ) = (1.0-OMEGA(lev))*MSRC_B  + OMEGA(lev)*EQB ; \
622                         \
623                         RAC(tcel,dNE) = (1.0-OMEGA(lev))*MSRC_NE + OMEGA(lev)*EQNE; \
624                         RAC(tcel,dNW) = (1.0-OMEGA(lev))*MSRC_NW + OMEGA(lev)*EQNW; \
625                         RAC(tcel,dSE) = (1.0-OMEGA(lev))*MSRC_SE + OMEGA(lev)*EQSE; \
626                         RAC(tcel,dSW) = (1.0-OMEGA(lev))*MSRC_SW + OMEGA(lev)*EQSW; \
627                         RAC(tcel,dNT) = (1.0-OMEGA(lev))*MSRC_NT + OMEGA(lev)*EQNT; \
628                         RAC(tcel,dNB) = (1.0-OMEGA(lev))*MSRC_NB + OMEGA(lev)*EQNB; \
629                         RAC(tcel,dST) = (1.0-OMEGA(lev))*MSRC_ST + OMEGA(lev)*EQST; \
630                         RAC(tcel,dSB) = (1.0-OMEGA(lev))*MSRC_SB + OMEGA(lev)*EQSB; \
631                         RAC(tcel,dET) = (1.0-OMEGA(lev))*MSRC_ET + OMEGA(lev)*EQET; \
632                         RAC(tcel,dEB) = (1.0-OMEGA(lev))*MSRC_EB + OMEGA(lev)*EQEB; \
633                         RAC(tcel,dWT) = (1.0-OMEGA(lev))*MSRC_WT + OMEGA(lev)*EQWT; \
634                         RAC(tcel,dWB) = (1.0-OMEGA(lev))*MSRC_WB + OMEGA(lev)*EQWB; 
635
636
637
638 #define OPTIMIZED_STREAMCOLLIDE_LES \
639                         /* only surrounded by fluid cells...!, so safe streaming here... */ \
640                         m[dC ] = CSRC_C ; \
641                         m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
642                         m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
643                         m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
644                         m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
645                         m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
646                         m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
647                         \
648                         rho = MSRC_C  + MSRC_N + MSRC_S  + MSRC_E + MSRC_W  + MSRC_T  \
649                                 + MSRC_B  + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT \
650                                 + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; \
651                         ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW \
652                                 + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB + mLevel[lev].gravity[0];  \
653                         uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW \
654                                 + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB + mLevel[lev].gravity[1];  \
655                         uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB \
656                                 + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB + mLevel[lev].gravity[2];  \
657                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
658                         COLL_CALCULATE_DFEQ(lcsmeq); \
659                         COLL_CALCULATE_NONEQTENSOR(lev, MSRC_) \
660                         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
661                         CSMOMEGA_STATS(lev,lcsmomega); \
662                         \
663                         RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
664                         RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ];  \
665                         RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ];  \
666                         RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
667                         RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ];  \
668                         RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ];  \
669                         RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
670                         \
671                         RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE];  \
672                         RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW];  \
673                         RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE];  \
674                         RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
675                         \
676                         RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT];  \
677                         RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB];  \
678                         RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST];  \
679                         RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
680                         \
681                         RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET];  \
682                         RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
683                         RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT];  \
684                         RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB];  \
685
686 #define OPTIMIZED_STREAMCOLLIDE_UNUSED \
687                         /* only surrounded by fluid cells...!, so safe streaming here... */ \
688                         rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T  \
689                                 + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
690                                 + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
691                         ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
692                                 + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB + mLevel[lev].gravity[0];  \
693                         uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
694                                 + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB + mLevel[lev].gravity[1];  \
695                         uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
696                                 + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB + mLevel[lev].gravity[2];  \
697                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
698                         COLL_CALCULATE_DFEQ(lcsmeq); \
699                         COLL_CALCULATE_NONEQTENSOR(lev, CSRC_) \
700                         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
701                         \
702                         RAC(tcel,dC ) = (1.0-lcsmomega)*CSRC_C  + lcsmomega*EQC ; \
703                         RAC(tcel,dN ) = (1.0-lcsmomega)*CSRC_N  + lcsmomega*lcsmeq[ dN ];  \
704                         RAC(tcel,dS ) = (1.0-lcsmomega)*CSRC_S  + lcsmomega*lcsmeq[ dS ];  \
705                         RAC(tcel,dE ) = (1.0-lcsmomega)*CSRC_E  + lcsmomega*lcsmeq[ dE ]; \
706                         RAC(tcel,dW ) = (1.0-lcsmomega)*CSRC_W  + lcsmomega*lcsmeq[ dW ];  \
707                         RAC(tcel,dT ) = (1.0-lcsmomega)*CSRC_T  + lcsmomega*lcsmeq[ dT ];  \
708                         RAC(tcel,dB ) = (1.0-lcsmomega)*CSRC_B  + lcsmomega*lcsmeq[ dB ]; \
709                         \
710                         RAC(tcel,dNE) = (1.0-lcsmomega)*CSRC_NE + lcsmomega*lcsmeq[ dNE];  \
711                         RAC(tcel,dNW) = (1.0-lcsmomega)*CSRC_NW + lcsmomega*lcsmeq[ dNW];  \
712                         RAC(tcel,dSE) = (1.0-lcsmomega)*CSRC_SE + lcsmomega*lcsmeq[ dSE];  \
713                         RAC(tcel,dSW) = (1.0-lcsmomega)*CSRC_SW + lcsmomega*lcsmeq[ dSW]; \
714                         \
715                         RAC(tcel,dNT) = (1.0-lcsmomega)*CSRC_NT + lcsmomega*lcsmeq[ dNT];  \
716                         RAC(tcel,dNB) = (1.0-lcsmomega)*CSRC_NB + lcsmomega*lcsmeq[ dNB];  \
717                         RAC(tcel,dST) = (1.0-lcsmomega)*CSRC_ST + lcsmomega*lcsmeq[ dST];  \
718                         RAC(tcel,dSB) = (1.0-lcsmomega)*CSRC_SB + lcsmomega*lcsmeq[ dSB]; \
719                         \
720                         RAC(tcel,dET) = (1.0-lcsmomega)*CSRC_ET + lcsmomega*lcsmeq[ dET];  \
721                         RAC(tcel,dEB) = (1.0-lcsmomega)*CSRC_EB + lcsmomega*lcsmeq[ dEB]; \
722                         RAC(tcel,dWT) = (1.0-lcsmomega)*CSRC_WT + lcsmomega*lcsmeq[ dWT];  \
723                         RAC(tcel,dWB) = (1.0-lcsmomega)*CSRC_WB + lcsmomega*lcsmeq[ dWB];  \
724
725 #define OPTIMIZED_STREAMCOLLIDE_NOLES \
726                         /* only surrounded by fluid cells...!, so safe streaming here... */ \
727                         rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T  \
728                                 + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
729                                 + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
730                         ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
731                                 + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB + mLevel[lev].gravity[0];  \
732                         uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
733                                 + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB + mLevel[lev].gravity[1];  \
734                         uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
735                                 + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB + mLevel[lev].gravity[2];  \
736                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
737                         RAC(tcel,dC ) = (1.0-OMEGA(lev))*CSRC_C  + OMEGA(lev)*EQC ; \
738                         RAC(tcel,dN ) = (1.0-OMEGA(lev))*CSRC_N  + OMEGA(lev)*EQN ;  \
739                         RAC(tcel,dS ) = (1.0-OMEGA(lev))*CSRC_S  + OMEGA(lev)*EQS ;  \
740                         RAC(tcel,dE ) = (1.0-OMEGA(lev))*CSRC_E  + OMEGA(lev)*EQE ; \
741                         RAC(tcel,dW ) = (1.0-OMEGA(lev))*CSRC_W  + OMEGA(lev)*EQW ;  \
742                         RAC(tcel,dT ) = (1.0-OMEGA(lev))*CSRC_T  + OMEGA(lev)*EQT ;  \
743                         RAC(tcel,dB ) = (1.0-OMEGA(lev))*CSRC_B  + OMEGA(lev)*EQB ; \
744                          \
745                         RAC(tcel,dNE) = (1.0-OMEGA(lev))*CSRC_NE + OMEGA(lev)*EQNE;  \
746                         RAC(tcel,dNW) = (1.0-OMEGA(lev))*CSRC_NW + OMEGA(lev)*EQNW;  \
747                         RAC(tcel,dSE) = (1.0-OMEGA(lev))*CSRC_SE + OMEGA(lev)*EQSE;  \
748                         RAC(tcel,dSW) = (1.0-OMEGA(lev))*CSRC_SW + OMEGA(lev)*EQSW; \
749                          \
750                         RAC(tcel,dNT) = (1.0-OMEGA(lev))*CSRC_NT + OMEGA(lev)*EQNT;  \
751                         RAC(tcel,dNB) = (1.0-OMEGA(lev))*CSRC_NB + OMEGA(lev)*EQNB;  \
752                         RAC(tcel,dST) = (1.0-OMEGA(lev))*CSRC_ST + OMEGA(lev)*EQST;  \
753                         RAC(tcel,dSB) = (1.0-OMEGA(lev))*CSRC_SB + OMEGA(lev)*EQSB; \
754                          \
755                         RAC(tcel,dET) = (1.0-OMEGA(lev))*CSRC_ET + OMEGA(lev)*EQET;  \
756                         RAC(tcel,dEB) = (1.0-OMEGA(lev))*CSRC_EB + OMEGA(lev)*EQEB; \
757                         RAC(tcel,dWT) = (1.0-OMEGA(lev))*CSRC_WT + OMEGA(lev)*EQWT;  \
758                         RAC(tcel,dWB) = (1.0-OMEGA(lev))*CSRC_WB + OMEGA(lev)*EQWB;  \
759
760
761 // debug version1
762 #define STREAMCHECK(ni,nj,nk,nl) 
763 #define COLLCHECK
764 #define OPTIMIZED_STREAMCOLLIDE_DEBUG \
765                         m[0] = RAC(ccel,0); \
766                         FORDF1 { /* df0 is set later on... */ \
767                                 if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); D::mPanic=1;  \
768                                 } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
769                                 STREAMCHECK(i+D::dfVecX[D::dfInv[l]], j+D::dfVecY[D::dfInv[l]],k+D::dfVecZ[D::dfInv[l]], l); \
770                         }   \
771                         rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
772                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
773                         D::collideArrays( m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago , &mDebugOmegaRet  ); \
774                         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
775                         FORDF0 { RAC(tcel,l) = m[l]; } \
776                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
777                         COLLCHECK;
778
779
780
781 // more debugging
782 /*DEBUG \
783                         m[0] = RAC(ccel,0); \
784                         FORDF1 { \
785                                 if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); D::mPanic=1;  \
786                                 } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
787                         }   \
788 f__printf(stderr,"QSDM at %d,%d,%d  lcsmqo=%25.15f, lcsmomega=%f \n", i,j,k, lcsmqo,lcsmomega ); \
789                         rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
790                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
791                         D::collideArrays( m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago  , &mDebugOmegaRet ); \
792                         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
793                         */
794 #if USE_LES==1
795 #define DEFAULT_COLLIDE DEFAULT_COLLIDE_LES
796 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_LES
797 #else 
798 #define DEFAULT_COLLIDE DEFAULT_COLLIDE_NOLES
799 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_NOLES
800 #endif
801
802 #endif
803
804 #define USQRMAXCHECK(Cusqr,Cux,Cuy,Cuz,  CmMaxVlen,CmMxvx,CmMxvy,CmMxvz) \
805                         if(Cusqr>CmMaxVlen) { \
806                                 CmMxvx = Cux; CmMxvy = Cuy; CmMxvz = Cuz; CmMaxVlen = Cusqr; \
807                         } /* stats */ 
808
809
810 // iso value defines
811 // border for marching cubes
812 #define ISOCORR 3
813
814
815 // helper for comparing floats with epsilon
816 #define GFX_FLOATNEQ(x,y) ( ABS((x)-(y)) > (VECTOR_EPSILON) )
817 #define LBM_FLOATNEQ(x,y) ( ABS((x)-(y)) > (10.0*LBM_EPSILON) )
818
819
820 // macros for loops over all DFs
821 #define FORDF0 for(int l= 0; l< LBM_DFNUM; ++l)
822 #define FORDF1 for(int l= 1; l< LBM_DFNUM; ++l)
823
824 /*****************************************************************************/
825 /*! cell access classes */
826 // WARNING - can be shadowed by class from lbmstdsolver.h 's
827 template<typename D>
828 class UniformFsgrCellIdentifier : 
829         public CellIdentifierInterface 
830 {
831         public:
832                 //! which grid level?
833                 int level;
834                 //! location in grid
835                 int x,y,z;
836
837                 //! reset constructor
838                 UniformFsgrCellIdentifier() :
839                         x(0), y(0), z(0) { };
840
841                 // implement CellIdentifierInterface
842                 virtual string getAsString() {
843                         std::ostringstream ret;
844                         ret <<"{ i"<<x<<",j"<<y;
845                         if(D::cDimension>2) ret<<",k"<<z;
846                         ret <<" }";
847                         return ret.str();
848                 }
849
850                 virtual bool equal(CellIdentifierInterface* other) {
851                         //UniformFsgrCellIdentifier<D> *cid = dynamic_cast<UniformFsgrCellIdentifier<D> *>( other );
852                         UniformFsgrCellIdentifier<D> *cid = (UniformFsgrCellIdentifier<D> *)( other );
853                         if(!cid) return false;
854                         if( x==cid->x && y==cid->y && z==cid->z && level==cid->level ) return true;
855                         return false;
856                 }
857 };
858
859 //! information needed for each level in the simulation
860 class FsgrLevelData {
861 public:
862         int id; // level number
863
864         //! node size on this level (geometric, in world coordinates, not simulation units!) 
865         LbmFloat nodeSize;
866         //! node size on this level in simulation units 
867         LbmFloat simCellSize;
868         //! quadtree node relaxation parameter 
869         LbmFloat omega;
870         //! size this level was advanced to 
871         LbmFloat time;
872         //! size of a single lbm step in time units on this level 
873         LbmFloat stepsize;
874         //! step count
875         int lsteps;
876         //! gravity force for this level
877         LbmVec gravity;
878         //! level array 
879         LbmFloat *mprsCells[2];
880         CellFlagType *mprsFlags[2];
881
882         //! smago params and precalculated values
883         LbmFloat lcsmago;
884         LbmFloat lcsmago_sqr;
885         LbmFloat lcnu;
886
887         // LES statistics per level
888         double avgOmega;
889         double avgOmegaCnt;
890
891         //! current set of dist funcs 
892         int setCurr;
893         //! target/other set of dist funcs 
894         int setOther;
895
896         //! mass&volume for this level
897         LbmFloat lmass;
898         LbmFloat lvolume;
899         LbmFloat lcellfactor;
900
901         //! local storage of mSizes
902         int lSizex, lSizey, lSizez;
903         int lOffsx, lOffsy, lOffsz;
904
905 };
906
907
908
909 /*****************************************************************************/
910 /*! class for solving a LBM problem */
911 template<class D>
912 class LbmFsgrSolver : 
913         public /*? virtual */ D // this means, the solver is a lbmData object and implements the lbmInterface
914 {
915
916         public:
917                 //! Constructor 
918                 LbmFsgrSolver();
919                 //! Destructor 
920                 virtual ~LbmFsgrSolver();
921                 //! id string of solver
922                 virtual string getIdString() { return string("FsgrSolver[") + D::getIdString(); }
923
924                 //! initilize variables fom attribute list 
925                 virtual void parseAttrList();
926                 //! Initialize omegas and forces on all levels (for init/timestep change)
927                 void initLevelOmegas();
928                 //! finish the init with config file values (allocate arrays...) 
929                 virtual bool initialize( ntlTree* /*tree*/, vector<ntlGeometryObject*>* /*objects*/ );
930
931 #if LBM_USE_GUI==1
932                 //! show simulation info (implement SimulationObject pure virtual func)
933                 virtual void debugDisplay(fluidDispSettings *set);
934 #endif
935                 
936                 
937                 // implement CellIterator<UniformFsgrCellIdentifier> interface
938                 typedef UniformFsgrCellIdentifier<typename D::LbmCellContents> stdCellId;
939                 virtual CellIdentifierInterface* getFirstCell( );
940                 virtual void advanceCell( CellIdentifierInterface* );
941                 virtual bool noEndCell( CellIdentifierInterface* );
942                 virtual void deleteCellIterator( CellIdentifierInterface** );
943                 virtual CellIdentifierInterface* getCellAt( ntlVec3Gfx pos );
944                 virtual int        getCellSet      ( CellIdentifierInterface* );
945                 virtual ntlVec3Gfx getCellOrigin   ( CellIdentifierInterface* );
946                 virtual ntlVec3Gfx getCellSize     ( CellIdentifierInterface* );
947                 virtual int        getCellLevel    ( CellIdentifierInterface* );
948                 virtual LbmFloat   getCellDensity  ( CellIdentifierInterface* ,int set);
949                 virtual LbmVec     getCellVelocity ( CellIdentifierInterface* ,int set);
950                 virtual LbmFloat   getCellDf       ( CellIdentifierInterface* ,int set, int dir);
951                 virtual LbmFloat   getCellMass     ( CellIdentifierInterface* ,int set);
952                 virtual LbmFloat   getCellFill     ( CellIdentifierInterface* ,int set);
953                 virtual CellFlagType getCellFlag   ( CellIdentifierInterface* ,int set);
954                 virtual LbmFloat   getEquilDf      ( int );
955                 virtual int        getDfNum        ( );
956                 // convert pointers
957                 stdCellId* convertBaseCidToStdCid( CellIdentifierInterface* basecid);
958
959                 //! perform geometry init (if switched on) 
960                 bool initGeometryFlags();
961                 //! init part for all freesurface testcases 
962                 void initFreeSurfaces();
963                 //! init density gradient if enabled
964                 void initStandingFluidGradient();
965
966                 /*! init a given cell with flag, density, mass and equilibrium dist. funcs */
967                 inline void initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass);
968                 void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel);
969
970                 /*! perform a single LBM step */
971                 virtual void step() { stepMain(); }
972                 void stepMain();
973                 void fineAdvance();
974                 void coarseAdvance(int lev);
975                 void coarseCalculateFluxareas(int lev);
976                 // coarsen a given level (returns true if sth. was changed)
977                 bool performCoarsening(int lev);
978                 bool performRefinement(int lev);
979                 //void oarseInterpolateToFineSpaceTime(int lev,LbmFloat t);
980                 void interpolateFineFromCoarse(int lev,LbmFloat t);
981                 void coarseRestrictFromFine(int lev);
982
983                 /*! init particle positions */
984                 virtual int initParticles(ParticleTracer *partt);
985                 /*! move all particles */
986                 virtual void advanceParticles(ParticleTracer *partt );
987
988
989                 /*! debug object display (used e.g. for preview surface) */
990                 virtual vector<ntlGeometryObject*> getDebugObjects();
991
992                 //! access the fillfrac field (for viz)
993                 float getFillFrac(int i, int j, int k) {
994                         return QCELL(mMaxRefine, i,j,k,mLevel[mMaxRefine].setOther, dFfrac);
995                 }
996
997                 //! retrieve the fillfrac field ready to display
998                 void getIsofieldWeighted(float *iso);
999                 void getIsofield(float *iso){ return getIsofieldWeighted(iso); }
1000                 //! for raytracing, preprocess
1001                 void prepareVisualization( void );
1002
1003                 // rt interface
1004                 void addDrop(bool active, float mx, float my);
1005                 void initDrop(float mx, float my);
1006                 void printCellStats();
1007                 int checkGfxEndTime(); // {return 9;};
1008                 //! get gfx geo setup id
1009                 int getGfxGeoSetup() { return mGfxGeoSetup; }
1010
1011                 /*! type for cells */
1012                 typedef typename D::LbmCell LbmCell;
1013                 
1014         protected:
1015
1016                 //! internal quick print function (for debugging) 
1017                 void printLbmCell(int level, int i, int j, int k,int set);
1018                 // debugging use CellIterator interface to mark cell
1019                 void debugMarkCellCall(int level, int vi,int vj,int vk);
1020
1021                 void mainLoop(int lev);
1022                 void adaptTimestep();
1023                 //! flag reinit step - always works on finest grid!
1024                 void reinitFlags( int workSet );
1025                 //! mass dist weights
1026                 LbmFloat getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l);
1027                 //! add point to mListNewInter list
1028                 inline void addToNewInterList( int ni, int nj, int nk );        
1029                 //! cell is interpolated from coarse level (inited into set, source sets are determined by t)
1030                 inline void interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet,bool markNbs);
1031
1032                 //! minimal and maximal z-coords (for 2D/3D loops)
1033                 int getForZMinBnd(int lev) { 
1034                         return 0; 
1035                 }
1036                 int getForZMaxBnd(int lev) { 
1037                         if(D::cDimension==2) return 1;
1038                         //return D::getSizeZ()-0; 
1039                         return mLevel[lev].lSizez -0;
1040                 }
1041                 int getForZMin1(int lev)   { 
1042                         if(D::cDimension==2) return 0;
1043                         return 1; 
1044                 }
1045                 int getForZMax1(int lev)   { 
1046                         if(D::cDimension==2) return 1;
1047                         //return D::getSizeZ()-1; 
1048                         return mLevel[lev].lSizez -1;
1049                 }
1050
1051
1052                 // member vars
1053
1054                 //! mass calculated during streaming step
1055                 LbmFloat mCurrentMass;
1056                 LbmFloat mCurrentVolume;
1057                 LbmFloat mInitialMass;
1058
1059                 //! count problematic cases, that occured so far...
1060                 int mNumProblems;
1061
1062                 // average mlsups, count how many so far...
1063                 double mAvgMLSUPS;
1064                 double mAvgMLSUPSCnt;
1065
1066                 //! Mcubes object for surface reconstruction 
1067                 IsoSurface *mpPreviewSurface;
1068                 int mLoopSubdivs;
1069                 float mSmoothSurface;
1070                 float mSmoothNormals;
1071                 
1072                 //! use time adaptivity? 
1073                 bool mTimeAdap;
1074
1075                 //! output surface preview? if >0 yes, and use as reduzed size 
1076                 int mOutputSurfacePreview;
1077                 LbmFloat mPreviewFactor;
1078                 //! fluid vol height
1079                 LbmFloat mFVHeight;
1080                 LbmFloat mFVArea;
1081                 bool mUpdateFVHeight;
1082
1083                 //! require some geo setup from the viz?
1084                 int mGfxGeoSetup;
1085                 //! force quit for gfx
1086                 LbmFloat mGfxEndTime;
1087                 //! smoother surface initialization?
1088                 int mInitSurfaceSmoothing;
1089
1090                 int mTimestepReduceLock;
1091                 int mTimeSwitchCounts;
1092                 //! total simulation time so far 
1093                 LbmFloat mSimulationTime;
1094                 //! smallest and largest step size so far 
1095                 LbmFloat mMinStepTime, mMaxStepTime;
1096                 //! track max. velocity
1097                 LbmFloat mMxvx, mMxvy, mMxvz, mMaxVlen;
1098
1099                 //! list of the cells to empty at the end of the step 
1100                 vector<LbmPoint> mListEmpty;
1101                 //! list of the cells to make fluid at the end of the step 
1102                 vector<LbmPoint> mListFull;
1103                 //! list of new interface cells to init
1104         vector<LbmPoint> mListNewInter;
1105                 //! class for handling redist weights in reinit flag function
1106                 class lbmFloatSet {
1107                         public:
1108                                 LbmFloat val[dTotalNum];
1109                 };
1110                 //! normalized vectors for all neighboring cell directions (for e.g. massdweight calc)
1111                 LbmVec mDvecNrm[27];
1112                 
1113                 
1114                 //! debugging
1115                 bool checkSymmetry(string idstring);
1116                 //! symmetric init?
1117                 //bool mStartSymm;
1118                 //! kepp track of max/min no. of filled cells
1119                 int mMaxNoCells, mMinNoCells;
1120                 long long int mAvgNumUsedCells;
1121
1122                 //! for interactive - how to drop drops?
1123                 int mDropMode;
1124                 LbmFloat mDropSize;
1125                 LbmVec mDropSpeed;
1126                 //! dropping variables
1127                 bool mDropping;
1128                 LbmFloat mDropX, mDropY;
1129                 LbmFloat mDropHeight;
1130
1131                 //! get isofield weights
1132                 int mIsoWeightMethod;
1133                 float mIsoWeight[27];
1134
1135                 // grid coarsening vars
1136                 
1137                 /*! vector for the data for each level */
1138 #               define MAX_LEV 5
1139                 FsgrLevelData mLevel[MAX_LEV];
1140
1141                 /*! minimal and maximal refinement levels */
1142                 int mMaxRefine;
1143
1144                 /*! df scale factors for level up/down */
1145                 LbmFloat mDfScaleUp, mDfScaleDown;
1146
1147                 /*! precomputed cell area values */
1148                 LbmFloat mFsgrCellArea[27];
1149
1150                 /*! LES C_smago paramter for finest grid */
1151                 float mInitialCsmago;
1152                 float mCsmagoRefineMultiplier;
1153                 /*! LES stats for non OPT3D */
1154                 LbmFloat mDebugOmegaRet;
1155
1156                 //! fluid stats
1157                 int mNumInterdCells;
1158                 int mNumInvIfCells;
1159                 int mNumInvIfTotal;
1160
1161                 //! debug function to disable standing f init
1162                 int mDisableStandingFluidInit;
1163                 //! debug function to force tadap syncing
1164                 int mForceTadapRefine;
1165
1166
1167 #               if FSGR_STRICT_DEBUG==1
1168 #               define STRICT_EXIT *((int *)0)=0;
1169                 int debLBMGI(int level, int ii,int ij,int ik, int is) {
1170                         if(level <  0){ errMsg("LbmStrict::debLBMGI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
1171                         if(level >  mMaxRefine){ errMsg("LbmStrict::debLBMGI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
1172                         
1173                         if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1174                         if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1175                         if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1176                         if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1177                         if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1178                         if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1179                         if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1180                         if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1181                         return _LBMGI(level, ii,ij,ik, is);
1182                 };
1183                 CellFlagType& debRFLAG(int level, int xx,int yy,int zz,int set){
1184                         return _RFLAG(level, xx,yy,zz,set);   
1185                 };
1186                 CellFlagType& debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir) {
1187                         if(dir<0)         { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1188                         // warning might access all spatial nbs
1189                         if(dir>D::cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1190                         return _RFLAG_NB(level, xx,yy,zz,set, dir);
1191                 };
1192                 CellFlagType& debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir) {
1193                         if(dir<0)         { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1194                         if(dir>D::cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1195                         return _RFLAG_NBINV(level, xx,yy,zz,set, dir);
1196                 };
1197                 int debLBMQI(int level, int ii,int ij,int ik, int is, int l) {
1198                         if(level <  0){ errMsg("LbmStrict::debLBMQI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
1199                         if(level >  mMaxRefine){ errMsg("LbmStrict::debLBMQI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
1200                         
1201                         if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1202                         if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1203                         if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1204                         if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1205                         if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1206                         if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1207                         if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1208                         if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1209                         if(l<0)        { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; }
1210                         if(l>D::cDfNum){  // dFfrac is an exception
1211                         if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } }
1212 #if COMPRESSGRIDS==1
1213                         //if((!D::mInitDone) && (is!=mLevel[level].setCurr)){ STRICT_EXIT; } // COMPRT debug
1214 #endif // COMPRESSGRIDS==1
1215                         return _LBMQI(level, ii,ij,ik, is, l);
1216                 };
1217                 LbmFloat& debQCELL(int level, int xx,int yy,int zz,int set,int l) {
1218                         //errMsg("LbmStrict","debQCELL debug: l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" l"<<l<<" index"<<LBMGI(level, xx,yy,zz,set)); 
1219                         return _QCELL(level, xx,yy,zz,set,l);
1220                 };
1221                 LbmFloat& debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l) {
1222                         if(dir<0)        { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1223                         if(dir>D::cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1224                         return _QCELL_NB(level, xx,yy,zz,set, dir,l);
1225                 };
1226                 LbmFloat& debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l) {
1227                         if(dir<0)        { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1228                         if(dir>D::cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1229                         return _QCELL_NBINV(level, xx,yy,zz,set, dir,l);
1230                 };
1231                 LbmFloat* debRACPNT(int level,  int ii,int ij,int ik, int is ) {
1232                         return _RACPNT(level, ii,ij,ik, is );
1233                 };
1234                 LbmFloat& debRAC(LbmFloat* s,int l) {
1235                         if(l<0)        { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; }
1236                         if(l>dTotalNum){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } 
1237                         //if(l>D::cDfNum){ // dFfrac is an exception 
1238                         //if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } }
1239                         return _RAC(s,l);
1240                 };
1241 #               endif // FSGR_STRICT_DEBUG==1
1242 };
1243
1244
1245
1246 /*****************************************************************************/
1247 // compilation settings
1248
1249
1250 // loops over _all_ cells (including boundary layer)
1251 #define FSGR_FORIJK_BOUNDS(leveli) \
1252         for(int k= getForZMinBnd(leveli); k< getForZMaxBnd(leveli); ++k) \
1253    for(int j=0;j<mLevel[leveli].lSizey-0;++j) \
1254     for(int i=0;i<mLevel[leveli].lSizex-0;++i) \
1255         
1256 // loops over _only inner_ cells 
1257 #define FSGR_FORIJK1(leveli) \
1258         for(int k= getForZMin1(leveli); k< getForZMax1(leveli); ++k) \
1259    for(int j=1;j<mLevel[leveli].lSizey-1;++j) \
1260     for(int i=1;i<mLevel[leveli].lSizex-1;++i) \
1261
1262
1263 /******************************************************************************
1264  * Lbm Constructor
1265  *****************************************************************************/
1266 template<class D>
1267 LbmFsgrSolver<D>::LbmFsgrSolver() :
1268         D(),
1269         mCurrentMass(0.0), mCurrentVolume(0.0),
1270         mNumProblems(0), 
1271         mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0),
1272         mpPreviewSurface(NULL), 
1273         mLoopSubdivs(0), mSmoothSurface(0.0), mSmoothNormals(0.0),
1274         mTimeAdap(false), 
1275         mOutputSurfacePreview(0), mPreviewFactor(0.25),
1276         mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
1277         mGfxGeoSetup(0), mGfxEndTime(-1.0), mInitSurfaceSmoothing(0),
1278         mTimestepReduceLock(0),
1279         mTimeSwitchCounts(0), 
1280         mSimulationTime(0.0),
1281         mMinStepTime(0.0), mMaxStepTime(0.0),
1282         mMaxNoCells(0), mMinNoCells(0), mAvgNumUsedCells(0),
1283         mDropMode(1), mDropSize(0.15), mDropSpeed(0.0),
1284         mDropping(false),
1285         mDropX(0.0), mDropY(0.0), mDropHeight(0.8),
1286         mIsoWeightMethod(2),
1287         mMaxRefine(1), 
1288         mDfScaleUp(-1.0), mDfScaleDown(-1.0),
1289         mInitialCsmago(0.04), mCsmagoRefineMultiplier(8.0), mDebugOmegaRet(0.0),
1290         mNumInvIfTotal(0),
1291         mDisableStandingFluidInit(0),
1292         mForceTadapRefine(-1)
1293 {
1294   /* not much to do here... */
1295         D::mpIso = new IsoSurface( D::mIsoValue, false );
1296
1297   /* init equilibrium dist. func */
1298   LbmFloat rho=1.0;
1299   FORDF0 {
1300                 D::dfEquil[l] = D::getCollideEq( l,rho,  0.0, 0.0, 0.0);
1301   }
1302
1303         int odm = 0;
1304         for(int m=0; m<D::cDimension; m++) { 
1305                 for(int l=0; l<D::cDfNum; l++) { 
1306                         D::lesCoeffDiag[m][l] = 
1307                         D::lesCoeffOffdiag[m][l] = 0.0;
1308                 }
1309         }
1310         for(int m=0; m<D::cDimension; m++) { 
1311                 for(int n=0; n<D::cDimension; n++) { 
1312                         //LbmFloat qadd = 0.0;
1313                         for(int l=1; l<D::cDfNum; l++) { 
1314                                 LbmFloat em;
1315                                 switch(m) {
1316                                         case 0: em = D::dfDvecX[l]; break;
1317                                         case 1: em = D::dfDvecY[l]; break;
1318                                         case 2: em = D::dfDvecZ[l]; break;
1319                                         default: em = -1.0; errMsg("SMAGO","err m="<<m); exit(1);
1320                                 }
1321                                 LbmFloat en;
1322                                 switch(n) {
1323                                         case 0: en = D::dfDvecX[l]; break;
1324                                         case 1: en = D::dfDvecY[l]; break;
1325                                         case 2: en = D::dfDvecZ[l]; break;
1326                                         default: en = -1.0; errMsg("SMAGO","err n="<<n); exit(1);
1327                                 }
1328                                 const LbmFloat coeff = em*en;
1329                                 if(m==n) {
1330                                         D::lesCoeffDiag[m][l] = coeff;
1331                                         // f__printf(stderr,"QSMDEB: CF DIAG m:%d l:%d = %f \n", m,l, coeff);
1332                                 } else {
1333                                         if(m>n) {
1334                                                 D::lesCoeffOffdiag[odm][l] = coeff;
1335                                                 // f__printf(stderr,"QSMDEB: CF OFFDIAG odm:%d l:%d = %f \n", odm,l, coeff);
1336                                         }
1337                                 }
1338                         }
1339
1340
1341                         if(m==n) {
1342                         } else {
1343                                 if(m>n) odm++;
1344                         }
1345                 }
1346         }
1347
1348         mDvecNrm[0] = LbmVec(0.0);
1349   FORDF1 {
1350                 mDvecNrm[l] = getNormalized( 
1351                         LbmVec(D::dfDvecX[D::dfInv[l]], D::dfDvecY[D::dfInv[l]], D::dfDvecZ[D::dfInv[l]] ) 
1352                         ) * -1.0; 
1353         }
1354
1355         addDrop(false,0,0);
1356 }
1357
1358 /*****************************************************************************/
1359 /* Destructor */
1360 /*****************************************************************************/
1361 template<class D>
1362 LbmFsgrSolver<D>::~LbmFsgrSolver()
1363 {
1364   if(!D::mInitDone){ debugOut("LbmFsgrSolver::LbmFsgrSolver : not inited...",0); return; }
1365
1366 #if COMPRESSGRIDS==1
1367         delete mLevel[mMaxRefine].mprsCells[1];
1368         mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
1369 #endif // COMPRESSGRIDS==1
1370
1371         for(int i=0; i<=mMaxRefine; i++) {
1372                 for(int s=0; s<2; s++) {
1373                         if(mLevel[i].mprsCells[s]) delete [] mLevel[i].mprsCells[s];
1374                         if(mLevel[i].mprsFlags[s]) delete [] mLevel[i].mprsFlags[s];
1375                 }
1376         }
1377         delete D::mpIso;
1378         if(mpPreviewSurface) delete mpPreviewSurface;
1379
1380         // always output performance estimate
1381         debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5);
1382   if(!D::mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10);
1383 }
1384
1385
1386
1387
1388 /******************************************************************************
1389  * initilize variables fom attribute list 
1390  *****************************************************************************/
1391 template<class D>
1392 void 
1393 LbmFsgrSolver<D>::parseAttrList()
1394 {
1395         LbmSolverInterface::parseStdAttrList();
1396
1397         string matIso("default");
1398         matIso = D::mpAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false );
1399         D::mpIso->setMaterialName( matIso );
1400         mOutputSurfacePreview = D::mpAttrs->readInt("surfacepreview", mOutputSurfacePreview, "SimulationLbm","mOutputSurfacePreview", false );
1401         mTimeAdap = D::mpAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false );
1402
1403         mIsoWeightMethod= D::mpAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false );
1404         mInitSurfaceSmoothing = D::mpAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false );
1405         mLoopSubdivs = D::mpAttrs->readInt("loopsubdivs", mLoopSubdivs, "SimulationLbm","mLoopSubdivs", false );
1406         mSmoothSurface = D::mpAttrs->readFloat("smoothsurface", mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
1407         mSmoothNormals = D::mpAttrs->readFloat("smoothnormals", mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
1408
1409         mInitialCsmago = D::mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
1410         mCsmagoRefineMultiplier = D::mpAttrs->readFloat("csmago_refinemultiplier", mCsmagoRefineMultiplier, "SimulationLbm","mCsmagoRefineMultiplier", false );
1411
1412         // refinement
1413         mMaxRefine  = D::mpAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", true);
1414         mDisableStandingFluidInit = D::mpAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
1415         mForceTadapRefine = D::mpAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
1416
1417         // demo mode settings
1418         mDropMode = D::mpAttrs->readInt("dropmode", mDropMode, "SimulationLbm","mDropMode", false );
1419         mDropSize = D::mpAttrs->readFloat("dropsize", mDropSize, "SimulationLbm","mDropSize", false );
1420         mDropHeight = D::mpAttrs->readFloat("dropheight", mDropHeight, "SimulationLbm","mDropHeight", false );
1421         mDropSpeed = vec2L( D::mpAttrs->readVec3d("dropspeed", ntlVec3d(0.0), "SimulationLbm","mDropSpeed", false ) );
1422         if( (mDropMode>2) || (mDropMode<-1) ) mDropMode=1;
1423         mGfxGeoSetup = D::mpAttrs->readInt("gfxgeosetup", mGfxGeoSetup, "SimulationLbm","mGfxGeoSetup", false );
1424         mGfxEndTime = D::mpAttrs->readFloat("gfxendtime", mGfxEndTime, "SimulationLbm","mGfxEndTime", false );
1425         mFVHeight = D::mpAttrs->readFloat("fvolheight", mFVHeight, "SimulationLbm","mFVHeight", false );
1426         mFVArea   = D::mpAttrs->readFloat("fvolarea", mFVArea, "SimulationLbm","mFArea", false );
1427
1428 }
1429
1430
1431 /******************************************************************************
1432  * Initialize omegas and forces on all levels (for init/timestep change)
1433  *****************************************************************************/
1434 template<class D>
1435 void 
1436 LbmFsgrSolver<D>::initLevelOmegas()
1437 {
1438         // no explicit settings
1439         D::mOmega = D::mpParam->calculateOmega();
1440         D::mGravity = vec2L( D::mpParam->calculateGravity() );
1441         D::mSurfaceTension = D::mpParam->calculateSurfaceTension(); // unused
1442
1443         if(mInitialCsmago<=0.0) {
1444                 if(OPT3D==1) {
1445                         errMsg("LbmFsgrSolver::initLevelOmegas","Csmago-LES = 0 not supported for optimized 3D version..."); 
1446                         exit(1);
1447                 }
1448         }
1449
1450         // use Tau instead of Omega for calculations
1451         int i = mMaxRefine;
1452         mLevel[i].omega    = D::mOmega;
1453         mLevel[i].stepsize = D::mpParam->getStepTime();
1454         mLevel[i].lcsmago = mInitialCsmago; //CSMAGO_INITIAL;
1455         mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
1456         mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
1457
1458         // init all sub levels
1459         for(int i=mMaxRefine-1; i>=0; i--) {
1460                 //mLevel[i].omega = 2.0 * (mLevel[i+1].omega-0.5) + 0.5;
1461                 double nomega = 0.5 * (  (1.0/(double)mLevel[i+1].omega) -0.5) + 0.5;
1462                 nomega                = 1.0/nomega;
1463                 mLevel[i].omega       = (LbmFloat)nomega;
1464                 mLevel[i].stepsize    = 2.0 * mLevel[i+1].stepsize;
1465                 mLevel[i].lcsmago     = mLevel[i+1].lcsmago*mCsmagoRefineMultiplier;
1466                 mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
1467                 mLevel[i].lcnu        = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
1468         }
1469         
1470         // for lbgk
1471         mLevel[ mMaxRefine ].gravity = D::mGravity / mLevel[ mMaxRefine ].omega;
1472         for(int i=mMaxRefine-1; i>=0; i--) {
1473                 // should be the same on all levels...
1474                 // for lbgk
1475                 mLevel[i].gravity = (mLevel[i+1].gravity * mLevel[i+1].omega) * 2.0 / mLevel[i].omega;
1476         }
1477
1478         // debug? invalidate old values...
1479         D::mGravity = -100.0;
1480         D::mOmega = -100.0;
1481
1482         for(int i=0; i<=mMaxRefine; i++) {
1483                 if(!D::mSilent) {
1484                         errMsg("LbmFsgrSolver", "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" offs:"<<mLevel[i].lOffsx<<","<<mLevel[i].lOffsy<<","<<mLevel[i].lOffsz 
1485                                         <<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", "
1486                                         <<" cmsagp:"<<mLevel[i].lcsmago<<", "
1487                                         << " ss"<<mLevel[i].stepsize<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize );
1488                 } else {
1489                         if(!D::mInitDone) {
1490                                 debMsgStd("LbmFsgrSolver", DM_MSG, "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" "
1491                                                 <<"omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity , 5);
1492                         }
1493                 }
1494         }
1495         if(mMaxRefine>0) {
1496                 mDfScaleUp   = (mLevel[0  ].stepsize/mLevel[0+1].stepsize)* (1.0/mLevel[0  ].omega-1.0)/ (1.0/mLevel[0+1].omega-1.0); // yu
1497                 mDfScaleDown = (mLevel[0+1].stepsize/mLevel[0  ].stepsize)* (1.0/mLevel[0+1].omega-1.0)/ (1.0/mLevel[0  ].omega-1.0); // yu
1498         }
1499 }
1500
1501
1502 /******************************************************************************
1503  * Init Solver (values should be read from config file)
1504  *****************************************************************************/
1505 template<class D>
1506 bool 
1507 LbmFsgrSolver<D>::initialize( ntlTree* /*tree*/, vector<ntlGeometryObject*>* /*objects*/ )
1508 {
1509   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... (Layout:"<<ALSTRING<<") ",1);
1510
1511         // fix size inits to force cubic cells and mult4 level dimensions
1512         const int debugGridsizeInit = 1;
1513         mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)D::mSizex;
1514         int maxGridSize = D::mSizex; // get max size
1515         if(D::mSizey>maxGridSize) maxGridSize = D::mSizey;
1516         if(D::mSizez>maxGridSize) maxGridSize = D::mSizez;
1517         LbmFloat maxGeoSize = (D::mvGeoEnd[0]-D::mvGeoStart[0]); // get max size
1518         if((D::mvGeoEnd[1]-D::mvGeoStart[1])>maxGridSize) maxGeoSize = (D::mvGeoEnd[1]-D::mvGeoStart[1]);
1519         if((D::mvGeoEnd[2]-D::mvGeoStart[2])>maxGridSize) maxGeoSize = (D::mvGeoEnd[2]-D::mvGeoStart[2]);
1520         // FIXME better divide max geo size by corresponding resolution rather than max? no prob for rx==ry==rz though
1521         LbmFloat cellSize = (maxGeoSize / (LbmFloat)maxGridSize);
1522   if(debugGridsizeInit) debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Start:"<<D::mvGeoStart<<" End:"<<D::mvGeoEnd<<" maxS:"<<maxGeoSize<<" maxG:"<<maxGridSize<<" cs:"<<cellSize, 10);
1523         // force grid sizes according to geom. size, rounded
1524         D::mSizex = (int) ((D::mvGeoEnd[0]-D::mvGeoStart[0]) / cellSize +0.5);
1525         D::mSizey = (int) ((D::mvGeoEnd[1]-D::mvGeoStart[1]) / cellSize +0.5);
1526         D::mSizez = (int) ((D::mvGeoEnd[2]-D::mvGeoStart[2]) / cellSize +0.5);
1527         // match refinement sizes, round downwards to multiple of 4
1528         int sizeMask = 0;
1529         int maskBits = mMaxRefine;
1530         if(PARALLEL==1) maskBits+=2;
1531         for(int i=0; i<maskBits; i++) { sizeMask |= (1<<i); }
1532         sizeMask = ~sizeMask;
1533   if(debugGridsizeInit) debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Size X:"<<D::mSizex<<" Y:"<<D::mSizey<<" Z:"<<D::mSizez<<" m"<<convertFlags2String(sizeMask) ,10);
1534         D::mSizex &= sizeMask;
1535         D::mSizey &= sizeMask;
1536         D::mSizez &= sizeMask;
1537         if(PARALLEL==1) {
1538                 // make (size+2)%4 == 0 , assumes MAX_THREADS=4
1539                 D::mSizex += 2;
1540                 D::mSizey += 2;
1541                 D::mSizez += 2;
1542         }
1543         // force geom size to match rounded grid sizes
1544         D::mvGeoEnd[0] = D::mvGeoStart[0] + cellSize*(LbmFloat)D::mSizex;
1545         D::mvGeoEnd[1] = D::mvGeoStart[1] + cellSize*(LbmFloat)D::mSizey;
1546         D::mvGeoEnd[2] = D::mvGeoStart[2] + cellSize*(LbmFloat)D::mSizez;
1547
1548   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Final domain size X:"<<D::mSizex<<" Y:"<<D::mSizey<<" Z:"<<D::mSizez<<
1549                         ", Domain: "<<D::mvGeoStart<<":"<<D::mvGeoEnd<<", "<<(D::mvGeoEnd-D::mvGeoStart) ,2);
1550   //debMsgStd("LbmFsgrSolver::initialize",DM_MSG, ,2);
1551         D::mpParam->setSize(D::mSizex, D::mSizey, D::mSizez);
1552
1553   //debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Size X:"<<D::mSizex<<" Y:"<<D::mSizey<<" Z:"<<D::mSizez ,2);
1554
1555 #if ELBEEM_BLENDER!=1
1556   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: "
1557                 <<"LBM_EPSILON="<<LBM_EPSILON       <<" "
1558                 <<"FSGR_STRICT_DEBUG="<<FSGR_STRICT_DEBUG       <<" "
1559                 <<"INTORDER="<<INTORDER        <<" "
1560                 <<"TIMEINTORDER="<<TIMEINTORDER        <<" "
1561                 <<"REFINEMENTBORDER="<<REFINEMENTBORDER        <<" "
1562                 <<"INTCFCOARSETEST="<<INTCFCOARSETEST<<" "
1563                 <<"OPT3D="<<OPT3D        <<" "
1564                 <<"COMPRESSGRIDS="<<COMPRESSGRIDS<<" "
1565                 <<"LS_FLUIDTHRESHOLD="<<LS_FLUIDTHRESHOLD        <<" "
1566                 <<"MASS_INVALID="<<MASS_INVALID        <<" "
1567                 <<"FSGR_LISTTRICK="<<FSGR_LISTTRICK            <<" "
1568                 <<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY          <<" "
1569                 <<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL           <<" "
1570                 <<"FSGR_MAGICNR="<<FSGR_MAGICNR              <<" " 
1571                 <<"USE_LES="<<USE_LES              <<" " 
1572                 ,10);
1573 #endif // ELBEEM_BLENDER!=1
1574
1575         // perform 2D corrections...
1576         if(D::cDimension == 2) D::mSizez = 1;
1577
1578         D::mpParam->setSimulationMaxSpeed(0.0);
1579         if(mFVHeight>0.0) D::mpParam->setFluidVolumeHeight(mFVHeight);
1580         D::mpParam->setTadapLevels( mMaxRefine+1 );
1581
1582         if(mForceTadapRefine>mMaxRefine) {
1583                 D::mpParam->setTadapLevels( mForceTadapRefine+1 );
1584                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<<mForceTadapRefine, 6);
1585         }
1586
1587         if(!D::mpParam->calculateAllMissingValues()) {
1588                 errMsg("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...");
1589                 exit(1);
1590         }
1591
1592
1593
1594         // init vectors
1595         if(mMaxRefine >= MAX_LEV) {
1596                 errMsg("LbmFsgrSolver::initializeLbmGridref"," error: Too many levels!");
1597                 exit(1);
1598         }
1599         for(int i=0; i<=mMaxRefine; i++) {
1600                 mLevel[i].id = i;
1601                 mLevel[i].nodeSize = 0.0; 
1602                 mLevel[i].simCellSize = 0.0; 
1603                 mLevel[i].omega = 0.0; 
1604                 mLevel[i].time = 0.0; 
1605                 mLevel[i].stepsize = 1.0;
1606                 mLevel[i].gravity = LbmVec(0.0); 
1607                 mLevel[i].mprsCells[0] = NULL;
1608                 mLevel[i].mprsCells[1] = NULL;
1609                 mLevel[i].mprsFlags[0] = NULL;
1610                 mLevel[i].mprsFlags[1] = NULL;
1611
1612                 mLevel[i].avgOmega = 0.0; 
1613                 mLevel[i].avgOmegaCnt = 0.0;
1614         }
1615
1616         // init sizes
1617         mLevel[mMaxRefine].lSizex = D::mSizex;
1618         mLevel[mMaxRefine].lSizey = D::mSizey;
1619         mLevel[mMaxRefine].lSizez = D::mSizez;
1620         for(int i=mMaxRefine-1; i>=0; i--) {
1621                 mLevel[i].lSizex = mLevel[i+1].lSizex/2;
1622                 mLevel[i].lSizey = mLevel[i+1].lSizey/2;
1623                 mLevel[i].lSizez = mLevel[i+1].lSizez/2;
1624                 if( 
1625                   ((mLevel[i].lSizex % 4) != 0) ||
1626                   ((mLevel[i].lSizey % 4) != 0) ||
1627                   ((mLevel[i].lSizez % 4) != 0) ) {
1628                         errMsg("LbmFsgrSolver","Init: error invalid sizes on level "<<i<<" "<<PRINT_VEC(mLevel[i].lSizex,mLevel[i].lSizey,mLevel[i].lSizez) );
1629                         exit(1);
1630                 }
1631         }
1632
1633         // estimate memory usage
1634         {
1635                 unsigned long int memCnt = 0;
1636                 unsigned long int rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);
1637                 memCnt += sizeof(CellFlagType) * (rcellSize/dTotalNum +4) *2;
1638 #if COMPRESSGRIDS==0
1639                 memCnt += sizeof(LbmFloat) * (rcellSize +4) *2;
1640 #else // COMPRESSGRIDS==0
1641                 unsigned long int compressOffset = (mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum*2);
1642                 memCnt += sizeof(LbmFloat) * (rcellSize+compressOffset +4);
1643 #endif // COMPRESSGRIDS==0
1644                 for(int i=mMaxRefine-1; i>=0; i--) {
1645                         rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
1646                         memCnt += sizeof(CellFlagType) * (rcellSize/dTotalNum +4) *2;
1647                         memCnt += sizeof(LbmFloat) * (rcellSize +4) *2;
1648                 }
1649                 double memd = memCnt;
1650                 char *sizeStr = "";
1651                 const double sfac = 1000.0;
1652                 if(memd>sfac){ memd /= sfac; sizeStr="KB"; }
1653                 if(memd>sfac){ memd /= sfac; sizeStr="MB"; }
1654                 if(memd>sfac){ memd /= sfac; sizeStr="GB"; }
1655                 if(memd>sfac){ memd /= sfac; sizeStr="TB"; }
1656                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Required Grid memory: "<< memd <<" "<< sizeStr<<" ",4);
1657         }
1658
1659         mLevel[ mMaxRefine ].nodeSize = ((D::mvGeoEnd[0]-D::mvGeoStart[0]) / (LbmFloat)(D::mSizex));
1660         mLevel[ mMaxRefine ].simCellSize = D::mpParam->getCellSize();
1661         mLevel[ mMaxRefine ].lcellfactor = 1.0;
1662         unsigned long int rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);
1663         // +4 for safety ?
1664         mLevel[ mMaxRefine ].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
1665         mLevel[ mMaxRefine ].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
1666
1667 #if COMPRESSGRIDS==0
1668         mLevel[ mMaxRefine ].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
1669         mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
1670 #else // COMPRESSGRIDS==0
1671         unsigned long int compressOffset = (mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum*2);
1672         mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +compressOffset +4 ];
1673         mLevel[ mMaxRefine ].mprsCells[0] = mLevel[ mMaxRefine ].mprsCells[1]+compressOffset;
1674 #endif // COMPRESSGRIDS==0
1675
1676         for(int i=mMaxRefine-1; i>=0; i--) {
1677                 mLevel[i].nodeSize = 2.0 * mLevel[i+1].nodeSize;
1678                 mLevel[i].simCellSize = 2.0 * mLevel[i+1].simCellSize;
1679                 LbmFloat dimFac = 8.0;
1680                 if(D::cDimension==2) dimFac = 4.0;
1681                 mLevel[i].lcellfactor = mLevel[i+1].lcellfactor * dimFac;
1682
1683                 if(D::cDimension==2){ mLevel[i].lSizez = 1; } // 2D
1684                 rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
1685                 mLevel[i].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
1686                 mLevel[i].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
1687                 mLevel[i].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
1688                 mLevel[i].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
1689         }
1690
1691         // init sizes for _all_ levels
1692         for(int i=mMaxRefine; i>=0; i--) {
1693                 mLevel[i].lOffsx = mLevel[i].lSizex;
1694                 mLevel[i].lOffsy = mLevel[i].lOffsx*mLevel[i].lSizey;
1695                 mLevel[i].lOffsz = mLevel[i].lOffsy*mLevel[i].lSizez;
1696         mLevel[i].setCurr  = 0;
1697         mLevel[i].setOther = 1;
1698         mLevel[i].lsteps = 0;
1699         mLevel[i].lmass = 0.0;
1700         mLevel[i].lvolume = 0.0;
1701         }
1702
1703         // calc omega, force for all levels
1704         initLevelOmegas();
1705         mMinStepTime = D::mpParam->getStepTime();
1706         mMaxStepTime = D::mpParam->getStepTime();
1707
1708         // init isosurf
1709         D::mpIso->setIsolevel( D::mIsoValue );
1710         D::mpIso->setLoopSubdivs( mLoopSubdivs );
1711         // approximate feature size with mesh resolution
1712         float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
1713         D::mpIso->setSmoothSurface( mSmoothSurface * featureSize );
1714         D::mpIso->setSmoothNormals( mSmoothNormals * featureSize );
1715
1716         // init iso weight values mIsoWeightMethod
1717         int wcnt = 0;
1718         float totw = 0.0;
1719         for(int ak=-1;ak<=1;ak++) 
1720                 for(int aj=-1;aj<=1;aj++) 
1721                         for(int ai=-1;ai<=1;ai++)  {
1722                                 switch(mIsoWeightMethod) {
1723                                 case 1: // light smoothing
1724                                         mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
1725                                         break;
1726                                 case 2: // very light smoothing
1727                                         mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
1728                                         mIsoWeight[wcnt] *= mIsoWeight[wcnt];
1729                                         break;
1730                                 case 3: // no smoothing
1731                                         if(ai==0 && aj==0 && ak==0) mIsoWeight[wcnt] = 1.0;
1732                                         else mIsoWeight[wcnt] = 0.0;
1733                                         break;
1734                                 default: // strong smoothing (=0)
1735                                         mIsoWeight[wcnt] = 1.0;
1736                                         break;
1737                                 }
1738                                 totw += mIsoWeight[wcnt];
1739                                 wcnt++;
1740                         }
1741         for(int i=0; i<27; i++) mIsoWeight[i] /= totw;
1742
1743         LbmVec isostart = vec2L(D::mvGeoStart);
1744         LbmVec isoend   = vec2L(D::mvGeoEnd);
1745         int twodOff = 0; // 2d slices
1746         if(D::cDimension==2) {
1747                 LbmFloat sn,se;
1748                 sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(D::mSizex+1.0))*0.5;
1749                 se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(D::mSizex+1.0))*0.5;
1750                 isostart[2] = sn;
1751                 isoend[2] = se;
1752                 twodOff = 2;
1753         }
1754         //errMsg(" SETISO ", " "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(D::mSizex+1.0))*0.5)<<" "<<(LbmFloat)(D::mSizex+1.0)<<" " );
1755         D::mpIso->setStart( vec2G(isostart) );
1756         D::mpIso->setEnd(   vec2G(isoend) );
1757         LbmVec isodist = isoend-isostart;
1758         D::mpIso->initializeIsosurface( D::mSizex+2, D::mSizey+2, D::mSizez+2+twodOff, vec2G(isodist) );
1759         for(int ak=0;ak<D::mSizez+2+twodOff;ak++) 
1760                 for(int aj=0;aj<D::mSizey+2;aj++) 
1761                         for(int ai=0;ai<D::mSizex+2;ai++) { *D::mpIso->getData(ai,aj,ak) = 0.0; }
1762
1763   /* init array (set all invalid first) */
1764         for(int lev=0; lev<=mMaxRefine; lev++) {
1765                 FSGR_FORIJK_BOUNDS(lev) {
1766                         initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0); 
1767                 }
1768         }
1769
1770         // init defaults
1771         mAvgNumUsedCells = 0;
1772         D::mFixMass= 0.0;
1773
1774   /* init boundaries */
1775   debugOut("LbmFsgrSolver::initialize : Boundary init...",10);
1776
1777
1778         // use the density init?
1779         initGeometryFlags();
1780         D::initGenericTestCases();
1781         
1782         // new - init noslip 1 everywhere...
1783         // half fill boundary cells?
1784
1785   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1786     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1787                         initEmptyCell(mMaxRefine, i,0,k, CFBnd, 0.0, BND_FILL); 
1788                         initEmptyCell(mMaxRefine, i,mLevel[mMaxRefine].lSizey-1,k, CFBnd, 0.0, BND_FILL); 
1789     }
1790
1791         if(D::cDimension == 3) {
1792                 // only for 3D
1793                 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
1794                         for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1795                                 initEmptyCell(mMaxRefine, i,j,0, CFBnd, 0.0, BND_FILL); 
1796                                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, CFBnd, 0.0, BND_FILL); 
1797                         }
1798         }
1799
1800   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1801     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1802                         initEmptyCell(mMaxRefine, 0,j,k, CFBnd, 0.0, BND_FILL); 
1803                         // for quads! tag +x border... CF4_EMPTY, CF4_BND, CF4_UNUSED
1804                         // do this last to prevent any overwriting...
1805                         //initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, CFBnd|CFBndMARK, 0.0, BND_FILL); 
1806                         initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, CFBnd, 0.0, BND_FILL); 
1807     }
1808
1809         /*for(int ii=0; ii<(int)pow(2.0,mMaxRefine)-1; ii++) {
1810                 errMsg("BNDTESTSYMM","set "<<mLevel[mMaxRefine].lSizex-2-ii );
1811                 for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1812                         for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1813                                 initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2-ii,j,k, CFBnd, 0.0, BND_FILL);  // SYMM!? 2D?
1814                         }
1815                 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
1816                         for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1817                                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, CFBnd, 0.0, BND_FILL);   // SYMM!? 3D?
1818                         }
1819         }
1820         // Symmetry tests */
1821
1822         // prepare interface cells
1823         initFreeSurfaces();
1824         D::mInitialMass = 0.0;
1825         initStandingFluidGradient();
1826
1827         // perform first step to init initial mass
1828         mInitialMass = 0.0;
1829         int inmCellCnt = 0;
1830         FSGR_FORIJK1(mMaxRefine) {
1831                 if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr),        CFFluid) ) {
1832                         mInitialMass += 1.0;
1833                         inmCellCnt ++;
1834                 } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) {
1835                         mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass);
1836                         inmCellCnt ++;
1837                 }
1838         }
1839         mCurrentVolume = mCurrentMass = mInitialMass;
1840
1841         ParamVec cspv = D::mpParam->calculateCellSize();
1842         if(D::cDimension==2) cspv[2] = 1.0;
1843         inmCellCnt = 1;
1844         double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0;
1845         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Initial Mass:"<<mInitialMass<<" normalized:"<<nrmMass, 3);
1846
1847         //mStartSymm = false;
1848 #if ELBEEM_BLENDER!=1
1849         if((D::cDimension==2)&&(D::mSizex<200)) {
1850                 if(!checkSymmetry("init")) {
1851                         errMsg("LbmFsgrSolver::initialize","Unsymmetric init...");
1852                 } else {
1853                         errMsg("LbmFsgrSolver::initialize","Symmetric init!");
1854                 }
1855         }
1856 #endif // ELBEEM_BLENDER!=1
1857         
1858
1859         // ----------------------------------------------------------------------
1860         // coarsen region
1861         myTime_t fsgrtstart = getTime(); 
1862         for(int lev=mMaxRefine-1; lev>=0; lev--) {
1863                 while(performCoarsening(lev)){
1864                         coarseRestrictFromFine(lev);
1865                         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8);
1866                 }
1867         }
1868         D::markedClearList();
1869         myTime_t fsgrtend = getTime(); 
1870         if(!D::mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< ((fsgrtend-fsgrtstart)/(double)1000.0)<<"s) " , 10 ); }
1871
1872         for(int l=0; l<D::cDirNum; l++) { 
1873                 LbmFloat area = 0.5 * 0.5 *0.5;
1874                 if(D::cDimension==2) area = 0.5 * 0.5;
1875
1876                 if(D::dfVecX[l]!=0) area *= 0.5;
1877                 if(D::dfVecY[l]!=0) area *= 0.5;
1878                 if(D::dfVecZ[l]!=0) area *= 0.5;
1879                 mFsgrCellArea[l] = area;
1880         } // l
1881         /*for(int lev=0; lev<mMaxRefine; lev++) {
1882         FSGR_FORIJK_BOUNDS(lev) {
1883                 if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) {
1884                         if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) {
1885                                 LbmFloat totArea = mFsgrCellArea[0]; // for l=0
1886                                 for(int l=1; l<D::cDirNum; l++) { 
1887                                         int ni=(2*i)+D::dfVecX[l], nj=(2*j)+D::dfVecY[l], nk=(2*k)+D::dfVecZ[l];
1888                                         if(RFLAG(lev+1, ni,nj,nk, mLevel[lev+1].setCurr)&
1889                                                         (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
1890                                                         //(CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
1891                                                         ) { 
1892                                                 //LbmFloat area = 0.25; if(D::dfVecX[l]!=0) area *= 0.5; if(D::dfVecY[l]!=0) area *= 0.5; if(D::dfVecZ[l]!=0) area *= 0.5;
1893                                                 totArea += mFsgrCellArea[l];
1894                                         }
1895                                 } // l
1896                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = totArea;
1897                         } else if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFEmpty) {
1898                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 1.0;
1899                         } else {
1900                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 0.0;
1901                         }
1902                         errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) );
1903                 }
1904         } } // */
1905
1906         // now really done...
1907   debugOut("LbmFsgrSolver::initialize : Init done ...",10);
1908         D::mInitDone = 1;
1909
1910         // make sure both sets are ok
1911         // copy from other to curr
1912         for(int lev=0; lev<=mMaxRefine; lev++) {
1913         FSGR_FORIJK_BOUNDS(lev) {
1914                 RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1915         } } // first copy flags */
1916 #if COMPRESSGRIDS==0
1917         /*for(int lev=0; lev<=mMaxRefine; lev++) {
1918         FSGR_FORIJK_BOUNDS(lev) {
1919                 // copy from other to curr
1920                 for(int l=0; l<D::cDfNum; l++) {
1921                         QCELL(lev, i,j,k,mLevel[lev].setOther, l) = QCELL(lev, i,j,k,mLevel[lev].setCurr, l);
1922                 }
1923                 QCELL(lev, i,j,k,mLevel[lev].setOther, dMass) = QCELL(lev, i,j,k,mLevel[lev].setCurr, dMass);
1924                 QCELL(lev, i,j,k,mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k,mLevel[lev].setCurr, dFfrac);
1925                 QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux);
1926         } } // COMPRT OFF */
1927 #endif // COMPRESSGRIDS==0
1928
1929
1930
1931         
1932         if(mOutputSurfacePreview) {
1933                 if(D::cDimension==2) {
1934                         errMsg("LbmFsgrSolver::init","No preview in 2D allowed!"); exit (1); 
1935                 }
1936
1937                 //int previewSize = mOutputSurfacePreview;
1938                 // same as normal one, but use reduced size
1939                 mpPreviewSurface = new IsoSurface( D::mIsoValue, false );
1940                 mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
1941                 mpPreviewSurface->setIsolevel( D::mIsoValue );
1942                 // usually dont display for rendering
1943                 mpPreviewSurface->setVisible( false );
1944
1945                 mpPreviewSurface->setStart( vec2G(isostart) );
1946                 mpPreviewSurface->setEnd(   vec2G(isoend) );
1947                 LbmVec isodist = isoend-isostart;
1948                 mpPreviewSurface->initializeIsosurface( (int)(mPreviewFactor*D::mSizex)+2, (int)(mPreviewFactor*D::mSizey)+2, (int)(mPreviewFactor*D::mSizez)+2, vec2G(isodist) );
1949                 //mpPreviewSurface->setName( D::getName() + "preview" );
1950                 mpPreviewSurface->setName( "preview" );
1951         
1952                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(mPreviewFactor*D::mSizex)<<","<<(mPreviewFactor*D::mSizey)<<","<<(mPreviewFactor*D::mSizez)<<" enabled",10);
1953         }
1954
1955 #if ELBEEM_BLENDER==1
1956         // make sure fill fracs are right for first surface generation
1957         stepMain();
1958 #endif // ELBEEM_BLENDER==1
1959
1960         // prepare once...
1961         prepareVisualization();
1962         // copy again for stats counting
1963         for(int lev=0; lev<=mMaxRefine; lev++) {
1964         FSGR_FORIJK_BOUNDS(lev) {
1965                 RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1966         } } // first copy flags */
1967
1968         /*{ int lev=mMaxRefine;
1969                 FSGR_FORIJK_BOUNDS(lev) { // COMPRT deb out
1970                 debMsgDirect("\n x="<<PRINT_IJK);
1971                 for(int l=0; l<D::cDfNum; l++) {
1972                         debMsgDirect(" df="<< QCELL(lev, i,j,k,mLevel[lev].setCurr, l) );
1973                 }
1974                 debMsgDirect(" m="<< QCELL(lev, i,j,k,mLevel[lev].setCurr, dMass) );
1975                 debMsgDirect(" f="<< QCELL(lev, i,j,k,mLevel[lev].setCurr, dFfrac) );
1976                 debMsgDirect(" x="<< QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) );
1977         } } // COMPRT ON */
1978         return true;
1979 }
1980
1981
1982
1983 /*****************************************************************************/
1984 /*! perform geometry init (if switched on) */
1985 /*****************************************************************************/
1986 template<class D>
1987 bool 
1988 LbmFsgrSolver<D>::initGeometryFlags() {
1989         if(!D::mPerformGeoInit) return false;
1990         int level = mMaxRefine;
1991         myTime_t geotimestart = getTime(); 
1992         ntlGeometryObject *pObj;
1993         // getCellSize (due to forced cubes, use x values)
1994         ntlVec3Gfx dvec( (D::mvGeoEnd[0]-D::mvGeoStart[0])/ ((LbmFloat)D::mSizex*2.0));
1995         // real cell size from now on...
1996         dvec *= 2.0; 
1997         ntlVec3Gfx nodesize = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
1998         dvec = nodesize;
1999         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< D::mGeoInitId <<") v"<<dvec,3);
2000
2001         /* set interface cells */
2002         D::initGeoTree(D::mGeoInitId);
2003         ntlVec3Gfx maxIniVel = vec2G( D::mpParam->calculateLattVelocityFromRw( vec2P(D::getGeoMaxInitialVelocity()) ));
2004         D::mpParam->setSimulationMaxSpeed( norm(maxIniVel) + norm(mLevel[level].gravity) );
2005         LbmFloat allowMax = D::mpParam->getTadapMaxSpeed();  // maximum allowed velocity
2006         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxIniVel <<", allowed Max="<<allowMax ,5);
2007         if(D::mpParam->getSimulationMaxSpeed() > allowMax) {
2008                 // similar to adaptTimestep();
2009                 LbmFloat nextmax = D::mpParam->getSimulationMaxSpeed();
2010                 LbmFloat newdt = D::mpParam->getStepTime() * (allowMax / nextmax); // newtr
2011                 debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< D::mpParam->getStepTime() <<" ",5);
2012                 D::mpParam->setDesiredStepTime( newdt );
2013                 D::mpParam->calculateAllMissingValues( D::mSilent );
2014                 maxIniVel = vec2G( D::mpParam->calculateLattVelocityFromRw( vec2P(D::getGeoMaxInitialVelocity()) ));
2015                 debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxIniVel,5);
2016         }
2017
2018
2019         ntlVec3Gfx pos,iniPos; // position of current cell
2020         LbmFloat rhomass = 0.0;
2021         int savedNodes = 0;
2022         int OId = -1;
2023         gfxReal distance;
2024         LbmVec nodevel(0.0);
2025
2026         // 2d display as rectangles
2027         if(D::cDimension==2) {
2028                 dvec[2] = 0.0; 
2029                 iniPos =(D::mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (D::mvGeoEnd[2]-D::mvGeoStart[2])*0.5 ))-(dvec*0.0);
2030                 //iniPos =(D::mvGeoStart + ntlVec3Gfx( 0.0 ))+dvec;
2031         } else {
2032                 iniPos =(D::mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
2033                 iniPos[2] = D::mvGeoStart[2] + dvec[2]*getForZMin1(level);
2034         }
2035
2036
2037         // first init boundary conditions
2038 #define GETPOS(i,j,k) \
2039                                                 ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
2040                                                 iniPos[1]+ dvec[1]*(gfxReal)(j), \
2041                                                 iniPos[2]+ dvec[2]*(gfxReal)(k) )
2042         //pos = iniPos;
2043         for(int k= getForZMin1(level); k< getForZMax1(level); ++k) {
2044                 //pos[1] = D::mvGeoStart[1] + dvec[1]*1.0;
2045                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
2046                         //pos[0] = D::mvGeoStart[0] + dvec[0]*1.0;
2047                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
2048
2049                                 //CellFlagType ntype = D::geoInitGetPointType( pos, nodesize, &pObj, distance );
2050                                 CellFlagType ntype = CFInvalid;
2051                                 if(D::geoInitCheckPointInside( GETPOS(i,j,k) , FGI_ALLBOUNDS, OId, distance)) {
2052                                         ntype = CFBnd;
2053                                         pObj = (*D::mpGiObjects)[OId];
2054                                 }
2055                                 //CellFlagType  convntype = ntype;
2056                                 if(ntype != CFInvalid) {
2057                                         // initDefaultCell
2058                                         nodevel = LbmVec(0.0);
2059                                         rhomass = BND_FILL;
2060                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, nodevel );
2061                                 }
2062
2063                                 // walk along x until hit for following inits
2064                                 if(distance<=-1.0) { distance = 100.0; }
2065                                 if(distance>0.0) {
2066                                         //distance += pos[0];
2067                                         //while(((pos[0]+dvec[0])<distance)&&(i+1<mLevel[level].lSizex-1)) {
2068                                         gfxReal dcnt=dvec[0];
2069                                         while((dcnt<distance)&&(i+1<mLevel[level].lSizex-1)) {
2070                                                 dcnt += dvec[0]; i++;
2071                                                 savedNodes++;
2072                                                 if(ntype != CFInvalid) {
2073                                                         // nodevel&rhomass are still inited from above
2074                                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, nodevel );
2075                                                 }
2076                                         }
2077                                 } //pos[0] += dvec[0];
2078                         } //pos[1] += dvec[1];
2079                 } //pos[2] += dvec[2];
2080         } // zmax
2081
2082
2083         // now init fluid layer
2084         //pos = iniPos;
2085         for(int k= getForZMin1(level); k< getForZMax1(level); ++k) {
2086                 //pos[1] = D::mvGeoStart[1] + dvec[1]*1.0;
2087                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
2088                         //pos[0] = D::mvGeoStart[0] + dvec[0]*1.0;
2089                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
2090                                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)&CFEmpty)) continue;
2091
2092                                 //CellFlagType ntype = D::geoInitGetPointType( pos, nodesize, &pObj, distance );
2093                                 CellFlagType ntype = CFInvalid;
2094                                 if(D::geoInitCheckPointInside( GETPOS(i,j,k) , FGI_FLUID, OId, distance)) {
2095                                         ntype = CFFluid;
2096                                         pObj = (*D::mpGiObjects)[OId];
2097                                 }
2098                                 //CellFlagType  convntype = ntype;
2099                                 if(ntype != CFInvalid) {
2100                                         // initDefaultCell
2101                                         nodevel = vec2L(D::mpParam->calculateLattVelocityFromRw( vec2P(pObj->getInitialVelocity() )));
2102                                         rhomass = 1.0;
2103                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, nodevel );
2104                                 }
2105
2106                                 // walk along x until hit for following inits
2107                                 if(distance<=-1.0) { distance = 100.0; }
2108                                 if(distance>0.0) {
2109                                         //distance += pos[0];
2110                                         //while(((pos[0]+dvec[0])<distance)&&(i+1<mLevel[level].lSizex-1)) {
2111                                                 //pos[0] += dvec[0]; i++;
2112                                         gfxReal dcnt=dvec[0];
2113                                         while((dcnt<distance)&&(i+1<mLevel[level].lSizex-1)) {
2114                                                 dcnt += dvec[0]; i++;
2115                                                 savedNodes++;
2116                                                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)&CFEmpty)) continue;
2117                                                 if(ntype != CFInvalid) {
2118                                                         // nodevel&rhomass are still inited from above
2119                                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, nodevel );
2120                                                 }
2121                                         }
2122                                 } //pos[0] += dvec[0];
2123                         } //pos[1] += dvec[1];
2124                 } //pos[2] += dvec[2];
2125         } // zmax
2126
2127         /*
2128         // now init fluid layer
2129         if(D::cDimension==2) {
2130                 dvec[2] = 0.0; 
2131                 pos =(D::mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (D::mvGeoEnd[2]-D::mvGeoStart[2])*0.5 ))+dvec;
2132                 pos[2] = D::mvGeoStart[2] + dvec[2]*getForZMin1(mMaxRefine);
2133         } else {
2134                 pos =(D::mvGeoStart + ntlVec3Gfx( 0.0, 0.0, 0.0 ))+dvec;
2135         }
2136         for(int k= getForZMin1(mMaxRefine); k< getForZMax1(mMaxRefine); ++k) {
2137                 pos[1] = D::mvGeoStart[1] + dvec[1]*1.0;
2138                 for(int j=1;j<mLevel[mMaxRefine].lSizey-1;j++) {
2139                         pos[0] = D::mvGeoStart[0] + dvec[0]*1.0;
2140                         for(int i=1;i<mLevel[mMaxRefine].lSizex-1;i++) {
2141                                 //debMsgInter("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Checking"<<PRINT_IJK,2,2000 );
2142                                 if(RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)&CFEmpty) {
2143                                         // overwrite these...
2144                                 } else {
2145                                         continue;
2146                                 }
2147
2148                                 CellFlagType ntype = D::geoInitGetPointType( pos, nodesize, &pObj, distance );
2149                                 CellFlagType  convntype = ntype;
2150                                 if(ntype != CFInvalid) {
2151                                         // initDefaultCell
2152                                         if(convntype & CFFluid) { nodevel = vec2L(D::mpParam->calculateLattVelocityFromRw( vec2P(pObj->getInitialVelocity() ))); }
2153                                         else { nodevel = LbmVec(0.0); }
2154
2155                                         if(convntype & CFFluid) { rhomass = 1.0; }
2156                                         else if(convntype & CFBnd) { rhomass = BND_FILL; }
2157                                         else { rhomass = 0.0; } // empty, inter
2158                                         //initEmptyCell(mMaxRefine, i,j,k, convntype, rhomass, rhomass );
2159                                         initVelocityCell(mMaxRefine, i,j,k, convntype, rhomass, rhomass, nodevel );
2160                                         //errMsg(" DT ",PRINT_IJK<<" "<<distance<<" "<<pos<<" i"<<i<<"  "<<ntype );
2161                                 }
2162
2163                                 // walk along x until hit for following inits
2164                                 //if(distance==-1.0) { distance = 100.0; }
2165                                 if(distance<=-1.0) { distance = 100.0; }
2166                                 if(distance>0.0) {
2167                                         distance += pos[0];
2168                                         //errMsg(" DS "," "<<distance<<" "<<pos<<" i"<<i<<"  type:"<<ntype<<" "<<convertFlags2String(ntype) );
2169                                         while(((pos[0]+dvec[0])<distance)&&(i+1<mLevel[mMaxRefine].lSizex-1)) {
2170                                                 if(RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)&CFEmpty) {
2171                                                         // overwrite these...
2172                                                 } else {
2173                                                         continue;
2174                                                 }
2175                                                 pos[0] += dvec[0]; i++;
2176                                                 savedNodes++;
2177                                                 //errMsg(" DT "," "<<distance<<" "<<pos<<" i"<<i<<"  "<<ntype );
2178                                                 if(ntype != CFInvalid) {
2179                                                         // initDefaultCell
2180                                                         // nodevel is still inited from above
2181                                                         if(convntype & CFFluid) { rhomass = 1.0; }
2182                                                         else if(convntype & CFBnd) { rhomass = BND_FILL; }
2183                                                         else { rhomass = 0.0; } // empty, inter
2184                                                         //initEmptyCell(mMaxRefine, i,j,k, convntype, rhomass, rhomass );
2185                                                         initVelocityCell(mMaxRefine, i,j,k, convntype, rhomass, rhomass, nodevel );
2186                                                 }
2187                                         }
2188                                 }
2189                                 pos[0] += dvec[0];
2190                         }
2191                         pos[1] += dvec[1];
2192                 }
2193                 pos[2] += dvec[2];
2194         } // */
2195
2196         D::freeGeoTree();
2197         myTime_t geotimeend = getTime(); 
2198         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< ((geotimeend-geotimestart)/(double)1000.0)<<"s) " , 10 ); 
2199         //errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez));
2200         return true;
2201 }
2202
2203 /*****************************************************************************/
2204 /* init part for all freesurface testcases */
2205 template<class D>
2206 void 
2207 LbmFsgrSolver<D>::initFreeSurfaces() {
2208         double interfaceFill = 0.45;   // filling level of interface cells
2209
2210         // set interface cells 
2211         FSGR_FORIJK1(mMaxRefine) {
2212
2213                 /*if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFEmpty )) {
2214                         int initInter = 0;
2215                         // check for neighboring fluid cells 
2216                         FORDF1 {
2217                                 if( TESTFLAG( RFLAG_NBINV(mMaxRefine, i, j, k,  mLevel[mMaxRefine].setCurr,l), CFFluid ) ) {
2218                                         initInter = 1;
2219                                 }
2220                         }
2221
2222                         if(initInter) {
2223                                 initEmptyCell(mMaxRefine, i,j,k, CFInter, 1.0, interfaceFill);
2224                         }
2225
2226                 } // empty cells  OLD */
2227
2228                 if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFFluid )) {
2229                         int initInter = 0; // check for neighboring empty cells 
2230                         FORDF1 {
2231                                 if( TESTFLAG( RFLAG_NBINV(mMaxRefine, i, j, k,  mLevel[mMaxRefine].setCurr,l), CFEmpty ) ) {
2232                                         initInter = 1;
2233                                 }
2234                         }
2235                         if(initInter) {
2236                                 QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr, dMass) = 
2237                                         //QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setOther, dMass) =  // COMPRT OFF
2238                                         interfaceFill;
2239                                 RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) = RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setOther) = CFInter;
2240                         }
2241                 }
2242         }
2243
2244         // remove invalid interface cells 
2245         FSGR_FORIJK1(mMaxRefine) {
2246                 // remove invalid interface cells 
2247                 if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) {
2248                         int delit = 0;
2249                         int NBs = 0; // neighbor flags or'ed 
2250                         int noEmptyNB = 1;
2251
2252                         FORDF1 {
2253                                 if( TESTFLAG( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l ), CFEmpty ) ) {
2254                                         noEmptyNB = 0;
2255                                 }
2256                                 NBs |= RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr, l);
2257                         }
2258
2259                         // remove cells with no fluid or interface neighbors
2260                         if((NBs & CFFluid)==0) { delit = 1; }
2261                         if((NBs & CFInter)==0) { delit = 1; }
2262
2263                         // remove cells with no empty neighbors
2264                         if(noEmptyNB) { delit = 2; }
2265
2266                         // now we can remove the cell 
2267                         if(delit==1) {
2268                                 initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0);
2269                         }
2270                         if(delit==2) {
2271                                 initEmptyCell(mMaxRefine, i,j,k, CFFluid, 1.0, 1.0);
2272                         }
2273                 } // interface 
2274         }
2275
2276         // another brute force init, make sure the fill values are right...
2277         // and make sure both sets are equal
2278         for(int lev=0; lev<=mMaxRefine; lev++) {
2279         FSGR_FORIJK_BOUNDS(lev) {
2280                 if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) { 
2281                         QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 
2282                                 //QCELL(lev, i,j,k,mLevel[mMaxRefine].setOther, dFfrac) =  // COMPRT OFF
2283                                 BND_FILL;
2284                         continue;
2285                 }
2286                 if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) { 
2287                         QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 
2288                                 //QCELL(lev, i,j,k,mLevel[mMaxRefine].setOther, dFfrac) =  // COMPRT OFF
2289                                 0.0;
2290                         continue;
2291                 }
2292
2293                 // copy from other to curr
2294                 //?? RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
2295                 //?? for(int l=0; l<D::cDfNum; l++) {
2296                         //?? QCELL(lev, i,j,k,mLevel[lev].setOther, l) = QCELL(lev, i,j,k,mLevel[lev].setCurr, l);
2297                 //?? }
2298                 //?? QCELL(lev, i,j,k,mLevel[lev].setOther, dMass) = QCELL(lev, i,j,k,mLevel[lev].setCurr, dMass);
2299                 //?? QCELL(lev, i,j,k,mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k,mLevel[lev].setCurr, dFfrac);
2300         } }
2301
2302         // ----------------------------------------------------------------------
2303         // smoother surface...
2304         if(mInitSurfaceSmoothing>0) {
2305                 debMsgStd("Surface Smoothing init", DM_MSG, "Performing "<<(mInitSurfaceSmoothing)<<" smoothing steps ",10);
2306 #if COMPRESSGRIDS==1
2307                 errMsg("NYI","COMPRESSGRIDS mInitSurfaceSmoothing"); exit(1);
2308 #endif // COMPRESSGRIDS==0
2309         }
2310         for(int s=0; s<mInitSurfaceSmoothing; s++) {
2311                 FSGR_FORIJK1(mMaxRefine) {
2312                         if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) {
2313                                 LbmFloat mass = 0.0;
2314                                 //LbmFloat nbdiv;
2315                                 FORDF0 {
2316                                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2317                                         if( RFLAG(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr) & CFFluid ){
2318                                                 mass += 1.0;
2319                                         }
2320                                         if( RFLAG(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr) & CFInter ){
2321                                                 mass += QCELL(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, dMass);
2322                                         }
2323                                         //nbdiv+=1.0;
2324                                 }
2325
2326                                 //errMsg(" I ", PRINT_IJK<<" m"<<mass );
2327                                 QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setOther, dMass) = (mass/19.0);
2328                                 QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setOther, dFfrac) = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setOther, dMass);
2329                         }
2330                 }
2331
2332                 mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr;
2333                 mLevel[mMaxRefine].setCurr ^= 1;
2334         }
2335
2336         // copy back...
2337         /*for(int lev=0; lev<=mMaxRefine; lev++) {
2338                 FSGR_FORIJK1(lev) {
2339                         if( TESTFLAG( RFLAG(lev, i,j,k, mLevel[lev].setCurr), CFInter) ) {
2340                                 QCELL(lev, i,j,k, mLevel[lev].setOther, dMass ) = QCELL(lev, i,j,k, mLevel[lev].setCurr, dMass);
2341                                 QCELL(lev, i,j,k, mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k, mLevel[lev].setCurr, dFfrac);
2342                                 QCELL(lev, i,j,k, mLevel[lev].setOther, dFlux) = QCELL(lev, i,j,k, mLevel[lev].setCurr, dFlux);
2343                         }
2344                 }
2345         } // COMPRT OFF */
2346
2347 }
2348
2349 /*****************************************************************************/
2350 /* init part for all freesurface testcases */
2351 template<class D>
2352 void 
2353 LbmFsgrSolver<D>::initStandingFluidGradient() {
2354         // ----------------------------------------------------------------------
2355         // standing fluid preinit
2356         const int debugStandingPreinit = 0;
2357         int haveStandingFluid = 0;
2358
2359 #define STANDFLAGCHECK(iindex) \
2360                                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || \
2361                                                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){  \
2362                                         if((iindex)>1) { haveStandingFluid=(iindex); } \
2363                                         j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=D::getForZMaxBnd(); \
2364                                         continue; \
2365                                 } 
2366         int gravIndex[3] = {0,0,0};
2367         int gravDir[3] = {1,1,1};
2368         int maxGravComp = 1; // by default y
2369         int gravComp1 = 0; // by default y
2370         int gravComp2 = 2; // by default y
2371         if( ABS(mLevel[mMaxRefine].gravity[0]) > ABS(mLevel[mMaxRefine].gravity[1]) ){ maxGravComp = 0; gravComp1=1; gravComp2=2; }
2372         if( ABS(mLevel[mMaxRefine].gravity[2]) > ABS(mLevel[mMaxRefine].gravity[0]) ){ maxGravComp = 2; gravComp1=0; gravComp2=1; }
2373
2374         int gravIMin[3] = { 0 , 0 , 0 };
2375         int gravIMax[3] = {
2376                 mLevel[mMaxRefine].lSizex + 0,
2377                 mLevel[mMaxRefine].lSizey + 0,
2378                 mLevel[mMaxRefine].lSizez + 0 };
2379         if(LBMDIM==2) gravIMax[2] = 1;
2380
2381         //int gravDir = 1;
2382         if( mLevel[mMaxRefine].gravity[maxGravComp] > 0.0 ) {
2383                 // swap directions
2384                 int i=maxGravComp;
2385                 int tmp = gravIMin[i];
2386                 gravIMin[i] = gravIMax[i] - 1;
2387                 gravIMax[i] = tmp - 1;
2388                 gravDir[i] = -1;
2389         }
2390 #define PRINTGDIRS \
2391         errMsg("Standing fp","X start="<<gravIMin[0]<<" end="<<gravIMax[0]<<" dir="<<gravDir[0] ); \
2392         errMsg("Standing fp","Y start="<<gravIMin[1]<<" end="<<gravIMax[1]<<" dir="<<gravDir[1] ); \
2393         errMsg("Standing fp","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] ); 
2394         // _PRINTGDIRS;
2395
2396         bool gravAbort = false;
2397 #define GRAVLOOP \
2398         gravAbort=false; \
2399         for(gravIndex[2]= gravIMin[2];     (gravIndex[2]!=gravIMax[2])&&(!gravAbort);  gravIndex[2] += gravDir[2]) \
2400                 for(gravIndex[1]= gravIMin[1];   (gravIndex[1]!=gravIMax[1])&&(!gravAbort);  gravIndex[1] += gravDir[1]) \
2401                         for(gravIndex[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort);  gravIndex[0] += gravDir[0]) 
2402
2403         GRAVLOOP {
2404                 int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
2405                 //if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) {debMsgStd("Standing fluid preinit", DM_MSG, "fluidheightinit check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 ); }
2406                 //STANDFLAGCHECK(gravIndex[maxGravComp]);
2407                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || 
2408                                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){  
2409                         int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp]));
2410                         if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] );
2411                         //if(gravIndex[maxGravComp]>1)  
2412                         if(fluidHeight>1) 
2413                         {
2414                                 haveStandingFluid = fluidHeight; //gravIndex[maxGravComp]; 
2415                                 gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp];
2416                         }
2417                         gravAbort = true; continue; 
2418                 } 
2419         } // GRAVLOOP
2420         // _PRINTGDIRS;
2421
2422         LbmFloat fluidHeight;
2423         //if(gravDir>0) { fluidHeight = (LbmFloat)haveStandingFluid;
2424         //} else { fluidHeight = (LbmFloat)haveStandingFluid; }
2425         fluidHeight = (LbmFloat)(ABS(gravIMax[maxGravComp]-gravIMin[maxGravComp]));
2426         if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "fheight="<<fluidHeight<<" min="<<PRINT_VEC(gravIMin[0],gravIMin[1],       gravIMin[2])<<" max="<<PRINT_VEC(gravIMax[0], gravIMax[1],gravIMax[2])<<
2427                         " mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10);
2428                                 
2429         if(mDisableStandingFluidInit) {
2430                 debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2);
2431                 haveStandingFluid=0;
2432         }
2433
2434         if(haveStandingFluid) {
2435                 int lev = mMaxRefine;
2436                 int rhoworkSet = mLevel[lev].setCurr;
2437                 myTime_t timestart = getTime(); // FIXME use user time here?
2438 #if OPT3D==true 
2439                 LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM], lcsmomega;
2440 #endif // OPT3D==true 
2441
2442                 GRAVLOOP {
2443                         int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
2444                         //debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 );
2445                         if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) ||
2446                                         ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){ 
2447                                 //gravAbort = true; 
2448                                 continue;
2449                         }
2450
2451                         LbmFloat rho = 1.0;
2452                         // 1/6 velocity from denisty gradient, 1/2 for delta of two cells
2453                         rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) * 
2454                                 (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega); 
2455                         if(debugStandingPreinit) 
2456                                 if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) { 
2457                                         errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK); 
2458                                 }
2459
2460                         if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) ||
2461                                         (RFLAG(lev, i,j,k, rhoworkSet) & CFInter) ) {
2462                                 FORDF0 { QCELL(lev, i,j,k, rhoworkSet, l) *= rho; }
2463                                 QCELL(lev, i,j,k, rhoworkSet, dMass) *= rho;
2464                         }
2465
2466                 } // GRAVLOOP
2467                 debMsgStd("Standing fluid preinit", DM_MSG, "Density gradient inited", 8);
2468                 
2469                 // copy flags and init , as no flags will be changed during grav init
2470                 CellFlagType nbflag[LBM_DFNUM], nbored; 
2471                 for(int k=D::getForZMinBnd();k<D::getForZMaxBnd();++k) {
2472                 for(int j=0;j<mLevel[lev].lSizey-0;++j) {
2473                 for(int i=0;i<mLevel[lev].lSizex-0;++i) {
2474                                 if( (RFLAG(lev, i,j,k,SRCS(lev)) & (CFFluid)) ) {
2475                                         nbored = 0;
2476                                         FORDF1 {
2477                                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
2478                                                 nbored |= nbflag[l];
2479                                         } 
2480                                         if(nbored&CFBnd) {
2481                                                 RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid);
2482                                         } else {
2483                                                 RFLAG(lev, i,j,k,SRCS(lev)) |= CFNoBndFluid;
2484                                         }
2485                                 }
2486                                 RFLAG(lev, i,j,k,TSET(lev)) = RFLAG(lev, i,j,k,SRCS(lev));
2487                 } } }
2488
2489                 int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
2490                 preinitSteps = (haveStandingFluid>>2); // not much use...?
2491                 //preinitSteps = 4; // DEBUG!!!!
2492                 //D::mInitDone = 1; // GRAVTEST
2493                 //preinitSteps = 0;
2494                 debMsgNnl("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
2495                 for(int s=0; s<preinitSteps; s++) {
2496                         int workSet = SRCS(lev); //mLevel[lev].setCurr;
2497                         int otherSet = TSET(lev); //mLevel[lev].setOther;
2498                         //debMsgDirect(".");
2499                         if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "s="<<s<<" curset="<<workSet<<" srcs"<<SRCS(lev), 10);
2500                         LbmFloat *ccel;
2501                         LbmFloat *tcel;
2502                         LbmFloat m[LBM_DFNUM];
2503
2504                 // grav loop not necessary here
2505 #define NBFLAG(l) (nbflag[(l)])
2506                 LbmFloat rho, ux,uy,uz, usqr; 
2507                 int kstart=D::getForZMinBnd(), kend=D::getForZMaxBnd();
2508 #if COMPRESSGRIDS==0
2509                 for(int k=kstart;k<kend;++k) {
2510 #else // COMPRESSGRIDS==0
2511                 int kdir = 1; // COMPRT ON
2512                 if(mLevel[lev].setCurr==1) {
2513                         kdir = -1;
2514                         int temp = kend;
2515                         kend = kstart-1;
2516                         kstart = temp-1;
2517                 } // COMPRT
2518                 for(int k=kstart;k!=kend;k+=kdir) {
2519
2520                 //errMsg("LbmFsgrSolver::mainLoop","k="<<k<<" ks="<<kstart<<" ke="<<kend<<" kdir="<<kdir ); // debug
2521 #endif // COMPRESSGRIDS==0
2522
2523                 for(int j=0;j<mLevel[lev].lSizey-0;++j) {
2524                 for(int i=0;i<mLevel[lev].lSizex-0;++i) {
2525                                 const CellFlagType currFlag = RFLAG(lev, i,j,k,workSet);
2526                                 if( (currFlag & (CFEmpty|CFBnd)) ) continue;
2527                                 ccel = RACPNT(lev, i,j,k,workSet); 
2528                                 tcel = RACPNT(lev, i,j,k,otherSet);
2529
2530                                 if( (currFlag & (CFInter)) ) {
2531                                         // copy all values
2532                                         for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
2533                                         continue;
2534                                 }
2535
2536                                 if( (currFlag & CFNoBndFluid)) {
2537                                         OPTIMIZED_STREAMCOLLIDE;
2538                                 } else {
2539                                         FORDF1 {
2540                                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
2541                                         } 
2542                                         DEFAULT_STREAM;
2543                                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
2544                                         DEFAULT_COLLIDE;
2545                                 }
2546                                 for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
2547                         } } } // GRAVLOOP
2548
2549                         mLevel[lev].setOther = mLevel[lev].setCurr;
2550                         mLevel[lev].setCurr ^= 1;
2551                 }
2552                 //D::mInitDone = 0; // GRAVTEST
2553                 // */
2554
2555                 myTime_t timeend = getTime();
2556                 //debMsgDirect(" done, "<<((timeend-timestart)/(double)1000.0)<<"s \n");
2557 #undef  NBFLAG
2558         }
2559 }
2560
2561
2562
2563 /*****************************************************************************/
2564 /* init a given cell with flag, density, mass and equilibrium dist. funcs */
2565 template<class D>
2566 void 
2567 LbmFsgrSolver<D>::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass) {
2568   /* init eq. dist funcs */
2569         LbmFloat *ecel;
2570         int workSet = mLevel[level].setCurr;
2571
2572         ecel = RACPNT(level, i,j,k, workSet);
2573         FORDF0 { RAC(ecel, l) = D::dfEquil[l] * rho; }
2574         RAC(ecel, dMass) = mass;
2575         RAC(ecel, dFfrac) = mass/rho;
2576         RAC(ecel, dFlux) = FLUX_INIT;
2577         RFLAG(level, i,j,k, workSet)= flag;
2578
2579   workSet ^= 1;
2580 #if COMPRESSGRIDS==0
2581         /*ecel = RACPNT(level, i,j,k, workSet); FIXME why doesnt this work with COMPRESSGRIDS off?
2582         FORDF0 { RAC(ecel, l) = D::dfEquil[l] * rho; }
2583         RAC(ecel, dMass) = mass;
2584         RAC(ecel, dFfrac) = mass/rho;
2585         RAC(ecel, dFlux) = FLUX_INIT; // COMPRT OFF */
2586 #endif // COMPRESSGRIDS==0
2587         RFLAG(level, i,j,k, workSet)= flag; 
2588         return;
2589 }
2590
2591 template<class D>
2592 void 
2593 LbmFsgrSolver<D>::initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel) {
2594         LbmFloat *ecel;
2595         int workSet = mLevel[level].setCurr;
2596
2597         ecel = RACPNT(level, i,j,k, workSet);
2598         FORDF0 { RAC(ecel, l) = D::getCollideEq(l, rho,vel[0],vel[1],vel[2]); }
2599         RAC(ecel, dMass) = mass;
2600         RAC(ecel, dFfrac) = mass/rho;
2601         RAC(ecel, dFlux) = FLUX_INIT;
2602         RFLAG(level, i,j,k, workSet) = flag;
2603
2604   workSet ^= 1;
2605 #if COMPRESSGRIDS==0
2606         /*ecel = RACPNT(level, i,j,k, workSet); FIXME why doesnt this work with COMPRESSGRIDS off?
2607         FORDF0 { RAC(ecel, l) = D::getCollideEq(l, rho,vel[0],vel[1],vel[2]); }
2608         RAC(ecel, dMass) = mass;
2609         RAC(ecel, dFfrac) = mass/rho;
2610         RAC(ecel, dFlux) = FLUX_INIT; // COMPRT OFF */
2611 #endif // COMPRESSGRIDS==0
2612         RFLAG(level, i,j,k, workSet) = flag;
2613         return;
2614 }
2615
2616 template<class D>
2617 bool 
2618 LbmFsgrSolver<D>::checkSymmetry(string idstring)
2619 {
2620         bool erro = false;
2621         bool symm = true;
2622         int msgs = 0;
2623         const int maxMsgs = 10;
2624         const bool markCells = false;
2625
2626         //for(int lev=0; lev<=mMaxRefine; lev++) {
2627         { int lev = mMaxRefine;
2628
2629         // no point if not symm.
2630         if( (mLevel[lev].lSizex==mLevel[lev].lSizey) && (mLevel[lev].lSizex==mLevel[lev].lSizez)) {
2631                 // ok
2632         } else {
2633                 return false;
2634         }
2635
2636         for(int s=0; s<2; s++) {
2637         FSGR_FORIJK1(lev) {
2638                 if(i<(mLevel[lev].lSizex/2)) {
2639                         int inb = (mLevel[lev].lSizey-1-i); 
2640
2641                         if(lev==mMaxRefine) inb -= 1;           // FSGR_SYMM_T
2642
2643                         if( RFLAG(lev, i,j,k,s) != RFLAG(lev, inb,j,k,s) ) { erro = true;
2644                                 if(D::cDimension==2) {
2645                                         if(msgs<maxMsgs) { msgs++;
2646                                                 errMsg("EFLAG", PRINT_IJK<<"s"<<s<<" flag "<<RFLAG(lev, i,j,k,s)<<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" flag "<<RFLAG(lev, inb,j,k,s) );
2647                                         }
2648                                 }
2649                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2650                                 symm = false;
2651                         }
2652                         if( LBM_FLOATNEQ(QCELL(lev, i,j,k,s, dMass), QCELL(lev, inb,j,k,s, dMass)) ) { erro = true;
2653                                 if(D::cDimension==2) {
2654                                         if(msgs<maxMsgs) { msgs++;
2655                                                 /*
2656                                                 debMsgDirect(" mass1 "<<QCELL(lev, i,j,k,s, dMass)<<" mass2 "<<QCELL(lev, inb,j,k,s, dMass) <<std::endl);
2657                                                 errMsg("EMASS", PRINT_IJK<<"s"<<s<<" mass "<<QCELL(lev, i,j,k,s, dMass)<<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" mass "<<QCELL(lev, inb,j,k,s, dMass) );
2658                                                 */
2659                                         }
2660                                 }
2661                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2662                                 symm = false;
2663                         }
2664
2665                         LbmFloat nbrho = QCELL(lev, i,j,k, s, dC);
2666                         FORDF1 { nbrho += QCELL(lev, i,j,k, s, l); }
2667                         LbmFloat otrho = QCELL(lev, inb,j,k, s, dC);
2668                         FORDF1 { otrho += QCELL(lev, inb,j,k, s, l); }
2669                         if( LBM_FLOATNEQ(nbrho, otrho) ) { erro = true;
2670                                 if(D::cDimension==2) {
2671                                         if(msgs<maxMsgs) { msgs++;
2672                                                 //debMsgDirect(" rho 1 "<<nbrho <<" rho 2 "<<otrho  <<std::endl);
2673                                                 errMsg("ERHO ", PRINT_IJK<<"s"<<s<<" rho  "<<nbrho <<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" rho  "<<otrho  );
2674                                         }
2675                                 }
2676                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2677                                 symm = false;
2678                         }
2679                 }
2680         } }
2681         } // lev
2682         LbmFloat maxdiv =0.0;
2683         if(erro) {
2684                 errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
2685                 //if(D::cDimension==2) D::mPanic = true; 
2686                 //return false;
2687         } else {
2688                 errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );
2689         }
2690         // all ok...
2691         return symm;
2692 }// */
2693
2694
2695 /*****************************************************************************/
2696 /*! debug object display */
2697 /*****************************************************************************/
2698 template<class D>
2699 vector<ntlGeometryObject*> LbmFsgrSolver<D>::getDebugObjects() { 
2700         vector<ntlGeometryObject*> debo; 
2701         if(mOutputSurfacePreview) {
2702                 debo.push_back( mpPreviewSurface );
2703         }
2704         return debo; 
2705 }
2706
2707 /*****************************************************************************/
2708 /*! perform a single LBM step */
2709 /*****************************************************************************/
2710
2711 template<class D>
2712 void 
2713 LbmFsgrSolver<D>::stepMain()
2714 {
2715 #if ELBEEM_BLENDER==1
2716                 // update gui display
2717                 SDL_mutexP(globalBakeLock);
2718                 if(globalBakeState<0) {
2719                         // this means abort... cause panic
2720                         D::mPanic = 1;
2721                         errMsg("LbmFsgrSolver::step","Got abort signal from GUI, causing panic, aborting...");
2722                 }
2723                 SDL_mutexV(globalBakeLock);
2724 #endif // ELBEEM_BLENDER==1
2725         D::markedClearList(); // DMC clearMarkedCellsList
2726         if(mDropping) {
2727                 initDrop(mDropX, mDropY);
2728         }
2729         if(mGfxGeoSetup==6) {
2730                 // xobs init hack
2731                 if(mSimulationTime<0.400) {
2732                         if((mSimulationTime>0.25) && (mSimulationTime<0.325)) {
2733                                 // stop shortly...
2734                                 mDropping = false;
2735                         } else {
2736                                 initDrop(0.0, 1.0);
2737                         }
2738                 } else {
2739                         mDropping=false;
2740                 }
2741         }
2742
2743         // safety check, counter reset
2744         D::mNumUsedCells = 0;
2745         mNumInterdCells = 0;
2746         mNumInvIfCells = 0;
2747
2748   //debugOutNnl("LbmFsgrSolver::step : "<<D::mStepCnt, 10);
2749   if(!D::mSilent){ debMsgNnl("LbmFsgrSolver::step", DM_MSG, D::mName<<" cnt:"<<D::mStepCnt<<"  ", 10); }
2750         //debMsgDirect(  "LbmFsgrSolver::step : "<<D::mStepCnt<<" ");
2751         myTime_t timestart = getTime();
2752         //myTime_t timestart = 0;
2753         //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 
2754
2755         // important - keep for tadap
2756         mCurrentMass = D::mFixMass; // reset here for next step
2757         mCurrentVolume = 0.0;
2758         
2759         //stats
2760         mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0;
2761
2762         //change to single step advance!
2763         int dsbits = D::mStepCnt ^ (D::mStepCnt-1);
2764         //errMsg("S"," step:"<<D::mStepCnt<<" s-1:"<<(D::mStepCnt-1)<<" xf:"<<convertFlags2String(dsbits));
2765         for(int lev=0; lev<=mMaxRefine; lev++) {
2766                 //if(! (D::mStepCnt&(1<<lev)) ) {
2767                 if( dsbits & (1<<(mMaxRefine-lev)) ) {
2768                         //errMsg("S"," l"<<lev);
2769
2770                         if(lev==mMaxRefine) {
2771                                 // always advance fine level...
2772                                 fineAdvance();
2773                                 //performRefinement(lev-1); // TEST here?
2774                         } else {
2775                                 performRefinement(lev); // TEST here?
2776                                 performCoarsening(lev); // TEST here?
2777                                 coarseRestrictFromFine(lev);
2778                                 coarseAdvance(lev);
2779                                 //performRefinement(lev-1); // TEST here?
2780                         }
2781 #if FSGR_OMEGA_DEBUG==1
2782                         errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) );
2783                         mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0;
2784 #endif // FSGR_OMEGA_DEBUG==1
2785
2786                         LbmFloat interTime = -10.0;
2787 #if TIMEINTORDER==1
2788                         interTime = 0.5;
2789                         if( D::mStepCnt & (1<<(mMaxRefine-lev)) ) interTime = 0.0;
2790                         // TEST influence... interTime = 0.0;
2791 #elif TIMEINTORDER==2
2792                         interTime = 1.0;
2793 #else // TIMEINTORDER==0
2794                         interTime = 0.0;
2795 #endif // TIMEINTORDER==1
2796
2797                         // test mix! , double 0.5 stablest?
2798 #if INTCFCOARSETEST==1
2799                         //if(lev!=mMaxRefine) { interpolateFineFromCoarse(lev,interTime); } // test!
2800 #else // INTCFCOARSETEST==1
2801                         interpolateFineFromCoarse(lev,interTime);
2802                         ooo
2803 #endif // INTCFCOARSETEST==1
2804                         // */
2805
2806                 }
2807                 mCurrentMass   += mLevel[lev].lmass;
2808                 mCurrentVolume += mLevel[lev].lvolume;
2809         }
2810
2811   // prepare next step
2812         D::mStepCnt++;
2813
2814
2815         // some dbugging output follows
2816         // calculate MLSUPS
2817         myTime_t timeend = getTime();
2818
2819         D::mNumUsedCells += mNumInterdCells; // count both types for MLSUPS
2820         mAvgNumUsedCells += D::mNumUsedCells;
2821         D::mMLSUPS = (D::mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000);
2822         if(D::mMLSUPS>10000){ D::mMLSUPS = -1; }
2823         else { mAvgMLSUPS += D::mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
2824         
2825         LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1(mMaxRefine))) / ((timeend-timestart)/(double)1000.0) ) / (1000000);
2826         if(totMLSUPS>10000) totMLSUPS = -1;
2827         mNumInvIfTotal += mNumInvIfCells; // debug
2828
2829   // do some formatting 
2830   if(!D::mSilent){ 
2831                 string sepStr(""); // DEBUG
2832                 debMsgDirect( 
2833                         "mlsups(curr:"<<D::mMLSUPS<<
2834                         " avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<<
2835                         " totcls:"<<(D::mNumUsedCells)<< sepStr<<
2836                         " avgcls:"<< (int)(mAvgNumUsedCells/(long long int)D::mStepCnt)<< sepStr<<
2837                         " intd:"<<mNumInterdCells<< sepStr<<
2838                         " invif:"<<mNumInvIfCells<< sepStr<<
2839                         " invift:"<<mNumInvIfTotal<< sepStr<<
2840                         " filled:"<<D::mNumFilledCells<<", emptied:"<<D::mNumEmptiedCells<< sepStr<<
2841                         " mMxv:"<<mMxvx<<","<<mMxvy<<","<<mMxvz<<", tscnts:"<<mTimeSwitchCounts<< sepStr<<
2842                         /*" rhoMax:"<<mRhoMax<<", rhoMin:"<<mRhoMin<<", vlenMax:"<<mMaxVlen<<", "*/
2843                         " probs:"<<mNumProblems<< sepStr<<
2844                         " simt:"<<mSimulationTime<< sepStr<<
2845                         " for '"<<D::mName<<"' " );
2846
2847                 //wrong?
2848                 debMsgDirect(", dccd:"<< mCurrentMass<<"/"<<mCurrentVolume<<"("<<D::mFixMass<<")" );
2849                 debMsgDirect(std::endl);
2850
2851                 // nicer output
2852                 debMsgDirect(std::endl); // 
2853                 //debMsgStd(" ",DM_MSG," ",10);
2854         } 
2855         // else {
2856                 //debMsgDirect(".");
2857                 //if((mStepCnt%10)==9) debMsgDirect("\n");
2858         //}
2859
2860         if(D::mStepCnt==1) {
2861                 mMinNoCells = mMaxNoCells = D::mNumUsedCells;
2862         } else {
2863                 if(D::mNumUsedCells>mMaxNoCells) mMaxNoCells = D::mNumUsedCells;
2864                 if(D::mNumUsedCells<mMinNoCells) mMinNoCells = D::mNumUsedCells;
2865         }
2866         
2867         // mass scale test
2868         if(mMaxRefine>0) {
2869                 LbmFloat mscale = mInitialMass/mCurrentMass;
2870
2871                 mscale = 1.0;
2872                 const LbmFloat dchh = 0.001;
2873                 if(mCurrentMass<mInitialMass) mscale = 1.0+dchh;
2874                 if(mCurrentMass>mInitialMass) mscale = 1.0-dchh;
2875
2876                 const int MSInter = 2;
2877                 static int mscount = 0;
2878                 if( (1) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){
2879                         // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37
2880                         errMsg("MDTDD","\n\n");
2881                         errMsg("MDTDD","FORCE RESCALE MASS! "
2882                                         <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
2883                                         <<" step:"<<D::mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
2884                                         );
2885                         errMsg("MDTDD","\n\n");
2886
2887                         mscount++;
2888                         for(int lev=mMaxRefine; lev>=0 ; lev--) {
2889                                 //for(int workSet = 0; workSet<=1; workSet++) {
2890                                 int wss = 0;
2891                                 int wse = 1;
2892 #if COMPRESSGRIDS==1
2893                                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
2894 #endif // COMPRESSGRIDS==1
2895                                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
2896
2897                                         FSGR_FORIJK1(lev) {
2898                                                 if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 
2899                                                         ) {
2900
2901                                                         FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; }
2902                                                         QCELL(lev, i,j,k,workSet, dMass) *= mscale;
2903                                                         QCELL(lev, i,j,k,workSet, dFfrac) *= mscale;
2904
2905                                                 } else {
2906                                                         continue;
2907                                                 }
2908                                         }
2909                                 }
2910                                 mLevel[lev].lmass *= mscale;
2911                         }
2912                 } 
2913
2914                 mCurrentMass *= mscale;
2915         }// if mass scale test */
2916
2917         // one of the last things to do - adapt timestep
2918         // was in fineAdvance before... 
2919         if(mTimeAdap) {
2920                 adaptTimestep();
2921         } // time adaptivity
2922
2923