mrtTemplates2D.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 
00025 /* Orestis Malaspinas contributed this code.
00026  */
00027 
00033 #ifndef MRT_TEMPLATES_2D_H
00034 #define MRT_TEMPLATES_2D_H
00035 
00036 #include "core/globalDefs.h"
00037 
00038 namespace plb {
00039 
00040 // Efficient specialization for D2Q9 lattice
00041 template<typename T> struct mrtTemplatesImpl<T, descriptors::MRTD2Q9DescriptorBase<T> > {
00042 
00043     typedef descriptors::D2Q9DescriptorBase<T> Descriptor; 
00044 
00046     static void computeEquilibrium( Array<T,Descriptor::q>& momentsEq,
00047                                     T rhoBar, Array<T,2> const& j, T jSqr )
00048     {
00049         T invRho = Descriptor::invRho(rhoBar);
00050         momentsEq[0] = rhoBar;
00051         momentsEq[1] = (T)3*jSqr*invRho-2*rhoBar;
00052         momentsEq[2] = -(T)3*jSqr*invRho+rhoBar;
00053         momentsEq[3] = j[0];
00054         momentsEq[4] = -j[0];
00055         momentsEq[5] = j[1];
00056         momentsEq[6] = -j[1];
00057         momentsEq[7] = (j[0]*j[0]-j[1]*j[1])*invRho;
00058         momentsEq[8] = j[1]*j[0]*invRho;
00059     }
00060     
00062     static void computeMoments(Array<T,Descriptor::q>& moments, const Array<T,Descriptor::q>& f)
00063     {
00064         T f1_f3 = f[1]-f[3];
00065         T f1pf3 = f[1]+f[3];
00066         
00067         T f2pf4 = f[2]+f[4];
00068         T f5_f7 = f[5]-f[7];
00069         T f5pf7 = f[5]+f[7];
00070         T f6pf8 = f[6]+f[8];
00071         
00072         moments[0] = f[0] + f1pf3 + f2pf4 + f5pf7 + f6pf8;
00073         moments[1] = -(T)4*f[0] + (T)2*(f1pf3 + f5pf7) - f2pf4 - f6pf8;
00074         moments[2] =  (T)4*f[0] + f1pf3 - (T)2*(f2pf4 + f6pf8) + f5pf7;
00075         moments[3] = -f1pf3 - f[2] + f5pf7 + f[6];
00076         moments[4] = -f1pf3 + (T)2*(f[2] - f[6]) + f5pf7;
00077         moments[5] =  f1_f3 - f[4] - f5_f7 + f[8];
00078         moments[6] =  f1_f3 + (T)2*(f[4] - f[8]) - f5_f7;
00079         moments[7] =  f[2] - f[4] + f[6] - f[8];
00080         moments[8] = -f1_f3 - f5_f7;
00081     }
00082     
00084     static T mrtCollision( Array<T,Descriptor::q>& f,
00085                            const T &rhoBar, const Array<T,2> & j,
00086                            T invM_S[9][9] )
00087     {
00088         Array<T,9> moments, momentsEq;
00089 
00090         computeMoments(moments,f);
00091         T jSqr = VectorTemplateImpl<T,2>::normSqr(j);
00092         
00093         computeEquilibrium(momentsEq,rhoBar,j,jSqr);
00094         
00095         T mom1 = moments[1] - momentsEq[1];
00096         T mom2 = moments[2] - momentsEq[2];
00097         T mom3 = moments[3] - momentsEq[3];
00098         T mom4 = moments[4] - momentsEq[4];
00099         T mom5 = moments[5] - momentsEq[5];
00100         T mom6 = moments[6] - momentsEq[6];
00101         T mom7 = moments[7] - momentsEq[7];
00102         T mom8 = moments[8] - momentsEq[8];
00103         
00104         f[0] -= invM_S[0][1]*mom1 + 
00105                    invM_S[0][2]*mom2;
00106         
00107         f[1] -= invM_S[1][1]*mom1 +
00108                    invM_S[1][2]*mom2 +
00109                    invM_S[1][3]*mom3 +
00110                    invM_S[1][4]*mom4 + 
00111                    invM_S[1][5]*mom5 + 
00112                    invM_S[1][6]*mom6 +
00113                    invM_S[1][8]*mom8;
00114         
00115         f[2] -= invM_S[2][1]*mom1 +
00116                    invM_S[2][2]*mom2 +
00117                    invM_S[2][3]*mom3 +
00118                    invM_S[2][4]*mom4 +
00119                    invM_S[2][7]*mom7;
00120         
00121         f[3] -= invM_S[3][1]*mom1 +
00122                    invM_S[3][2]*mom2 +
00123                    invM_S[3][3]*mom3 +
00124                    invM_S[3][4]*mom4 +
00125                    invM_S[3][5]*mom5 +
00126                    invM_S[3][6]*mom6 +
00127                    invM_S[3][8]*mom8;
00128         
00129         f[4] -= invM_S[4][1]*mom1 +
00130                    invM_S[4][2]*mom2 +
00131                    invM_S[4][5]*mom5 +
00132                    invM_S[4][6]*mom6 +
00133                    invM_S[4][7]*mom7;
00134         
00135         f[5] -= invM_S[5][1]*mom1 +
00136                    invM_S[5][2]*mom2 +
00137                    invM_S[5][3]*mom3 +
00138                    invM_S[5][4]*mom4 +
00139                    invM_S[5][5]*mom5 +
00140                    invM_S[5][6]*mom6 +
00141                    invM_S[5][8]*mom8;
00142         
00143         f[6] -= invM_S[6][1]*mom1 +
00144                    invM_S[6][2]*mom2 +
00145                    invM_S[6][3]*mom3 +
00146                    invM_S[6][4]*mom4 +
00147                    invM_S[6][7]*mom7;
00148         
00149         f[7] -= invM_S[7][1]*mom1 +
00150                    invM_S[7][2]*mom2 +
00151                    invM_S[7][3]*mom3 +
00152                    invM_S[7][4]*mom4 +
00153                    invM_S[7][5]*mom5 +
00154                    invM_S[7][6]*mom6 +
00155                    invM_S[7][8]*mom8;
00156         
00157         f[8] -= invM_S[8][1]*mom1 +
00158                    invM_S[8][2]*mom2 +
00159                    invM_S[8][5]*mom5 +
00160                    invM_S[8][6]*mom6 +
00161                    invM_S[8][7]*mom7;
00162 
00163         return jSqr;
00164     }
00165     
00166     static void addGuoForce( Array<T,Descriptor::q>& f, const Array<T,Descriptor::d>& force,
00167                              Array<T,Descriptor::d> const& u,
00168                              T invM_S[Descriptor::q][Descriptor::q], T amplitude )
00169     {
00170         Array<T,Descriptor::q> forcing;
00171         T g_u   = force[0] * u[0] + force[1] * u[1];
00172         T g_u_a = force[0] * u[0] - force[1] * u[1];            
00173         
00174         T mom1 = -(T)3 * g_u;
00175         T mom2 =  (T)3 * g_u;
00176         T mom3 = -(T)0.5 * force[0];
00177         T mom4 = -mom3;
00178         T mom5 = -(T)0.5 * force[1];
00179         T mom6 = -mom5;
00180         T mom7 = -g_u_a;
00181         T mom8 = -(T)0.5 * (force[0]*u[1] + force[1]*u[0]);
00182         
00183         f[0] += invM_S[0][1]*mom1 + 
00184                    invM_S[0][2]*mom2;
00185         
00186         f[1] += invM_S[1][1]*mom1 +
00187                    invM_S[1][2]*mom2 +
00188                    invM_S[1][3]*mom3 +
00189                    invM_S[1][4]*mom4 + 
00190                    invM_S[1][5]*mom5 + 
00191                    invM_S[1][6]*mom6 +
00192                    invM_S[1][8]*mom8;
00193         
00194         f[2] += invM_S[2][1]*mom1 +
00195                    invM_S[2][2]*mom2 +
00196                    invM_S[2][3]*mom3 +
00197                    invM_S[2][4]*mom4 +
00198                    invM_S[2][7]*mom7;
00199         
00200         f[3] += invM_S[3][1]*mom1 +
00201                    invM_S[3][2]*mom2 +
00202                    invM_S[3][3]*mom3 +
00203                    invM_S[3][4]*mom4 +
00204                    invM_S[3][5]*mom5 +
00205                    invM_S[3][6]*mom6 +
00206                    invM_S[3][8]*mom8;
00207         
00208         f[4] += invM_S[4][1]*mom1 +
00209                    invM_S[4][2]*mom2 +
00210                    invM_S[4][5]*mom5 +
00211                    invM_S[4][6]*mom6 +
00212                    invM_S[4][7]*mom7;
00213         
00214         f[5] += invM_S[5][1]*mom1 +
00215                    invM_S[5][2]*mom2 +
00216                    invM_S[5][3]*mom3 +
00217                    invM_S[5][4]*mom4 +
00218                    invM_S[5][5]*mom5 +
00219                    invM_S[5][6]*mom6 +
00220                    invM_S[5][8]*mom8;
00221         
00222         f[6] += invM_S[6][1]*mom1 +
00223                    invM_S[6][2]*mom2 +
00224                    invM_S[6][3]*mom3 +
00225                    invM_S[6][4]*mom4 +
00226                    invM_S[6][7]*mom7;
00227         
00228         f[7] += invM_S[7][1]*mom1 +
00229                    invM_S[7][2]*mom2 +
00230                    invM_S[7][3]*mom3 +
00231                    invM_S[7][4]*mom4 +
00232                    invM_S[7][5]*mom5 +
00233                    invM_S[7][6]*mom6 +
00234                    invM_S[7][8]*mom8;
00235         
00236         f[8] += invM_S[8][1]*mom1 +
00237                    invM_S[8][2]*mom2 +
00238                    invM_S[8][5]*mom5 +
00239                    invM_S[8][6]*mom6 +
00240                    invM_S[8][7]*mom7;
00241                    
00242         f[0] += -(T)2/(T)3*g_u;
00243         f[1] += (-force[0]+force[1]+(T)2*g_u-(T)3*g_u)/(T)12;
00244         f[2] += (-force[0]+2*force[0]*u[0]-force[1]*u[1])/(T)3;
00245         f[3] += (-force[0]-force[1]+(T)2*g_u+(T)3*g_u) / (T)12;
00246         f[4] += -(force[1]+force[0]*u[0] - (T)2*force[1]*u[1]) / (T)3;
00247         f[5] += (force[0]-force[1]+(T)2*g_u-(T)3*g_u) / (T)12;
00248         f[6] += (force[0]+(T)2*force[0]*u[0]-force[1]*u[1]) / (T)3;
00249         f[7] += (force[0]+force[1]+(T)2*g_u+(T)3*g_u) / (T)12;
00250         f[8] += -(-force[1]+force[0]*u[0]-2*force[1]*u[1])/(T)3;
00251     }
00252     
00254     static T mrtCollisionWithForce( Array<T,Descriptor::q>& f,
00255                                     const T &rhoBar, const Array<T,Descriptor::d> & u,
00256                                     T invM_S[Descriptor::q][Descriptor::q], 
00257                                     const Array<T,Descriptor::d> &force, T amplitude) 
00258     {
00259         Array<T,Descriptor::d> j = Descriptor::fullRho(rhoBar)*u;
00260         T jSqr = mrtCollision( f, rhoBar, j, invM_S );
00261         addGuoForce( f, force, u, invM_S, amplitude );
00262         
00263         return jSqr;
00264     }
00265     
00267     static void computeQuasiIncEquilibrium( Array<T,Descriptor::q>& momentsEq,
00268                                     T rhoBar, Array<T,2> const& j, T jSqr )
00269     {
00270         momentsEq[0] = rhoBar;
00271         momentsEq[1] = (T)3*jSqr-2*rhoBar;
00272         momentsEq[2] = -(T)3*jSqr+rhoBar;
00273         momentsEq[3] = j[0];
00274         momentsEq[4] = -j[0];
00275         momentsEq[5] = j[1];
00276         momentsEq[6] = -j[1];
00277         momentsEq[7] = (j[0]*j[0]-j[1]*j[1]);
00278         momentsEq[8] = j[1]*j[0];
00279     }
00280     
00282     static T quasiIncMrtCollision( Array<T,Descriptor::q>& f,
00283                                    const T &rhoBar, const Array<T,2> & j,
00284                                    T invM_S[9][9] )
00285     {
00286         Array<T,9> moments, momentsEq;
00287 
00288         computeMoments(moments,f);
00289         T jSqr = VectorTemplateImpl<T,2>::normSqr(j);
00290         
00291         computeQuasiIncEquilibrium(momentsEq,rhoBar,j,jSqr);
00292         
00293         T mom1 = moments[1] - momentsEq[1];
00294         T mom2 = moments[2] - momentsEq[2];
00295         T mom3 = moments[3] - momentsEq[3];
00296         T mom4 = moments[4] - momentsEq[4];
00297         T mom5 = moments[5] - momentsEq[5];
00298         T mom6 = moments[6] - momentsEq[6];
00299         T mom7 = moments[7] - momentsEq[7];
00300         T mom8 = moments[8] - momentsEq[8];
00301         
00302         f[0] -= invM_S[0][1]*mom1 + 
00303                    invM_S[0][2]*mom2;
00304         
00305         f[1] -= invM_S[1][1]*mom1 +
00306                    invM_S[1][2]*mom2 +
00307                    invM_S[1][4]*mom4 + 
00308                    invM_S[1][6]*mom6 +
00309                    invM_S[1][8]*mom8;
00310         
00311         f[2] -= invM_S[2][1]*mom1 +
00312                    invM_S[2][2]*mom2 +
00313                    invM_S[2][4]*mom4 +
00314                    invM_S[2][7]*mom7;
00315         
00316         f[3] -= invM_S[3][1]*mom1 +
00317                    invM_S[3][2]*mom2 +
00318                    invM_S[3][4]*mom4 +
00319                    invM_S[3][6]*mom6 +
00320                    invM_S[3][8]*mom8;
00321         
00322         f[4] -= invM_S[4][1]*mom1 +
00323                    invM_S[4][2]*mom2 +
00324                    invM_S[4][6]*mom6 +
00325                    invM_S[4][7]*mom7;
00326         
00327         f[5] -= invM_S[5][1]*mom1 +
00328                    invM_S[5][2]*mom2 +
00329                    invM_S[5][4]*mom4 +
00330                    invM_S[5][6]*mom6 +
00331                    invM_S[5][8]*mom8;
00332         
00333         f[6] -= invM_S[6][1]*mom1 +
00334                    invM_S[6][2]*mom2 +
00335                    invM_S[6][4]*mom4 +
00336                    invM_S[6][7]*mom7;
00337         
00338         f[7] -= invM_S[7][1]*mom1 +
00339                    invM_S[7][2]*mom2 +
00340                    invM_S[7][4]*mom4 +
00341                    invM_S[7][6]*mom6 +
00342                    invM_S[7][8]*mom8;
00343         
00344         f[8] -= invM_S[8][1]*mom1 +
00345                    invM_S[8][2]*mom2 +
00346                    invM_S[8][6]*mom6 +
00347                    invM_S[8][7]*mom7;
00348 
00349         return jSqr;
00350     }
00351 
00352 };
00353 
00354 
00355 }  // namespace plb
00356 
00357 #endif