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