summaryrefslogtreecommitdiff
path: root/src/dynamics/lbHelpersD2Q9.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/dynamics/lbHelpersD2Q9.h')
-rw-r--r--src/dynamics/lbHelpersD2Q9.h343
1 files changed, 343 insertions, 0 deletions
diff --git a/src/dynamics/lbHelpersD2Q9.h b/src/dynamics/lbHelpersD2Q9.h
new file mode 100644
index 0000000..62e6a4f
--- /dev/null
+++ b/src/dynamics/lbHelpersD2Q9.h
@@ -0,0 +1,343 @@
+/* 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
+ * <http://www.openlb.net/>
+ *
+ * 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
+ * Template specializations for some computationally intensive LB
+ * functions of the header file lbHelpers.h, for the D2Q9 grid.
+ */
+
+#ifndef LB_HELPERS_D2Q9_H
+#define LB_HELPERS_D2Q9_H
+
+namespace olb {
+
+// Efficient specialization for D2Q9 base lattice
+template<typename T, typename... FIELDS>
+struct lbDynamicsHelpers<T, descriptors::D2Q9<FIELDS...> > {
+ using SpecializedCellBase = CellBase<T,descriptors::D2Q9<FIELDS...>>;
+
+ using SpecializedDescriptor = descriptors::D2Q9<FIELDS...>;
+
+ static T equilibrium(int iPop, T rho, const T u[2], T uSqr)
+ {
+ typedef descriptors::D2Q9<> L;
+ T c_u = descriptors::c<L>(iPop,0)*u[0] + descriptors::c<L>(iPop,1)*u[1];
+ return rho * descriptors::t<T,L>(iPop) * (
+ 1. + 3.*c_u + 4.5*c_u*c_u - 1.5*uSqr )
+ - descriptors::t<T,L>(iPop);
+ }
+
+ static T equilibriumFirstOrder(int iPop, T rho, const T u[2])
+ {
+ typedef descriptors::D2Q9<> L;
+ T c_u = descriptors::c<L>(iPop,0) * u[0] + descriptors::c<L>(iPop,1) * u[1];
+
+ return rho * descriptors::t<T,L>(iPop) * ( ( T )1 + c_u * descriptors::invCs2<T,L>() ) - descriptors::t<T,L>(iPop);
+ }
+
+ static T incEquilibrium(int iPop, const T j[2], const T jSqr, const T pressure)
+ {
+ typedef descriptors::D2Q9<> L;
+ T c_j = descriptors::c<L>(iPop,0)*j[0] + descriptors::c<L>(iPop,1)*j[1];
+ return descriptors::t<T,L>(iPop) * (
+ 3.*pressure + 3.*c_j + 4.5*c_j*c_j - 1.5*jSqr )
+ - descriptors::t<T,L>(iPop);
+ }
+
+ static void computeFneq(SpecializedCellBase const& cell, T fNeq[9], T rho, const T u[2])
+ {
+ const T uSqr = u[0]*u[0] + u[1]*u[1];
+ for (int iPop=0; iPop < 9; ++iPop) {
+ fNeq[iPop] = cell[iPop] - equilibrium(iPop, rho, u, uSqr);
+ }
+ }
+
+ static T bgkCollision(SpecializedCellBase& cell, T const& rho, const T u[2], T const& omega)
+ {
+ T uxSqr = u[0]*u[0];
+ T uySqr = u[1]*u[1];
+
+ T ux_ = (T)3 * u[0];
+ T uy_ = (T)3 * u[1];
+
+ T uxSqr_ = (T)3 * uxSqr;
+ T uySqr_ = (T)3 * uySqr;
+ T uxSqr__ = (T)3/(T)2 * uxSqr;
+ T uySqr__ = (T)3/(T)2 * uySqr;
+ T uSqr_ = uxSqr__ + uySqr__;
+
+ T uxPySqr_ = (T)9/(T)2 * (u[0]+u[1])*(u[0]+u[1]);
+ T uxMySqr_ = (T)9/(T)2 * (u[0]-u[1])*(u[0]-u[1]);
+
+ T rho_ = (T)4/(T)9 * rho;
+ T cf_ = (T)4/(T)9 * (rho-(T)1);
+
+ cell[0] *= (T)1-omega;
+ cell[0] += omega*(cf_ + rho_*(- uxSqr__ - uySqr__));
+
+ rho_ = (T)1/(T)9 * rho;
+ cf_ = (T)1/(T)9 * (rho-(T)1);
+
+ cell[6] *= (T)1-omega;
+ cell[6] += omega*(cf_ + rho_*(ux_ + uxSqr_ - uySqr__));
+ cell[8] *= (T)1-omega;
+ cell[8] += omega*(cf_ + rho_*(uy_ + uySqr_ - uxSqr__));
+ cell[2] *= (T)1-omega;
+ cell[2] += omega*(cf_ + rho_*(-ux_ + uxSqr_ - uySqr__));
+ cell[4] *= (T)1-omega;
+ cell[4] += omega*(cf_ + rho_*(-uy_ + uySqr_ - uxSqr__));
+
+ rho_ = (T)1/(T)36 * rho;
+ cf_ = (T)1/(T)36 * (rho-(T)1);
+
+ cell[7] *= (T)1-omega;
+ cell[7] += omega*(cf_ + rho_*(ux_ + uy_ + uxPySqr_ - uSqr_));
+ cell[1] *= (T)1-omega;
+ cell[1] += omega*(cf_ + rho_*(-ux_ + uy_ + uxMySqr_ - uSqr_));
+ cell[3] *= (T)1-omega;
+ cell[3] += omega*(cf_ + rho_*(-ux_ - uy_ + uxPySqr_ - uSqr_));
+ cell[5] *= (T)1-omega;
+ cell[5] += omega*(cf_ + rho_*(ux_ - uy_ + uxMySqr_ - uSqr_));
+
+ return uxSqr + uySqr;
+ }
+
+ static T incBgkCollision(SpecializedCellBase& cell, T pressure, const T j[2], T omega)
+ {
+ const T jSqr = util::normSqr<T,descriptors::D2Q9<>::d>(j);
+ for (int iPop=0; iPop < descriptors::D2Q9<>::q; ++iPop) {
+ cell[iPop] *= (T)1-omega;
+ cell[iPop] += omega * lbHelpers<T,SpecializedDescriptor>::incEquilibrium (
+ iPop, j, jSqr, pressure );
+ }
+ return jSqr;
+ }
+
+ static T constRhoBgkCollision(SpecializedCellBase& cell, T rho, const T u[2], T ratioRho, T omega)
+ {
+ const T uSqr = util::normSqr<T,descriptors::D2Q9<>::d>(u);
+ for (int iPop=0; iPop < descriptors::D2Q9<>::q; ++iPop) {
+ T feq = lbHelpers<T,SpecializedDescriptor>::equilibrium(iPop, rho, u, uSqr );
+ cell[iPop] = ratioRho*(feq+descriptors::t<T,SpecializedDescriptor>(iPop))
+ -descriptors::t<T,SpecializedDescriptor>(iPop) +
+ ((T)1-omega)*(cell[iPop]-feq);
+ }
+ return uSqr;
+ }
+
+
+ static void partial_rho(SpecializedCellBase const& cell,
+ T& lineX_P1, T& lineX_0, T& lineX_M1, T& lineY_P1, T& lineY_M1)
+ {
+ lineX_P1 = cell[5] + cell[6] + cell[7];
+ lineX_0 = cell[0] + cell[4] + cell[8];
+ lineX_M1 = cell[1] + cell[2] + cell[3];
+
+ lineY_P1 = cell[7] + cell[8] + cell[1];
+ lineY_M1 = cell[3] + cell[4] + cell[5];
+ }
+
+ static T computeRho(SpecializedCellBase const& cell)
+ {
+ T rho = cell[0] + cell[1] + cell[2] + cell[3] + cell[4]
+ + cell[5] + cell[6] + cell[7] + cell[8] + (T)1;
+ return rho;
+ }
+
+ static void computeRhoU(SpecializedCellBase const& cell, T& rho, T u[2])
+ {
+ T lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1;
+ partial_rho(cell, lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1);
+
+ rho = lineX_P1 + lineX_0 + lineX_M1 + (T)1;
+ T invRho= 1./rho;
+ u[0] = (lineX_P1 - lineX_M1)*invRho;
+ u[1] = (lineY_P1 - lineY_M1)*invRho;
+ }
+
+ static void computeRhoJ(SpecializedCellBase const& cell, T& rho, T j[2])
+ {
+ T lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1;
+ partial_rho(cell, lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1);
+
+ rho = lineX_P1 + lineX_0 + lineX_M1 + (T)1;
+ j[0] = (lineX_P1 - lineX_M1);
+ j[1] = (lineY_P1 - lineY_M1);
+ }
+
+ static void computeJ(SpecializedCellBase const& cell, T j[2] )
+ {
+ T lineX_P1, lineX_M1, lineY_P1, lineY_M1;
+
+ lineX_P1 = cell[5] + cell[6] + cell[7];
+ lineX_M1 = cell[1] + cell[2] + cell[3];
+ lineY_P1 = cell[7] + cell[8] + cell[1];
+ lineY_M1 = cell[3] + cell[4] + cell[5];
+
+ j[0] = (lineX_P1 - lineX_M1);
+ j[1] = (lineY_P1 - lineY_M1);
+ }
+
+ static void computeStress(SpecializedCellBase const& cell, T rho, const T u[2], T pi[3])
+ {
+ typedef descriptors::D2Q9<> L;
+ // Workaround for Intel(r) compiler 9.1;
+ // "using namespace util::tensorIndices2D" is not sufficient
+ using util::tensorIndices2D::xx;
+ using util::tensorIndices2D::yy;
+ using util::tensorIndices2D::xy;
+
+ T lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1;
+ partial_rho(cell, lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1);
+
+ pi[xx] = lineX_P1+lineX_M1 - 1./descriptors::invCs2<T,L>()*(rho-(T)1) - rho*u[0]*u[0];
+ pi[yy] = lineY_P1+lineY_M1 - 1./descriptors::invCs2<T,L>()*(rho-(T)1) - rho*u[1]*u[1];
+ pi[xy] = -cell[1] + cell[3] - cell[5] + cell[7] - rho*u[0]*u[1];
+ }
+
+ static void computeAllMomenta(SpecializedCellBase const& cell, T& rho, T u[2], T pi[3] )
+ {
+ typedef descriptors::D2Q9<> L;
+ // Workaround for Intel(r) compiler 9.1;
+ // "using namespace util::tensorIndices2D" is not sufficient
+ using util::tensorIndices2D::xx;
+ using util::tensorIndices2D::yy;
+ using util::tensorIndices2D::xy;
+
+ T lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1;
+ partial_rho(cell, lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1);
+
+ rho = lineX_P1 + lineX_0 + lineX_M1 + (T)1;
+
+ T rhoU0 = (lineX_P1 - lineX_M1);
+ T rhoU1 = (lineY_P1 - lineY_M1);
+ u[0] = rhoU0/rho;
+ u[1] = rhoU1/rho;
+
+ pi[xx] = lineX_P1 + lineX_M1 - 1./descriptors::invCs2<T,L>()*(rho-(T)1) - rhoU0*u[0];
+ pi[yy] = lineY_P1 + lineY_M1 - 1./descriptors::invCs2<T,L>()*(rho-(T)1) - rhoU1*u[1];
+ pi[xy] = -cell[1] + cell[3] - cell[5] + cell[7] - rhoU0*u[1];
+ }
+
+ static void modifyVelocity(SpecializedCellBase const& cell, const T newU[2])
+ {
+ T rho, oldU[2];
+ computeRhoU(cell, rho, oldU);
+ const T oldUSqr = util::normSqr<T,2>(oldU);
+ const T newUSqr = util::normSqr<T,2>(newU);
+ for (int iPop=0; iPop<9; ++iPop) {
+ cell[iPop] = cell[iPop]
+ - equilibrium(iPop, rho, oldU, oldUSqr)
+ + equilibrium(iPop, rho, newU, newUSqr);
+ }
+ }
+
+}; //struct lbHelpers<D2Q9>
+
+// Efficient specialization for D2Q9 lattice with force
+template<typename T>
+struct lbExternalHelpers<T, descriptors::D2Q9<descriptors::FORCE>> {
+
+ static void addExternalForce(
+ Cell<T,descriptors::D2Q9<descriptors::FORCE>>& cell,
+ const T u[descriptors::D2Q9<descriptors::FORCE>::d], T omega, T amplitude)
+ {
+ const T* force = cell.template getFieldPointer<descriptors::FORCE>();
+ const T mu = amplitude*((T)1-omega/(T)2);
+
+ cell[0] += mu *(T)4/(T)3 *( force[0] * (- u[0] ) +
+ force[1] * ( - u[1] ) );
+ cell[1] += mu *(T)1/(T)12 *( force[0] * ( 2*u[0] - 3*u[1] - 1) +
+ force[1] * (-3*u[0] + 2*u[1] + 1) );
+ cell[2] += mu *(T)1/(T)3 *( force[0] * ( 2*u[0] - 1) +
+ force[1] * ( - u[1] ) );
+ cell[3] += mu *(T)1/(T)12 *( force[0] * ( 2*u[0] + 3*u[1] - 1) +
+ force[1] * ( 3*u[0] + 2*u[1] - 1) );
+ cell[4] += mu *(T)1/(T)3 *( force[0] * (- u[0] ) +
+ force[1] * ( + 2*u[1] - 1) );
+ cell[5] += mu *(T)1/(T)12 *( force[0] * ( 2*u[0] - 3*u[1] + 1) +
+ force[1] * (-3*u[0] + 2*u[1] - 1) );
+ cell[6] += mu *(T)1/(T)3 *( force[0] * ( 2*u[0] + 1) +
+ force[1] * ( - u[1] ) );
+ cell[7] += mu *(T)1/(T)12 *( force[0] * ( 2*u[0] + 3*u[1] + 1) +
+ force[1] * ( 3*u[0] + 2*u[1] + 1) );
+ cell[8] += mu *(T)1/(T)3 *( force[0] * (- u[0] ) +
+ force[1] * ( + 2*u[1] + 1) );
+ }
+};
+
+// Efficient specialization for D2Q9 lattice and for forced D2Q9 lattice
+// (operations applying to the whole lattice)
+
+template<typename T>
+struct lbLatticeHelpers<T, descriptors::D2Q9<>> {
+
+ static void swapAndStreamCell (
+ Cell<T,descriptors::D2Q9<>> **grid,
+ int iX, int iY, int nX, int nY, int iPop, T& fTmp )
+ {
+ fTmp = grid[iX][iY][iPop];
+ grid[iX][iY][iPop] = grid[iX][iY][iPop+4];
+ grid[iX][iY][iPop+4] = grid[nX][nY][iPop];
+ grid[nX][nY][iPop] = fTmp;
+ }
+
+ static void swapAndStream2D (
+ Cell<T,descriptors::D2Q9<>> **grid, int iX, int iY )
+ {
+ T fTmp;
+ swapAndStreamCell(grid, iX, iY, iX-1, iY+1, 1, fTmp);
+ swapAndStreamCell(grid, iX, iY, iX-1, iY, 2, fTmp);
+ swapAndStreamCell(grid, iX, iY, iX-1, iY-1, 3, fTmp);
+ swapAndStreamCell(grid, iX, iY, iX, iY-1, 4, fTmp);
+ }
+};
+
+template<typename T>
+struct lbLatticeHelpers<T, descriptors::D2Q9<descriptors::FORCE>> {
+
+ static void swapAndStreamCell (
+ Cell<T,descriptors::D2Q9<descriptors::FORCE>> **grid,
+ int iX, int iY, int nX, int nY, int iPop, T& fTmp )
+ {
+ fTmp = grid[iX][iY][iPop];
+ grid[iX][iY][iPop] = grid[iX][iY][iPop+4];
+ grid[iX][iY][iPop+4] = grid[nX][nY][iPop];
+ grid[nX][nY][iPop] = fTmp;
+ }
+
+ static void swapAndStream2D (
+ Cell<T,descriptors::D2Q9<descriptors::FORCE>> **grid, int iX, int iY )
+ {
+ T fTmp;
+ swapAndStreamCell(grid, iX, iY, iX-1, iY+1, 1, fTmp);
+ swapAndStreamCell(grid, iX, iY, iX-1, iY, 2, fTmp);
+ swapAndStreamCell(grid, iX, iY, iX-1, iY-1, 3, fTmp);
+ swapAndStreamCell(grid, iX, iY, iX, iY-1, 4, fTmp);
+ }
+
+};
+
+} // namespace olb
+
+#endif