From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 24 Jun 2019 14:43:36 +0200 Subject: Initialize at openlb-1-3 --- src/dynamics/lbHelpers.h | 678 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 678 insertions(+) create mode 100644 src/dynamics/lbHelpers.h (limited to 'src/dynamics/lbHelpers.h') diff --git a/src/dynamics/lbHelpers.h b/src/dynamics/lbHelpers.h new file mode 100644 index 0000000..e0d9ccd --- /dev/null +++ b/src/dynamics/lbHelpers.h @@ -0,0 +1,678 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2006, 2007 Jonas Latt + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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 "core/cell.h" +#include "core/util.h" + + +namespace olb { + + +// Forward declarations +template struct lbDynamicsHelpers; +template struct lbExternalHelpers; +template struct lbLatticeHelpers; + +/// This structure forwards the calls to the appropriate helper class +template +struct lbHelpers { + + static T equilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], const T uSqr) + { + return lbDynamicsHelpers + ::equilibrium(iPop, rho, u, uSqr); + } + + static T equilibriumFirstOrder(int iPop, T rho, const T u[DESCRIPTOR::d]) + { + return lbDynamicsHelpers + ::equilibriumFirstOrder(iPop, rho, u); + } + + static T incEquilibrium(int iPop, const T j[DESCRIPTOR::d], const T jSqr, const T pressure) + { + return lbDynamicsHelpers + ::incEquilibrium(iPop, j, jSqr, pressure); + } + + static void computeFneq ( Cell const& cell, + T fNeq[DESCRIPTOR::q], T rho, const T u[DESCRIPTOR::d] ) + { + lbDynamicsHelpers::computeFneq(cell, fNeq, rho, u); + } + + static T bgkCollision(Cell& cell, T const& rho, const T u[DESCRIPTOR::d], T const& omega) + { + return lbDynamicsHelpers + ::bgkCollision(cell, rho, u, omega); + } + + static T bgkCollision(Cell& cell, T const& rho, const T u[DESCRIPTOR::d], const T omega[DESCRIPTOR::q]) + { + const T uSqr = util::normSqr(u); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] *= (T)1-omega[iPop]; + cell[iPop] += omega[iPop] * lbDynamicsHelpers::equilibrium (iPop, rho, u, uSqr ); + } + return uSqr; + } + + static T rlbCollision( Cell& cell, T rho, const T u[DESCRIPTOR::d], T omega ) + { + return lbDynamicsHelpers + ::rlbCollision( cell, rho, u, omega ); + } + + static T mrtCollision( Cell& cell, T const& rho, const T u[DESCRIPTOR::d], T invM_S[DESCRIPTOR::q][DESCRIPTOR::q] ) { + + return lbDynamicsHelpers + ::mrtCollision(cell, rho, u, invM_S ); + } + + static T incBgkCollision(Cell& cell, T pressure, const T j[DESCRIPTOR::d], T omega) + { + return lbDynamicsHelpers + ::incBgkCollision(cell, pressure, j, omega); + } + + static T constRhoBgkCollision(Cell& cell, + T rho, const T u[DESCRIPTOR::d], T ratioRho, T omega) + { + return lbDynamicsHelpers + ::constRhoBgkCollision(cell, rho, u, ratioRho, omega); + } + + static T computeRho(Cell const& cell) + { + return lbDynamicsHelpers + ::computeRho(cell); + } + + static void computeJ(Cell const& cell, T j[DESCRIPTOR::d] ) + { + lbDynamicsHelpers + ::computeJ(cell, j); + } + + static void computeRhoU(Cell const& cell, T& rho, T u[DESCRIPTOR::d]) + { + lbDynamicsHelpers + ::computeRhoU(cell, rho, u); + } + + static void computeFeq(Cell const& cell, T fEq[DESCRIPTOR::q]) + { + T rho{}; + T u[2] {}; + computeRhoU(cell, rho, u); + const T uSqr = u[0]*u[0] + u[1]*u[1]; + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fEq[iPop] = equilibrium(iPop, rho, u, uSqr); + } + } + + static void computeFneq(Cell const& cell, T fNeq[DESCRIPTOR::q]) + { + T rho{}; + T u[2] {}; + computeRhoU(cell, rho, u); + computeFneq(cell, fNeq, rho, u); + } + + static void computeRhoJ(Cell const& cell, T& rho, T j[DESCRIPTOR::d]) + { + lbDynamicsHelpers + ::computeRhoJ(cell, rho, j); + } + + static void computeStress(Cell const& cell, T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) + { + lbDynamicsHelpers + ::computeStress(cell, rho, u, pi); + } + + static void computeAllMomenta(Cell const& cell, T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) + { + lbDynamicsHelpers + ::computeAllMomenta(cell, rho, u, pi); + } + + static T computePiNeqNormSqr(Cell const& cell) + { + //return lbDynamicsHelpers + //::computePiNeqNormSqr(cell); + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + computeAllMomenta(cell, rho, u, pi); + T PiNeqNormSqr = pi[0]*pi[0] + 2.*pi[1]*pi[1] + pi[2]*pi[2]; + if (util::TensorVal::n == 6) { + PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.*pi[4]*pi[4] +pi[5]*pi[5]; + } + return PiNeqNormSqr; + } + + static T computeForcedPiNeqNormSqr(Cell const& cell) + { + //return lbDynamicsHelpers + //::computeForcedPiNeqNormSqr(cell); + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + computeAllMomenta(cell, rho, u, pi); + //Creation of body force tensor (rho/2.)*(G_alpha*U_beta + U_alpha*G_Beta) + T ForceTensor[util::TensorVal::n]; + int iPi = 0; + for (int Alpha=0; Alpha()[Alpha]*u[Beta] + u[Alpha]*cell.template getFieldPointer()[Beta]); + ++iPi; + } + } + // Creation of second-order moment off-equilibrium tensor + for (int iPi=0; iPi < util::TensorVal::n; ++iPi){ + pi[iPi] += ForceTensor[iPi]; + } + T PiNeqNormSqr = pi[0]*pi[0] + 2.*pi[1]*pi[1] + pi[2]*pi[2]; + if (util::TensorVal::n == 6) { + PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.*pi[4]*pi[4] +pi[5]*pi[5]; + } + return PiNeqNormSqr; + } + + static void modifyVelocity(Cell& cell, const T newU[DESCRIPTOR::d]) + { + lbDynamicsHelpers + ::modifyVelocity(cell, newU); + } + + static void addExternalForce(Cell& cell, const T u[DESCRIPTOR::d], T omega, T amplitude=(T)1) + { + lbExternalHelpers::addExternalForce(cell, u, omega, amplitude); + } + + static void swapAndStream2D(Cell **grid, int iX, int iY) + { + lbLatticeHelpers::swapAndStream2D(grid, iX, iY); + } + + static void swapAndStream3D(Cell ***grid, int iX, int iY, int iZ) + { + lbLatticeHelpers::swapAndStream3D(grid, iX, iY, iZ); + } + +}; // struct lbHelpers + + +/// All helper functions are inside this structure +template +struct lbDynamicsHelpers { + /// Computation of equilibrium distribution, second order in u + static T equilibrium(int iPop, T rho, const T u[DESCRIPTORBASE::d], const T uSqr) + { + T c_u = T(); + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + c_u += descriptors::c(iPop,iD)*u[iD]; + } + return rho * descriptors::t(iPop) * ( (T)1 + descriptors::invCs2() * c_u + + descriptors::invCs2() * descriptors::invCs2() * (T)0.5 * c_u *c_u + - descriptors::invCs2() * (T)0.5 * uSqr ) + - descriptors::t(iPop); + } + + /// Computation of equilibrium distribution, first order in u + static T equilibriumFirstOrder(int iPop, T rho, const T u[DESCRIPTORBASE::d]) + { + T c_u = T(); + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + c_u += descriptors::c(iPop,iD)*u[iD]; + } + return rho * descriptors::t(iPop) * ( (T)1 + c_u * descriptors::invCs2() ) - descriptors::t(iPop); + } + + static T incEquilibrium( int iPop, const T j[DESCRIPTORBASE::d], + const T jSqr, const T pressure ) + { + T c_j = T(); + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + c_j += descriptors::c(iPop,iD)*j[iD]; + } + + return descriptors::t(iPop) * ( descriptors::invCs2() * pressure + + descriptors::invCs2() * c_j + + descriptors::invCs2() * descriptors::invCs2()/(T)2 * c_j*c_j - + descriptors::invCs2()/(T)2 * jSqr + ) - descriptors::t(iPop); + } + + /// Computation of non-equilibrium distribution + static void computeFneq(CellBase const& cell, T fNeq[DESCRIPTORBASE::q], T rho, const T u[DESCRIPTORBASE::d]) + { + const T uSqr = util::normSqr(u); + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + fNeq[iPop] = cell[iPop] - equilibrium(iPop, rho, u, uSqr); + } + } + + /// BGK collision step + static T bgkCollision(CellBase& cell, T const& rho, const T u[DESCRIPTORBASE::d], T const& omega) + { + const T uSqr = util::normSqr(u); + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + cell[iPop] *= (T)1-omega; + cell[iPop] += omega * lbDynamicsHelpers::equilibrium(iPop, rho, u, uSqr ); + } + return uSqr; + } + + /// Incompressible BGK collision step + static T incBgkCollision(CellBase& cell, T pressure, const T j[DESCRIPTORBASE::d], T omega) + { + const T jSqr = util::normSqr(j); + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + cell[iPop] *= (T)1-omega; + cell[iPop] += omega * lbDynamicsHelpers::incEquilibrium ( + iPop, j, jSqr, pressure ); + } + return jSqr; + } + + /// BGK collision step with density correction + static T constRhoBgkCollision(CellBase& cell, T rho, const T u[DESCRIPTORBASE::d], T ratioRho, T omega) + { + const T uSqr = util::normSqr(u); + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + T feq = lbDynamicsHelpers::equilibrium(iPop, rho, u, uSqr ); + cell[iPop] = + ratioRho*(feq+descriptors::t(iPop))-descriptors::t(iPop) + + ((T)1-omega)*(cell[iPop]-feq); + } + return uSqr; + } + + /// RLB advection diffusion collision step + static T rlbCollision(CellBase& cell, T rho, const T u[DESCRIPTORBASE::d], T omega ) + { + const T uSqr = util::normSqr( u ); + // First-order moment for the regularization + T j1[DESCRIPTORBASE::d]; + for ( int iD = 0; iD < DESCRIPTORBASE::d; ++iD ) { + j1[iD] = T(); + } + + T fEq[DESCRIPTORBASE::q]; + for ( int iPop = 0; iPop < DESCRIPTORBASE::q; ++iPop ) { + fEq[iPop] = lbDynamicsHelpers::equilibriumFirstOrder( iPop, rho, u ); + for ( int iD = 0; iD < DESCRIPTORBASE::d; ++iD ) { + j1[iD] += descriptors::c(iPop,iD) * ( cell[iPop] - fEq[iPop] ); + } + } + + // Collision step + for ( int iPop = 0; iPop < DESCRIPTORBASE::q; ++iPop ) { + T fNeq = T(); + for ( int iD = 0; iD < DESCRIPTORBASE::d; ++iD ) { + fNeq += descriptors::c(iPop,iD) * j1[iD]; + } + fNeq *= descriptors::t(iPop) * descriptors::invCs2(); + cell[iPop] = fEq[iPop] + ( (T)1 - omega ) * fNeq; + } + return uSqr; + } + + +///// Computation of all equilibrium distribution (in momenta space) + static void computeMomentaEquilibrium( T momentaEq[DESCRIPTORBASE::q], T rho, const T u[DESCRIPTORBASE::d], T uSqr ) + { + for (int iPop = 0; iPop < DESCRIPTORBASE::q; ++iPop) { + momentaEq[iPop] = T(); + for (int jPop = 0; jPop < DESCRIPTORBASE::q; ++jPop) { + momentaEq[iPop] += DESCRIPTORBASE::M[iPop][jPop] * + (lbDynamicsHelpers::equilibrium(jPop,rho,u,uSqr) + DESCRIPTORBASE::t[jPop]); + } + } + } + + static void computeMomenta(T momenta[DESCRIPTORBASE::q], CellBase& cell) + { + for (int iPop = 0; iPop < DESCRIPTORBASE::q; ++iPop) { + momenta[iPop] = T(); + for (int jPop = 0; jPop < DESCRIPTORBASE::q; ++jPop) { + momenta[iPop] += DESCRIPTORBASE::M[iPop][jPop] * + (cell[jPop] + DESCRIPTORBASE::t[jPop]); + } + } + } + + static T mrtCollision( CellBase& cell, T const& rho, const T u[DESCRIPTORBASE::d], T invM_S[DESCRIPTORBASE::q][DESCRIPTORBASE::q] ) + { + //// Implemented in advectionDiffusionMRTlbHelpers2D.h and advectionDiffusionMRTlbHelpers3D.h + T uSqr = util::normSqr(u); + T momenta[DESCRIPTORBASE::q]; + T momentaEq[DESCRIPTORBASE::q]; + + computeMomenta(momenta, cell); + computeMomentaEquilibrium(momentaEq, rho, u, uSqr); + +// std::cout << "momenta = "; +// for (int i=0; i < DESCRIPTORBASE::q; ++i) { +// std::cout << momenta[i] << ", "; +// } +// std::cout << std::endl; + +// std::cout << "momentaEq = "; +// for (int i=0; i < DESCRIPTORBASE::q; ++i) { +// std::cout << momentaEq[i] << ", "; +// } +// std::cout << std::endl; + + for (int iPop = 0; iPop < DESCRIPTORBASE::q; ++iPop) { + T collisionTerm = T(); + for (int jPop = 0; jPop < DESCRIPTORBASE::q; ++jPop) { + collisionTerm += invM_S[iPop][jPop] * (momenta[jPop] - momentaEq[jPop]); + } + cell[iPop] -= collisionTerm; + } + return uSqr; + } + + /// Computation of density + static T computeRho(CellBase const& cell) + { + T rho = T(); + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + rho += cell[iPop]; + } + rho += (T)1; + return rho; + } + + /// Computation of momentum + static void computeJ(CellBase const& cell, T j[DESCRIPTORBASE::d]) + { + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + j[iD] = T(); + } + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + j[iD] += cell[iPop]*descriptors::c(iPop,iD); + } + } + } + + /// Computation of hydrodynamic variables + static void computeRhoU(CellBase const& cell, T& rho, T u[DESCRIPTORBASE::d]) + { + rho = T(); + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + u[iD] = T(); + } + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + rho += cell[iPop]; + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + u[iD] += cell[iPop]*descriptors::c(iPop,iD); + } + } + rho += (T)1; + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + u[iD] /= rho; + } + } + + /// Computation of hydrodynamic variables + static void computeRhoJ(CellBase const& cell, T& rho, T j[DESCRIPTORBASE::d]) + { + rho = T(); + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + j[iD] = T(); + } + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + rho += cell[iPop]; + for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) { + j[iD] += cell[iPop]*descriptors::c(iPop,iD); + } + } + rho += (T)1; + } + + /// Computation of stress tensor + static void computeStress(CellBase const& cell, T rho, const T u[DESCRIPTORBASE::d], + T pi[util::TensorVal::n] ) + { + int iPi = 0; + for (int iAlpha=0; iAlpha < DESCRIPTORBASE::d; ++iAlpha) { + for (int iBeta=iAlpha; iBeta < DESCRIPTORBASE::d; ++iBeta) { + pi[iPi] = T(); + for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) { + pi[iPi] += descriptors::c(iPop,iAlpha)* + descriptors::c(iPop,iBeta) * cell[iPop]; + } + // stripe off equilibrium contribution + pi[iPi] -= rho*u[iAlpha]*u[iBeta]; + if (iAlpha==iBeta) { + pi[iPi] -= 1./descriptors::invCs2()*(rho-(T)1); + } + ++iPi; + } + } + } + + /// Computation of all hydrodynamic variables + static void computeAllMomenta(CellBase const& cell, T& rho, T u[DESCRIPTORBASE::d], + T pi[util::TensorVal::n] ) + { + computeRhoU(cell, rho, u); + computeStress(cell, rho, u, pi); + } + + /* + /// Computes squared norm of non-equilibrium part of 2nd momentum + static T computePiNeqNormSqr(CellBase const& cell) + { + T rho, u[DESCRIPTORBASE::d], pi[util::TensorVal::n]; + computeAllMomenta(cell, rho, u, pi); + T PiNeqNormSqr = pi[0]*pi[0] + 2.*pi[1]*pi[1] + pi[2]*pi[2]; + if (util::TensorVal::n == 6) { + PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.*pi[4]*pi[4] +pi[5]*pi[5]; + } + return PiNeqNormSqr; + } + + /// Computes squared norm of forced non-equilibrium part of 2nd momentum + static T computeForcedPiNeqNormSqr(CellBase const& cell) + { + T rho, u[DESCRIPTORBASE::d], pi[util::TensorVal::n]; + computeAllMomenta(cell, rho, u, pi); + //Creation of body force tensor (rho/2.)*(G_alpha*U_beta + U_alpha*G_Beta) + T ForceTensor[util::TensorVal::n]; + int iPi = 0; + for (int Alpha=0; Alpha()][Alpha]*u[Beta] + u[Alpha]*cell[DESCRIPTORBASE::template index()][Beta]); + ++iPi; + } + } + // Creation of second-order moment off-equilibrium tensor + for (int iPi=0; iPi < util::TensorVal::n; ++iPi){ + pi[iPi] += ForceTensor[iPi]; + } + T PiNeqNormSqr = pi[0]*pi[0] + 2.*pi[1]*pi[1] + pi[2]*pi[2]; + if (util::TensorVal::n == 6) { + PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.*pi[4]*pi[4] +pi[5]*pi[5]; + } + return PiNeqNormSqr; + } + */ + + static void modifyVelocity(CellBase& cell, const T newU[DESCRIPTORBASE::d]) + { + T rho, oldU[DESCRIPTORBASE::d]; + computeRhoU(cell, rho, oldU); + const T oldUSqr = util::normSqr(oldU); + const T newUSqr = util::normSqr(newU); + for (int iPop=0; iPop +struct lbExternalHelpers { + /// Add a force term after BGK collision + static void addExternalForce(Cell& cell, const T u[DESCRIPTOR::d], T omega, T amplitude) + { + const T* force = cell.template getFieldPointer(); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + T c_u = T(); + for (int iD=0; iD < DESCRIPTOR::d; ++iD) { + c_u += descriptors::c(iPop,iD)*u[iD]; + } + c_u *= descriptors::invCs2()*descriptors::invCs2(); + T forceTerm = T(); + for (int iD=0; iD < DESCRIPTOR::d; ++iD) { + forceTerm += + ( ((T)descriptors::c(iPop,iD)-u[iD]) * descriptors::invCs2() + + c_u * descriptors::c(iPop,iD) + ) + * force[iD]; + } + forceTerm *= descriptors::t(iPop); + forceTerm *= T(1) - omega/T(2); + forceTerm *= amplitude; + cell[iPop] += forceTerm; + } + } +}; // struct externalFieldHelpers + +/// Helper functions with full-lattice access +template +struct lbLatticeHelpers { + /// Swap ("bounce-back") values of a cell (2D), and apply streaming step + static void swapAndStream2D(Cell **grid, int iX, int iY) + { + const int half = DESCRIPTOR::q/2; + for (int iPop=1; iPop<=half; ++iPop) { + int nextX = iX + descriptors::c(iPop,0); + int nextY = iY + descriptors::c(iPop,1); + T fTmp = grid[iX][iY][iPop]; + grid[iX][iY][iPop] = grid[iX][iY][iPop+half]; + grid[iX][iY][iPop+half] = grid[nextX][nextY][iPop]; + grid[nextX][nextY][iPop] = fTmp; + } + } + + /// Swap ("bounce-back") values of a cell (3D), and apply streaming step + static void swapAndStream3D(Cell ***grid, + int iX, int iY, int iZ) + { + const int half = DESCRIPTOR::q/2; + for (int iPop=1; iPop<=half; ++iPop) { + int nextX = iX + descriptors::c(iPop,0); + int nextY = iY + descriptors::c(iPop,1); + int nextZ = iZ + descriptors::c(iPop,2); + T fTmp = grid[iX][iY][iZ][iPop]; + grid[iX][iY][iZ][iPop] = grid[iX][iY][iZ][iPop+half]; + grid[iX][iY][iZ][iPop+half] = grid[nextX][nextY][nextZ][iPop]; + grid[nextX][nextY][nextZ][iPop] = fTmp; + } + } +}; + +/// All boundary helper functions are inside this structure +template +struct BoundaryHelpers { + static void computeStress ( + Cell const& cell, T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) + { + typedef DESCRIPTOR L; + const T uSqr = util::normSqr(u); + + std::vector const& onWallIndices = util::subIndex(); + std::vector const& normalIndices = util::subIndex(); + + T fNeq[DESCRIPTOR::q]; + for (unsigned fIndex=0; fIndex::equilibrium(iPop, rho, u, uSqr); + } + } + for (unsigned fIndex=0; fIndex::equilibrium(iPop, rho, u, uSqr); + } + + int iPi = 0; + for (int iAlpha=0; iAlpha(iPop,iAlpha)*descriptors::c(iPop,iBeta)*fNeq[iPop]; + } + for (unsigned fIndex=0; fIndex(iPop,iAlpha)*descriptors::c(iPop,iBeta)* + fNeq[iPop]; + } + ++iPi; + } + } + } + +}; // struct boundaryHelpers + +} // namespace olb + +// The specialized code is directly included. That is because we never want +// it to be precompiled so that in both the precompiled and the +// "include-everything" version, the compiler can apply all the +// optimizations it wants. +#include "lbHelpersD2Q5.h" +#include "lbHelpersD2Q9.h" +#include "lbHelpersD3Q7.h" +#include "lbHelpersD3Q15.h" +#include "lbHelpersD3Q19.h" +#include "lbHelpersD3Q27.h" + +#endif -- cgit v1.2.3