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