lbhelpers.h

open lattice boltzmann project www.openlb.org

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2006, 2007 Jonas Latt
 *  Address: Rue General Dufour 24,  1211 Geneva 4, Switzerland 
 *  E-mail: [e-mail]
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public 
 *  License along with this program; if not, write to the Free 
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/

/** \file
 * Helper functions for the implementation of LB dynamics. This file is all
 * about efficiency. The generic template code is specialized for commonly
 * used Lattices, so that a maximum performance can be taken out of each
 * case.
 */
#ifndef LB_HELPERS_H
#define LB_HELPERS_H

#include "latticeDescriptors.h"
#include "cell.h"
#include "util.h"


namespace olb {


// Forward declarations
template<typename T, class Descriptor> struct lbDynamicsHelpers;
template<typename T, template<typename U> class Lattice> struct lbExternalHelpers;
template<typename T, template<typename U> class Lattice> struct lbLatticeHelpers;

/// This structure forwards the calls to the appropriate helper class
template<typename T, template<typename U> class Lattice>
struct lbHelpers {

    static T equilibrium(int iPop, T rho, const T u[Lattice<T>::d], const T uSqr) {
        return lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::equilibrium(iPop, rho, u, uSqr);
    }

    static T incEquilibrium(int iPop, const T j[Lattice<T>::d], const T jSqr, const T pressure) {
        return lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::incEquilibrium(iPop, j, jSqr, pressure);
    }

    static void computeFneq ( Cell<T,Lattice> const& cell,
                              T fNeq[Lattice<T>::q], T rho, const T u[Lattice<T>::d] )
    {
        lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::computeFneq(&cell[0], fNeq, rho, u);
    }

    static T bgkCollision(Cell<T,Lattice>& cell, T rho, const T u[Lattice<T>::d], T omega)
    {
        return lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::bgkCollision(&cell[0], rho, u, omega);
    }

    static T incBgkCollision(Cell<T,Lattice>& cell, T pressure, const T j[Lattice<T>::d], T omega)
    {
        return lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::incBgkCollision(&cell[0], pressure, j, omega);
    }

    static T constRhoBgkCollision(Cell<T,Lattice>& cell,
                                  T rho, const T u[Lattice<T>::d], T ratioRho, T omega)
    {
        return lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
                   ::constRhoBgkCollision(&cell[0], rho, u, ratioRho, omega);
    }

    static T computeRho(Cell<T,Lattice> const& cell) {
        return lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::computeRho(&cell[0]);
    }

    static void computeJ(Cell<T,Lattice> const& cell, T j[Lattice<T>::d] ) {
        lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::computeJ(&cell[0], j);
    }

    static void computeRhoU(Cell<T,Lattice> const& cell, T& rho, T u[Lattice<T>::d])
    {
        lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::computeRhoU(&cell[0], rho, u);
    }

    static void computeStress(Cell<T,Lattice> const& cell, T rho, const T u[Lattice<T>::d],
                              T pi[util::TensorVal<Lattice<T> >::n] )
    {
        lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::computeStress(&cell[0], rho, u, pi);
    }
    
    static void computeAllMomenta(Cell<T,Lattice> const& cell, T& rho, T u[Lattice<T>::d],
                                  T pi[util::TensorVal<Lattice<T> >::n] )
    {
        lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::computeAllMomenta(&cell[0], rho, u, pi);
    }

    static void modifyVelocity(Cell<T,Lattice>& cell, const T newU[Lattice<T>::d]) {
        lbDynamicsHelpers<T,typename Lattice<T>::BaseDescriptor>
            ::modifyVelocity(&cell[0], newU);
    }

    static void addExternalForce(Cell<T,Lattice>& cell, const T u[Lattice<T>::d], T omega, T amplitude=(T)1)
    {
        lbExternalHelpers<T,Lattice>::addExternalForce(cell, u, omega, amplitude);
    }

