287b73c77b2920146866cc532ea5786fc75cb134
[blender-staging.git] / intern / elbeem / intern / solver_relax.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-2006 Nils Thuerey
7  *
8  * Combined 2D/3D Lattice Boltzmann relaxation macros
9  *
10  *****************************************************************************/
11
12 #if FSGR_STRICT_DEBUG==1
13 #define CAUSE_PANIC { this->mPanic=1; /* *((int*)(0x0)) = 1; crash*/ }
14 #else // FSGR_STRICT_DEBUG==1
15 #define CAUSE_PANIC { this->mPanic=1; } /*set flag*/
16 #endif // FSGR_STRICT_DEBUG==1
17
18 #if LBM_INCLUDE_TESTSOLVERS!=1
19
20 #define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
21         ux += (grav)[0]; \
22         uy += (grav)[1]; \
23         uz += (grav)[2]; 
24
25 #define TEST_IF_CHECK 
26
27 #else // LBM_INCLUDE_TESTSOLVERS!=1
28 // defined in test.h
29
30 #define NEWDIRVELMOTEST 0
31 #if NEWDIRVELMOTEST==1
32 // off for non testing
33 #undef PRECOLLIDE_MODS
34 #define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
35         ux += (grav)[0]; \
36         uy += (grav)[1]; \
37         uz += (grav)[2];  \
38   { \
39                 int lev = mMaxRefine, nomb=0; \
40                 LbmFloat bcnt = 0.,nux=0.,nuy=0.,nuz=0.; \
41                 for(int l=1; l<this->cDfNum; l++) {  \
42                         if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { \
43                                 if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBndMoving) { \
44                                         nux += QCELL_NB(lev, i,j,k,SRCS(lev),l, dMass); \
45                                         nuy += QCELL_NB(lev, i,j,k,SRCS(lev),l, dFfrac); \
46                                         bcnt += 1.; \
47                                 }       else { \
48                                         nomb++; \
49                                 } \
50                         }  \
51                 } \
52                 if((bcnt>0.)&&(nomb==0)) { \
53                         ux = nux/bcnt; \
54                         uy = nuy/bcnt; \
55                         uz = nuz/bcnt;  \
56                 } \
57         } 
58 #else // NEWDIRVELMOTEST==1
59 // off for non testing
60 #endif // NEWDIRVELMOTEST==1
61
62 #endif // LBM_INCLUDE_TESTSOLVERS!=1
63
64         
65 /******************************************************************************
66  * normal relaxation
67  *****************************************************************************/
68
69 // standard arrays
70 #define CSRC_C    RAC(ccel                                , dC )
71 #define CSRC_E    RAC(ccel + (-1)             *(dTotalNum), dE )
72 #define CSRC_W    RAC(ccel + (+1)             *(dTotalNum), dW )
73 #define CSRC_N    RAC(ccel + (-mLevel[lev].lOffsx)        *(dTotalNum), dN )
74 #define CSRC_S    RAC(ccel + (+mLevel[lev].lOffsx)        *(dTotalNum), dS )
75 #define CSRC_NE   RAC(ccel + (-mLevel[lev].lOffsx-1)      *(dTotalNum), dNE)
76 #define CSRC_NW   RAC(ccel + (-mLevel[lev].lOffsx+1)      *(dTotalNum), dNW)
77 #define CSRC_SE   RAC(ccel + (+mLevel[lev].lOffsx-1)      *(dTotalNum), dSE)
78 #define CSRC_SW   RAC(ccel + (+mLevel[lev].lOffsx+1)      *(dTotalNum), dSW)
79 #define CSRC_T    RAC(ccel + (-mLevel[lev].lOffsy)        *(dTotalNum), dT )
80 #define CSRC_B    RAC(ccel + (+mLevel[lev].lOffsy)        *(dTotalNum), dB )
81 #define CSRC_ET   RAC(ccel + (-mLevel[lev].lOffsy-1)      *(dTotalNum), dET)
82 #define CSRC_EB   RAC(ccel + (+mLevel[lev].lOffsy-1)      *(dTotalNum), dEB)
83 #define CSRC_WT   RAC(ccel + (-mLevel[lev].lOffsy+1)      *(dTotalNum), dWT)
84 #define CSRC_WB   RAC(ccel + (+mLevel[lev].lOffsy+1)      *(dTotalNum), dWB)
85 #define CSRC_NT   RAC(ccel + (-mLevel[lev].lOffsy-mLevel[lev].lOffsx) *(dTotalNum), dNT)
86 #define CSRC_NB   RAC(ccel + (+mLevel[lev].lOffsy-mLevel[lev].lOffsx) *(dTotalNum), dNB)
87 #define CSRC_ST   RAC(ccel + (-mLevel[lev].lOffsy+mLevel[lev].lOffsx) *(dTotalNum), dST)
88 #define CSRC_SB   RAC(ccel + (+mLevel[lev].lOffsy+mLevel[lev].lOffsx) *(dTotalNum), dSB)
89
90 #define XSRC_C(x)    RAC(ccel + (x)                 *dTotalNum, dC )
91 #define XSRC_E(x)    RAC(ccel + ((x)-1)             *dTotalNum, dE )
92 #define XSRC_W(x)    RAC(ccel + ((x)+1)             *dTotalNum, dW )
93 #define XSRC_N(x)    RAC(ccel + ((x)-mLevel[lev].lOffsx)        *dTotalNum, dN )
94 #define XSRC_S(x)    RAC(ccel + ((x)+mLevel[lev].lOffsx)        *dTotalNum, dS )
95 #define XSRC_NE(x)   RAC(ccel + ((x)-mLevel[lev].lOffsx-1)      *dTotalNum, dNE)
96 #define XSRC_NW(x)   RAC(ccel + ((x)-mLevel[lev].lOffsx+1)      *dTotalNum, dNW)
97 #define XSRC_SE(x)   RAC(ccel + ((x)+mLevel[lev].lOffsx-1)      *dTotalNum, dSE)
98 #define XSRC_SW(x)   RAC(ccel + ((x)+mLevel[lev].lOffsx+1)      *dTotalNum, dSW)
99 #define XSRC_T(x)    RAC(ccel + ((x)-mLevel[lev].lOffsy)        *dTotalNum, dT )
100 #define XSRC_B(x)    RAC(ccel + ((x)+mLevel[lev].lOffsy)        *dTotalNum, dB )
101 #define XSRC_ET(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy-1)      *dTotalNum, dET)
102 #define XSRC_EB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy-1)      *dTotalNum, dEB)
103 #define XSRC_WT(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy+1)      *dTotalNum, dWT)
104 #define XSRC_WB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy+1)      *dTotalNum, dWB)
105 #define XSRC_NT(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy-mLevel[lev].lOffsx) *dTotalNum, dNT)
106 #define XSRC_NB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy-mLevel[lev].lOffsx) *dTotalNum, dNB)
107 #define XSRC_ST(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy+mLevel[lev].lOffsx) *dTotalNum, dST)
108 #define XSRC_SB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy+mLevel[lev].lOffsx) *dTotalNum, dSB)
109
110
111
112 #define OMEGA(l) mLevel[(l)].omega
113
114 #define EQC (  DFL1*(rho - usqr))
115 #define EQN (  DFL2*(rho + uy*(4.5*uy + 3.0) - usqr))
116 #define EQS (  DFL2*(rho + uy*(4.5*uy - 3.0) - usqr))
117 #define EQE (  DFL2*(rho + ux*(4.5*ux + 3.0) - usqr))
118 #define EQW (  DFL2*(rho + ux*(4.5*ux - 3.0) - usqr))
119 #define EQT (  DFL2*(rho + uz*(4.5*uz + 3.0) - usqr))
120 #define EQB (  DFL2*(rho + uz*(4.5*uz - 3.0) - usqr))
121                     
122 #define EQNE ( DFL3*(rho + (+ux+uy)*(4.5*(+ux+uy) + 3.0) - usqr))
123 #define EQNW ( DFL3*(rho + (-ux+uy)*(4.5*(-ux+uy) + 3.0) - usqr))
124 #define EQSE ( DFL3*(rho + (+ux-uy)*(4.5*(+ux-uy) + 3.0) - usqr))
125 #define EQSW ( DFL3*(rho + (-ux-uy)*(4.5*(-ux-uy) + 3.0) - usqr))
126 #define EQNT ( DFL3*(rho + (+uy+uz)*(4.5*(+uy+uz) + 3.0) - usqr))
127 #define EQNB ( DFL3*(rho + (+uy-uz)*(4.5*(+uy-uz) + 3.0) - usqr))
128 #define EQST ( DFL3*(rho + (-uy+uz)*(4.5*(-uy+uz) + 3.0) - usqr))
129 #define EQSB ( DFL3*(rho + (-uy-uz)*(4.5*(-uy-uz) + 3.0) - usqr))
130 #define EQET ( DFL3*(rho + (+ux+uz)*(4.5*(+ux+uz) + 3.0) - usqr))
131 #define EQEB ( DFL3*(rho + (+ux-uz)*(4.5*(+ux-uz) + 3.0) - usqr))
132 #define EQWT ( DFL3*(rho + (-ux+uz)*(4.5*(-ux+uz) + 3.0) - usqr))
133 #define EQWB ( DFL3*(rho + (-ux-uz)*(4.5*(-ux-uz) + 3.0) - usqr))
134
135
136 // this is a bit ugly, but necessary for the CSRC_ access...
137 #define MSRC_C    m[dC ]
138 #define MSRC_N    m[dN ]
139 #define MSRC_S    m[dS ]
140 #define MSRC_E    m[dE ]
141 #define MSRC_W    m[dW ]
142 #define MSRC_T    m[dT ]
143 #define MSRC_B    m[dB ]
144 #define MSRC_NE   m[dNE]
145 #define MSRC_NW   m[dNW]
146 #define MSRC_SE   m[dSE]
147 #define MSRC_SW   m[dSW]
148 #define MSRC_NT   m[dNT]
149 #define MSRC_NB   m[dNB]
150 #define MSRC_ST   m[dST]
151 #define MSRC_SB   m[dSB]
152 #define MSRC_ET   m[dET]
153 #define MSRC_EB   m[dEB]
154 #define MSRC_WT   m[dWT]
155 #define MSRC_WB   m[dWB]
156
157 // this is a bit ugly, but necessary for the ccel local access...
158 #define CCEL_C    RAC(ccel, dC )
159 #define CCEL_N    RAC(ccel, dN )
160 #define CCEL_S    RAC(ccel, dS )
161 #define CCEL_E    RAC(ccel, dE )
162 #define CCEL_W    RAC(ccel, dW )
163 #define CCEL_T    RAC(ccel, dT )
164 #define CCEL_B    RAC(ccel, dB )
165 #define CCEL_NE   RAC(ccel, dNE)
166 #define CCEL_NW   RAC(ccel, dNW)
167 #define CCEL_SE   RAC(ccel, dSE)
168 #define CCEL_SW   RAC(ccel, dSW)
169 #define CCEL_NT   RAC(ccel, dNT)
170 #define CCEL_NB   RAC(ccel, dNB)
171 #define CCEL_ST   RAC(ccel, dST)
172 #define CCEL_SB   RAC(ccel, dSB)
173 #define CCEL_ET   RAC(ccel, dET)
174 #define CCEL_EB   RAC(ccel, dEB)
175 #define CCEL_WT   RAC(ccel, dWT)
176 #define CCEL_WB   RAC(ccel, dWB)
177 // for coarse to fine interpol access
178 #define CCELG_C(f)    (RAC(ccel, dC )*mGaussw[(f)])
179 #define CCELG_N(f)    (RAC(ccel, dN )*mGaussw[(f)])
180 #define CCELG_S(f)    (RAC(ccel, dS )*mGaussw[(f)])
181 #define CCELG_E(f)    (RAC(ccel, dE )*mGaussw[(f)])
182 #define CCELG_W(f)    (RAC(ccel, dW )*mGaussw[(f)])
183 #define CCELG_T(f)    (RAC(ccel, dT )*mGaussw[(f)])
184 #define CCELG_B(f)    (RAC(ccel, dB )*mGaussw[(f)])
185 #define CCELG_NE(f)   (RAC(ccel, dNE)*mGaussw[(f)])
186 #define CCELG_NW(f)   (RAC(ccel, dNW)*mGaussw[(f)])
187 #define CCELG_SE(f)   (RAC(ccel, dSE)*mGaussw[(f)])
188 #define CCELG_SW(f)   (RAC(ccel, dSW)*mGaussw[(f)])
189 #define CCELG_NT(f)   (RAC(ccel, dNT)*mGaussw[(f)])
190 #define CCELG_NB(f)   (RAC(ccel, dNB)*mGaussw[(f)])
191 #define CCELG_ST(f)   (RAC(ccel, dST)*mGaussw[(f)])
192 #define CCELG_SB(f)   (RAC(ccel, dSB)*mGaussw[(f)])
193 #define CCELG_ET(f)   (RAC(ccel, dET)*mGaussw[(f)])
194 #define CCELG_EB(f)   (RAC(ccel, dEB)*mGaussw[(f)])
195 #define CCELG_WT(f)   (RAC(ccel, dWT)*mGaussw[(f)])
196 #define CCELG_WB(f)   (RAC(ccel, dWB)*mGaussw[(f)])
197
198
199 #if PARALLEL==1
200 #define CSMOMEGA_STATS(dlev, domega) 
201 #else // PARALLEL==1
202 #if FSGR_OMEGA_DEBUG==1
203 #define CSMOMEGA_STATS(dlev, domega) \
204         mLevel[dlev].avgOmega += domega; mLevel[dlev].avgOmegaCnt+=1.0; 
205 #else // FSGR_OMEGA_DEBUG==1
206 #define CSMOMEGA_STATS(dlev, domega) 
207 #endif // FSGR_OMEGA_DEBUG==1
208 #endif // PARALLEL==1
209
210
211 // used for main loops and grav init
212 // source set
213 #define SRCS(l) mLevel[(l)].setCurr
214 // target set
215 #define TSET(l) mLevel[(l)].setOther
216
217 // handle mov. obj 
218 #if FSGR_STRICT_DEBUG==1
219
220 #define  LBMDS_ADDMOV(linv,l)  \
221          \
222         if((nbflag[linv]&CFBndMoving)&&(!(nbflag[l]&CFBnd))){ \
223          \
224         LbmFloat dte=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)-(mSimulationTime+this->mpParam->getTimestep()); \
225         if( ABS(dte)< 1e-15 ) { \
226         m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
227         } else { \
228         const int sdx = i+this->dfVecX[linv], sdy = j+this->dfVecY[linv], sdz = k+this->dfVecZ[linv]; \
229          \
230         errMsg("INVALID_MOV_OBJ_TIME"," at "<<PRINT_IJK<<" from l"<<l<<" "<<PRINT_VEC(sdx,sdy,sdz)<<" t="<<(mSimulationTime+this->mpParam->getTimestep())<<" ct="<<QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)<<" dte="<<dte); \
231         debugMarkCell(lev,sdx,sdy,sdz); \
232         } \
233         } \
234
235
236
237 #else // FSGR_STRICT_DEBUG==1
238
239 #define  LBMDS_ADDMOV(linv,l)  \
240          \
241          \
242         if((nbflag[linv]&CFBndMoving)&&(!(nbflag[l]&CFBnd))){ \
243          \
244         m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
245         } \
246
247
248
249 #endif // !FSGR_STRICT_DEBUG==1
250
251 // treatment of freeslip reflection
252 // used both for OPT and nonOPT
253 #define  DEFAULT_STREAM_FREESLIP(l,invl,mnbf)  \
254          \
255         int nb1 = 0, nb2 = 0; \
256         LbmFloat newval = 0.0; \
257         const int dx = this->dfVecX[invl], dy = this->dfVecY[invl], dz = this->dfVecZ[invl]; \
258          \
259          \
260          \
261         const LbmFloat movadd = ( \
262         ((nbflag[invl]&CFBndMoving)&&(!(nbflag[l]&CFBnd))) ? \
263         (QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l)) : 0.); \
264          \
265         if(dz==0) { \
266         nb1 = !(RFLAG(lev, i,   j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
267         nb2 = !(RFLAG(lev, i+dx,j,   k, SRCS(lev))&(CFFluid|CFInter)); \
268         if((nb1)&&(!nb2)) { \
269          \
270         newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
271         } else \
272         if((!nb1)&&(nb2)) { \
273          \
274         newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
275         } else { \
276          \
277         newval = RAC(ccel, this->dfInv[l] ) +movadd /* */; \
278         } \
279         } else \
280         if(dy==0) { \
281         nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
282         nb2 = !(RFLAG(lev, i+dx,j,k, SRCS(lev))&(CFFluid|CFInter)); \
283         if((nb1)&&(!nb2)) { \
284          \
285         newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
286         } else \
287         if((!nb1)&&(nb2)) { \
288          \
289         newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
290         } else { \
291          \
292         newval = RAC(ccel, this->dfInv[l] )  +movadd /* */; \
293         } \
294          \
295         } else \
296          \
297         { \
298          \
299         nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
300         nb2 = !(RFLAG(lev, i,j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
301         if((nb1)&&(!nb2)) { \
302          \
303         newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
304         } else \
305         if((!nb1)&&(nb2)) { \
306          \
307         newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
308         } else { \
309          \
310         newval = RAC(ccel, this->dfInv[l] )  +movadd /* */; \
311         } \
312         } \
313          \
314         if(mnbf & CFBndPartslip) { \
315         const LbmFloat partv = mObjectPartslips[(int)(mnbf>>24)]; \
316          \
317         m[l] = (RAC(ccel, this->dfInv[l] )  +movadd /* d *(1./1.) */ ) * partv + newval * (1.0-partv); \
318         } else { \
319         m[l] = newval; \
320         } \
321          \
322
323
324
325
326 // complete default stream&collide, 2d/3d
327 /* read distribution funtions of adjacent cells = sweep step */ 
328 #if OPT3D==0 
329
330 #if FSGR_STRICT_DEBUG==1
331 #define MARKCELLCHECK \
332         debugMarkCell(lev,i,j,k); CAUSE_PANIC;
333 #define STREAMCHECK(id,ni,nj,nk,nl) \
334         if((!(m[nl] > -1.0) && (m[nl]<1.0)) ) {\
335                 errMsg("STREAMCHECK","ID"<<id<<" Invalid streamed DF nl"<<nl<<" value:"<<m[nl]<<" at "<<PRINT_IJK<<" from "<<PRINT_VEC(ni,nj,nk)<<" nl"<<(nl)<<\
336                                 " nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther)  ); \
337                 /*FORDF0{ errMsg("STREAMCHECK"," at "<<PRINT_IJK<<" df "<<l<<"="<<m[l] ); } */ \
338                 MARKCELLCHECK; \
339                 m[nl] = dfEquil[nl]; /* REPAIR */ \
340         }
341 #define COLLCHECK \
342         if( (rho>2.0) || (rho<-1.0) || (ABS(ux)>1.0) || (ABS(uy)>1.0) |(ABS(uz)>1.0) ) {\
343                 errMsg("COLLCHECK","Invalid collision values r:"<<rho<<" u:"PRINT_VEC(ux,uy,uz)<<" at? "<<PRINT_IJK ); \
344                 /*FORDF0{ errMsg("COLLCHECK"," at? "<<PRINT_IJK<<" df "<<l<<"="<<m[l] ); }*/ \
345                 rho=ux=uy=uz= 0.; /* REPAIR */ \
346                 MARKCELLCHECK; \
347         }
348 #else
349 #define STREAMCHECK(id, ni,nj,nk,nl) 
350 #define COLLCHECK
351 #endif
352
353 // careful ux,uy,uz need to be inited before!
354 #define  DEFAULT_STREAM  \
355         m[dC] = RAC(ccel,dC); \
356         STREAMCHECK(1, i,j,k, dC); \
357         FORDF1 { \
358         CellFlagType nbf = nbflag[ this->dfInv[l] ]; \
359         if(nbf & CFBnd) { \
360         if(nbf & CFBndNoslip) { \
361          \
362         m[l] = RAC(ccel, this->dfInv[l] ); \
363         LBMDS_ADDMOV(this->dfInv[l],l); \
364         STREAMCHECK(2, i,j,k, l); \
365         } else if(nbf & (CFBndFreeslip|CFBndPartslip)) { \
366          \
367         if(l<=LBMDIM*2) { \
368         m[l] = RAC(ccel, this->dfInv[l] ); STREAMCHECK(3, i,j,k, l); \
369         LBMDS_ADDMOV(this->dfInv[l],l); \
370         } else { \
371         const int inv_l = this->dfInv[l]; \
372         DEFAULT_STREAM_FREESLIP(l,inv_l,nbf); \
373         } \
374          \
375         } \
376         else { \
377         errMsg("LbmFsgrSolver","Invalid Bnd type at "<<PRINT_IJK<<" f"<<convertCellFlagType2String(nbf)<<",nbdir"<<this->dfInv[l] ); \
378         } \
379         } else { \
380         m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
381         if(RFLAG(lev, i,j,k, mLevel[lev].setCurr)&CFFluid) { \
382         if(!(nbf&(CFFluid|CFInter)) ) { \
383         int ni=i+this->dfVecX[this->dfInv[l]], nj=j+this->dfVecY[this->dfInv[l]], nk=k+this->dfVecZ[this->dfInv[l]]; \
384         errMsg("STREAMCHECK"," Invalid nbflag, streamed DF l"<<l<<" value:"<<m[l]<<" at "<<PRINT_IJK<<" from "<< \
385         PRINT_VEC(ni,nj,nk) <<" l"<<(l)<< \
386         " nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther)  ); \
387          \
388          \
389         } } \
390         STREAMCHECK(4, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
391         } \
392         } \
393
394
395
396
397 // careful ux,uy,uz need to be inited before!
398 #define  DEFAULT_COLLIDEG(grav)  \
399         this->collideArrays(lev, i,j,k, m, rho,ux,uy,uz, OMEGA(lev), grav, mLevel[lev].lcsmago, &mDebugOmegaRet, &lcsmqo ); \
400         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
401         FORDF0 { RAC(tcel,l) = m[l]; } \
402         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
403         COLLCHECK; \
404
405
406
407 #define  OPTIMIZED_STREAMCOLLIDE  \
408         m[0] = RAC(ccel,0); \
409         FORDF1 { \
410          \
411         if(RFLAG_NBINV(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); CAUSE_PANIC; \
412         } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
413         STREAMCHECK(8, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
414         } \
415         rho=m[0]; \
416         DEFAULT_COLLIDEG(mLevel[lev].gravity) \
417
418
419
420 #define  OPTIMIZED_STREAMCOLLIDE___UNUSED  \
421          \
422         this->collideArrays(lev, i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].gravity, mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo   ); \
423         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
424         FORDF0 { RAC(tcel,l) = m[l]; } \
425         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
426         COLLCHECK; \
427
428
429
430 #else  // 3D, opt OPT3D==true
431
432
433 // default stream opt3d add moving bc val
434 #define  DEFAULT_STREAM  \
435         m[dC] = RAC(ccel,dC); \
436          \
437         if((!nbored & CFBnd)) { \
438          \
439         m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
440         m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
441         m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
442         m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
443         m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
444         m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
445         } else { \
446          \
447         if(nbflag[dS ]&CFBnd) { m[dN ] = RAC(ccel,dS ); LBMDS_ADDMOV(dS ,dN ); } else { m[dN ] = CSRC_N ; } \
448         if(nbflag[dN ]&CFBnd) { m[dS ] = RAC(ccel,dN ); LBMDS_ADDMOV(dN ,dS ); } else { m[dS ] = CSRC_S ; } \
449         if(nbflag[dW ]&CFBnd) { m[dE ] = RAC(ccel,dW ); LBMDS_ADDMOV(dW ,dE ); } else { m[dE ] = CSRC_E ; } \
450         if(nbflag[dE ]&CFBnd) { m[dW ] = RAC(ccel,dE ); LBMDS_ADDMOV(dE ,dW ); } else { m[dW ] = CSRC_W ; } \
451         if(nbflag[dB ]&CFBnd) { m[dT ] = RAC(ccel,dB ); LBMDS_ADDMOV(dB ,dT ); } else { m[dT ] = CSRC_T ; } \
452         if(nbflag[dT ]&CFBnd) { m[dB ] = RAC(ccel,dT ); LBMDS_ADDMOV(dT ,dB ); } else { m[dB ] = CSRC_B ; } \
453          \
454          \
455         if(nbflag[dSW]&CFBnd) { if(nbflag[dSW]&CFBndNoslip){ m[dNE] = RAC(ccel,dSW); LBMDS_ADDMOV(dSW,dNE); }else{ DEFAULT_STREAM_FREESLIP(dNE,dSW,nbflag[dSW]);} } else { m[dNE] = CSRC_NE; } \
456         if(nbflag[dSE]&CFBnd) { if(nbflag[dSE]&CFBndNoslip){ m[dNW] = RAC(ccel,dSE); LBMDS_ADDMOV(dSE,dNW); }else{ DEFAULT_STREAM_FREESLIP(dNW,dSE,nbflag[dSE]);} } else { m[dNW] = CSRC_NW; } \
457         if(nbflag[dNW]&CFBnd) { if(nbflag[dNW]&CFBndNoslip){ m[dSE] = RAC(ccel,dNW); LBMDS_ADDMOV(dNW,dSE); }else{ DEFAULT_STREAM_FREESLIP(dSE,dNW,nbflag[dNW]);} } else { m[dSE] = CSRC_SE; } \
458         if(nbflag[dNE]&CFBnd) { if(nbflag[dNE]&CFBndNoslip){ m[dSW] = RAC(ccel,dNE); LBMDS_ADDMOV(dNE,dSW); }else{ DEFAULT_STREAM_FREESLIP(dSW,dNE,nbflag[dNE]);} } else { m[dSW] = CSRC_SW; } \
459         if(nbflag[dSB]&CFBnd) { if(nbflag[dSB]&CFBndNoslip){ m[dNT] = RAC(ccel,dSB); LBMDS_ADDMOV(dSB,dNT); }else{ DEFAULT_STREAM_FREESLIP(dNT,dSB,nbflag[dSB]);} } else { m[dNT] = CSRC_NT; } \
460         if(nbflag[dST]&CFBnd) { if(nbflag[dST]&CFBndNoslip){ m[dNB] = RAC(ccel,dST); LBMDS_ADDMOV(dST,dNB); }else{ DEFAULT_STREAM_FREESLIP(dNB,dST,nbflag[dST]);} } else { m[dNB] = CSRC_NB; } \
461         if(nbflag[dNB]&CFBnd) { if(nbflag[dNB]&CFBndNoslip){ m[dST] = RAC(ccel,dNB); LBMDS_ADDMOV(dNB,dST); }else{ DEFAULT_STREAM_FREESLIP(dST,dNB,nbflag[dNB]);} } else { m[dST] = CSRC_ST; } \
462         if(nbflag[dNT]&CFBnd) { if(nbflag[dNT]&CFBndNoslip){ m[dSB] = RAC(ccel,dNT); LBMDS_ADDMOV(dNT,dSB); }else{ DEFAULT_STREAM_FREESLIP(dSB,dNT,nbflag[dNT]);} } else { m[dSB] = CSRC_SB; } \
463         if(nbflag[dWB]&CFBnd) { if(nbflag[dWB]&CFBndNoslip){ m[dET] = RAC(ccel,dWB); LBMDS_ADDMOV(dWB,dET); }else{ DEFAULT_STREAM_FREESLIP(dET,dWB,nbflag[dWB]);} } else { m[dET] = CSRC_ET; } \
464         if(nbflag[dWT]&CFBnd) { if(nbflag[dWT]&CFBndNoslip){ m[dEB] = RAC(ccel,dWT); LBMDS_ADDMOV(dWT,dEB); }else{ DEFAULT_STREAM_FREESLIP(dEB,dWT,nbflag[dWT]);} } else { m[dEB] = CSRC_EB; } \
465         if(nbflag[dEB]&CFBnd) { if(nbflag[dEB]&CFBndNoslip){ m[dWT] = RAC(ccel,dEB); LBMDS_ADDMOV(dEB,dWT); }else{ DEFAULT_STREAM_FREESLIP(dWT,dEB,nbflag[dEB]);} } else { m[dWT] = CSRC_WT; } \
466         if(nbflag[dET]&CFBnd) { if(nbflag[dET]&CFBndNoslip){ m[dWB] = RAC(ccel,dET); LBMDS_ADDMOV(dET,dWB); }else{ DEFAULT_STREAM_FREESLIP(dWB,dET,nbflag[dET]);} } else { m[dWB] = CSRC_WB; } \
467         } \
468
469
470
471
472
473 #define  COLL_CALCULATE_DFEQ(dstarray)  \
474         dstarray[dN ] = EQN ; dstarray[dS ] = EQS ; \
475         dstarray[dE ] = EQE ; dstarray[dW ] = EQW ; \
476         dstarray[dT ] = EQT ; dstarray[dB ] = EQB ; \
477         dstarray[dNE] = EQNE; dstarray[dNW] = EQNW; dstarray[dSE] = EQSE; dstarray[dSW] = EQSW; \
478         dstarray[dNT] = EQNT; dstarray[dNB] = EQNB; dstarray[dST] = EQST; dstarray[dSB] = EQSB; \
479         dstarray[dET] = EQET; dstarray[dEB] = EQEB; dstarray[dWT] = EQWT; dstarray[dWB] = EQWB; \
480
481
482
483 #define  COLL_CALCULATE_NONEQTENSOR(csolev, srcArray )  \
484         lcsmqadd  = (srcArray##NE - lcsmeq[ dNE ]); \
485         lcsmqadd -= (srcArray##NW - lcsmeq[ dNW ]); \
486         lcsmqadd -= (srcArray##SE - lcsmeq[ dSE ]); \
487         lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
488         lcsmqo = (lcsmqadd*    lcsmqadd); \
489         lcsmqadd  = (srcArray##ET - lcsmeq[  dET ]); \
490         lcsmqadd -= (srcArray##EB - lcsmeq[  dEB ]); \
491         lcsmqadd -= (srcArray##WT - lcsmeq[  dWT ]); \
492         lcsmqadd += (srcArray##WB - lcsmeq[  dWB ]); \
493         lcsmqo += (lcsmqadd*    lcsmqadd); \
494         lcsmqadd  = (srcArray##NT - lcsmeq[  dNT ]); \
495         lcsmqadd -= (srcArray##NB - lcsmeq[  dNB ]); \
496         lcsmqadd -= (srcArray##ST - lcsmeq[  dST ]); \
497         lcsmqadd += (srcArray##SB - lcsmeq[  dSB ]); \
498         lcsmqo += (lcsmqadd*    lcsmqadd); \
499         lcsmqo *= 2.0; \
500         lcsmqadd  = (srcArray##E  -  lcsmeq[ dE  ]); \
501         lcsmqadd += (srcArray##W  -  lcsmeq[ dW  ]); \
502         lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
503         lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
504         lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
505         lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
506         lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
507         lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
508         lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
509         lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
510         lcsmqo += (lcsmqadd*    lcsmqadd); \
511         lcsmqadd  = (srcArray##N  -  lcsmeq[ dN  ]); \
512         lcsmqadd += (srcArray##S  -  lcsmeq[ dS  ]); \
513         lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
514         lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
515         lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
516         lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
517         lcsmqadd += (srcArray##NT  - lcsmeq[ dNT ]); \
518         lcsmqadd += (srcArray##NB  - lcsmeq[ dNB ]); \
519         lcsmqadd += (srcArray##ST  - lcsmeq[ dST ]); \
520         lcsmqadd += (srcArray##SB  - lcsmeq[ dSB ]); \
521         lcsmqo += (lcsmqadd*    lcsmqadd); \
522         lcsmqadd  = (srcArray##T  -  lcsmeq[ dT  ]); \
523         lcsmqadd += (srcArray##B  -  lcsmeq[ dB  ]); \
524         lcsmqadd += (srcArray##NT -  lcsmeq[ dNT ]); \
525         lcsmqadd += (srcArray##NB -  lcsmeq[ dNB ]); \
526         lcsmqadd += (srcArray##ST -  lcsmeq[ dST ]); \
527         lcsmqadd += (srcArray##SB -  lcsmeq[ dSB ]); \
528         lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
529         lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
530         lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
531         lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
532         lcsmqo += (lcsmqadd*    lcsmqadd); \
533         lcsmqo = sqrt(lcsmqo); \
534
535
536
537 //                      COLL_CALCULATE_CSMOMEGAVAL(csolev, lcsmomega); 
538
539 // careful - need lcsmqo 
540 #define  COLL_CALCULATE_CSMOMEGAVAL(csolev, dstomega )  \
541         dstomega =  1.0/ \
542         ( 3.0*( mLevel[(csolev)].lcnu+mLevel[(csolev)].lcsmago_sqr*( \
543         -mLevel[(csolev)].lcnu + sqrt( mLevel[(csolev)].lcnu*mLevel[(csolev)].lcnu + 18.0*mLevel[(csolev)].lcsmago_sqr* lcsmqo ) \
544         / (6.0*mLevel[(csolev)].lcsmago_sqr)) \
545         ) +0.5 ); \
546
547
548
549 #define  DEFAULT_COLLIDE_LES(grav)  \
550         rho = + MSRC_C  + MSRC_N \
551         + MSRC_S  + MSRC_E \
552         + MSRC_W  + MSRC_T \
553         + MSRC_B  + MSRC_NE \
554         + MSRC_NW + MSRC_SE \
555         + MSRC_SW + MSRC_NT \
556         + MSRC_NB + MSRC_ST \
557         + MSRC_SB + MSRC_ET \
558         + MSRC_EB + MSRC_WT \
559         + MSRC_WB; \
560          \
561         ux = MSRC_E - MSRC_W \
562         + MSRC_NE - MSRC_NW \
563         + MSRC_SE - MSRC_SW \
564         + MSRC_ET + MSRC_EB \
565         - MSRC_WT - MSRC_WB ; \
566          \
567         uy = MSRC_N - MSRC_S \
568         + MSRC_NE + MSRC_NW \
569         - MSRC_SE - MSRC_SW \
570         + MSRC_NT + MSRC_NB \
571         - MSRC_ST - MSRC_SB ; \
572          \
573         uz = MSRC_T - MSRC_B \
574         + MSRC_NT - MSRC_NB \
575         + MSRC_ST - MSRC_SB \
576         + MSRC_ET - MSRC_EB \
577         + MSRC_WT - MSRC_WB ; \
578         PRECOLLIDE_MODS(rho,ux,uy,uz, grav); \
579         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
580         COLL_CALCULATE_DFEQ(lcsmeq); \
581         COLL_CALCULATE_NONEQTENSOR(lev, MSRC_); \
582         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
583         CSMOMEGA_STATS(lev,lcsmomega); \
584          \
585         RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
586          \
587         RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ]; \
588         RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ]; \
589         RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
590         RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ]; \
591         RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ]; \
592         RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
593          \
594         RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
595         RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
596         RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
597         RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
598         RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
599         RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
600         RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
601         RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
602         RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
603         RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
604         RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
605         RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \
606
607
608
609 #define  DEFAULT_COLLIDE_NOLES(grav)  \
610         rho = + MSRC_C  + MSRC_N \
611         + MSRC_S  + MSRC_E \
612         + MSRC_W  + MSRC_T \
613         + MSRC_B  + MSRC_NE \
614         + MSRC_NW + MSRC_SE \
615         + MSRC_SW + MSRC_NT \
616         + MSRC_NB + MSRC_ST \
617         + MSRC_SB + MSRC_ET \
618         + MSRC_EB + MSRC_WT \
619         + MSRC_WB; \
620          \
621         ux = MSRC_E - MSRC_W \
622         + MSRC_NE - MSRC_NW \
623         + MSRC_SE - MSRC_SW \
624         + MSRC_ET + MSRC_EB \
625         - MSRC_WT - MSRC_WB ; \
626          \
627         uy = MSRC_N - MSRC_S \
628         + MSRC_NE + MSRC_NW \
629         - MSRC_SE - MSRC_SW \
630         + MSRC_NT + MSRC_NB \
631         - MSRC_ST - MSRC_SB ; \
632          \
633         uz = MSRC_T - MSRC_B \
634         + MSRC_NT - MSRC_NB \
635         + MSRC_ST - MSRC_SB \
636         + MSRC_ET - MSRC_EB \
637         + MSRC_WT - MSRC_WB ; \
638         PRECOLLIDE_MODS(rho, ux,uy,uz, grav); \
639         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
640          \
641         RAC(tcel,dC ) = (1.0-OMEGA(lev))*MSRC_C  + OMEGA(lev)*EQC ; \
642          \
643         RAC(tcel,dN ) = (1.0-OMEGA(lev))*MSRC_N  + OMEGA(lev)*EQN ; \
644         RAC(tcel,dS ) = (1.0-OMEGA(lev))*MSRC_S  + OMEGA(lev)*EQS ; \
645         RAC(tcel,dE ) = (1.0-OMEGA(lev))*MSRC_E  + OMEGA(lev)*EQE ; \
646         RAC(tcel,dW ) = (1.0-OMEGA(lev))*MSRC_W  + OMEGA(lev)*EQW ; \
647         RAC(tcel,dT ) = (1.0-OMEGA(lev))*MSRC_T  + OMEGA(lev)*EQT ; \
648         RAC(tcel,dB ) = (1.0-OMEGA(lev))*MSRC_B  + OMEGA(lev)*EQB ; \
649          \
650         RAC(tcel,dNE) = (1.0-OMEGA(lev))*MSRC_NE + OMEGA(lev)*EQNE; \
651         RAC(tcel,dNW) = (1.0-OMEGA(lev))*MSRC_NW + OMEGA(lev)*EQNW; \
652         RAC(tcel,dSE) = (1.0-OMEGA(lev))*MSRC_SE + OMEGA(lev)*EQSE; \
653         RAC(tcel,dSW) = (1.0-OMEGA(lev))*MSRC_SW + OMEGA(lev)*EQSW; \
654         RAC(tcel,dNT) = (1.0-OMEGA(lev))*MSRC_NT + OMEGA(lev)*EQNT; \
655         RAC(tcel,dNB) = (1.0-OMEGA(lev))*MSRC_NB + OMEGA(lev)*EQNB; \
656         RAC(tcel,dST) = (1.0-OMEGA(lev))*MSRC_ST + OMEGA(lev)*EQST; \
657         RAC(tcel,dSB) = (1.0-OMEGA(lev))*MSRC_SB + OMEGA(lev)*EQSB; \
658         RAC(tcel,dET) = (1.0-OMEGA(lev))*MSRC_ET + OMEGA(lev)*EQET; \
659         RAC(tcel,dEB) = (1.0-OMEGA(lev))*MSRC_EB + OMEGA(lev)*EQEB; \
660         RAC(tcel,dWT) = (1.0-OMEGA(lev))*MSRC_WT + OMEGA(lev)*EQWT; \
661         RAC(tcel,dWB) = (1.0-OMEGA(lev))*MSRC_WB + OMEGA(lev)*EQWB; \
662
663
664
665
666
667 #define  OPTIMIZED_STREAMCOLLIDE_LES  \
668          \
669         m[dC ] = CSRC_C ; \
670         m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
671         m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
672         m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
673         m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
674         m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
675         m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
676          \
677         rho = MSRC_C  + MSRC_N + MSRC_S  + MSRC_E + MSRC_W  + MSRC_T \
678         + MSRC_B  + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT \
679         + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; \
680         ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW \
681         + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB; \
682         uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW \
683         + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB; \
684         uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB \
685         + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB; \
686         PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
687         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
688         COLL_CALCULATE_DFEQ(lcsmeq); \
689         COLL_CALCULATE_NONEQTENSOR(lev, MSRC_) \
690         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
691         CSMOMEGA_STATS(lev,lcsmomega); \
692          \
693         RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
694         RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ]; \
695         RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ]; \
696         RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
697         RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ]; \
698         RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ]; \
699         RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
700          \
701         RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
702         RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
703         RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
704         RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
705          \
706         RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
707         RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
708         RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
709         RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
710          \
711         RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
712         RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
713         RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
714         RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \
715
716
717
718 #define  OPTIMIZED_STREAMCOLLIDE_UNUSED  \
719          \
720         rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T \
721         + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
722         + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
723         ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
724         + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \
725         uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
726         + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \
727         uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
728         + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \
729         PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
730         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
731         COLL_CALCULATE_DFEQ(lcsmeq); \
732         COLL_CALCULATE_NONEQTENSOR(lev, CSRC_) \
733         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
734          \
735         RAC(tcel,dC ) = (1.0-lcsmomega)*CSRC_C  + lcsmomega*EQC ; \
736         RAC(tcel,dN ) = (1.0-lcsmomega)*CSRC_N  + lcsmomega*lcsmeq[ dN ]; \
737         RAC(tcel,dS ) = (1.0-lcsmomega)*CSRC_S  + lcsmomega*lcsmeq[ dS ]; \
738         RAC(tcel,dE ) = (1.0-lcsmomega)*CSRC_E  + lcsmomega*lcsmeq[ dE ]; \
739         RAC(tcel,dW ) = (1.0-lcsmomega)*CSRC_W  + lcsmomega*lcsmeq[ dW ]; \
740         RAC(tcel,dT ) = (1.0-lcsmomega)*CSRC_T  + lcsmomega*lcsmeq[ dT ]; \
741         RAC(tcel,dB ) = (1.0-lcsmomega)*CSRC_B  + lcsmomega*lcsmeq[ dB ]; \
742          \
743         RAC(tcel,dNE) = (1.0-lcsmomega)*CSRC_NE + lcsmomega*lcsmeq[ dNE]; \
744         RAC(tcel,dNW) = (1.0-lcsmomega)*CSRC_NW + lcsmomega*lcsmeq[ dNW]; \
745         RAC(tcel,dSE) = (1.0-lcsmomega)*CSRC_SE + lcsmomega*lcsmeq[ dSE]; \
746         RAC(tcel,dSW) = (1.0-lcsmomega)*CSRC_SW + lcsmomega*lcsmeq[ dSW]; \
747          \
748         RAC(tcel,dNT) = (1.0-lcsmomega)*CSRC_NT + lcsmomega*lcsmeq[ dNT]; \
749         RAC(tcel,dNB) = (1.0-lcsmomega)*CSRC_NB + lcsmomega*lcsmeq[ dNB]; \
750         RAC(tcel,dST) = (1.0-lcsmomega)*CSRC_ST + lcsmomega*lcsmeq[ dST]; \
751         RAC(tcel,dSB) = (1.0-lcsmomega)*CSRC_SB + lcsmomega*lcsmeq[ dSB]; \
752          \
753         RAC(tcel,dET) = (1.0-lcsmomega)*CSRC_ET + lcsmomega*lcsmeq[ dET]; \
754         RAC(tcel,dEB) = (1.0-lcsmomega)*CSRC_EB + lcsmomega*lcsmeq[ dEB]; \
755         RAC(tcel,dWT) = (1.0-lcsmomega)*CSRC_WT + lcsmomega*lcsmeq[ dWT]; \
756         RAC(tcel,dWB) = (1.0-lcsmomega)*CSRC_WB + lcsmomega*lcsmeq[ dWB]; \
757
758
759
760 #define  OPTIMIZED_STREAMCOLLIDE_NOLES  \
761          \
762         rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T \
763         + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
764         + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
765         ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
766         + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \
767         uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
768         + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \
769         uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
770         + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \
771         PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
772         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
773         RAC(tcel,dC ) = (1.0-OMEGA(lev))*CSRC_C  + OMEGA(lev)*EQC ; \
774         RAC(tcel,dN ) = (1.0-OMEGA(lev))*CSRC_N  + OMEGA(lev)*EQN ; \
775         RAC(tcel,dS ) = (1.0-OMEGA(lev))*CSRC_S  + OMEGA(lev)*EQS ; \
776         RAC(tcel,dE ) = (1.0-OMEGA(lev))*CSRC_E  + OMEGA(lev)*EQE ; \
777         RAC(tcel,dW ) = (1.0-OMEGA(lev))*CSRC_W  + OMEGA(lev)*EQW ; \
778         RAC(tcel,dT ) = (1.0-OMEGA(lev))*CSRC_T  + OMEGA(lev)*EQT ; \
779         RAC(tcel,dB ) = (1.0-OMEGA(lev))*CSRC_B  + OMEGA(lev)*EQB ; \
780          \
781         RAC(tcel,dNE) = (1.0-OMEGA(lev))*CSRC_NE + OMEGA(lev)*EQNE; \
782         RAC(tcel,dNW) = (1.0-OMEGA(lev))*CSRC_NW + OMEGA(lev)*EQNW; \
783         RAC(tcel,dSE) = (1.0-OMEGA(lev))*CSRC_SE + OMEGA(lev)*EQSE; \
784         RAC(tcel,dSW) = (1.0-OMEGA(lev))*CSRC_SW + OMEGA(lev)*EQSW; \
785          \
786         RAC(tcel,dNT) = (1.0-OMEGA(lev))*CSRC_NT + OMEGA(lev)*EQNT; \
787         RAC(tcel,dNB) = (1.0-OMEGA(lev))*CSRC_NB + OMEGA(lev)*EQNB; \
788         RAC(tcel,dST) = (1.0-OMEGA(lev))*CSRC_ST + OMEGA(lev)*EQST; \
789         RAC(tcel,dSB) = (1.0-OMEGA(lev))*CSRC_SB + OMEGA(lev)*EQSB; \
790          \
791         RAC(tcel,dET) = (1.0-OMEGA(lev))*CSRC_ET + OMEGA(lev)*EQET; \
792         RAC(tcel,dEB) = (1.0-OMEGA(lev))*CSRC_EB + OMEGA(lev)*EQEB; \
793         RAC(tcel,dWT) = (1.0-OMEGA(lev))*CSRC_WT + OMEGA(lev)*EQWT; \
794         RAC(tcel,dWB) = (1.0-OMEGA(lev))*CSRC_WB + OMEGA(lev)*EQWB; \
795
796
797
798
799
800 // LES switching for OPT3D
801 #if USE_LES==1
802 #define DEFAULT_COLLIDEG(grav) DEFAULT_COLLIDE_LES(grav)
803 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_LES
804 #else 
805 #define DEFAULT_COLLIDEG(grav) DEFAULT_COLLIDE_NOLES(grav)
806 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_NOLES
807 #endif
808
809 #endif  // 3D, opt OPT3D==true
810
811 #define USQRMAXCHECK(Cusqr,Cux,Cuy,Cuz,  CmMaxVlen,CmMxvx,CmMxvy,CmMxvz) \
812                         if(Cusqr>CmMaxVlen) { \
813                                 CmMxvx = Cux; CmMxvy = Cuy; CmMxvz = Cuz; CmMaxVlen = Cusqr; \
814                         } /* stats */ 
815
816
817
818 /******************************************************************************
819  * interpolateCellFromCoarse macros
820  *****************************************************************************/
821
822
823 // WOXDY_N = Weight Order X Dimension Y _ number N
824 #define WO1D1   ( 1.0/ 2.0)
825 #define WO1D2   ( 1.0/ 4.0)
826 #define WO1D3   ( 1.0/ 8.0)
827
828 #define WO2D1_1 (-1.0/16.0)
829 #define WO2D1_9 ( 9.0/16.0)
830
831 #define WO2D2_11 (WO2D1_1 * WO2D1_1)
832 #define WO2D2_19 (WO2D1_9 * WO2D1_1)
833 #define WO2D2_91 (WO2D1_9 * WO2D1_1)
834 #define WO2D2_99 (WO2D1_9 * WO2D1_9)
835
836 #define WO2D3_111 (WO2D1_1 * WO2D1_1 * WO2D1_1)
837 #define WO2D3_191 (WO2D1_9 * WO2D1_1 * WO2D1_1)
838 #define WO2D3_911 (WO2D1_9 * WO2D1_1 * WO2D1_1)
839 #define WO2D3_991 (WO2D1_9 * WO2D1_9 * WO2D1_1)
840 #define WO2D3_119 (WO2D1_1 * WO2D1_1 * WO2D1_9)
841 #define WO2D3_199 (WO2D1_9 * WO2D1_1 * WO2D1_9)
842 #define WO2D3_919 (WO2D1_9 * WO2D1_1 * WO2D1_9)
843 #define WO2D3_999 (WO2D1_9 * WO2D1_9 * WO2D1_9)
844
845 #if FSGR_STRICT_DEBUG==1
846 #define ADD_INT_DFSCHECK(alev, ai,aj,ak, at, afac, l) \
847                                 if(     (((1.0-(at))>0.0) && (!(QCELL((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr , l) > -1.0 ))) || \
848                                                 (((    (at))>0.0) && (!(QCELL((alev), (ai),(aj),(ak),mLevel[(alev)].setOther, l) > -1.0 ))) ){ \
849                                         errMsg("INVDFSCHECK", " l"<<(alev)<<" "<<PRINT_VEC((ai),(aj),(ak))<<" fc:"<<RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr )<<" fo:"<<RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setOther )<<" dfl"<<l ); \
850                                         debugMarkCell((alev), (ai),(aj),(ak));\
851                                         CAUSE_PANIC; \
852                                 }
853                                 // end ADD_INT_DFSCHECK
854 #define ADD_INT_FLAGCHECK(alev, ai,aj,ak, at, afac) \
855                                 if(     (((1.0-(at))>0.0) && (!(RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr )&(CFInter|CFFluid|CFGrCoarseInited) ))) || \
856                                                 (((    (at))>0.0) && (!(RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setOther)&(CFInter|CFFluid|CFGrCoarseInited) ))) ){ \
857                                         errMsg("INVFLAGCINTCHECK", " l"<<(alev)<<" at:"<<(at)<<" "<<PRINT_VEC((ai),(aj),(ak))<<\
858                                                         " fc:"<<   convertCellFlagType2String(RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr  )) <<\
859                                                         " fold:"<< convertCellFlagType2String(RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setOther )) ); \
860                                         debugMarkCell((alev), (ai),(aj),(ak));\
861                                         CAUSE_PANIC; \
862                                 }
863                                 // end ADD_INT_DFSCHECK
864                                 
865 #define INTUNUTCHECK(ix,iy,iz) \
866                                 if(     (RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) != (CFFluid|CFGrFromCoarse)) ){\
867                                         errMsg("INTFLAGUNU_CHECK", PRINT_VEC(i,j,k)<<" child not unused at l"<<(lev+1)<<" "<<PRINT_VEC((ix),(iy),(iz))<<" flag: "<<  RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) ); \
868                                         debugMarkCell((lev+1), (ix),(iy),(iz));\
869                                         CAUSE_PANIC; \
870                                 }\
871                                 RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) |= CFGrCoarseInited; \
872                                 // INTUNUTCHECK 
873 #define INTSTRICTCHECK(ix,iy,iz,caseId) \
874                                 if(     QCELL(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr, l) <= 0.0 ){\
875                                         errMsg("INVDFCCELLCHECK", "caseId:"<<caseId<<" "<<PRINT_VEC(i,j,k)<<" child inter at "<<PRINT_VEC((ix),(iy),(iz))<<" invalid df "<<l<<" = "<< QCELL(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr, l) ); \
876                                         debugMarkCell((lev+1), (ix),(iy),(iz));\
877                                         CAUSE_PANIC; \
878                                 }\
879                                 // INTSTRICTCHECK
880
881 #else// FSGR_STRICT_DEBUG==1
882 #define ADD_INT_FLAGCHECK(alev, ai,aj,ak, at, afac) 
883 #define ADD_INT_DFSCHECK(alev, ai,aj,ak, at, afac, l) 
884 #define INTSTRICTCHECK(x,y,z,caseId) 
885 #define INTUNUTCHECK(ix,iy,iz) 
886 #endif// FSGR_STRICT_DEBUG==1
887
888
889 #if FSGR_STRICT_DEBUG==1
890 #define INTDEBOUT \
891                 { /*LbmFloat rho,ux,uy,uz;*/ \
892                         rho = ux=uy=uz=0.0; \
893                         FORDF0{ LbmFloat m = QCELL(lev,i,j,k, dstSet, l); \
894                                 rho += m; ux  += (this->dfDvecX[l]*m); uy  += (this->dfDvecY[l]*m); uz  += (this->dfDvecZ[l]*m);  \
895                                 if(ABS(m)>1.0) { errMsg("interpolateCellFromCoarse", "ICFC_DFCHECK cell  "<<PRINT_IJK<<" m"<<l<<":"<< m );CAUSE_PANIC;}\
896                                 /*errMsg("interpolateCellFromCoarse", " cell "<<PRINT_IJK<<" df"<<l<<":"<<m );*/ \
897                         }  \
898                         /*if(this->mPanic) { errMsg("interpolateCellFromCoarse", "ICFC_DFOUT cell  "<<PRINT_IJK<<" rho:"<<rho<<" u:"<<PRINT_VEC(ux,uy,uz)<<" b"<<PRINT_VEC(betx,bety,betz) ); }*/ \
899                         if(markNbs) errMsg("interpolateCellFromCoarse", " cell "<<PRINT_IJK<<" rho:"<<rho<<" u:"<<PRINT_VEC(ux,uy,uz)<<" b"<<PRINT_VEC(betx,bety,betz) );  \
900                         /*errMsg("interpolateCellFromCoarse", "ICFC_DFDEBUG cell  "<<PRINT_IJK<<" rho:"<<rho<<" u:"<<PRINT_VEC(ux,uy,uz)<<" b"<<PRINT_VEC(betx,bety,betz) ); */\
901                 } \
902                 /* both cases are ok to interpolate */  \
903                 if( (!(RFLAG(lev,i,j,k, dstSet) & CFGrFromCoarse)) &&   \
904                                 (!(RFLAG(lev,i,j,k, dstSet) & CFUnused)) ) {    \
905                         /* might also have CFGrCoarseInited (shouldnt be a problem here)*/      \
906                         errMsg("interpolateCellFromCoarse", "CHECK cell not CFGrFromCoarse? "<<PRINT_IJK<<" flag:"<< RFLAG(lev,i,j,k, dstSet)<<" fstr:"<<convertCellFlagType2String(  RFLAG(lev,i,j,k, dstSet) ));      \
907                         /* FIXME check this warning...? return; this can happen !? */   \
908                         /*CAUSE_PANIC;*/        \
909                 }       \
910                 // end INTDEBOUT
911 #else // FSGR_STRICT_DEBUG==1
912 #define INTDEBOUT 
913 #endif // FSGR_STRICT_DEBUG==1
914
915         
916 // t=0.0 -> only current
917 // t=0.5 -> mix
918 // t=1.0 -> only other
919 #if OPT3D==0 
920 #define ADD_INT_DFS(alev, ai,aj,ak, at, afac) \
921                                                 ADD_INT_FLAGCHECK(alev, ai,aj,ak, at, afac); \
922                                                 FORDF0{ \
923                                                         LbmFloat df = ( \
924                                                                         QCELL((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr , l)*(1.0-(at)) + \
925                                                                         QCELL((alev), (ai),(aj),(ak),mLevel[(alev)].setOther, l)*(    (at)) \
926                                                                         ) ; \
927                                                         ADD_INT_DFSCHECK(alev, ai,aj,ak, at, afac, l); \
928                                                         df *= (afac); \
929                                                         rho += df;  \
930                                                         ux  += (this->dfDvecX[l]*df);  \
931                                                         uy  += (this->dfDvecY[l]*df);   \
932                                                         uz  += (this->dfDvecZ[l]*df);   \
933                                                         intDf[l] += df; \
934                                                 } 
935 // write interpolated dfs back to cell (correct non-eq. parts)
936 #define IDF_WRITEBACK_ \
937                 FORDF0{ \
938                         LbmFloat eq = getCollideEq(l, rho,ux,uy,uz);\
939                         QCELL(lev,i,j,k, dstSet, l) = (eq+ (intDf[l]-eq)*mDfScaleDown);\
940                 } \
941                 /* check that all values are ok */ \
942                 INTDEBOUT
943 #define IDF_WRITEBACK \
944                 LbmFloat omegaDst, omegaSrc;\
945                 /* smago new */ \
946                 LbmFloat feq[LBM_DFNUM]; \
947                 LbmFloat dfScale = mDfScaleDown; \
948                 FORDF0{ \
949                         feq[l] = getCollideEq(l, rho,ux,uy,uz); \
950                 } \
951                 if(mLevel[lev  ].lcsmago>0.0) {\
952                         LbmFloat Qo = this->getLesNoneqTensorCoeff(intDf,feq); \
953                         omegaDst  = this->getLesOmega(mLevel[lev+0].omega,mLevel[lev+0].lcsmago,Qo); \
954                         omegaSrc = this->getLesOmega(mLevel[lev-1].omega,mLevel[lev-1].lcsmago,Qo); \
955                 } else {\
956                         omegaDst = mLevel[lev+0].omega; \
957                         omegaSrc = mLevel[lev-1].omega;\
958                 } \
959                  \
960                 dfScale   = (mLevel[lev+0].timestep/mLevel[lev-1].timestep)* (1.0/omegaDst-1.0)/ (1.0/omegaSrc-1.0);  \
961                 FORDF0{ \
962                         /*errMsg("SMAGO"," org"<<mDfScaleDown<<" n"<<dfScale<<" qc"<< QCELL(lev,i,j,k, dstSet, l)<<" idf"<<intDf[l]<<" eq"<<feq[l] ); */ \
963                         QCELL(lev,i,j,k, dstSet, l) = (feq[l]+ (intDf[l]-feq[l])*dfScale);\
964                 } \
965                 /* check that all values are ok */ \
966                 INTDEBOUT
967
968 #else //OPT3D==0 
969
970 #define ADDALLVALS \
971         addVal = addDfFacT * RAC(addfcel , dC ); \
972                                                 intDf[dC ] += addVal; rho += addVal; \
973         addVal  = addDfFacT * RAC(addfcel , dN ); \
974                      uy+=addVal;               intDf[dN ] += addVal; rho += addVal; \
975         addVal  = addDfFacT * RAC(addfcel , dS ); \
976                      uy-=addVal;               intDf[dS ] += addVal; rho += addVal; \
977         addVal  = addDfFacT * RAC(addfcel , dE ); \
978         ux+=addVal;                            intDf[dE ] += addVal; rho += addVal; \
979         addVal  = addDfFacT * RAC(addfcel , dW ); \
980         ux-=addVal;                            intDf[dW ] += addVal; rho += addVal; \
981         addVal  = addDfFacT * RAC(addfcel , dT ); \
982                                   uz+=addVal;  intDf[dT ] += addVal; rho += addVal; \
983         addVal  = addDfFacT * RAC(addfcel , dB ); \
984                                   uz-=addVal;  intDf[dB ] += addVal; rho += addVal; \
985         addVal  = addDfFacT * RAC(addfcel , dNE); \
986         ux+=addVal; uy+=addVal;               intDf[dNE] += addVal; rho += addVal; \
987         addVal  = addDfFacT * RAC(addfcel , dNW); \
988         ux-=addVal; uy+=addVal;               intDf[dNW] += addVal; rho += addVal; \
989         addVal  = addDfFacT * RAC(addfcel , dSE); \
990         ux+=addVal; uy-=addVal;               intDf[dSE] += addVal; rho += addVal; \
991         addVal  = addDfFacT * RAC(addfcel , dSW); \
992         ux-=addVal; uy-=addVal;               intDf[dSW] += addVal; rho += addVal; \
993         addVal  = addDfFacT * RAC(addfcel , dNT); \
994                      uy+=addVal; uz+=addVal;  intDf[dNT] += addVal; rho += addVal; \
995         addVal  = addDfFacT * RAC(addfcel , dNB); \
996                      uy+=addVal; uz-=addVal;  intDf[dNB] += addVal; rho += addVal; \
997         addVal  = addDfFacT * RAC(addfcel , dST); \
998                      uy-=addVal; uz+=addVal;  intDf[dST] += addVal; rho += addVal; \
999         addVal  = addDfFacT * RAC(addfcel , dSB); \
1000                      uy-=addVal; uz-=addVal;  intDf[dSB] += addVal; rho += addVal; \
1001         addVal  = addDfFacT * RAC(addfcel , dET); \
1002         ux+=addVal;              uz+=addVal;  intDf[dET] += addVal; rho += addVal; \
1003         addVal  = addDfFacT * RAC(addfcel , dEB); \
1004         ux+=addVal;              uz-=addVal;  intDf[dEB] += addVal; rho += addVal; \
1005         addVal  = addDfFacT * RAC(addfcel , dWT); \
1006         ux-=addVal;              uz+=addVal;  intDf[dWT] += addVal; rho += addVal; \
1007         addVal  = addDfFacT * RAC(addfcel , dWB); \
1008         ux-=addVal;              uz-=addVal;  intDf[dWB] += addVal; rho += addVal; 
1009
1010 #define ADD_INT_DFS(alev, ai,aj,ak, at, afac) \
1011         addDfFacT = at*afac; \
1012         addfcel = RACPNT((alev), (ai),(aj),(ak),mLevel[(alev)].setOther); \
1013         ADDALLVALS\
1014         addDfFacT = (1.0-at)*afac; \
1015         addfcel = RACPNT((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr); \
1016         ADDALLVALS
1017
1018 // also ugly...
1019 #define INTDF_C    intDf[dC ]
1020 #define INTDF_N    intDf[dN ]
1021 #define INTDF_S    intDf[dS ]
1022 #define INTDF_E    intDf[dE ]
1023 #define INTDF_W    intDf[dW ]
1024 #define INTDF_T    intDf[dT ]
1025 #define INTDF_B    intDf[dB ]
1026 #define INTDF_NE   intDf[dNE]
1027 #define INTDF_NW   intDf[dNW]
1028 #define INTDF_SE   intDf[dSE]
1029 #define INTDF_SW   intDf[dSW]
1030 #define INTDF_NT   intDf[dNT]
1031 #define INTDF_NB   intDf[dNB]
1032 #define INTDF_ST   intDf[dST]
1033 #define INTDF_SB   intDf[dSB]
1034 #define INTDF_ET   intDf[dET]
1035 #define INTDF_EB   intDf[dEB]
1036 #define INTDF_WT   intDf[dWT]
1037 #define INTDF_WB   intDf[dWB]
1038
1039
1040 // write interpolated dfs back to cell (correct non-eq. parts)
1041 #define IDF_WRITEBACK_LES \
1042                 dstcell = RACPNT(lev, i,j,k,dstSet); \
1043                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
1044                 \
1045                 lcsmeq[dC] = EQC ; \
1046                 COLL_CALCULATE_DFEQ(lcsmeq); \
1047                 COLL_CALCULATE_NONEQTENSOR(lev, INTDF_ )\
1048                 COLL_CALCULATE_CSMOMEGAVAL(lev+0, lcsmDstOmega); \
1049                 COLL_CALCULATE_CSMOMEGAVAL(lev-1, lcsmSrcOmega); \
1050                 \
1051                 lcsmdfscale   = (mLevel[lev+0].timestep/mLevel[lev-1].timestep)* (1.0/lcsmDstOmega-1.0)/ (1.0/lcsmSrcOmega-1.0);  \
1052                 RAC(dstcell, dC ) = (lcsmeq[dC ] + (intDf[dC ]-lcsmeq[dC ] )*lcsmdfscale);\
1053                 RAC(dstcell, dN ) = (lcsmeq[dN ] + (intDf[dN ]-lcsmeq[dN ] )*lcsmdfscale);\
1054                 RAC(dstcell, dS ) = (lcsmeq[dS ] + (intDf[dS ]-lcsmeq[dS ] )*lcsmdfscale);\
1055                 RAC(dstcell, dE ) = (lcsmeq[dE ] + (intDf[dE ]-lcsmeq[dE ] )*lcsmdfscale);\
1056                 RAC(dstcell, dW ) = (lcsmeq[dW ] + (intDf[dW ]-lcsmeq[dW ] )*lcsmdfscale);\
1057                 RAC(dstcell, dT ) = (lcsmeq[dT ] + (intDf[dT ]-lcsmeq[dT ] )*lcsmdfscale);\
1058                 RAC(dstcell, dB ) = (lcsmeq[dB ] + (intDf[dB ]-lcsmeq[dB ] )*lcsmdfscale);\
1059                 RAC(dstcell, dNE) = (lcsmeq[dNE] + (intDf[dNE]-lcsmeq[dNE] )*lcsmdfscale);\
1060                 RAC(dstcell, dNW) = (lcsmeq[dNW] + (intDf[dNW]-lcsmeq[dNW] )*lcsmdfscale);\
1061                 RAC(dstcell, dSE) = (lcsmeq[dSE] + (intDf[dSE]-lcsmeq[dSE] )*lcsmdfscale);\
1062                 RAC(dstcell, dSW) = (lcsmeq[dSW] + (intDf[dSW]-lcsmeq[dSW] )*lcsmdfscale);\
1063                 RAC(dstcell, dNT) = (lcsmeq[dNT] + (intDf[dNT]-lcsmeq[dNT] )*lcsmdfscale);\
1064                 RAC(dstcell, dNB) = (lcsmeq[dNB] + (intDf[dNB]-lcsmeq[dNB] )*lcsmdfscale);\
1065                 RAC(dstcell, dST) = (lcsmeq[dST] + (intDf[dST]-lcsmeq[dST] )*lcsmdfscale);\
1066                 RAC(dstcell, dSB) = (lcsmeq[dSB] + (intDf[dSB]-lcsmeq[dSB] )*lcsmdfscale);\
1067                 RAC(dstcell, dET) = (lcsmeq[dET] + (intDf[dET]-lcsmeq[dET] )*lcsmdfscale);\
1068                 RAC(dstcell, dEB) = (lcsmeq[dEB] + (intDf[dEB]-lcsmeq[dEB] )*lcsmdfscale);\
1069                 RAC(dstcell, dWT) = (lcsmeq[dWT] + (intDf[dWT]-lcsmeq[dWT] )*lcsmdfscale);\
1070                 RAC(dstcell, dWB) = (lcsmeq[dWB] + (intDf[dWB]-lcsmeq[dWB] )*lcsmdfscale);\
1071                 /* IDF_WRITEBACK optimized */
1072
1073 #define IDF_WRITEBACK_NOLES \
1074                 dstcell = RACPNT(lev, i,j,k,dstSet); \
1075                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
1076                 \
1077                 RAC(dstcell, dC ) = (EQC  + (intDf[dC ]-EQC  )*mDfScaleDown);\
1078                 RAC(dstcell, dN ) = (EQN  + (intDf[dN ]-EQN  )*mDfScaleDown);\
1079                 RAC(dstcell, dS ) = (EQS  + (intDf[dS ]-EQS  )*mDfScaleDown);\
1080                 /*old*/ RAC(dstcell, dE ) = (EQE  + (intDf[dE ]-EQE  )*mDfScaleDown);\
1081                 RAC(dstcell, dW ) = (EQW  + (intDf[dW ]-EQW  )*mDfScaleDown);\
1082                 RAC(dstcell, dT ) = (EQT  + (intDf[dT ]-EQT  )*mDfScaleDown);\
1083                 RAC(dstcell, dB ) = (EQB  + (intDf[dB ]-EQB  )*mDfScaleDown);\
1084                 /*old*/ RAC(dstcell, dNE) = (EQNE + (intDf[dNE]-EQNE )*mDfScaleDown);\
1085                 RAC(dstcell, dNW) = (EQNW + (intDf[dNW]-EQNW )*mDfScaleDown);\
1086                 RAC(dstcell, dSE) = (EQSE + (intDf[dSE]-EQSE )*mDfScaleDown);\
1087                 RAC(dstcell, dSW) = (EQSW + (intDf[dSW]-EQSW )*mDfScaleDown);\
1088                 RAC(dstcell, dNT) = (EQNT + (intDf[dNT]-EQNT )*mDfScaleDown);\
1089                 RAC(dstcell, dNB) = (EQNB + (intDf[dNB]-EQNB )*mDfScaleDown);\
1090                 RAC(dstcell, dST) = (EQST + (intDf[dST]-EQST )*mDfScaleDown);\
1091                 RAC(dstcell, dSB) = (EQSB + (intDf[dSB]-EQSB )*mDfScaleDown);\
1092                 RAC(dstcell, dET) = (EQET + (intDf[dET]-EQET )*mDfScaleDown);\
1093                 /*old*/ RAC(dstcell, dEB) = (EQEB + (intDf[dEB]-EQEB )*mDfScaleDown);\
1094                 RAC(dstcell, dWT) = (EQWT + (intDf[dWT]-EQWT )*mDfScaleDown);\
1095                 RAC(dstcell, dWB) = (EQWB + (intDf[dWB]-EQWB )*mDfScaleDown);\
1096                 /* IDF_WRITEBACK optimized */
1097
1098 #if USE_LES==1
1099 #define IDF_WRITEBACK IDF_WRITEBACK_LES
1100 #else 
1101 #define IDF_WRITEBACK IDF_WRITEBACK_NOLES
1102 #endif
1103
1104 #endif// OPT3D==0 
1105
1106
1107
1108 /******************************************************************************/
1109 /*! relaxation LES functions */
1110 /******************************************************************************/
1111
1112
1113 inline LbmFloat LbmFsgrSolver::getLesNoneqTensorCoeff(
1114                 LbmFloat df[],                          
1115                 LbmFloat feq[] ) {
1116         LbmFloat Qo = 0.0;
1117         for(int m=0; m< ((LBMDIM*LBMDIM)-LBMDIM)/2 ; m++) { 
1118                 LbmFloat qadd = 0.0;
1119                 for(int l=1; l<this->cDfNum; l++) { 
1120                         if(this->lesCoeffOffdiag[m][l]==0.0) continue;
1121                         qadd += this->lesCoeffOffdiag[m][l]*(df[l]-feq[l]);
1122                 }
1123                 Qo += (qadd*qadd);
1124         }
1125         Qo *= 2.0; // off diag twice
1126         for(int m=0; m<LBMDIM; m++) { 
1127                 LbmFloat qadd = 0.0;
1128                 for(int l=1; l<this->cDfNum; l++) { 
1129                         if(this->lesCoeffDiag[m][l]==0.0) continue;
1130                         qadd += this->lesCoeffDiag[m][l]*(df[l]-feq[l]);
1131                 }
1132                 Qo += (qadd*qadd);
1133         }
1134         Qo = sqrt(Qo);
1135         return Qo;
1136 };
1137
1138 inline LbmFloat LbmFsgrSolver::getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo) {
1139         const LbmFloat tau = 1.0/omega;
1140         const LbmFloat nu = (2.0*tau-1.0) * (1.0/6.0);
1141         const LbmFloat C = csmago;
1142         const LbmFloat Csqr = C*C;
1143         LbmFloat S = -nu + sqrt( nu*nu + 18.0*Csqr*Qo ) / (6.0*Csqr);
1144         return( 1.0/( 3.0*( nu+Csqr*S ) +0.5 ) );
1145 }
1146
1147 #define DEBUG_CALCPRINTCELL(str,df) {\
1148                 LbmFloat prho=df[0], pux=0., puy=0., puz=0.; \
1149                 for(int dfl=1; dfl<this->cDfNum; dfl++) { \
1150                         prho += df[dfl];  \
1151                         pux  += (this->dfDvecX[dfl]*df[dfl]);  \
1152                         puy  += (this->dfDvecY[dfl]*df[dfl]);  \
1153                         puz  += (this->dfDvecZ[dfl]*df[dfl]);  \
1154                 } \
1155                 errMsg("DEBUG_CALCPRINTCELL",">"<<str<<" rho="<<prho<<" vel="<<ntlVec3Gfx(pux,puy,puz) ); \
1156         } /* END DEBUG_CALCPRINTCELL */ 
1157
1158 // "normal" collision
1159 inline void LbmFsgrSolver::collideArrays(
1160                 int lev, int i, int j, int k, // position - more for debugging
1161                 LbmFloat df[],                          
1162                 LbmFloat &outrho, // out only!
1163                 // velocity modifiers (returns actual velocity!)
1164                 LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
1165                 LbmFloat omega, 
1166                 LbmVec gravity,
1167                 LbmFloat csmago, 
1168                 LbmFloat *newOmegaRet, LbmFloat *newQoRet
1169         ) {
1170         int l;
1171         LbmFloat rho=df[0]; 
1172         LbmFloat ux = 0; //mux;
1173         LbmFloat uy = 0; //muy;
1174         LbmFloat uz = 0; //muz; 
1175         LbmFloat feq[19];
1176         LbmFloat omegaNew;
1177         LbmFloat Qo = 0.0;
1178
1179         for(l=1; l<this->cDfNum; l++) { 
1180                 rho += df[l]; 
1181                 ux  += (this->dfDvecX[l]*df[l]); 
1182                 uy  += (this->dfDvecY[l]*df[l]);  
1183                 uz  += (this->dfDvecZ[l]*df[l]);  
1184         }  
1185
1186
1187         PRECOLLIDE_MODS(rho,ux,uy,uz, gravity);
1188         for(l=0; l<this->cDfNum; l++) { 
1189                 feq[l] = getCollideEq(l,rho,ux,uy,uz); 
1190         }
1191
1192         if(csmago>0.0) {
1193                 Qo = getLesNoneqTensorCoeff(df,feq);
1194                 omegaNew = getLesOmega(omega,csmago,Qo);
1195         } else {
1196                 omegaNew = omega; // smago off...
1197         }
1198         if(newOmegaRet) *newOmegaRet = omegaNew; // return value for stats
1199         if(newQoRet)    *newQoRet = Qo; // return value of non-eq. stress tensor
1200
1201         for(l=0; l<this->cDfNum; l++) { 
1202                 df[l] = (1.0-omegaNew ) * df[l] + omegaNew * feq[l]; 
1203         }  
1204         //if((i==16)&&(j==10)) DEBUG_CALCPRINTCELL( "2dcoll "<<PRINT_IJK, df);
1205
1206         mux = ux;
1207         muy = uy;
1208         muz = uz;
1209         outrho = rho;
1210
1211         lev=i=j=k; // debug, remove warnings
1212 };
1213