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