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/porousBGKdynamics.hh | 671 ++++++++++++++++++++++++++++++++++++++
1 file changed, 671 insertions(+)
create mode 100644 src/dynamics/porousBGKdynamics.hh
(limited to 'src/dynamics/porousBGKdynamics.hh')
diff --git a/src/dynamics/porousBGKdynamics.hh b/src/dynamics/porousBGKdynamics.hh
new file mode 100644
index 0000000..0da06c1
--- /dev/null
+++ b/src/dynamics/porousBGKdynamics.hh
@@ -0,0 +1,671 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2016 Thomas Henn, Mathias J. Krause, 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
+ * BGK Dynamics for porous -- generic implementation.
+ */
+#ifndef POROUS_BGK_DYNAMICS_HH
+#define POROUS_BGK_DYNAMICS_HH
+
+#include "porousBGKdynamics.h"
+#include "core/cell.h"
+#include "dynamics.h"
+#include "core/util.h"
+#include "lbHelpers.h"
+#include "math.h"
+
+namespace olb {
+
+////////////////////// Class PorousBGKdynamics //////////////////////////
+
+template
+PorousBGKdynamics::PorousBGKdynamics (
+ T omega_, Momenta& momenta_)
+ : BGKdynamics(omega_,momenta_),
+ omega(omega_)
+{ }
+
+template
+void PorousBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T* porosity = cell.template getFieldPointer();
+ for (int i=0; i::bgkCollision(cell, rho, u, omega);
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T PorousBGKdynamics::getOmega() const
+{
+ return omega;
+}
+
+template
+void PorousBGKdynamics::setOmega(T omega_)
+{
+ omega = omega_;
+}
+
+
+//////////////////// Class ExtendedPorousBGKdynamics ////////////////////
+
+template
+ExtendedPorousBGKdynamics::ExtendedPorousBGKdynamics (
+ T omega_, Momenta& momenta_)
+ : BGKdynamics(omega_,momenta_),
+ omega(omega_)
+{
+}
+
+template
+void ExtendedPorousBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T* porosity = cell.template getFieldPointer();
+ T* localVelocity = cell.template getFieldPointer();
+
+ cell.template setField(u);
+
+ for (int i=0; i::bgkCollision(cell, rho, u, omega);
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T ExtendedPorousBGKdynamics::getOmega() const
+{
+ return omega;
+}
+
+template
+void ExtendedPorousBGKdynamics::setOmega(T omega_)
+{
+ omega = omega_;
+}
+
+//////////////////// Class SubgridParticleBGKdynamics ////////////////////
+
+template
+SubgridParticleBGKdynamics::SubgridParticleBGKdynamics (
+ T omega_, Momenta& momenta_)
+ : BGKdynamics(omega_,momenta_),
+ omega(omega_)
+{
+ _fieldTmp[0] = T();
+ _fieldTmp[1] = T();
+ _fieldTmp[2] = T();
+ _fieldTmp[3] = T();
+}
+
+template
+void SubgridParticleBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T* porosity = cell.template getFieldPointer();
+ T* extVelocity = cell.template getFieldPointer();
+// if (porosity[0] != 0) {
+// cout << "extVelocity: " << extVelocity[0] << " " << extVelocity[1] << " " << extVelocity[2] << " " << std::endl;
+// cout << "porosity: " << porosity[0] << std::endl;
+// }
+ for (int i=0; i::bgkCollision(cell, rho, u, omega);
+
+ statistics.incrementStats(rho, uSqr);
+ cell.template setField(0);
+ cell.template setField(0);
+ cell.template setField(0);
+}
+
+template
+T SubgridParticleBGKdynamics::getOmega() const
+{
+ return omega;
+}
+
+template
+void SubgridParticleBGKdynamics::setOmega(T omega_)
+{
+ omega = omega_;
+}
+
+//////////////////// Class PorousParticleBGKdynamics ////////////////////
+
+template
+PorousParticleBGKdynamics::PorousParticleBGKdynamics (
+ T omega_, Momenta& momenta_)
+ : BGKdynamics(omega_,momenta_),
+ omega(omega_)
+{}
+
+template
+void PorousParticleBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T tmpMomentumLoss[DESCRIPTOR::d] = { };
+
+ T* const velNumerator = cell.template getFieldPointer();
+ T* const external = cell.template getFieldPointer(); // TODO: Remove implicit layout requirements in favor of descriptor fields
+ T* const velDenominator = cell.template getFieldPointer();
+
+// use HLBM with Shan-Chen forcing as in previous versions
+#ifdef HLBM_FORCING_SHANCHEN
+ if (*velDenominator > std::numeric_limits::epsilon()) {
+ effectiveVelocity::calculate(velNumerator, u);
+ const T tmp_uSqr = util::normSqr(u);
+ for(int tmp_i=0; tmp_i::equilibrium(tmp_iPop, rho, u, tmp_uSqr)-cell[tmp_iPop]);
+ }
+ }
+ }
+ T uSqr = lbHelpers::bgkCollision(cell, rho, u, omega);
+
+// use Kuperstokh forcing by default
+#else
+ T uSqr = lbHelpers::bgkCollision(cell, rho, u, omega);
+ T uPlus[DESCRIPTOR::d] = { };
+ T diff[DESCRIPTOR::q] = {};
+ if (*velDenominator > std::numeric_limits::epsilon()) {
+ for(int iDim=0; iDim::calculate(external, uPlus);
+ const T uPlusSqr = util::normSqr(uPlus);
+ for(int tmp_iPop=0; tmp_iPop::equilibrium(tmp_iPop, rho, uPlus, uPlusSqr)
+ - lbHelpers::equilibrium(tmp_iPop, rho, u, uSqr);
+ cell[tmp_iPop] += diff[tmp_iPop];
+ for(int iDim=0; iDim(tmp_iPop,iDim) * diff[tmp_iPop];
+ }
+ }
+ // alternate option to calculate the force
+ //for(int iDim=0; iDim
+template
+struct PorousParticleBGKdynamics::effectiveVelocity {
+ static void calculate(T* pExternal, T* pVelocity) {
+ for (int i=0; i() - DESCRIPTOR::q])
+ * (pExternal[DESCRIPTOR::template index() - DESCRIPTOR::q +i] / pExternal[DESCRIPTOR::template index() - DESCRIPTOR::q] - pVelocity[i]);
+ }
+ pExternal[DESCRIPTOR::template index() - DESCRIPTOR::q] = 1.;
+ pExternal[DESCRIPTOR::template index() - DESCRIPTOR::q] = 0.;
+ }
+};
+
+// TODO: Update for meta descriptor
+template
+template
+struct PorousParticleBGKdynamics::effectiveVelocity {
+ static void calculate(T* pExternal, T* pVelocity) {
+ for (int i=0; i() - Lattice::q]) * pVelocity[i];
+ }
+ }
+};
+
+template
+T PorousParticleBGKdynamics::getOmega() const
+{
+ return omega;
+}
+
+template
+void PorousParticleBGKdynamics::setOmega(T omega_)
+{
+ omega = omega_;
+}
+
+//////////////////// Class KrauseHBGKdynamics ////////////////////
+
+template
+KrauseHBGKdynamics::KrauseHBGKdynamics (T omega_,
+ Momenta& momenta_, T smagoConst_, T dx_, T dt_ )
+ : BGKdynamics(omega_,momenta_), omega(omega_), smagoConst(smagoConst_),
+ preFactor(computePreFactor(omega_,smagoConst_) )
+{
+ _fieldTmp[0] = T(1);
+ _fieldTmp[1] = T();
+ _fieldTmp[2] = T();
+ _fieldTmp[3] = T();
+}
+
+template
+void KrauseHBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d];
+ T newOmega[DESCRIPTOR::q];
+ this->_momenta.computeRhoU(cell, rho, u);
+ computeOmega(this->getOmega(), cell, preFactor, rho, u, newOmega);
+
+ T vel_denom = *cell.template getFieldPointer();
+ if (vel_denom > std::numeric_limits::epsilon()) {
+ T porosity = *cell.template getFieldPointer(); // prod(1-smoothInd)
+ T* vel_num = cell.template getFieldPointer();
+ porosity = 1.-porosity; // 1-prod(1-smoothInd)
+ for (int i=0; i::bgkCollision(cell, rho, u, newOmega);
+ statistics.incrementStats(rho, uSqr);
+
+ cell.template setField(_fieldTmp[0]);
+ cell.template setField({_fieldTmp[1], _fieldTmp[2]});
+ cell.template setField(_fieldTmp[3]);
+}
+
+template
+T KrauseHBGKdynamics::getOmega() const
+{
+ return omega;
+}
+
+template
+void KrauseHBGKdynamics::setOmega(T omega)
+{
+ this->setOmega(omega);
+ preFactor = computePreFactor(omega, smagoConst);
+}
+
+template
+T KrauseHBGKdynamics::computePreFactor(T omega, T smagoConst)
+{
+ return (T)smagoConst*smagoConst*descriptors::invCs2()*descriptors::invCs2()*2*sqrt(2);
+}
+
+
+template
+void KrauseHBGKdynamics::computeOmega(T omega0, Cell& cell, T preFactor, T rho,
+ T u[DESCRIPTOR::d], T newOmega[DESCRIPTOR::q])
+{
+ T uSqr = u[0]*u[0];
+ for (int iDim=0; iDim::equilibrium(iPop, rho, u, uSqr));
+ /// Turbulent realaxation time
+ T tau_turb = 0.5*(sqrt(tau_mol*tau_mol+(preFactor*fNeq))-tau_mol);
+ /// Effective realaxation time
+ tau_eff = tau_mol + tau_turb;
+ newOmega[iPop] = 1./tau_eff;
+ }
+}
+
+
+//////////////////// Class ParticlePorousBGKdynamics ////////////////////
+/*
+template
+ParticlePorousBGKdynamics::ParticlePorousBGKdynamics (
+ T omega_, Momenta& momenta_)
+ : BGKdynamics(omega_,momenta_),
+ omega(omega_)
+{ }
+
+template
+void ParticlePorousBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T* porosity = cell.template getFieldPointer();
+ T* localVelocity = cell.template getFieldPointer();
+ for (int i=0; i::bgkCollision(cell, rho, u, omega);
+ statistics.incrementStats(rho, uSqr);
+
+// *cell.template getFieldPointer() = 1;
+// *cell.template getFieldPointer() = 0.;
+// *(cell.template getFieldPointer()+1) = 0.;
+}
+
+template
+T ParticlePorousBGKdynamics::getOmega() const {
+ return omega;
+}
+
+template
+void ParticlePorousBGKdynamics::setOmega(T omega_) {
+ omega = omega_;
+}
+*/
+
+//////////////////// Class SmallParticleBGKdynamics ////////////////////
+
+template
+SmallParticleBGKdynamics::SmallParticleBGKdynamics (
+ T omega_, Momenta& momenta_)
+ : BGKdynamics(omega_,momenta_),
+ omega(omega_)
+{ }
+
+template
+void SmallParticleBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T* porosity = cell.template getFieldPointer();
+ T* localVelocity = cell.template getFieldPointer();
+
+ //cout << porosity[0] << " " << localVelocity[0]<< " " << localVelocity[1]<< " " << localVelocity[2]<::bgkCollision(cell, rho, u, omega);
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T SmallParticleBGKdynamics::getOmega() const
+{
+ return omega;
+}
+
+template
+void SmallParticleBGKdynamics::setOmega(T omega_)
+{
+ omega = omega_;
+}
+
+////////////////////// Class PSMBGKdynamics //////////////////////////
+
+/** \param omega relaxation parameter, related to the dynamic viscosity
+ * \param momenta a Momenta object to know how to compute velocity momenta
+ */
+template
+PSMBGKdynamics::PSMBGKdynamics (
+ T omega_, Momenta& momenta_, int mode_ )
+ : BGKdynamics(omega_, momenta_),
+ omega(omega_), paramA(1. / omega_ - 0.5)
+{
+ mode = (Mode) mode_;
+}
+
+template
+void PSMBGKdynamics::computeU (Cell const& cell, T u[DESCRIPTOR::d] ) const
+{
+ T rho;
+ this->_momenta.computeRhoU(cell, rho, u);
+// T epsilon = 1. - *(cell.template getFieldPointer());
+// // speed up paramB
+// T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
+// // speed up paramC
+// T paramC = (1. - paramB);
+// for (int iVel=0; iVel()[iVel];
+// }
+}
+
+template
+void PSMBGKdynamics::computeRhoU (Cell const& cell, T& rho, T u[DESCRIPTOR::d] ) const
+{
+ this->_momenta.computeRhoU(cell, rho, u);
+// T epsilon = 1. - *(cell.template getFieldPointer());
+// // speed up paramB
+// T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
+// // speed up paramC
+// T paramC = (1. - paramB);
+// for (int iVel=0; iVel()[iVel];
+// }
+}
+
+template
+void PSMBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d], uSqr;
+ T epsilon = 1. - *(cell.template getFieldPointer());
+
+ this->_momenta.computeRhoU(cell, rho, u);
+ // velocity at the boundary
+ T u_s[DESCRIPTOR::d];
+ for (int i = 0; i < DESCRIPTOR::d; i++) {
+ u_s[i] = (cell.template getFieldPointer())[i];
+ }
+
+ if (epsilon < 1e-5) {
+ uSqr = lbHelpers::bgkCollision(cell, rho, u, omega);
+ } else {
+ // speed up paramB
+ T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
+ // speed up paramC
+ T paramC = (1. - paramB);
+
+ T omega_s[DESCRIPTOR::q];
+ T cell_tmp[DESCRIPTOR::q];
+
+ const T uSqr_s = util::normSqr(u_s);
+
+ uSqr = util::normSqr(u);
+ for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+ cell_tmp[iPop] = cell[iPop];
+ switch(mode){
+ case M2: omega_s[iPop] = (lbDynamicsHelpers::equilibrium(iPop, rho, u_s, uSqr_s ) - cell[iPop])
+ + (1 - omega) * (cell[iPop] - lbDynamicsHelpers::equilibrium(iPop, rho, u, uSqr )); break;
+ case M3: omega_s[iPop] = (cell[descriptors::opposite(iPop)] - lbDynamicsHelpers::equilibrium(descriptors::opposite(iPop), rho, u_s, uSqr_s ))
+ - (cell[iPop] - lbDynamicsHelpers::equilibrium(iPop, rho, u_s, uSqr_s ));
+ }
+
+ }
+
+ uSqr = lbHelpers::bgkCollision(cell, rho, u, omega);
+
+ for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+ cell[iPop] = cell_tmp[iPop] + paramC * (cell[iPop] - cell_tmp[iPop]);
+ cell[iPop] += paramB * omega_s[iPop];
+ }
+ for (int iVel=0; iVel(u);
+ }
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T PSMBGKdynamics::getOmega() const
+{
+ return omega;
+}
+
+template
+void PSMBGKdynamics::setOmega(T omega_)
+{
+ paramA = (1. / omega_ - 0.5);
+ omega = omega_;
+
+}
+
+////////////////////// Class ForcedPSMBGKdynamics //////////////////////////
+
+/** \param omega relaxation parameter, related to the dynamic viscosity
+ * \param momenta a Momenta object to know how to compute velocity momenta
+ */
+template
+ForcedPSMBGKdynamics::ForcedPSMBGKdynamics (
+ T omega_, Momenta& momenta_, int mode_ )
+ : ForcedBGKdynamics(omega_, momenta_),
+ omega(omega_), paramA(1. / omega_ - 0.5)
+{
+ mode = (Mode) mode_;
+}
+
+
+template
+void ForcedPSMBGKdynamics::computeU (Cell const& cell, T u[DESCRIPTOR::d] ) const
+{
+ T rho;
+ this->_momenta.computeRhoU(cell, rho, u);
+ T epsilon = 1. - *(cell.template getFieldPointer());
+ // speed up paramB
+ T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
+ // speed up paramC
+ T paramC = (1. - paramB);
+ for (int iVel=0; iVel()[iVel] / (T)2.) +
+ paramB * cell.template getFieldPointer()[iVel];
+ }
+}
+
+template
+void ForcedPSMBGKdynamics::computeRhoU (Cell const& cell, T& rho, T u[DESCRIPTOR::d] ) const
+{
+ this->_momenta.computeRhoU(cell, rho, u);
+ T epsilon = 1. - *(cell.template getFieldPointer());
+ // speed up paramB
+ T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
+ // speed up paramC
+ T paramC = (1. - paramB);
+ for (int iVel=0; iVel()[iVel] / (T)2.) +
+ paramB * cell.template getFieldPointer()[iVel];
+ }
+}
+
+template
+void ForcedPSMBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d], uSqr;
+ T epsilon = 1. - *(cell.template getFieldPointer());
+
+ this->_momenta.computeRhoU(cell, rho, u);
+
+ T* force = cell.template getFieldPointer();
+ for (int iVel=0; iVel())[i];
+ }
+
+ if (epsilon < 1e-5) {
+ uSqr = lbHelpers::bgkCollision(cell, rho, u, omega);
+ lbHelpers::addExternalForce(cell, u, omega, rho);
+ } else {
+ // speed up paramB
+ T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
+ // speed up paramC
+ T paramC = (1. - paramB);
+
+ T omega_s[DESCRIPTOR::q];
+ T cell_tmp[DESCRIPTOR::q];
+
+ const T uSqr_s = util::normSqr(u_s);
+
+ for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+ cell_tmp[iPop] = cell[iPop];
+ switch(mode){
+ case M2: omega_s[iPop] = (lbDynamicsHelpers::equilibrium(iPop, rho, u_s, uSqr_s ) - cell[iPop])
+ + (1 - omega) * (cell[iPop] - lbDynamicsHelpers::equilibrium(iPop, rho, u, uSqr )); break;
+ case M3: omega_s[iPop] = (cell[descriptors::opposite(iPop)] - lbDynamicsHelpers::equilibrium(descriptors::opposite(iPop), rho, u_s, uSqr_s ))
+ - (cell[iPop] - lbDynamicsHelpers::equilibrium(iPop, rho, u_s, uSqr_s ));
+ }
+ }
+
+ uSqr = lbHelpers::bgkCollision(cell, rho, u, omega);
+ lbHelpers::addExternalForce(cell, u, omega, rho);
+
+ for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+ cell[iPop] = cell_tmp[iPop] + paramC * (cell[iPop] - cell_tmp[iPop]);
+ cell[iPop] += paramB * omega_s[iPop];
+ }
+ for (int iVel=0; iVel(u);
+ }
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T ForcedPSMBGKdynamics::getOmega() const
+{
+ return omega;
+}
+
+template
+void ForcedPSMBGKdynamics::setOmega(T omega_)
+{
+ paramA = (1. / omega_ - 0.5);
+ omega = omega_;
+
+}
+
+
+} // olb
+
+#endif
--
cgit v1.2.3