momentTemplates.h

00001 /* This file is part of the Palabos library.
00002  *
00003  * Copyright (C) 2011 FlowKit Sarl
00004  * Avenue de Chailly 23
00005  * 1012 Lausanne, Switzerland
00006  * E-mail contact: contact@flowkit.com
00007  *
00008  * The most recent release of Palabos can be downloaded at 
00009  * <http://www.palabos.org/>
00010  *
00011  * The library Palabos is free software: you can redistribute it and/or
00012  * modify it under the terms of the GNU Affero General Public License as
00013  * published by the Free Software Foundation, either version 3 of the
00014  * License, or (at your option) any later version.
00015  *
00016  * The library is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU Affero General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU Affero General Public License
00022  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00023 */
00024 
00031 #ifndef MOMENT_TEMPLATES_H
00032 #define MOMENT_TEMPLATES_H
00033 
00034 #include "core/globalDefs.h"
00035 #include "core/cell.h"
00036 #include "core/util.h"
00037 #include "latticeBoltzmann/geometricOperationTemplates.h"
00038 #include "latticeBoltzmann/roundOffPolicy.h"
00039 
00040 namespace plb {
00041 
00042 template<typename T, class Descriptor> struct momentTemplatesImpl;
00043 
00044 // This structure forwards the calls to the appropriate helper class
00045 template<typename T, template<typename U> class Descriptor>
00046 struct momentTemplates {
00047 
00048 static T get_rhoBar(Cell<T,Descriptor> const& cell) {
00049     return momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00050         ::get_rhoBar(cell.getRawPopulations());
00051 }
00052 
00053 static void get_j(Cell<T,Descriptor> const& cell, Array<T,Descriptor<T>::d>& j ) {
00054     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00055         ::get_j(cell.getRawPopulations(), j);
00056 }
00057 
00058 static T get_eBar(Cell<T,Descriptor> const& cell) {
00059     return momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00060         ::get_eBar(cell.getRawPopulations());
00061 }
00062 
00063 static void get_rhoBar_j(Cell<T,Descriptor> const& cell, T& rhoBar, Array<T,Descriptor<T>::d>& j ) {
00064     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00065         ::get_rhoBar_j(cell.getRawPopulations(), rhoBar, j);
00066 }
00067 
00068 static void get_rhoBar_j_thetaBar(Cell<T,Descriptor> const& cell, T& rhoBar, Array<T,Descriptor<T>::d>& j, T &thetaBar ) {
00069     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00070     ::get_rhoBar_j_thetaBar(cell.getRawPopulations(), rhoBar, j, thetaBar);
00071 }
00072 
00073 static T compute_rho(Cell<T,Descriptor> const& cell) {
00074     return momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00075         ::compute_rho(cell.getRawPopulations());
00076 }
00077 
00079 
00083 static void compute_uLb(Cell<T,Descriptor> const& cell, Array<T,Descriptor<T>::d>& uLb ) {
00084     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00085         ::compute_uLb(cell.getRawPopulations(), uLb);
00086 }
00087 
00088 static void compute_rho_uLb(Cell<T,Descriptor> const& cell, T& rho, Array<T,Descriptor<T>::d>& uLb ) {
00089     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00090         ::compute_rho_uLb(cell.getRawPopulations(), rho, uLb);
00091 }
00092 
00093 static T compute_e(Cell<T,Descriptor> const& cell) {
00094     return momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00095         ::compute_e(cell.getRawPopulations());
00096 }
00097 
00098 static T compute_rhoThetaBar(Cell<T,Descriptor> const& cell, T rhoBar, T jSqr) {
00099     return momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00100         ::compute_rhoThetaBar(cell.getRawPopulations(), rhoBar, jSqr);
00101 }
00102 
00103 static void compute_rho_rhoThetaBar(Cell<T,Descriptor> const& cell, T& rho, T& rhoThetaBar) {
00104     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00105         ::compute_rho_rhoThetaBar(cell.getRawPopulations(), rho, rhoThetaBar);
00106 }
00107 
00108 static T compute_theta(Cell<T,Descriptor> const& cell, T rhoBar, T jSqr) {
00109     return momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00110         ::compute_theta(cell.getRawPopulations(), rhoBar, jSqr);
00111 }
00112 
00113 static T compute_rhoEpsilon(Cell<T,Descriptor> const& cell, T rhoBar, T jSqr) {
00114     return momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00115         ::compute_rhoEpsilon(cell.getRawPopulations(), rhoBar, jSqr);
00116 }
00117 
00118 static void compute_PiNeq(Cell<T,Descriptor> const& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j,
00119                           Array<T,SymmetricTensor<T,Descriptor>::n>& PiNeq, bool incompr=false )
00120 {
00121     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00122         ::compute_PiNeq(cell.getRawPopulations(), rhoBar, j, PiNeq, incompr);
00123 }
00124 
00125 static void compute_thermal_PiNeq(Cell<T,Descriptor> const& cell, T rhoBar, T thetaBar,
00126                                   Array<T,Descriptor<T>::d> const& j,
00127                                   Array<T,SymmetricTensor<T,Descriptor>::n>& PiNeq )
00128 {
00129     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00130         ::compute_thermal_PiNeq(cell.getRawPopulations(), rhoBar, thetaBar, j, PiNeq);
00131 }
00132 
00133 static void compute_rhoBar_j_PiNeq(Cell<T,Descriptor> const& cell, T& rhoBar, Array<T,Descriptor<T>::d>& j,
00134                                    Array<T,SymmetricTensor<T,Descriptor>::n>& PiNeq, bool incompr=false )
00135 {
00136     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00137         ::compute_rhoBar_j_PiNeq(cell.getRawPopulations(), rhoBar, j, PiNeq, incompr);
00138 }
00139 
00140 static void compute_rhoBar_thetaBar_j_PiNeq(Cell<T,Descriptor> const& cell, T& rhoBar, T& thetaBar,
00141                                             Array<T,Descriptor<T>::d> & j,
00142                                             Array<T,SymmetricTensor<T,Descriptor>::n>& PiNeq )
00143 {
00144     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00145         ::compute_rhoBar_thetaBar_j_PiNeq(cell.getRawPopulations(), rhoBar, thetaBar, j, PiNeq);
00146 }
00147                                                  
00148 static void compute_rhoBar_thetaBar_j_PiNeq_qNeq(Cell<T,Descriptor> const& cell, T& rhoBar, T& thetaBar,
00149                                                  Array<T,Descriptor<T>::d> & j,
00150                                                  Array<T,SymmetricTensor<T,Descriptor>::n>& PiNeq,
00151                                                  Array<T,SymmetricRankThreeTensor<T,Descriptor>::n>& qNeq)
00152 {
00153     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00154         ::compute_rhoBar_thetaBar_j_PiNeq_qNeq(cell.getRawPopulations(), rhoBar, thetaBar, j, PiNeq, qNeq);
00155 }
00156 
00157 
00159 
00162 static void compute_P(Cell<T,Descriptor> const& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j,
00163                       Array<T,SymmetricTensor<T,Descriptor>::n>& P)
00164 {
00165     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00166         ::compute_P(cell.getRawPopulations(), rhoBar, j, P);
00167 }
00168 
00169 static void modifyJ(T& cell, Array<T,Descriptor<T>::d> const& newJ) {
00170     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00171         ::modifyVelocity(cell.getRawPopulations(), newJ);
00172 }
00173 
00174 static void compute_Qneq(Cell<T,Descriptor> const& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j,
00175                          T thetaBar,
00176                          Array<T,SymmetricRankThreeTensor<T,Descriptor>::n>& qNeq )
00177 {
00178     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00179         ::compute_Qneq(cell.getRawPopulations(), rhoBar, j, thetaBar, qNeq);
00180 }
00181 
00182 static void compute_heat_flux(Cell<T,Descriptor> const& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j,
00183                               T thetaBar,
00184                               Array<T,Descriptor<T>::d>& q )
00185 {
00186     momentTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
00187         ::compute_heat_flux(cell.getRawPopulations(), rhoBar, j, thetaBar, q);
00188 }
00189 
00190 };  // struct momentTemplates
00191 
00192 
00193 // This structure forwards the calls to the appropriate helper class
00194 template<typename T, class Descriptor>
00195 struct momentTemplatesImpl {
00196 
00197 static T get_rhoBar(Array<T,Descriptor::q> const& f) {
00198     T rhoBar = f[0];
00199     for (plint iPop=1; iPop < Descriptor::q; ++iPop) {
00200         rhoBar += f[iPop];
00201     }
00202     return rhoBar;
00203 }
00204 
00205 static void get_j(Array<T,Descriptor::q> const& f, Array<T,Descriptor::d>& j ) {
00206     for (int iD=0; iD < Descriptor::d; ++iD) {
00207         j[iD] = f[0]*Descriptor::c[0][iD];
00208     }
00209     for (plint iPop=1; iPop < Descriptor::q; ++iPop) {
00210         for (int iD=0; iD < Descriptor::d; ++iD) {
00211             j[iD] += f[iPop]*Descriptor::c[iPop][iD];
00212         }
00213     }
00214 }
00215 
00216 static T get_eBar(Array<T,Descriptor::q> const& f) {
00217     T eBar = f[0] * Descriptor::cNormSqr[0];
00218     for (plint iPop=1; iPop < Descriptor::q; ++iPop) {
00219         eBar += f[iPop] * Descriptor::cNormSqr[iPop];
00220     }
00221     return eBar;
00222 }
00223 
00224 static void get_rhoBar_j(Array<T,Descriptor::q> const& f, T& rhoBar, Array<T,Descriptor::d>& j ) {
00225     rhoBar = get_rhoBar(f);
00226     get_j(f, j);
00227 }
00228 
00229 static void get_rhoBar_j_thetaBar(Array<T,Descriptor::q> const& f, T& rhoBar, Array<T,Descriptor::d>& j, T &thetaBar ) {
00230     get_rhoBar_j(f,rhoBar,j);
00231     T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
00232     T invRho = Descriptor::invRho(rhoBar);
00233     thetaBar = compute_rhoThetaBar(f, rhoBar, jSqr) * invRho;
00234 }
00235 
00236 
00237 static T compute_rho(Array<T,Descriptor::q> const& f) {
00238     return Descriptor::fullRho(get_rhoBar(f));
00239 }
00240 
00241 static void compute_uLb(Array<T,Descriptor::q> const& f, Array<T,Descriptor::d>& uLb ) {
00242     get_j(f, uLb);
00243     T invRho = Descriptor::invRho(get_rhoBar(f));
00244     for (int iD=0; iD < Descriptor::d; ++iD) {
00245         uLb[iD] *= invRho;
00246     }
00247 }
00248 
00249 static void compute_rho_uLb(Array<T,Descriptor::q> const& f, T& rho, Array<T,Descriptor::d>& uLb ) {
00250     get_j(f, uLb);
00251     T rhoBar = get_rhoBar(f);
00252     T invRho = Descriptor::invRho(rhoBar);
00253     rho = Descriptor::fullRho(rhoBar);
00254     for (int iD=0; iD < Descriptor::d; ++iD) {
00255         uLb[iD] *= invRho;
00256     }
00257 }
00258 
00259 static T compute_e(Array<T,Descriptor::q> const& f) {
00260     return get_eBar(f) + Descriptor::SkordosFactor() * Descriptor::d * Descriptor::cs2;
00261 }
00262 
00263 static T compute_rhoThetaBar(Array<T,Descriptor::q> const& f, T rhoBar, T jSqr) {
00264     T invRho = Descriptor::invRho(rhoBar);
00265     return Descriptor::invCs2 * Descriptor::invD * (get_eBar(f) - invRho*jSqr) - rhoBar;
00266 }
00267 
00268 static void compute_rho_rhoThetaBar(Array<T,Descriptor::q> const& f, T& rho, T& rhoThetaBar) {
00269     T rhoBar, j[Descriptor::d];
00270     get_rhoBar_j(f, rhoBar, j);
00271     T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
00272     rho = Descriptor::fullRho(rhoBar);
00273     rhoThetaBar = compute_rhoThetaBar(f, rhoBar, jSqr);
00274 }
00275 
00276 static T compute_theta(Array<T,Descriptor::q> const& f, T rhoBar, T jSqr) {
00277     T invRho = Descriptor::invRho(rhoBar);
00278     T e = compute_e(f);
00279     return invRho * Descriptor::invD * Descriptor::invCs2 * (e - invRho*jSqr);
00280 }
00281 
00282 static T compute_rhoEpsilon(Array<T,Descriptor::q> const& f, T rhoBar, T jSqr) {
00283     T invRho = Descriptor::invRho(rhoBar);
00284     T e = compute_e(f);
00285     return (e - invRho*jSqr) / (T)2;
00286 }
00287 
00288 static void compute_PiNeq(Array<T,Descriptor::q> const& f, T rhoBar, Array<T,Descriptor::d> const& j,
00289                           Array<T,SymmetricTensorImpl<T,Descriptor::d>::n>& PiNeq, bool incompr )
00290 {
00291     T invRho = incompr ? 1. : Descriptor::invRho(rhoBar);
00292     int iPi = 0;
00293     for (int iAlpha=0; iAlpha < Descriptor::d; ++iAlpha) {
00294         int iDiagonal = iPi;
00295         for (int iBeta=iAlpha; iBeta < Descriptor::d; ++iBeta) {
00296             PiNeq[iPi] = Descriptor::c[0][iAlpha]*
00297                          Descriptor::c[0][iBeta] * f[0];
00298             for (plint iPop=1; iPop < Descriptor::q; ++iPop) {
00299                 PiNeq[iPi] += Descriptor::c[iPop][iAlpha]*
00300                               Descriptor::c[iPop][iBeta] * f[iPop];
00301             }
00302             // Stripe off relative velocity
00303             PiNeq[iPi] -= invRho*j[iAlpha]*j[iBeta];
00304             ++iPi;
00305         }
00306         // Stripe off diagonal term
00307         PiNeq[iDiagonal] -= Descriptor::cs2 * rhoBar;
00308     }
00309 }
00310 
00311 static void compute_thermal_PiNeq(Array<T,Descriptor::q> const& f, T rhoBar, T thetaBar,
00312                                   Array<T,Descriptor::d> const& j,
00313                                   Array<T,SymmetricTensorImpl<T,Descriptor::d>::n>& PiNeq )
00314 {
00315     // rhoTheta_bar == rho*theta - 1 
00316     T rhoTheta_bar = rhoBar * thetaBar + rhoBar + Descriptor::SkordosFactor() * thetaBar;
00317     
00318     T invRho = Descriptor::invRho(rhoBar);
00319     int iPi = 0;
00320     for (int iAlpha=0; iAlpha < Descriptor::d; ++iAlpha) {
00321         int iDiagonal = iPi;
00322         for (int iBeta=iAlpha; iBeta < Descriptor::d; ++iBeta) {
00323             PiNeq[iPi] = Descriptor::c[0][iAlpha]*
00324                          Descriptor::c[0][iBeta] * f[0];
00325             for (plint iPop=1; iPop < Descriptor::q; ++iPop) {
00326                 PiNeq[iPi] += Descriptor::c[iPop][iAlpha]*
00327                               Descriptor::c[iPop][iBeta] * f[iPop];
00328             }
00329             // Stripe off relative velocity
00330             PiNeq[iPi] -= invRho*j[iAlpha]*j[iBeta];
00331             ++iPi;
00332         }
00333         // Stripe off diagonal term
00334         PiNeq[iDiagonal] -= Descriptor::cs2 * rhoTheta_bar;
00335     }
00336 }
00337 
00338 static void compute_rhoBar_j_PiNeq(Array<T,Descriptor::q> const& f, T& rhoBar, Array<T,Descriptor::d>& j,
00339                                    Array<T,SymmetricTensorImpl<T,Descriptor::d>::n>& PiNeq, bool incompr )
00340 {
00341     get_rhoBar_j(f, rhoBar, j);
00342     compute_PiNeq(f, rhoBar, j, PiNeq, incompr);
00343 }
00344 
00345 static void compute_rhoBar_thetaBar_j_PiNeq(Array<T,Descriptor::q> const& f, T& rhoBar, T& thetaBar,
00346                                             Array<T,Descriptor::d> & j,
00347                                             Array<T,SymmetricTensorImpl<T,Descriptor::d>::n>& PiNeq )
00348 {
00349     get_rhoBar_j_thetaBar(f,rhoBar,j, thetaBar);
00350     compute_thermal_PiNeq(f, rhoBar, thetaBar, j, PiNeq);
00351 }
00352 
00353 static void compute_rhoBar_thetaBar_j_PiNeq_qNeq(Array<T,Descriptor::q> const& f, T& rhoBar, T& thetaBar,
00354                                                  Array<T,Descriptor::d> & j,
00355                                                  Array<T,SymmetricTensorImpl<T,Descriptor::d>::n>& PiNeq,
00356                                                  Array<T,SymmetricRankThreeTensorImpl<T,Descriptor::d>::n>& qNeq)
00357 {
00358     compute_rhoBar_thetaBar_j_PiNeq(f, rhoBar, thetaBar, j, PiNeq);
00359     compute_Qneq(f, rhoBar, j, thetaBar, qNeq );
00360 }
00361 
00362 
00363 static void compute_P(Array<T,Descriptor::q> const& f, T rhoBar, Array<T,Descriptor::d> const& j,
00364                       Array<T,SymmetricTensorImpl<T,Descriptor::d>::n>& P )
00365 {
00366     T invRho = Descriptor::invRho(rhoBar);
00367     plint iP = 0;
00368     for (int iAlpha=0; iAlpha < Descriptor::d; ++iAlpha) {
00369         for (int iBeta=iAlpha; iBeta < Descriptor::d; ++iBeta) {
00370             P[iP] = Descriptor::c[0][iAlpha]*
00371                     Descriptor::c[0][iBeta] * f[0];
00372             for (plint iPop=1; iPop < Descriptor::q; ++iPop) {
00373                 P[iP] += Descriptor::c[iPop][iAlpha]*
00374                          Descriptor::c[iPop][iBeta] * f[iPop];
00375             }
00376             // Stripe off relative velocity
00377             P[iP] -= invRho*j[iAlpha]*j[iBeta];
00378             ++iP;
00379         }
00380     }
00381 }
00382 
00383 static void modifyJ(T* f, Array<T,Descriptor::d> const& newJ) {
00384     T rhoBar;
00385     Array<T,Descriptor::d> oldJ;
00386     get_rhoBar_j(f, rhoBar, oldJ);
00387     T invRho = Descriptor::invRho(rhoBar);
00388     const T oldJSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(oldJ);
00389     const T newJSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(newJ);
00390     for (plint iPop=0; iPop<Descriptor::q; ++iPop) {
00391         f[iPop] = f[iPop]
00392                          - bgk_ma2_equilibrium(iPop, rhoBar, invRho, oldJ, oldJSqr)
00393                          + bgk_ma2_equilibrium(iPop, rhoBar, invRho, newJ, newJSqr);
00394     }
00395 }
00396 
00397 static void compute_Qneq(Array<T,Descriptor::q> const& f, T rhoBar, Array<T,Descriptor::d> const& j,
00398                          T thetaBar,
00399                          Array<T,SymmetricRankThreeTensorImpl<T,Descriptor::d>::n>& qNeq )
00400 {
00401     typedef Descriptor L;
00402     T invRho = L::invRho(rhoBar);
00403     T factor = L::cs2+thetaBar;
00404     plint iQ = 0;
00405     for (plint iA = 0; iA < L::d; ++iA) {
00406         for (plint iB = iA; iB < L::d; ++iB) {
00407             for (plint iC = iB; iC < L::d; ++iC) {
00408                 qNeq[iQ] = - j[iA]*j[iB]*j[iC]*invRho*invRho;
00409                 for (plint iPop = 0; iPop < L::q; ++iPop) {
00410                     qNeq[iQ] = L::c[iPop][iA]*L::c[iPop][iB]*L::c[iPop][iC]*f[iPop];
00411                 }
00412                 if (iA == iB && iB == iC) {
00413                     qNeq[iQ] -= (T)3 * factor * j[iA];
00414                 }
00415                 else if (iA == iB && iB != iC) {
00416                     qNeq[iQ] -= factor * j[iC];
00417                 }
00418                 else if (iA == iC && iC != iB) {
00419                     qNeq[iQ] -= factor * j[iB];
00420                 }
00421                 else if (iB == iC && iC != iA) {
00422                     qNeq[iQ] -= factor * j[iA];
00423                 }
00424                 
00425                 ++iQ;
00426             }
00427         }
00428     }
00429     
00430 }
00431 
00432 static void compute_heat_flux(Array<T,Descriptor::q> const& f, T rhoBar, Array<T,Descriptor::d> const& j,
00433                               T thetaBar,
00434                               Array<T,Descriptor::d>& q)
00435 {
00436     
00437     typedef Descriptor L;
00438     Array<T,SymmetricRankThreeTensorImpl<T,Descriptor::d>::n> qNeq;
00439     compute_Qneq(f, rhoBar, j,thetaBar,qNeq );
00440     SymmetricRankThreeTensorImpl<T,L::d>::contractLastTwoIndexes(qNeq,q);
00441 }
00442 
00443 };  // struct momentTemplatesImpl
00444 
00445 }  // namespace plb
00446 
00447 #include "latticeBoltzmann/momentTemplates2D.h"
00448 #include "latticeBoltzmann/momentTemplates3D.h"
00449 
00450 #endif  // MOMENT_TEMPLATES_H