externalForceTemplates3D.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 /* Main author: Orestis Malaspinas
00026  */
00027 
00031 #ifndef EXTERNAL_FORCE_TEMPLATES_3D_H
00032 #define EXTERNAL_FORCE_TEMPLATES_3D_H
00033 
00034 namespace plb {
00035 
00036 template<typename T>
00037 struct externalForceTemplatesImpl<T, descriptors::ForcedD3Q19Descriptor<T> > 
00038 {
00039 
00040 static void addNaiveForce(
00041                 Array<T,descriptors::ForcedD3Q19Descriptor<T>::q>& f,
00042                 T* externalScalars,
00043                 Array<T,descriptors::ForcedD3Q19Descriptor<T>::d> const& u, T omega, T amplitude )
00044 {
00045     static const int forceBeginsAt
00046         = descriptors::ForcedD3Q19Descriptor<T>::ExternalField::forceBeginsAt;
00047     T* force = externalScalars + forceBeginsAt;
00048     T& fx = force[0];
00049     T& fy = force[1];
00050     T& fz = force[2];
00051 
00052     f[1] +=  (-fx       )* (T)1/(T)6;
00053     f[2] +=  (   -fy    )* (T)1/(T)6;
00054     f[3] +=  (      -fz )* (T)1/(T)6;
00055     f[4] +=  (-fx-fy    )* (T)1/(T)12;
00056     f[5] +=  (-fx+fy    )* (T)1/(T)12;
00057     f[6] +=  (-fx   -fz )* (T)1/(T)12;
00058     f[7] +=  (-fx   +fz )* (T)1/(T)12;
00059     f[8] +=  (   -fy-fz )* (T)1/(T)12;
00060     f[9] +=  (   -fy+fz )* (T)1/(T)12;
00061                                    
00062     f[10] += ( fx       )* (T)1/(T)6;
00063     f[11] += (    fy    )* (T)1/(T)6;
00064     f[12] += (       fz )* (T)1/(T)6;
00065     f[13] += ( fx+fy    )* (T)1/(T)12;
00066     f[14] += ( fx-fy    )* (T)1/(T)12;
00067     f[15] += ( fx   +fz )* (T)1/(T)12;
00068     f[16] += ( fx   -fz )* (T)1/(T)12;
00069     f[17] += (    fy+fz )* (T)1/(T)12;
00070     f[18] += (    fy-fz )* (T)1/(T)12;
00071 
00072 }
00073 
00074 static void addGuoForce(
00075                 Array<T,descriptors::ForcedD3Q19Descriptor<T>::q>& f,
00076                 T* externalScalars,
00077                 Array<T,descriptors::ForcedD3Q19Descriptor<T>::d> const& u, T omega, T amplitude )
00078 {
00079     static const int forceBeginsAt
00080         = descriptors::ForcedD3Q19Descriptor<T>::ExternalField::forceBeginsAt;
00081     T* force = externalScalars + forceBeginsAt;
00082     T mu = amplitude*((T)1-omega/(T)2);
00083     
00084     static const T oneOver6 = (T)1/(T)6;
00085     static const T oneOver12 = (T)1/(T)12;
00086     
00087     f[0] += -mu*(force[0]*u[0]+force[1]*u[1]+force[2]*u[2]);
00088     
00089     f[1] += oneOver6*mu*(force[0]*(-(T)1+2*u[0])-force[1]*u[1]-force[2]*u[2]);
00090     
00091     f[2] += -oneOver6*mu*(force[0]*u[0]+force[1]*((T)1-2*u[1])+force[2]*u[2]);
00092     
00093     f[3] += -oneOver6*mu*(force[0]*u[0]
00094                            + force[1]*u[1]
00095                            + force[2]*((T)1-2*u[2]));
00096     
00097     f[4] += oneOver12*mu*(force[0]*(-(T)1+2*u[0]+3*u[1])
00098                            + force[1]*(-(T)1+2*u[1]+3*u[0])
00099                            - force[2]*u[2]);
00100     
00101     f[5] += oneOver12*mu*( force[0]*(-(T)1+2*u[0]-3*u[1])
00102                             + force[1]*((T)1+2*u[1]-3*u[0])
00103                             - force[2]*u[2]);
00104     
00105     f[6] += oneOver12*mu*(force[0]*(-(T)1+2*u[0]+3*u[2])
00106                            - force[1]*u[1]
00107                            + force[2]*(-(T)1+2*u[2]+3*u[0]));
00108     
00109     f[7] += oneOver12*mu*(force[0]*(-(T)1+2*u[0]-3*u[2])
00110                            - force[1]*u[1]
00111                            + force[2]*((T)1+2*u[2]-3*u[0]));
00112     
00113     f[8] += -oneOver12*mu*(force[0]*u[0]
00114                             + force[1]*((T)1-2*u[1]-3*u[2])
00115                             + force[2]*((T)1-2*u[2]-3*u[1]));
00116     
00117     f[9] += -oneOver12*mu*(force[0]*u[0]
00118                             + force[1]*((T)1-2*u[1]+3*u[2])
00119                             + force[2]*(-(T)1-2*u[2]+3*u[1]));
00120     
00121     f[10] += oneOver6*mu*(force[0]*((T)1+2*u[0])
00122                             -force[1]*u[1]
00123                             -force[2]*u[2]);
00124     
00125     f[11] += -oneOver6*mu*(force[0]*u[0]
00126                             +force[1]*(-(T)1-2*u[1])
00127                             +force[2]*u[2]);
00128     
00129     f[12] += -oneOver6*mu*(force[0]*u[0]
00130                              +force[1]*u[1]
00131                              +force[2]*(-(T)1-2*u[2]));
00132     
00133     f[13] += oneOver12*mu*(force[0]*((T)1+2*u[0]+3*u[1])
00134                              +force[1]*((T)1+2*u[1]+3*u[0])
00135                              -force[2]*u[2]);
00136     
00137     f[14] += oneOver12*mu*(force[0]*((T)1+2*u[0]-3*u[1])
00138                              +force[1]*(-(T)1+2*u[1]-3*u[0])
00139                              -force[2]*u[2]);
00140     
00141     f[15] += oneOver12*mu*(force[0]*((T)1+2*u[0]+3*u[2])
00142                              -force[1]*u[1]
00143                              +force[2]*((T)1+2*u[2]+3*u[0]));
00144     
00145     f[16] += oneOver12*mu*(force[0]*((T)1+2*u[0]-3*u[2])
00146                              -force[1]*u[1]
00147                              +force[2]*(-(T)1+2*u[2]-3*u[0]));
00148     
00149     f[17] += -oneOver12*mu*(force[0]*u[0]
00150                               +force[1]*(-(T)1-2*u[1]-3*u[2])
00151                               +force[2]*(-(T)1-2*u[2]-3*u[1]));
00152     
00153     f[18] += -oneOver12*mu*(force[0]*u[0]
00154                               +force[1]*(-(T)1-2*u[1]+3*u[2])
00155                               +force[2]*((T)1-2*u[2]+3*u[1]));
00156 }
00157 
00158 static T heForcedBGKCollision (
00159         Array<T,descriptors::ForcedD3Q19Descriptor<T>::q>& f, T* externalScalars,
00160         T rhoBar, Array<T,descriptors::ForcedD3Q19Descriptor<T>::d> const& j, T omega )
00161 {
00162     static const int forceBeginsAt
00163         = descriptors::ForcedD3Q19Descriptor<T>::ExternalField::forceBeginsAt;
00164     T* force = externalScalars + forceBeginsAt;
00165     T invRho = descriptors::ForcedD3Q19Descriptor<T>::invRho(rhoBar);
00166     const T jSqr = VectorTemplateImpl<T,descriptors::ForcedD3Q19Descriptor<T>::d>::normSqr(j);
00167         
00168         Array<T,descriptors::ForcedD3Q19Descriptor<T>::d> uLB(j[0]*invRho, j[1]*invRho, j[2]*invRho);
00169         const T uSqrLB = VectorTemplateImpl<T,descriptors::ForcedD3Q19Descriptor<T>::d>::normSqr(uLB);
00170     
00171     dynamicsTemplatesImpl<T,descriptors::D3Q19DescriptorBase<T> >::bgk_ma2_collision(f, rhoBar, j, omega);
00172         
00173     for (plint iPop=0; iPop < descriptors::ForcedD3Q19Descriptor<T>::q; ++iPop) {
00174         T ug = (descriptors::ForcedD3Q19Descriptor<T>::c[iPop][0]-uLB[0])*force[0] +
00175                (descriptors::ForcedD3Q19Descriptor<T>::c[iPop][1]-uLB[1])*force[1] +
00176                (descriptors::ForcedD3Q19Descriptor<T>::c[iPop][2]-uLB[2])*force[2];
00177                 
00178                 f[iPop] += descriptors::ForcedD3Q19Descriptor<T>::invCs2 * ug * dynamicsTemplatesImpl<T,descriptors::ForcedD3Q19Descriptor<T> >
00179                         ::bgk_ma2_equilibrium( iPop, (T)1, (T)1, uLB, uSqrLB );
00180     }
00181     T uSqr = util::sqr(uLB[0] + 0.5*force[0]) +
00182                          util::sqr(uLB[1] + 0.5*force[1]) +
00183                          util::sqr(uLB[2] + 0.5*force[2]);
00184     return uSqr;
00185 }
00186 
00187 };
00188 
00189 }  // namespace plb
00190 
00191 #endif