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