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/dynamics.hh | 2026 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 2026 insertions(+) create mode 100644 src/dynamics/dynamics.hh (limited to 'src/dynamics/dynamics.hh') diff --git a/src/dynamics/dynamics.hh b/src/dynamics/dynamics.hh new file mode 100644 index 0000000..12b77bf --- /dev/null +++ b/src/dynamics/dynamics.hh @@ -0,0 +1,2026 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2006-2015 Jonas Latt, Mathias J. Krause + * 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 + * A collection of dynamics classes (e.g. BGK) with which a Cell object + * can be instantiated -- generic implementation. + */ +#ifndef LB_DYNAMICS_HH +#define LB_DYNAMICS_HH + +#include +#include +#include "dynamics.h" +#include "core/cell.h" +#include "lbHelpers.h" +#include "firstOrderLbHelpers.h" +#include "d3q13Helpers.h" + +namespace olb { + +////////////////////// Class Dynamics //////////////////////// + +template +T Dynamics::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return lbHelpers::equilibrium(iPop, rho, u, uSqr); +} + +template +void Dynamics::iniEquilibrium(Cell& cell, T rho, const T u[DESCRIPTOR::d]) +{ + T uSqr = util::normSqr(u); + for (int iPop=0; iPop +void Dynamics::setBoundaryIntersection(int iPop, T distance) +{ } + +template +bool Dynamics::getBoundaryIntersection(int iPop, T point[DESCRIPTOR::d]) +{ + return 0; +} + +template +void Dynamics::defineRho(int iPop, T rho) +{ } + +template +void Dynamics::defineU(const T u[DESCRIPTOR::d]) +{ } + +template +void Dynamics::defineU(int iPop, const T u[DESCRIPTOR::d]) +{ } + +template +T Dynamics::getVelocityCoefficient(int iPop) +{ + return 0; +} + +////////////////////// Class BasicDynamics //////////////////////// + +template +BasicDynamics::BasicDynamics(Momenta& momenta) + : _momenta(momenta) +{ } + +template +T BasicDynamics::computeRho(Cell const& cell) const +{ + return _momenta.computeRho(cell); +} + +template +void BasicDynamics::computeU ( + Cell const& cell, + T u[DESCRIPTOR::d]) const +{ + _momenta.computeU(cell, u); +} + +template +void BasicDynamics::computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d]) const +{ + _momenta.computeJ(cell, j); +} + +template +void BasicDynamics::computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + _momenta.computeStress(cell, rho, u, pi); +} + +template +void BasicDynamics::computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const +{ + _momenta.computeRhoU(cell, rho, u); +} + +template +void BasicDynamics::computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + this->computeRhoU(cell, rho, u); + this->computeStress(cell, rho, u, pi); +} + +template +void BasicDynamics::defineRho(Cell& cell, T rho) +{ + _momenta.defineRho(cell, rho); +} + +template +void BasicDynamics::defineU ( + Cell& cell, + const T u[DESCRIPTOR::d]) +{ + _momenta.defineU(cell, u); +} + +template +void BasicDynamics::defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) +{ + _momenta.defineRhoU(cell, rho, u); +} + +template +void BasicDynamics::defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) +{ + _momenta.defineAllMomenta(cell, rho, u, pi); +} + + +////////////////////// Class BGKdynamics ////////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template +BGKdynamics::BGKdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta), + _omega(omega) +{ } + +template +void BGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T uSqr = lbHelpers::bgkCollision(cell, rho, u, _omega); + statistics.incrementStats(rho, uSqr); +} + +template +T BGKdynamics::getOmega() const +{ + return _omega; +} + +template +void BGKdynamics::setOmega(T omega) +{ + _omega = omega; +} + + +////////////////////// Class TRTdynamics ////////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template +TRTdynamics::TRTdynamics ( + T omega, Momenta& momenta, T magicParameter ) + : BasicDynamics(momenta), + _omega(omega), _omega2(1/(magicParameter/(1/omega-0.5)+0.5)), _magicParameter(magicParameter) +{ } + +template +void TRTdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + // T uSqr = lbHelpers::bgkCollision(cell, rho, u, _omega); + const T uSqr = util::normSqr(u); + T fPlus[DESCRIPTOR::q], fMinus[DESCRIPTOR::q]; + T fEq[DESCRIPTOR::q], fEqPlus[DESCRIPTOR::q], fEqMinus[DESCRIPTOR::q]; + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fPlus[iPop] = 0.5 * ( cell[iPop] + cell[descriptors::opposite(iPop)] ); + fMinus[iPop] = 0.5 * ( cell[iPop] - cell[descriptors::opposite(iPop)] ); + fEq[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr); + } + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fEqPlus[iPop] = 0.5 * ( fEq[iPop] + fEq[descriptors::opposite(iPop)] ); + fEqMinus[iPop] = 0.5 * ( fEq[iPop] - fEq[descriptors::opposite(iPop)] ); + } + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] -= _omega * (fPlus[iPop] - fEqPlus[iPop]) + _omega2 * (fMinus[iPop] - fEqMinus[iPop]); + } + statistics.incrementStats(rho, uSqr); +} + +template +T TRTdynamics::getOmega() const +{ + return _omega; +} + +template +void TRTdynamics::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class ConstRhoBGKdynamics ////////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template +ConstRhoBGKdynamics::ConstRhoBGKdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta), + _omega(omega) +{ } + +template +void ConstRhoBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + + T deltaRho = (T)1 - (statistics).getAverageRho(); + T ratioRho = (T)1 + deltaRho/rho; + + T uSqr = lbHelpers::constRhoBgkCollision ( + cell, rho, u, ratioRho, _omega ); + statistics.incrementStats(rho+deltaRho, uSqr); +} + +template +T ConstRhoBGKdynamics::getOmega() const +{ + return _omega; +} + +template +void ConstRhoBGKdynamics::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class IncBGKdynamics ////////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template +IncBGKdynamics::IncBGKdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta), _omega(omega) +{ } + +template +void IncBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho = this->_momenta.computeRho(cell); + T p = rho / descriptors::invCs2(); + T j[DESCRIPTOR::d]; + this->_momenta.computeJ(cell, j); + T uSqr = lbHelpers::incBgkCollision(cell, p, j, _omega); + statistics.incrementStats(rho, uSqr); +} + +template +T IncBGKdynamics::getOmega() const +{ + return _omega; +} + +template +void IncBGKdynamics::setOmega(T omega) +{ + _omega = omega; +} + + + +////////////////////// Class RLBdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template +RLBdynamics::RLBdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta), + _omega(omega) +{ } + +template +void RLBdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell, rho, u, pi); + T uSqr = rlbHelpers::rlbCollision(cell, rho, u, pi, _omega); + statistics.incrementStats(rho, uSqr); +} + +template +T RLBdynamics::getOmega() const +{ + return _omega; +} + +template +void RLBdynamics::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class CombinedRLBdynamics ///////////////////////// + +template +CombinedRLBdynamics::CombinedRLBdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta), + _boundaryDynamics(omega, momenta) +{ } + +template +T CombinedRLBdynamics:: +computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return _boundaryDynamics.computeEquilibrium(iPop, rho, u, uSqr); +} + +template +void CombinedRLBdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + typedef DESCRIPTOR L; + + T rho, u[L::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell,rho,u,pi); + + T uSqr = util::normSqr(u); + + for (int iPop = 0; iPop < L::q; ++iPop) { + cell[iPop] = computeEquilibrium(iPop,rho,u,uSqr) + + firstOrderLbHelpers::fromPiToFneq(iPop, pi); + } + + _boundaryDynamics.collide(cell, statistics); +} + +template +T CombinedRLBdynamics::getOmega() const +{ + return _boundaryDynamics.getOmega(); +} + +template +void CombinedRLBdynamics::setOmega(T omega) +{ + _boundaryDynamics.setOmega(omega); +} + + +////////////////////// Class ForcedBGKdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta Momenta object to know how to compute velocity momenta + */ +template +ForcedBGKdynamics::ForcedBGKdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta), _omega(omega) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides() ); +} + +template +void ForcedBGKdynamics::computeU (Cell const& cell, T u[DESCRIPTOR::d] ) const +{ + T rho; + this->_momenta.computeRhoU(cell, rho, u); + const T* force = cell.template getFieldPointer(); + for (int iVel=0; iVel < DESCRIPTOR::d; ++iVel) { + u[iVel] += force[iVel] / (T)2.; + } +} + +template +void ForcedBGKdynamics::computeRhoU (Cell const& cell, T& rho, T u[DESCRIPTOR::d] ) const +{ + this->_momenta.computeRhoU(cell, rho, u); + const T* force = cell.template getFieldPointer(); + for (int iVel=0; iVel +void ForcedBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + const T* force = cell.template getFieldPointer(); + for (int iVel=0; iVel::bgkCollision(cell, rho, u, _omega); + lbHelpers::addExternalForce(cell, u, _omega, rho); + statistics.incrementStats(rho, uSqr); +} + +template +T ForcedBGKdynamics::getOmega() const +{ + return _omega; +} + +template +void ForcedBGKdynamics::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class ForcedKupershtokhBGKdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta Momenta object to know how to compute velocity momenta + */ +template +ForcedKupershtokhBGKdynamics::ForcedKupershtokhBGKdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta), _omega(omega) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides() ); +} + +template +void ForcedKupershtokhBGKdynamics::computeU (Cell const& cell, T u[DESCRIPTOR::d] ) const +{ + T rho; + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel()[iVel] / (T)2.; + } +} + +template +void ForcedKupershtokhBGKdynamics::computeRhoU (Cell const& cell, T& rho, T u[DESCRIPTOR::d] ) const +{ + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel()[iVel] / (T)2.; + } +} + + +template +void ForcedKupershtokhBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d], uPlusDeltaU[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T* force = cell.template getFieldPointer(); + T uSqr = lbHelpers::bgkCollision(cell, rho, u, _omega); + + for (int iVel=0; iVel(uPlusDeltaU); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] += lbHelpers::equilibrium(iPop, rho, uPlusDeltaU, uPlusDeltaUSqr) + - lbHelpers::equilibrium(iPop, rho, u, uSqr); + } + + statistics.incrementStats(rho, uSqr); +} + +template +T ForcedKupershtokhBGKdynamics::getOmega() const +{ + return _omega; +} + +template +void ForcedKupershtokhBGKdynamics::setOmega(T omega) +{ + _omega = omega; +} + + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta Momenta object to know how to compute velocity momenta + * \param sink counterpart of a source term + */ +template +PoissonDynamics::PoissonDynamics ( + T omega, Momenta& momenta , T sink) + : BasicDynamics(momenta), _omega(omega), _sink(sink) +{ +} + +template +T PoissonDynamics::computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const +{ + return descriptors::t(iPop) * rho - descriptors::t(iPop); +} + + +template +T PoissonDynamics::computeRho(Cell const& cell) const +{ + return lbHelpers::computeRho(cell); +} + +template +void PoissonDynamics::computeU(Cell const& cell, T u[DESCRIPTOR::d]) const +{ + for ( int iDim = 0; iDim < DESCRIPTOR::d; iDim++ ) { + u[iDim] = T(); + } +} + +template +void PoissonDynamics::computeJ(Cell const& cell, T j[DESCRIPTOR::d]) const +{ + lbHelpers::computeJ(cell, j); +} + +template +void PoissonDynamics::computeStress(Cell const& cell, T rho, + const T u[DESCRIPTOR::d], T pi[util::TensorVal::n] ) const +{ + for ( int iDim = 0; iDim < DESCRIPTOR::d; iDim++ ) { + pi[iDim] = T(); + } +} + +template +void PoissonDynamics::computeRhoU( Cell const& cell, T& rho, T u[DESCRIPTOR::d]) const +{ + rho = computeRho(cell); + computeU(cell, u); +} + +template +void PoissonDynamics::computeAllMomenta( Cell const& cell, T &rho, + T u[DESCRIPTOR::q], T pi[util::TensorVal::n] ) const +{ + rho = computeRho(cell); + computeU(cell, u); + computeStress(cell, rho, u, pi); +} + + +template +void PoissonDynamics::collide( Cell& cell, LatticeStatistics& statistics ) +{ + T rho = computeRho(cell); + + for ( int iPop = 0; iPop < DESCRIPTOR::q; ++iPop ) { + cell[iPop] = (1 - _omega) * (cell[iPop] + descriptors::t(iPop)) + + _omega * descriptors::t(iPop) * rho + - _sink*(cell[iPop] + descriptors::t(iPop)) - descriptors::t(iPop); + } + +// // add spherical harmonic definition f_i = 3*pi/4 *density + pi/4 *flux \cdot e_i +// T density_post_collision = lbDynamicsHelpers::computeRho(cell); +// T flux_post_collision[DESCRIPTOR::q]; +// lbDynamicsHelpers::computeJ(cell,flux_post_collision); +// T s = flux_post_collision[0]*descriptors::t(0)+ +// flux_post_collision[1]*descriptors::t(1)+ +// flux_post_collision[2]*descriptors::t(2); +// std::cout << "s " << s << std::endl; +// for ( int iPop = 0; iPop < DESCRIPTOR::q; ++iPop ) { +// cell[iPop] = 3*M_PI/4. *density_post_collision +// + M_PI/4. *s; +// } + statistics.incrementStats(rho, 0); +} + +template +T PoissonDynamics::getOmega() const +{ + return _omega; +} + +template +void PoissonDynamics::setOmega(T omega) +{ + _omega = omega; +} + + +template +P1Dynamics::P1Dynamics ( + T omega, Momenta& momenta , T absorption, T scattering) + : BasicDynamics(momenta), _omega(omega), _absorption(absorption), _scattering(scattering) +{ +} + +template +T P1Dynamics::computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const +{ + return descriptors::t(iPop) * rho - descriptors::t(iPop); +} + + +template +T P1Dynamics::computeRho(Cell const& cell) const +{ + return lbHelpers::computeRho(cell); +} + +template +void P1Dynamics::computeU(Cell const& cell, T u[DESCRIPTOR::d]) const +{ + for ( int iDim = 0; iDim < DESCRIPTOR::d; iDim++ ) { + u[iDim] = T(); + } +} + +template +void P1Dynamics::computeJ(Cell const& cell, T j[DESCRIPTOR::d]) const +{ + lbHelpers::computeJ(cell, j); +} + +template +void P1Dynamics::computeStress(Cell const& cell, T rho, + const T u[DESCRIPTOR::d], T pi[util::TensorVal::n] ) const +{ + for ( int iDim = 0; iDim < DESCRIPTOR::d; iDim++ ) { + pi[iDim] = T(); + } +} + +template +void P1Dynamics::computeRhoU( Cell const& cell, T& rho, T u[DESCRIPTOR::d]) const +{ + rho = computeRho(cell); + computeU(cell, u); +} + +template +void P1Dynamics::computeAllMomenta( Cell const& cell, T &rho, + T u[DESCRIPTOR::q], T pi[util::TensorVal::n] ) const +{ + rho = computeRho(cell); + computeU(cell, u); + computeStress(cell, rho, u, pi); +} + + +template +void P1Dynamics::collide( Cell& cell, LatticeStatistics& statistics ) +{ + T rho = computeRho(cell); + + // compute equilibrium f^eq_i eq. (5.32) from DOI:10.1002/9780470177013 + + // ----move to computeEq + // give me a std::array, please + std::array cellShifted; + for (int iPop = 0; iPop (iPop); + } + double moment0 = std::accumulate(cellShifted.begin(),cellShifted.end(), T(0)); + std::array moment1; + moment1.fill( T(0) ); + // sum_j v_j f_j + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + for (int iDim = 0; iDim < DESCRIPTOR::d; ++iDim) { + moment1[iDim] += descriptors::c(iPop,iDim)*cellShifted[iPop]; + } + } + // do scalar product in eq. sum_j v_j f_j * v_i + std::array moment1_v; + moment1_v.fill( T(0) ); + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + // scalar product + for (int iDim = 0; iDim < DESCRIPTOR::d; ++iDim) { + moment1_v[iPop] += moment1[iDim]*descriptors::c(iPop,iDim); + } + } + std::array fEq; + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + fEq[iPop] = descriptors::t(iPop) *(moment0 +moment1_v[iPop]); + //fEq[iPop] = 1./(4*M_PI) *moment0 +3./(4*M_PI) * moment1_v[iPop] -descriptors::t(iPop); + } + + + std::array relaxCell; + for (int iPop = 0; iPop (iPop)*cellShifted[iPop]; + } + + for (int iPop = 0; iPop (iPop); + } + + statistics.incrementStats(rho, 0); +} + +template +T P1Dynamics::getOmega() const +{ + return _omega; +} + +template +void P1Dynamics::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class ResettingForcedBGKdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template +ResettingForcedBGKdynamics::ResettingForcedBGKdynamics ( + T omega, Momenta& momenta ) + : ForcedBGKdynamics(omega, momenta) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides() ); +} + +template +void ResettingForcedBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T* force = cell.template getFieldPointer(); + if ( !util::nearZero(force[0]) || !util::nearZero(force[1]) || !util::nearZero(force[2]) ) // TODO: unnecessary?? + for (int iVel=0; iVel::bgkCollision(cell, rho, u, this->_omega); + lbHelpers::addExternalForce(cell, u, this->_omega, rho); + statistics.incrementStats(rho, uSqr); + + force[0] = _frc[0]; + force[1] = _frc[1]; + force[2] = _frc[2]; +// force[0] = 0.; +// force[1] = 0.; +// force[2] = 0.; +} + +////////////////////// Class ForcedShanChenBGKdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template +ForcedShanChenBGKdynamics::ForcedShanChenBGKdynamics ( + T omega, Momenta& momenta ) + : ForcedBGKdynamics(omega, momenta ) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides() ); +} + +template +void ForcedShanChenBGKdynamics::computeU (Cell const& cell, T u[DESCRIPTOR::d] ) const +{ + T rho; + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel()[iVel] / (T)2.; + } +} + +template +void ForcedShanChenBGKdynamics::computeRhoU (Cell const& cell, T& rho, T u[DESCRIPTOR::d] ) const +{ + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel()[iVel] / (T)2.; + } +} + +template +void ForcedShanChenBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T* force = cell.template getFieldPointer(); + for (int iVel=0; iVelgetOmega(); + } + T uSqr = lbHelpers::bgkCollision(cell, rho, u, this->getOmega() ); + uSqr=0.; + for (int iVel=0; iVelgetOmega(); + uSqr += u[iVel]*u[iVel]; + } + statistics.incrementStats(rho, uSqr); +} + +////////////////////// Class ForcedTRTdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta Momenta object to know how to compute velocity momenta + */ +template +ForcedTRTdynamics::ForcedTRTdynamics ( + T omega, Momenta& momenta, T magicParameter ) + : BasicDynamics(momenta), + _omega(omega), _omega2(1/(magicParameter/(1/omega-0.5)+0.5)), _magicParameter(magicParameter) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides() ); +} + +template +void ForcedTRTdynamics::computeU (Cell const& cell, T u[DESCRIPTOR::d] ) const +{ + T rho; + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel()[iVel] / (T)2.; + } +} + +template +void ForcedTRTdynamics::computeRhoU (Cell const& cell, T& rho, T u[DESCRIPTOR::d] ) const +{ + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel()[iVel] / (T)2.; + } +} + + +template +void ForcedTRTdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T* force = cell.template getFieldPointer(); + for (int iVel=0; iVel(u); + T fPlus[DESCRIPTOR::q], fMinus[DESCRIPTOR::q]; + T fEq[DESCRIPTOR::q], fEqPlus[DESCRIPTOR::q], fEqMinus[DESCRIPTOR::q]; + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fPlus[iPop] = 0.5 * ( cell[iPop] + cell[descriptors::opposite(iPop)] ); + fMinus[iPop] = 0.5 * ( cell[iPop] - cell[descriptors::opposite(iPop)] ); + fEq[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr); + } + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fEqPlus[iPop] = 0.5 * ( fEq[iPop] + fEq[descriptors::opposite(iPop)] ); + fEqMinus[iPop] = 0.5 * ( fEq[iPop] - fEq[descriptors::opposite(iPop)] ); + } + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] -= _omega * (fPlus[iPop] - fEqPlus[iPop]) + _omega2 * (fMinus[iPop] - fEqMinus[iPop]); + } + lbHelpers::addExternalForce(cell, u, _omega, rho); + statistics.incrementStats(rho, uSqr); +} + +template +T ForcedTRTdynamics::getOmega() const +{ + return _omega; +} + +template +void ForcedTRTdynamics::setOmega(T omega) +{ + _omega = omega; +} +////////////////////// Class D3Q13dynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template +D3Q13dynamics::D3Q13dynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta) +{ + setOmega(omega); +} + +template +T D3Q13dynamics::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + // To get at the equilibrium, execute collision with relaxation parameters 1 + Cell tmp; + for (int pop=0; pop(pop); + } + d3q13Helpers::collision(tmp, rho, u, (T)1, (T)1); + return tmp[iPop]; +} + +template +void D3Q13dynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T uSqr = d3q13Helpers::collision ( + cell, rho, u, lambda_nu, lambda_nu_prime ); + statistics.incrementStats(rho, uSqr); +} + +template +T D3Q13dynamics::getOmega() const +{ + return (T)4 / ( (T)3/lambda_nu + (T)1/(T)2 ); +} + +template +void D3Q13dynamics::setOmega(T omega) +{ + lambda_nu = (T)3 / ( (T)4/omega - (T)1/(T)2 ); + lambda_nu_prime = (T)3 / ( (T)2/omega + (T)1/(T)2 ); +} + +////////////////////// Class Momenta ////////////////////////////// + +template +void Momenta::computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const +{ + rho = this->computeRho(cell); + this->computeU(cell, u); + +} + +template +void Momenta::computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + this->computeRhoU(cell, rho, u); + this->computeStress(cell, rho, u, pi); +} + +template +void Momenta::defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) +{ + this->defineRho(cell, rho); + this->defineU(cell, u); + +} + +////////////////////// Class BulkMomenta ////////////////////////// + +template +T BulkMomenta::computeRho(Cell const& cell) const +{ + return lbHelpers::computeRho(cell); +} + +template +void BulkMomenta::computeU(Cell const& cell, T u[DESCRIPTOR::d]) const +{ + T dummyRho; + lbHelpers::computeRhoU(cell, dummyRho, u); +} + +template +void BulkMomenta::computeJ(Cell const& cell, T j[DESCRIPTOR::d]) const +{ + lbHelpers::computeJ(cell, j); +} + +template +void BulkMomenta::computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + lbHelpers::computeStress(cell, rho, u, pi); +} + +template +void BulkMomenta::computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d] ) const +{ + lbHelpers::computeRhoU(cell, rho,u); +} + +template +void BulkMomenta::computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + lbHelpers::computeAllMomenta(cell, rho, u, pi); +} + +template +void BulkMomenta::defineRho(Cell& cell, T rho) +{ + T oldRho, u[DESCRIPTOR::d]; + computeRhoU(cell, oldRho, u); + T uSqr = util::normSqr(u); + T fNeq[DESCRIPTOR::q]; + lbHelpers::computeFneq(cell, fNeq, oldRho, u); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + + fNeq[iPop]; + } +} + +template +void BulkMomenta::defineU ( + Cell& cell, + const T u[DESCRIPTOR::d]) +{ + T rho, oldU[DESCRIPTOR::d]; + computeRhoU(cell, rho, oldU); + T uSqr = util::normSqr(u); + T fNeq[DESCRIPTOR::q]; + lbHelpers::computeFneq(cell, fNeq, rho, oldU); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + + fNeq[iPop]; + } + +} + +template +void BulkMomenta::defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) +{ + T oldRho, oldU[DESCRIPTOR::d]; + computeRhoU(cell, oldRho, oldU); + T uSqr = util::normSqr(u); + T fNeq[DESCRIPTOR::q]; + lbHelpers::computeFneq(cell, fNeq, oldRho, oldU); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + + fNeq[iPop]; + } +} + +template +void BulkMomenta::defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) +{ + T uSqr = util::normSqr(u); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + + firstOrderLbHelpers::fromPiToFneq(iPop, pi); + } +} + +////////////////////// Class ExternalVelocityMomenta ////////////////////////// + +template +T ExternalVelocityMomenta::computeRho(Cell const& cell) const +{ + return lbHelpers::computeRho(cell); +} + +template +void ExternalVelocityMomenta::computeU(Cell const& cell, T u[DESCRIPTOR::d]) const +{ + const T* uExt = cell.template getFieldPointer(); + for (int iD=0; iD < DESCRIPTOR::d; ++iD) { + u[iD] = uExt[iD]; + } +} + +template +void ExternalVelocityMomenta::computeJ(Cell const& cell, T j[DESCRIPTOR::d]) const +{ + const auto rho = computeRho(cell); + const T* uExt = cell.template getFieldPointer(); + for (int iD=0; iD < DESCRIPTOR::d; ++iD) { + j[iD] = uExt[iD] * rho; + } +} + +template +void ExternalVelocityMomenta::computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + lbHelpers::computeStress(cell, rho, u, pi); +} + +template +void ExternalVelocityMomenta::computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d] ) const +{ + rho = computeRho(cell); + computeU(cell,u); +} + +template +void ExternalVelocityMomenta::computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + computeRhoU(cell, rho,u); + computeStress(cell, rho, u, pi); +} + +template +void ExternalVelocityMomenta::defineRho(Cell& cell, T rho) +{ + T oldRho, u[DESCRIPTOR::d]; + computeRhoU(cell, oldRho, u); + T uSqr = util::normSqr(u); + T fNeq[DESCRIPTOR::q]; + lbHelpers::computeFneq(cell, fNeq, oldRho, u); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + + fNeq[iPop]; + } +} + +template +void ExternalVelocityMomenta::defineU ( + Cell& cell, + const T u[DESCRIPTOR::d]) +{ + T* const uExt = cell.template getFieldPointer(); + for (int iD=0; iD +void ExternalVelocityMomenta::defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) +{ + defineRho(cell, rho); + defineU(cell, u); +} + +template +void ExternalVelocityMomenta::defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) +{ + defineU(cell, u); + T uSqr = util::normSqr(u); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + + firstOrderLbHelpers::fromPiToFneq(iPop, pi); + } +} + +////////////////////// Class BounceBack /////////////////////////// + +template +BounceBack::BounceBack() +{ + _rhoFixed=false; +} + +template +BounceBack::BounceBack(T rho) + :_rho(rho) +{ + _rhoFixed=true; +} + +template +void BounceBack::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + // !do not touch element 0! + for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) { + std::swap(cell[iPop], cell[descriptors::opposite(iPop)]); + } +} + +template +T BounceBack::computeRho(Cell const& cell) const +{ + + if (_rhoFixed) { + return _rho; + } + return lbHelpers::computeRho(cell); +} + +template +void BounceBack::computeU ( + Cell const& cell, + T u[DESCRIPTOR::d]) const +{ + for (int iD=0; iD +void BounceBack::computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d]) const +{ + for (int iD=0; iD +void BounceBack::computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + for (int iPi=0; iPi::n; ++iPi) { + pi[iPi] = T();//std::numeric_limits::signaling_NaN(); + } +} + +template +void BounceBack::computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const +{ + rho = computeRho(cell); + computeU(cell, u); +} + +template +void BounceBack::computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + computeRhoU(cell, rho, u); + computeStress(cell, rho, u, pi); +} + +template +void BounceBack::defineRho(Cell& cell, T rho) +{ } + +template +void BounceBack::defineU ( + Cell& cell, + const T u[DESCRIPTOR::d]) +{ } + +template +void BounceBack::defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) +{ } + +template +void BounceBack::defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) +{ } + +template +T BounceBack::getOmega() const +{ + return T();//std::numeric_limits::signaling_NaN(); +} + +template +void BounceBack::setOmega(T omega) +{ } + + +////////////////////// Class BounceBackVelocity /////////////////////////// + +template +BounceBackVelocity::BounceBackVelocity(const T u[DESCRIPTOR::d]) +{ + _rhoFixed=false; + for (int iD=0; iD +BounceBackVelocity::BounceBackVelocity(const T rho, const T u[DESCRIPTOR::d]) + :_rho(rho) +{ + _rhoFixed=true; + for (int iD=0; iD +void BounceBackVelocity::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) { + std::swap(cell[iPop], cell[iPop+DESCRIPTOR::q/2]); + } + for (int iPop=1; iPop < DESCRIPTOR::q; ++iPop) { + for (int iD=0; iD(iPop,iD)*descriptors::t(iPop)*2*descriptors::invCs2(); + } + } +} + +template +T BounceBackVelocity::computeRho(Cell const& cell) const +{ + if (_rhoFixed) { + return _rho; + } + return lbHelpers::computeRho(cell); +} + +template +void BounceBackVelocity::computeU ( + Cell const& cell, + T u[DESCRIPTOR::d]) const +{ + for (int iD=0; iD +void BounceBackVelocity::computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d]) const +{ + for (int iD=0; iD +void BounceBackVelocity::computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + for (int iPi=0; iPi::n; ++iPi) { + pi[iPi] = T();//std::numeric_limits::signaling_NaN(); + } +} + +template +void BounceBackVelocity::computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const +{ + rho = computeRho(cell); + computeU(cell, u); +} + +template +void BounceBackVelocity::computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + computeRhoU(cell, rho, u); + computeStress(cell, rho, u, pi); +} + +template +void BounceBackVelocity::defineRho(Cell& cell, T rho) +{ } + +template +void BounceBackVelocity::defineU ( + Cell& cell, + const T u[DESCRIPTOR::d]) +{ + for (int iD=0; iD +void BounceBackVelocity::defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) +{ + defineRho(cell,rho); + defineU(cell,u); +} + +template +void BounceBackVelocity::defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) +{ } + +template +T BounceBackVelocity::getOmega() const +{ + return T();//std::numeric_limits::signaling_NaN(); +} + +template +void BounceBackVelocity::setOmega(T omega) +{ } + +////////////////////// Class BounceBackAnti /////////////////////////// + +template +BounceBackAnti::BounceBackAnti() +{ + _rhoFixed = false; + _rho = T(1); +} + +template +BounceBackAnti::BounceBackAnti(const T rho) + :_rho(rho) +{ + _rhoFixed = true; +} + +template +void BounceBackAnti::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + /* + for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) { + std::swap(cell[iPop], cell[iPop+DESCRIPTOR::q/2]); + } + for (int iPop=1; iPop < DESCRIPTOR::q; ++iPop) { + if (descriptors::c(iPop,0) == -1) + cell[iPop] = -cell[descriptors::opposite(iPop)] + (computeRho(cell) - T(1))*descriptors::t(iPop)*2; + } + */ + //T rho, u[DESCRIPTOR::d]; + //computeRhoU(cell, rho, u); + //T uSqr = lbHelpers::bgkCollision(cell, rho, u, 1.78571); + //statistics.incrementStats(rho, uSqr); + +} + +template +T BounceBackAnti::computeRho(Cell const& cell) const +{ + + if (_rhoFixed) { + return _rho; + } + return lbHelpers::computeRho(cell); +} + +template +void BounceBackAnti::computeU ( + Cell const& cell, + T u[DESCRIPTOR::d]) const +{ + for (int iD=0; iD +void BounceBackAnti::computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d]) const +{ + computeU(cell, j); + for (int iD=0; iD +void BounceBackAnti::computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + for (int iPi=0; iPi::n; ++iPi) { + pi[iPi] = T();//std::numeric_limits::signaling_NaN(); + } +} + +template +void BounceBackAnti::computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const +{ + rho = computeRho(cell); + computeU(cell, u); +} + +template +void BounceBackAnti::computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + computeRhoU(cell, rho, u); + computeStress(cell, rho, u, pi); +} + +template +void BounceBackAnti::defineRho(Cell& cell, T rho) +{ + _rho = rho; +} + +template +void BounceBackAnti::defineU ( + Cell& cell, + const T u[DESCRIPTOR::d]) +{ +} + +template +void BounceBackAnti::defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) +{ + defineRho(cell,rho); + defineU(cell,u); +} + +template +void BounceBackAnti::defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) +{ } + +template +T BounceBackAnti::getOmega() const +{ + return T();//std::numeric_limits::signaling_NaN(); +} + +template +void BounceBackAnti::setOmega(T omega) +{ } + + +////////////////////// Class partialBounceBack /////////////////////////// + +template +PartialBounceBack::PartialBounceBack(T rf) : _rf(rf) +{ +} + +template +T PartialBounceBack::computeEquilibrium ( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const +{ + return lbHelpers::equilibriumFirstOrder( iPop, rho, u ); +} + +template +void PartialBounceBack::collide( Cell& cell, LatticeStatistics& statistics ) +{ + for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) { + std::swap(cell[iPop], cell[iPop+DESCRIPTOR::q/2]); + } + for (int iPop=1; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = (_rf -1) * (cell[iPop] + descriptors::t(iPop)) - descriptors::t(iPop); + } +} + + + +////////////////////// Class NoDynamics /////////////////////////// + +template +NoDynamics::NoDynamics(T rho) :_rho(rho) +{ +} + +template +T NoDynamics::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return T(); +} + +template +void NoDynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ } + +template +T NoDynamics::computeRho(Cell const& cell) const +{ + return _rho; +} + +template +void NoDynamics::computeU ( + Cell const& cell, + T u[DESCRIPTOR::d]) const +{ + for (int iD=0; iD +void NoDynamics::computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d]) const +{ + for (int iD=0; iD +void NoDynamics::computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + for (int iPi=0; iPi::n; ++iPi) { + pi[iPi] = T();//std::numeric_limits::signaling_NaN(); + } +} + +template +void NoDynamics::computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const +{ + rho = computeRho(cell); + computeU(cell, u); +} + +template +void NoDynamics::computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const +{ + computeRhoU(cell, rho, u); + computeStress(cell, rho, u, pi); +} + +template +void NoDynamics::defineRho(Cell& cell, T rho) +{ } + +template +void NoDynamics::defineU ( + Cell& cell, + const T u[DESCRIPTOR::d]) +{ } + +template +void NoDynamics::defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) +{ } + +template +void NoDynamics::defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) +{ } + +template +T NoDynamics::getOmega() const +{ + return T();//std::numeric_limits::signaling_NaN(); +} + +template +void NoDynamics::setOmega(T omega) +{ } + +////////////////////// Class offDynamics /////////////////////////// + +template +OffDynamics::OffDynamics(const T _location[DESCRIPTOR::d]) +{ + typedef DESCRIPTOR L; + for (int iD = 0; iD < L::d; iD++) { + location[iD] = _location[iD]; + } + for (int iPop = 0; iPop < L::q; iPop++) { + distances[iPop] = -1; + velocityCoefficient[iPop] = 0; + for (int iD = 0; iD < L::d; iD++) { + boundaryIntersection[iPop][iD] = _location[iD]; + _u[iPop][iD] = T(); + } + } + _rho=T(1); +} + +template +OffDynamics::OffDynamics(const T _location[DESCRIPTOR::d], T _distances[DESCRIPTOR::q]) +{ + typedef DESCRIPTOR L; + for (int iD = 0; iD < L::d; iD++) { + location[iD] = _location[iD]; + } + for (int iPop = 0; iPop < L::q; iPop++) { + distances[iPop] = _distances[iPop]; + velocityCoefficient[iPop] = 0; + for (int iD = 0; iD < L::d; iD++) { + boundaryIntersection[iPop][iD] = _location[iD] - _distances[iPop]*descriptors::c(iPop,iD); + _u[iPop][iD] = T(); + } + } + _rho=T(1); +} + +template +T OffDynamics::computeRho(Cell const& cell) const +{ + /*typedef DESCRIPTOR L; + T rhoTmp = T(); + T counter = T(); + int counter2 = int(); + for (int iPop = 0; iPop < L::q; iPop++) { + if (distances[iPop] != -1) { + rhoTmp += (cell[iPop] + descriptors::t(iPop))*descriptors::t(iPop); + counter += descriptors::t(iPop); + counter2++; + } + } + //if (rhoTmp/counter + 1<0.1999) std::cout << rhoTmp/counter2 + 1 <1.001) std::cout << rhoTmp/counter2 + 1 < +void OffDynamics::computeU(Cell const& cell, T u[DESCRIPTOR::d] ) const +{ + typedef DESCRIPTOR L; + for (int iD = 0; iD < L::d; iD++) { + u[iD] = T(); + } + int counter = 0; + for (int iPop = 0; iPop < L::q; iPop++) { + if ( !util::nearZero(distances[iPop]+1) ) { + for (int iD = 0; iD < L::d; iD++) { + u[iD] += _u[iPop][iD]; + } + counter++; + } + } + if (counter!=0) { + for (int iD = 0; iD < L::d; iD++) { + u[iD] /= counter; + } + } + return; +} + +template +void OffDynamics::setBoundaryIntersection(int iPop, T distance) +{ + /// direction points from the fluid node into the solid domain + /// distance is the distance from the fluid node to the solid wall + typedef DESCRIPTOR L; + distances[iPop] = distance; + for (int iD = 0; iD < L::d; iD++) { + boundaryIntersection[iPop][iD] = location[iD] - distance*descriptors::c(iPop,iD); + } +} + +template +bool OffDynamics::getBoundaryIntersection(int iPop, T intersection[DESCRIPTOR::d]) +{ + typedef DESCRIPTOR L; + if ( !util::nearZero(distances[iPop]+1) ) { + for (int iD = 0; iD < L::d; iD++) { + intersection[iD] = boundaryIntersection[iPop][iD]; + } + return true; + } + return false; +} + +template +void OffDynamics::defineRho(Cell& cell, T rho) +{ + _rho=rho; +} + +template +void OffDynamics::defineRho(int iPop, T rho) +{ + _rho=rho; +} + +template +void OffDynamics::defineU ( + Cell& cell, + const T u[DESCRIPTOR::d]) +{ + defineU(u); +} + +template +void OffDynamics::defineU(const T u[DESCRIPTOR::d]) +{ + typedef DESCRIPTOR L; + for (int iPop = 0; iPop < L::q; iPop++) { + if ( !util::nearZero(distances[iPop]+1) ) { + defineU(iPop, u); + } + } +} + +/// Bouzidi velocity boundary condition formulas for the Coefficients: +/** 2* invCs2*weight*(c,u) for dist < 1/2 + * 1/dist*invCs2*weight*(c,u) for dist >= 1/2 + */ + +template +void OffDynamics::defineU( + int iPop, const T u[DESCRIPTOR::d]) +{ + OLB_PRECONDITION(distances[iPop] != -1) + typedef DESCRIPTOR L; + velocityCoefficient[iPop] = 0; + // scalar product of c(iPop) and u + for (int sum = 0; sum < L::d; sum++) { // +/- problem because of first stream than postprocess + velocityCoefficient[iPop] -= descriptors::c(iPop,sum)*u[sum]; + } + // compute summand for boundary condition + velocityCoefficient[iPop] *= 2*descriptors::invCs2() * descriptors::t(iPop); + + for (int iD = 0; iD < L::d; iD++) { + _u[iPop][iD] = u[iD]; + } +} + +template +T OffDynamics::getVelocityCoefficient(int iPop) +{ + return velocityCoefficient[iPop]; +} + +////////////////////// Class ZeroDistributionDynamics /////////////////////////// + +template +ZeroDistributionDynamics::ZeroDistributionDynamics() +{ +} + +template +void ZeroDistributionDynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = -descriptors::t(iPop); + } +} + +template +T ZeroDistributionDynamics::computeRho(Cell const& cell) const +{ + return lbHelpers::computeRho(cell); +} + +/////////////// Singletons ////////////////////////////////// + +namespace instances { + +template +BulkMomenta& getBulkMomenta() +{ + static BulkMomenta bulkMomentaSingleton; + return bulkMomentaSingleton; +} + +template +ExternalVelocityMomenta& getExternalVelocityMomenta() +{ + static ExternalVelocityMomenta externalVelocityMomentaSingleton; + return externalVelocityMomentaSingleton; +} + +template +BounceBack& getBounceBack() +{ + static BounceBack bounceBackSingleton; + return bounceBackSingleton; +} + +template +PartialBounceBack& getPartialBounceBack(const double rf) +{ + static PartialBounceBack partialBounceBackSingleton(rf); + return partialBounceBackSingleton; +} + + +template +BounceBackVelocity& getBounceBackVelocity(const double rho, const double u[DESCRIPTOR::d]) +{ + static BounceBackVelocity bounceBackSingleton(rho,u); + return bounceBackSingleton; +} + +template +BounceBackAnti& getBounceBackAnti(const double rho) +{ + static BounceBackAnti bounceBackSingleton(rho); + return bounceBackSingleton; +} + +template +NoDynamics& getNoDynamics(T rho) +{ + static NoDynamics noDynamicsSingleton(rho); + return noDynamicsSingleton; +} + +template +ZeroDistributionDynamics& getZeroDistributionDynamics() +{ + static ZeroDistributionDynamics zeroSingleton; + return zeroSingleton; +} + +} + +} + +#endif -- cgit v1.2.3