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/lbHelpersD2Q5.h | 222 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 222 insertions(+) create mode 100644 src/dynamics/lbHelpersD2Q5.h (limited to 'src/dynamics/lbHelpersD2Q5.h') diff --git a/src/dynamics/lbHelpersD2Q5.h b/src/dynamics/lbHelpersD2Q5.h new file mode 100644 index 0000000..04593d1 --- /dev/null +++ b/src/dynamics/lbHelpersD2Q5.h @@ -0,0 +1,222 @@ +/* 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 + * Template specializations for some computationally intensive LB + * functions of the header file lbHelpers.h, for the D2Q5 grid. + */ + +#ifndef LB_HELPERS_D2Q5_H +#define LB_HELPERS_D2Q5_H + +#include "mrtLatticeDescriptors.h" + +namespace olb { + +// Efficient specialization for D2Q5 lattice +template +struct lbDynamicsHelpers > { + using SpecializedCellBase = CellBase>; + + static T equilibrium( int iPop, T rho, const T u[2], T uSqr ) + { + typedef descriptors::D2Q5<> L; + T c_u = descriptors::c(iPop,0)*u[0] + descriptors::c(iPop,1)*u[1]; + return rho * descriptors::t(iPop) * ( + 1. + 3.*c_u + 4.5*c_u*c_u - 1.5*uSqr ) + - descriptors::t(iPop); + } + + /// equilibrium distribution + static T equilibriumFirstOrder( int iPop, T rho, const T u[2] ) + { + typedef descriptors::D2Q5<> L; + T c_u = descriptors::c(iPop,0) * u[0] + descriptors::c(iPop,1) * u[1]; + + return rho * descriptors::t(iPop) * ( ( T )1 + c_u * descriptors::invCs2() ) - descriptors::t(iPop); + } + + /// RLB advection diffusion collision step + static T rlbCollision( SpecializedCellBase& cell, T rho, const T u[2], T omega ) + { + typedef descriptors::D2Q5<> L; + const T uSqr = u[0] * u[0] + u[1] * u[1]; + + const T Cs2 = ( T )1 / descriptors::invCs2(); + + T rho_1 = rho - ( T )1; + cell[0] = ( ( T )1 - ( T )2 * Cs2 ) * rho_1; //f[0]=(1-2c_s^2)(rho-1) + + const T omega_ = ( T )1 - omega; + + const T f1_3 = ( T )0.5 * omega_ * ( cell[1] - cell[3] ); + const T f2_4 = ( T )0.5 * omega_ * ( cell[2] - cell[4] ); + + rho_1 *= ( T )0.5 * Cs2; + + const T ux_ = ( T )0.5 * omega * rho * u[0]; + + cell[1] = rho_1 + f1_3 - ux_; //f[1]=1/2*(c_s^2(rho-1)+(1-omega)*(f[1]-f[3])-omega*rho*u[x]) + cell[3] = rho_1 - f1_3 + ux_; //f[3]=1/2*(c_s^2(rho-1)-(1-omega)*(f[1]-f[3])+omega*rho*u[x]) + + const T uy_ = ( T )0.5 * omega * rho * u[1]; + + cell[2] = rho_1 + f2_4 - uy_; //f[2]=1/2*(c_s^2(rho-1)+(1-omega)*(f[2]-f[4])-omega*rho*u[y]) + cell[4] = rho_1 - f2_4 + uy_; //f[4]=1/2*(c_s^2(rho-1)-(1-omega)*(f[2]-f[4])+omega*rho*u[y]) + + return uSqr; + } + + // BGK advection diffusion collision step + static T bgkCollision( SpecializedCellBase& cell, T rho, const T u[2], T omega ) + { + typedef descriptors::D2Q5<> L; + + const T Cs2 = ( T )1 / descriptors::invCs2(); + const T uSqr = u[0] * u[0] + u[1] * u[1]; + + const T omega_ = ( T )1 - omega; + const T omega_2 = ( T )0.5 * omega; + T rho_ = ( rho - ( T )1 ); + + cell[0] = omega_ * cell[0] + omega * ( ( T )1 - ( T )2 * Cs2 ) * rho_; + + const T jx = rho * u[0]; + const T jy = rho * u[1]; + + rho_ *= Cs2; + cell[1] = omega_ * cell[1] + omega_2 * ( rho_ - jx ); + cell[2] = omega_ * cell[2] + omega_2 * ( rho_ - jy ); + cell[3] = omega_ * cell[3] + omega_2 * ( rho_ + jx ); + cell[4] = omega_ * cell[4] + omega_2 * ( rho_ + jy ); + + return uSqr; + } + + static void computeMomentaEquilibrium( T momentaEq[5], T rho, const T u[2], T uSqr ) + { + momentaEq[0] = rho; + momentaEq[1] = rho*u[0]; + momentaEq[2] = rho*u[1]; + momentaEq[3] = -rho*(T)2/(T)3; + momentaEq[4] = T(); + } + + static void computeMomenta( T momenta[5], SpecializedCellBase const& cell ) + { + momenta[0] = cell[0] + cell[1] + cell[2] + cell[3] + cell[4]; + momenta[1] = -cell[1] + cell[3]; + momenta[2] = -cell[2] + cell[4]; + momenta[3] = -(T)4*cell[0] + cell[1] + cell[2] + cell[3] + cell[4] - (T)2 / (T)3; + momenta[4] = cell[1] - cell[2] + cell[3] - cell[4]; + } + + static T mrtCollision( SpecializedCellBase& cell, T rho, const T u[2], T invM_S[5][5] ) + { + T uSqr = util::normSqr(u); + T momenta[5]; + T momentaEq[5]; + + computeMomenta(momenta, cell); + computeMomentaEquilibrium(momentaEq, rho, u, uSqr); + + T mom1 = momenta[1] - momentaEq[1]; + T mom2 = momenta[2] - momentaEq[2]; + T mom3 = momenta[3] - momentaEq[3]; + T mom4 = momenta[4] - momentaEq[4]; + + cell[0] -= invM_S[0][1]*mom1 + invM_S[0][2]*mom2 + + invM_S[0][3]*mom3 + invM_S[0][4]*mom4; + + cell[1] -= invM_S[1][1]*mom1 + invM_S[1][2]*mom2 + + invM_S[1][3]*mom3 + invM_S[1][4]*mom4; + + cell[2] -= invM_S[2][1]*mom1 + invM_S[2][2]*mom2 + + invM_S[2][3]*mom3 + invM_S[2][4]*mom4; + + cell[3] -= invM_S[3][1]*mom1 + invM_S[3][2]*mom2 + + invM_S[3][3]*mom3 + invM_S[3][4]*mom4; + + cell[4] -= invM_S[4][1]*mom1 + invM_S[4][2]*mom2 + + invM_S[4][3]*mom3 + invM_S[4][4]*mom4; + + return uSqr; + } + + static void computeFneq ( SpecializedCellBase const& cell, T fNeq[5], T rho, const T u[2] ) + { + for (int iPop=0; iPop < 5; ++iPop) { + // printf("WARNING: First order equilibrium function is used for the calculation of the non-equilibrium parts!\n"); + fNeq[iPop] = cell[iPop] - equilibriumFirstOrder(iPop, rho, u); + } + } + + static T computeRho(SpecializedCellBase const& cell) + { + T rho = cell[0] + cell[1] + cell[2] + cell[3] + cell[4] + (T)1; + return rho; + } + + static void computeRhoU(SpecializedCellBase const& cell, T& rho, T u[2]) + { + rho = cell[0] + cell[1] + cell[2] + cell[3] + cell[4] + (T)1; + T invRho= 1./rho; + u[0] = (cell[3] - cell[1])*invRho; + u[1] = (cell[4] - cell[2])*invRho; + } + + static void computeRhoJ(SpecializedCellBase const& cell, T& rho, T j[2]) + { + rho = cell[0] + cell[1] + cell[2] + cell[3] + cell[4] + (T)1; + + j[0] = (cell[3] - cell[1]); + j[1] = (cell[4] - cell[2]); + } + + static void computeJ(SpecializedCellBase const& cell, T j[2]) + { + j[0] = (cell[3] - cell[1]); + j[1] = (cell[4] - cell[2]); + } + + static void computeStress(SpecializedCellBase const& cell, T rho, const T u[2], T pi[6]) + { + // printf("ERROR: Stress tensor not defined in D2Q5!\n"); + } + + static void computeAllMomenta(SpecializedCellBase const& cell, T& rho, T u[2], T pi[6]) + { + rho = cell[0] + cell[1] + cell[2] + cell[3] + cell[4] + (T)1; + T invRho= 1./rho; + u[0] = (cell[3] - cell[1])*invRho; + u[1] = (cell[4] - cell[2])*invRho; + assert(false); + } + +}; + + +} // namespace olb + +#endif + -- cgit v1.2.3