offEquilibriumTemplates.h

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OFF_EQUILIBRIUM_TEMPLATES_H
#define OFF_EQUILIBRIUM_TEMPLATES_H

#include "core/globalDefs.h"
#include "core/cell.h"
#include "core/util.h"
#include "hermitePolynomialsTemplates.h"
#include "geometricOperationTemplates.h"

namespace plb {

template<typename T, class Descriptor> struct offEquilibriumTemplatesImpl;

template<typename T, template<typename U> class Descriptor>
struct offEquilibriumTemplates {


static T fromPiToFneq(plint iPop, Array<T,SymmetricTensor<T,Descriptor>::n> const& pi) {
    return offEquilibriumTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::fromPiToFneq(iPop, pi);
}

static T fromPiAndQtoFneq(plint iPop, Array<T,SymmetricTensor<T,Descriptor>::n> const& pi, 
                          Array<T,SymmetricRankThreeTensor<T,Descriptor>::n> const& q) {
    return offEquilibriumTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::fromPiAndQtoFneq(iPop, pi, q);
}


static T fromStrainToFneq(plint iPop, Array<T,SymmetricTensor<T,Descriptor>::n> const& S, T density, T omega) {
    return offEquilibriumTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::fromStrainToFneq(iPop, S, density, omega);
}

};  // struct offEquilibriumTemplates

template<typename T, class Descriptor>
struct offEquilibriumTemplatesImpl {

static T fromPiToFneq (
    plint iPop, Array<T,SymmetricTensorImpl<T,Descriptor::d>::n> const& PiNeq )
{
    typedef Descriptor L;
    T fNeq = T();
    plint iPi = 0;
    // Iterate only over superior triangle + diagonal, and add
    // the elements under the diagonal by symmetry
    for (int iAlpha=0; iAlpha<L::d; ++iAlpha) {
        // Treat diagonal term first
        fNeq += PiNeq[iPi] * (L::c[iPop][iAlpha]*L::c[iPop][iAlpha]
                              - L::cs2);
        ++iPi;
        // Then, treat off-diagonal terms
        for (int iBeta=iAlpha+1; iBeta<L::d; ++iBeta) {
            // Multiply off-diagonal elements by 2 because
            // the Q tensor is symmetric
            fNeq += PiNeq[iPi]
                      * (T)2 * L::c[iPop][iAlpha]*L::c[iPop][iBeta];
            ++iPi;
        }
    }
    fNeq *= L::t[iPop] * L::invCs2 * L::invCs2 / (T)2;
    return fNeq;
}

static T fromPiAndQtoFneq (
    plint iPop, Array<T,SymmetricTensorImpl<T,Descriptor::d>::n> const& PiNeq,
    Array<T,SymmetricRankThreeTensorImpl<T,Descriptor::d>::n> const& Q)
{
    typedef Descriptor L;
    
    Array<T,SymmetricTensorImpl<T,Descriptor::d>::n> H2 = 
        HermiteTemplateImpl<T,Descriptor,Descriptor::d>::order2(iPop);
    Array<T,SymmetricRankThreeTensorImpl<T,Descriptor::d>::n> H3 = 
        HermiteTemplateImpl<T,Descriptor,Descriptor::d>::order3(iPop);
        
    T factor = L::t[iPop] * L::invCs2 * L::invCs2 / (T)2;
    
    T fNeqPi = factor * SymmetricTensorImpl<T,Descriptor::d>::contractIndexes(H2,PiNeq);
    T fNeqQ = factor * L::invCs2 / (T)3 * 
                SymmetricRankThreeTensorImpl<T,Descriptor::d>::contractIndexes(H3,Q);
    
    return fNeqPi + fNeqQ;
}


static T fromStrainToFneq (
    plint iPop, Array<T,SymmetricTensorImpl<T,Descriptor::d>::n> const& S, T density, T omega)
{
    typedef Descriptor L;
    T fNeq = fromPiToFneq(iPop,S) * (-(T)2 * density * L::cs2 / omega);
    return fNeq;
}

};  // struct offEquilibriumTemplates

}  // namespace plb

#include "latticeBoltzmann/offEquilibriumTemplates2D.h"
#include "latticeBoltzmann/offEquilibriumTemplates3D.h"

#endif