    static void swapAndStream2D(Cell<T,Lattice> **grid, int iX, int iY)
    {
        lbLatticeHelpers<T,Lattice>::swapAndStream2D(grid, iX, iY);
    }

    static void swapAndStream3D(Cell<T,Lattice> ***grid, int iX, int iY, int iZ)
    {
        lbLatticeHelpers<T,Lattice>::swapAndStream3D(grid, iX, iY, iZ);
    }

};  // struct lbHelpers


/// All helper functions are inside this structure
template<typename T, class Descriptor>
struct lbDynamicsHelpers {
    /// Computation of equilibrium distribution
    static T equilibrium(int iPop, T rho, const T u[Descriptor::d], const T uSqr) {
        T c_u = T();
        for (int iD=0; iD < Descriptor::d; ++iD) {
           c_u += Descriptor::c[iPop][iD]*u[iD];
        }
        return rho * Descriptor::t[iPop] * (
               (T)1 + Descriptor::invCs2 * c_u +
               Descriptor::invCs2 * Descriptor::invCs2/(T)2 * c_u*c_u -
               Descriptor::invCs2/(T)2 * uSqr
           ) - Descriptor::t[iPop];
    }

    static T incEquilibrium( int iPop, const T j[Descriptor::d],
                             const T jSqr, const T pressure )
    {
        T c_j = T();
        for (int iD=0; iD < Descriptor::d; ++iD) {
           c_j += Descriptor::c[iPop][iD]*j[iD];
        }
        T rho = (T)1 + pressure*Descriptor::invCs2;
        return Descriptor::t[iPop] * ( rho +
               Descriptor::invCs2 * c_j +
               Descriptor::invCs2 * Descriptor::invCs2/(T)2 * c_j*c_j -
               Descriptor::invCs2/(T)2 * jSqr
           ) - Descriptor::t[iPop];
    }

    static void computeFneq(T const* cell, T fNeq[Descriptor::q], T rho, const T u[Descriptor::d]) {
        const T uSqr = util::normSqr<T,Descriptor::d>(u);
        for (int iPop=0; iPop < Descriptor::q; ++iPop) {
            fNeq[iPop] = cell[iPop] - equilibrium(iPop, rho, u, uSqr);
        }
    }

    /// BGK collision step
    static T bgkCollision(T* cell, T rho, const T u[Descriptor::d], T omega) {
        const T uSqr = util::normSqr<T,Descriptor::d>(u);
        for (int iPop=0; iPop < Descriptor::q; ++iPop) {
            cell[iPop] *= (T)1-omega;
            cell[iPop] += omega * lbDynamicsHelpers<T,Descriptor>::equilibrium (
                                    iPop, rho, u, uSqr );
        }
        return uSqr;
    }

    /// Incompressible BGK collision step
    static T incBgkCollision(T* cell, T pressure, const T j[Descriptor::d], T omega) {
        const T jSqr = util::normSqr<T,Descriptor::d>(j);
        for (int iPop=0; iPop < Descriptor::q; ++iPop) {
            cell[iPop] *= (T)1-omega;
            cell[iPop] += omega * lbDynamicsHelpers<T,Descriptor>::incEquilibrium (
                              iPop, j, jSqr, pressure );
        }
        return jSqr;
    }

    /// BGK collision step with density correction
    static T constRhoBgkCollision(T* cell, T rho, const T u[Descriptor::d], T ratioRho, T omega) {
        const T uSqr = util::normSqr<T,Descriptor::d>(u);
        for (int iPop=0; iPop < Descriptor::q; ++iPop) {
            T feq = lbDynamicsHelpers<T,Descriptor>::equilibrium(iPop, rho, u, uSqr );
            cell[iPop] =
              ratioRho*(feq+Descriptor::t[iPop])-Descriptor::t[iPop] +
              ((T)1-omega)*(cell[iPop]-feq);
        }
        return uSqr;
    }

    /// Computation of density
    static T computeRho(T const* cell) {