entropicLbTemplates2D.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 
00032 #ifndef ENTROPIC_LB_HELPERS_2D_H
00033 #define ENTROPIC_LB_HELPERS_2D_H
00034 
00035 #include "core/globalDefs.h"
00036 #include "latticeBoltzmann/nearestNeighborLattices2D.h"
00037 #include <cmath>
00038 
00039 namespace plb {
00040 
00041 template<typename T>
00042 struct entropicLbTemplates<T, descriptors::D2Q9Descriptor>
00043 {
00044 
00045     typedef descriptors::D2Q9Descriptor<T> Descriptor;
00048     static T equilibrium( plint iPop, T rho, Array<T,2> const& u)
00049     {
00050         T c_u = Descriptor::c[iPop][0]*u[0] + Descriptor::c[iPop][1]*u[1];
00051         T c_u2 = c_u*c_u;
00052         T c_u3 = c_u2*c_u;
00053         T c_u4 = c_u3*c_u;
00054         T c_u5 = c_u4*c_u;
00055         T c_u6 = c_u5*c_u;
00056         T c_u7 = c_u6*c_u;
00057         
00058         T uSqr = u[0]*u[0] + u[1]*u[1];
00059         T uSqr2 = uSqr*uSqr;
00060         T uSqr3 = uSqr2*uSqr;
00061         
00062         T powUx = u[0]*u[0]*u[0]*u[0]*u[0]; // u_x^5
00063         T powUy = u[1]*u[1]*u[1]*u[1]*u[1]; // u_y^5
00064         
00065         T C = Descriptor::c[iPop][0] * powUx + Descriptor::c[iPop][1] * powUy;
00066         
00067         powUx *= u[0]; // u_x^6
00068         powUy *= u[1]; // u_y^6
00069         
00070         T E = powUx + powUy;
00071         
00072         powUx *= u[0]; // u_x^7
00073         powUy *= u[1]; // u_y^7
00074         
00075         T F = Descriptor::c[iPop][0] * powUx + Descriptor::c[iPop][1] * powUy;
00076         
00077         return Descriptor::t[iPop] * rho * 
00078               ((T)1
00079                 + c_u*(C*(T)81/(T)20 + uSqr2*(T)27/(T)8 - uSqr*(T)9/(T)2
00080                 - E*(T)81/(T)24 - uSqr3*(T)81/(T)48 + (T)3)
00081                 
00082                 + c_u2*(uSqr2*(T)81/(T)16 - uSqr*(T)27/(T)4 
00083                 + C*(T)243/(T)40 + (T)9/(T)2)
00084                 
00085                 + c_u3*(uSqr2*(T)243/(T)48 - uSqr*(T)81/(T)12 + (T)27/(T)6)
00086                 
00087                 - c_u4*uSqr*(T)243/(T)48
00088                 + c_u4*(T)81/(T)24
00089                 
00090                 - c_u5*uSqr*(T)729/(T)240
00091                 + c_u5*(T)243/(T)120
00092                 
00093                 + c_u6*(T)729/(T)720
00094                 
00095                 + c_u7*(T)2187/(T)5040
00096                 
00097                 - C*uSqr*(T)81/(T)40
00098                 
00099                 + C*(T)27/(T)20 - uSqr3*(T)27/(T)48 - E*(T)27/(T)24
00100                 - F*(T)81/(T)56 - uSqr*(T)3/(T)2 + uSqr2*(T)9/(T)8
00101                 )
00102                 - Descriptor::SkordosFactor() * Descriptor::t[iPop];
00103     }
00104     
00105 };
00106 
00107 }
00108 
00109 #endif