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/smagorinskyBGKdynamics.hh | 1361 ++++++++++++++++++++++++++++++++
1 file changed, 1361 insertions(+)
create mode 100644 src/dynamics/smagorinskyBGKdynamics.hh
(limited to 'src/dynamics/smagorinskyBGKdynamics.hh')
diff --git a/src/dynamics/smagorinskyBGKdynamics.hh b/src/dynamics/smagorinskyBGKdynamics.hh
new file mode 100644
index 0000000..4814fa4
--- /dev/null
+++ b/src/dynamics/smagorinskyBGKdynamics.hh
@@ -0,0 +1,1361 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2012-2015 Mathias J. Krause, Jonas Latt, Patrick Nathen
+ * 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 with adjusted omega -- generic implementation.
+ */
+#ifndef SMAGORINSKY_BGK_DYNAMICS_HH
+#define SMAGORINSKY_BGK_DYNAMICS_HH
+
+#include "smagorinskyBGKdynamics.h"
+#include "core/cell.h"
+#include "core/util.h"
+#include "lbHelpers.h"
+#include "math.h"
+
+#include // For shear kalman Smagorinsky - Populations
+
+
+namespace olb {
+
+/// Smagorinsky Dynamics
+template
+SmagorinskyDynamics::SmagorinskyDynamics(T smagoConst_)
+ : smagoConst(smagoConst_), preFactor(computePreFactor())
+{ }
+
+template
+T SmagorinskyDynamics::computePreFactor()
+{
+ return (T)smagoConst*smagoConst*descriptors::invCs2()*descriptors::invCs2()*2*sqrt(2);
+}
+
+template
+T SmagorinskyDynamics::getPreFactor()
+{
+ return preFactor;
+}
+
+template
+T SmagorinskyDynamics::getSmagoConst()
+{
+ return smagoConst;
+}
+
+///////////////////////// ADM BGK /////////////////////////////
+
+/*template
+ADMBGKdynamics::ADMBGKdynamics(T omega_, Momenta& momenta_ )
+ : BGKdynamics(omega_,momenta_), omega(omega_)
+{ }
+
+template
+void ADMBGKdynamics::collide(Cell& cell, LatticeStatistics& statistics )
+{
+ T uSqr = lbHelpers::bgkCollision(cell, *cell., cell[velocityBeginsAt], omega);
+ statistics.incrementStats(*cell[rhoIsAt], uSqr);
+}*/
+
+///////////////////////// ForcedADM BGK /////////////////////////////
+
+template
+ForcedADMBGKdynamics::ForcedADMBGKdynamics (
+ T omega_, Momenta& momenta_ )
+ : BGKdynamics(omega_,momenta_),
+ omega(omega_)
+{ }
+
+template
+void ForcedADMBGKdynamics::collide (
+ Cell& cell,
+ LatticeStatistics& statistics )
+{
+ OstreamManager clout(std::cout,"Forced ADM collide:");
+ T rho, u[DESCRIPTOR::d], utst[DESCRIPTOR::d];
+
+// this->momenta.computeAllMomenta(cell, rho, utst, pi);
+
+ T* rho_fil = cell.template getFieldPointer();
+ T* u_filX = cell.template getFieldPointer();
+ T* u_filY = cell.template getFieldPointer();
+ T* u_filZ = cell.template getFieldPointer();
+
+ u[0] = *u_filX;/// *rho_fil;
+ u[1] = *u_filY;/// *rho_fil;
+ u[2] = *u_filZ;/// *rho_fil;
+
+ T* force = cell.template getFieldPointer();
+ for (int iVel=0; iVel::bgkCollision(cell, *rho_fil, u, omega);
+
+ lbHelpers::addExternalForce(cell, u, omega);
+
+ statistics.incrementStats(rho, uSqr);
+}
+
+///////////////////// Class ConSmagorinskyBGKdynamics //////////////////////////
+
+/////////////////// Consistent Smagorinsky BGK --> Malaspinas/Sagaut //////////////
+
+template
+ConSmagorinskyBGKdynamics::ConSmagorinskyBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyBGKdynamics(omega_, momenta_, smagoConst_)
+{ }
+
+template
+T ConSmagorinskyBGKdynamics::computeEffectiveOmega (Cell& cell)
+{
+
+ T H[util::TensorVal::n];
+ T conSmagoR[DESCRIPTOR::q];
+ T S[util::TensorVal::n];
+ T tau_mol = 1./this->getOmega();
+ T cs2 = 1./descriptors::invCs2();
+ T smagoConst = this->getSmagoConst();
+
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(2*PiNeqNormSqr);
+
+ //Strain rate Tensor
+ if (PiNeqNorm != 0) {
+ T Phi = (-0.5*(-rho*tau_mol*cs2+sqrt(rho*rho*tau_mol*tau_mol*cs2*cs2+2.0*(smagoConst*smagoConst)*rho*PiNeqNorm))/(smagoConst*smagoConst*rho*PiNeqNorm));
+ for (int n = 0; n < util::TensorVal::n; ++n) {
+ S[n] = Phi*pi[n];
+ }
+ } else {
+ for (int n = 0; n < util::TensorVal::n; ++n) {
+ S[n] = 0;
+ }
+ }
+
+ //Strain rate Tensor Norm
+ T SNormSqr = S[0]*S[0] + 2.0*S[1]*S[1] + S[2]*S[2];
+ if (util::TensorVal::n == 6) {
+ SNormSqr += S[2]*S[2] + S[3]*S[3] + 2.0*S[4]*S[4] + S[5]*S[5];
+ }
+ T SNorm = sqrt(2*SNormSqr);
+
+ //consistent Samagorinsky additional R term
+ for (int q = 0; q < DESCRIPTOR::q; ++q) {
+ T t = descriptors::t(q); //lattice weights
+
+ //Hermite-Polynom H = c*c-cs^2*kronDelta
+ H[0] = descriptors::c(q,0)*descriptors::c(q,0)-cs2;
+ H[1] = descriptors::c(q,0)*descriptors::c(q,1);
+ H[2] = descriptors::c(q,1)*descriptors::c(q,1)-cs2;//2D
+ if (util::TensorVal::n == 6) {
+ H[2] = descriptors::c(q,0)*descriptors::c(q,2);//3D
+ H[3] = descriptors::c(q,1)*descriptors::c(q,1)-cs2;
+ H[4] = descriptors::c(q,1)*descriptors::c(q,2);
+ H[5] = descriptors::c(q,2)*descriptors::c(q,2)-cs2;
+ }
+
+ //contraction or scalar product H*S
+ T contractHS = H[0]*S[0] + 2.0*H[1]*S[1] + H[2]*S[2];
+ if (util::TensorVal::n == 6) {
+ contractHS += H[2]*S[2] + H[3]*S[3] + 2.0*H[4]*S[4] + H[5]*S[5];
+ }
+
+ //additional term
+ conSmagoR[q] = t*this->getPreFactor()*SNorm*contractHS;
+ }
+ return conSmagoR[0];
+}
+
+/////////////////// Consistent Strain Smagorinsky BGK --> Malaspinas/Sagaut //////////////
+
+template
+ConStrainSmagorinskyBGKdynamics::ConStrainSmagorinskyBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyBGKdynamics(omega_, momenta_, smagoConst_)
+{ }
+
+template
+T ConStrainSmagorinskyBGKdynamics::computeEffectiveOmega(Cell& cell)
+{
+ T S[util::TensorVal::n];
+ T cs2 = 1./descriptors::invCs2();
+ T tau_mol = 1./this->getOmega();
+ T smagoConst_ = this->getSmagoConst();
+
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.0*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(2*PiNeqNormSqr);
+
+
+ //Strain Tensor
+ if ( !util::nearZero(PiNeqNorm) ) {
+ T Phi = (-0.5*(-rho*tau_mol*cs2+sqrt(rho*rho*tau_mol*tau_mol*cs2*cs2+2.0*(smagoConst_*smagoConst_)*rho*PiNeqNorm))/(smagoConst_*smagoConst_*rho*PiNeqNorm));
+ for (int n = 0; n < util::TensorVal::n; ++n) {
+ S[n] = Phi*pi[n];
+ }
+ } else {
+ for (int n = 0; n < util::TensorVal::n; ++n) {
+ S[n] = 0;
+ }
+ }
+
+ //Strain Tensor Norm
+ T SNormSqr = S[0]*S[0] + 2.0*S[1]*S[1] + S[2]*S[2];
+ if (util::TensorVal::n == 6) {
+ SNormSqr += S[2]*S[2] + S[3]*S[3] + 2.0*S[4]*S[4] + S[5]*S[5];
+ }
+ T SNorm = sqrt(2*SNormSqr);
+
+ /// Turbulent realaxation time
+ T tau_turb = pow(smagoConst_,2)*SNorm/cs2;
+ /// Effective realaxation time
+ T tau_eff = tau_mol+tau_turb;
+ T omega_new= 1./tau_eff;
+ return omega_new;
+}
+
+///////////////////////// DYNAMIC SMAGO BGK /////////////////////////////
+template
+DynSmagorinskyBGKdynamics::DynSmagorinskyBGKdynamics (
+ T omega_, Momenta& momenta_)
+ : SmagorinskyBGKdynamics(omega_, momenta_, T(0.))
+{ }
+
+template
+T DynSmagorinskyBGKdynamics::computeEffectiveOmega(Cell& cell)
+{
+ // computation of the relaxation time
+ T v_t = 0;
+ T* dynSmago = cell.template getFieldPointer();
+
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(PiNeqNormSqr);
+ //v_t = *dynSmago*dx*dx*PiNeqNorm;
+ v_t = *dynSmago*PiNeqNorm;
+ T tau_t = 3*v_t;
+ T tau_0 = 1/this->getOmega();
+ T omega_new = 1/(tau_t+tau_0);
+ return omega_new;
+}
+
+///////////////////SHEAR IMPROVED SMAGORINSKY//////////////////////////
+template
+ShearSmagorinskyBGKdynamics::ShearSmagorinskyBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyBGKdynamics(omega_, momenta_, smagoConst_)
+{ }
+
+template
+void ShearSmagorinskyBGKdynamics::collide(Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T newOmega = computeEffectiveOmega(cell,statistics.getTime());
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T uSqr = lbHelpers::bgkCollision(cell, rho, u, newOmega);
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T ShearSmagorinskyBGKdynamics::getEffectiveOmega(Cell& cell)
+{
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ // computation of the relaxation time
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(PiNeqNormSqr);
+
+ T* avShear = cell.template getFieldPointer();
+ T tau_0 = 1./this->getOmega();
+ T PiNeqNorm_SISM = PiNeqNorm - *avShear;
+ T tau_t = 0.5*(sqrt(tau_0*tau_0+(this->getPreFactor()*PiNeqNorm_SISM/rho))-tau_0);
+
+ T omega_new = 1./(tau_t+tau_0);
+
+ return omega_new;
+}
+
+template
+T ShearSmagorinskyBGKdynamics::computeEffectiveOmega(Cell& cell, int iT)
+{
+ OstreamManager clout(std::cout,"shearImprovedCollide");
+
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ // computation of the relaxation time
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(PiNeqNormSqr);
+
+ T* avShear = cell.template getFieldPointer();
+ *avShear = (*avShear*iT + PiNeqNorm)/(iT+1);
+ T tau_0 = 1./this->getOmega();
+ T PiNeqNorm_SISM = PiNeqNorm - *avShear;
+ T tau_t = 0.5*(sqrt(tau_0*tau_0+(this->getPreFactor()*PiNeqNorm_SISM/rho))-tau_0);
+
+ T omega_new = 1./(tau_t+tau_0);
+
+ return omega_new;
+}
+
+///////////////////////// FORCED SHEAR SMAGO BGK /////////////////////////////
+template
+ShearSmagorinskyForcedBGKdynamics::ShearSmagorinskyForcedBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyForcedBGKdynamics(omega_, momenta_, smagoConst_)
+{ }
+
+template
+void ShearSmagorinskyForcedBGKdynamics::collide(Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T newOmega = computeEffectiveOmega(cell, statistics.getTime());
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T* force = cell.template getFieldPointer();
+ for (int iVel=0; iVel::bgkCollision(cell, rho, u, newOmega);
+ lbHelpers::addExternalForce(cell, u, newOmega, rho);
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T ShearSmagorinskyForcedBGKdynamics::getEffectiveOmega(Cell& cell)
+{
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ // computation of the relaxation time
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(PiNeqNormSqr);
+
+ T* avShear = cell.template getFieldPointer();
+ T tau_0 = 1./this->getOmega();
+ T PiNeqNorm_SISM = PiNeqNorm - *avShear;
+ T tau_t = 0.5*(sqrt(tau_0*tau_0+(this->getPreFactor()*PiNeqNorm_SISM/rho))-tau_0);
+
+ T omega_new = 1./(tau_t+tau_0);
+
+ return omega_new;
+}
+
+template
+T ShearSmagorinskyForcedBGKdynamics::computeEffectiveOmega(Cell& cell, int iT)
+{
+ OstreamManager clout(std::cout,"shearImprovedCollide");
+
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ //Creation of body force tensor (rho/2.)*(G_alpha*U_beta + U_alpha*G_Beta)
+ T ForceTensor[util::TensorVal::n];
+ int iPi = 0;
+ for (int Alpha=0; Alpha()[Alpha]*u[Beta] + u[Alpha]*cell.template getFieldPointer()[Beta]);
+ ++iPi;
+ }
+ }
+ // Creation of second-order moment off-equilibrium tensor
+ for (int iPi=0; iPi < util::TensorVal::n; ++iPi){
+ pi[iPi] += ForceTensor[iPi];
+ }
+ // computation of the relaxation time
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(PiNeqNormSqr);
+
+ T* avShear = cell.template getFieldPointer();
+ *avShear = (*avShear*iT+PiNeqNorm)/(iT+1);
+
+ T tau_0 = 1./this->getOmega();
+ T PiNeqNorm_SISM = PiNeqNorm - *avShear;
+ T tau_t = 0.5*(sqrt(tau_0*tau_0+(this->getPreFactor()*PiNeqNorm_SISM/rho))-tau_0);
+
+ T omega_new = 1./(tau_t+tau_0);
+
+ return omega_new;
+}
+
+////////////////////// Class SmagorinskyBGKdynamics //////////////////////////
+/** \param vs2_ speed of sound
+ * \param momenta_ a Momenta object to know how to compute velocity momenta
+ * \param momenta_ a Momenta object to know how to compute velocity momenta
+ */
+template
+SmagorinskyBGKdynamics::SmagorinskyBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyDynamics(smagoConst_),
+ BGKdynamics(omega_,momenta_)
+{ }
+
+template
+void SmagorinskyBGKdynamics::collide(Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T newOmega = computeEffectiveOmega(cell);
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T uSqr = lbHelpers::bgkCollision(cell, rho, u, newOmega);
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T SmagorinskyBGKdynamics::getEffectiveOmega(Cell& cell)
+{
+ T newOmega = computeEffectiveOmega(cell);
+ return newOmega;
+}
+
+template
+T SmagorinskyBGKdynamics::computeEffectiveOmega(Cell& cell)
+{
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(PiNeqNormSqr);
+ /// Molecular realaxation time
+ T tau_mol = 1. /this->getOmega();
+ /// Turbulent realaxation time
+ T tau_turb = 0.5*(sqrt(tau_mol*tau_mol + this->getPreFactor()/rho*PiNeqNorm) - tau_mol);
+ /// Effective realaxation time
+ T tau_eff = tau_mol+tau_turb;
+ T omega_new= 1./tau_eff;
+ return omega_new;
+
+}
+
+///////////////////////// FORCED SMAGO BGK /////////////////////////////
+template
+SmagorinskyForcedBGKdynamics::SmagorinskyForcedBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyDynamics(smagoConst_),
+ ForcedBGKdynamics(omega_,momenta_)
+{ }
+
+template
+void SmagorinskyForcedBGKdynamics::collide(Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T newOmega = computeEffectiveOmega(cell);
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T* force = cell.template getFieldPointer();
+ for (int iVel=0; iVel::bgkCollision(cell, rho, u, newOmega);
+ lbHelpers::addExternalForce(cell, u, newOmega, rho);
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T SmagorinskyForcedBGKdynamics::getEffectiveOmega(Cell& cell)
+{
+ T newOmega = computeEffectiveOmega(cell);
+ return newOmega;
+}
+
+template
+T SmagorinskyForcedBGKdynamics::computeEffectiveOmega(Cell& cell)
+{
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ //Creation of body force tensor (rho/2.)*(G_alpha*U_beta + U_alpha*G_Beta)
+ T ForceTensor[util::TensorVal::n];
+ int iPi = 0;
+ for (int Alpha=0; Alpha()[Alpha]*u[Beta] + u[Alpha]*cell.template getFieldPointer()[Beta]);
+ ++iPi;
+ }
+ }
+ // Creation of second-order moment off-equilibrium tensor
+ for (int iPi=0; iPi < util::TensorVal::n; ++iPi){
+ pi[iPi] += ForceTensor[iPi];
+ }
+ T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
+ if (util::TensorVal::n == 6) {
+ PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
+ }
+ T PiNeqNorm = sqrt(PiNeqNormSqr);
+ /// Molecular realaxation time
+ T tau_mol = 1. /this->getOmega();
+ /// Turbulent realaxation time
+ T tau_turb = 0.5*(sqrt(tau_mol*tau_mol + this->getPreFactor()/rho*PiNeqNorm) - tau_mol);
+ /// Effective realaxation time
+ T tau_eff = tau_mol+tau_turb;
+ T omega_new= 1./tau_eff;
+ return omega_new;
+}
+
+///////////////////////// External TAU EFF LES BGK /////////////////////////////
+template
+ExternalTauEffLESBGKdynamics::ExternalTauEffLESBGKdynamics(T omega_, Momenta& momenta_, T smagoConst_)
+ : SmagorinskyBGKdynamics(omega_, momenta_, smagoConst_)
+{ }
+
+template
+void ExternalTauEffLESBGKdynamics::collide(Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T newTau = *(cell.template getFieldPointer());
+ T newOmega = 1./newTau;
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T uSqr = lbHelpers::bgkCollision(cell, rho, u, newOmega);
+ statistics.incrementStats(rho, uSqr);
+}
+
+///////////////////////// External TAU EFF FORCED LES BGK /////////////////////////////
+template
+ExternalTauEffLESForcedBGKdynamics::ExternalTauEffLESForcedBGKdynamics(T omega_, Momenta& momenta_, T smagoConst_)
+ : SmagorinskyForcedBGKdynamics(omega_, momenta_, smagoConst_)
+{ }
+
+template
+void ExternalTauEffLESForcedBGKdynamics::collide(Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T newTau = *(cell.template getFieldPointer());
+ T newOmega = 1./newTau;
+ T rho, u[DESCRIPTOR::d];
+ this->_momenta.computeRhoU(cell, rho, u);
+ T* force = cell.template getFieldPointer();
+ for (int iVel=0; iVel::bgkCollision(cell, rho, u, newOmega);
+ lbHelpers::addExternalForce(cell, u, newOmega, rho);
+ statistics.incrementStats(rho, uSqr);
+}
+
+///////////////////////// FORCED Linear Velocity SMAGO BGK /////////////////////////////
+template
+SmagorinskyLinearVelocityForcedBGKdynamics::SmagorinskyLinearVelocityForcedBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyForcedBGKdynamics(omega_, momenta_, smagoConst_)
+{ }
+
+template
+void SmagorinskyLinearVelocityForcedBGKdynamics::collide(Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_momenta.computeAllMomenta(cell, rho, u, pi);
+ T newOmega = computeEffectiveOmega(cell);
+ T* force = cell.template getFieldPointer();
+ int nDim = DESCRIPTOR::d;
+ T forceSave[nDim];
+ // adds a+Bu to force, where
+ // d=2: a1=v[0], a2=v[1], B11=v[2], B12=v[3], B21=v[4], B22=v[5]
+ // d=2: a1=v[0], a2=v[1], a3=v[2], B11=v[3], B12=v[4], B13=v[5], B21=v[6], B22=v[7], B23=v[8], B31=v[9], B32=v[10], B33=v[11]
+ T* v = cell.template getFieldPointer();
+ for (int iDim=0; iDim::bgkCollision(cell, rho, u, newOmega);
+ lbHelpers::addExternalForce(cell, u, newOmega, rho);
+ statistics.incrementStats(rho, uSqr);
+ // Writing back to froce fector
+ for (int iVel=0; iVel
+KrauseBGKdynamics::KrauseBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyBGKdynamics(omega_, momenta_, smagoConst_),
+ preFactor(computePreFactor() )
+{ }
+
+template
+void KrauseBGKdynamics::collide(Cell& cell,
+ LatticeStatistics& statistics )
+{
+ T rho, u[DESCRIPTOR::d];
+ T newOmega[DESCRIPTOR::q];
+ this->_momenta.computeRhoU(cell, rho, u);
+ computeEffectiveOmega(this->getOmega(), cell, preFactor, rho, u, newOmega);
+ T uSqr = lbHelpers::bgkCollision(cell, rho, u, newOmega);
+ statistics.incrementStats(rho, uSqr);
+}
+
+template
+T KrauseBGKdynamics::getEffectiveOmega(Cell& cell)
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::n];
+ T newOmega[DESCRIPTOR::q];
+ this->_momenta.computeAllMomenta(cell, rho, uTemp, pi);
+ computeEffectiveOmega(this->getOmega(), cell, preFactor, rho, uTemp, newOmega);
+ T newOmega_average = 0.;
+ for (int iPop=0; iPop
+T KrauseBGKdynamics::computePreFactor()
+{
+ return (T)this->getSmagoConst()*this->getSmagoConst()*3*descriptors::invCs2()*descriptors::invCs2()*2*sqrt(2);
+}
+
+template
+void KrauseBGKdynamics::computeEffectiveOmega(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_/rho*fNeq) - tau_mol);
+ /// Effective realaxation time
+ T tau_eff = tau_mol + tau_turb;
+ newOmega[iPop] = 1./tau_eff;
+ }
+}
+
+////////////////////// Class WALEBGKdynamics //////////////////////////
+/** \param vs2_ speed of sound
+ * \param momenta_ a Momenta object to know how to compute velocity momenta
+ * \param momenta_ a Momenta object to know how to compute velocity momenta
+ */
+template
+WALEBGKdynamics::WALEBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyBGKdynamics(omega_, momenta_, smagoConst_)
+{
+ this->preFactor = this->getSmagoConst()*this->getSmagoConst();
+}
+
+template
+T WALEBGKdynamics::computePreFactor()
+{
+ return (T)this->getSmagoConst()*this->getSmagoConst();
+}
+
+template
+T WALEBGKdynamics::computeEffectiveOmega(Cell& cell_)
+{
+ // velocity gradient tensor
+ T g[3][3];
+ for ( int i = 0; i < 3; i++) {
+ for ( int j = 0; j < 3; j++) {
+ g[i][j] = *(cell_.template getFieldPointer()+(i*3 + j));
+ }
+ }
+ // strain rate tensor
+ T s[3][3];
+ for ( int i = 0; i < 3; i++) {
+ for ( int j = 0; j < 3; j++) {
+ s[i][j] = (g[i][j] + g[j][i]) / 2.;
+ }
+ }
+ // traceless symmetric part of the square of the velocity gradient tensor
+ T G[3][3];
+ for ( int i = 0; i < 3; i++) {
+ for ( int j = 0; j < 3; j++) {
+ G[i][j] = 0.;
+ }
+ }
+
+ for ( int i = 0; i < 3; i++) {
+ for ( int j = 0; j < 3; j++) {
+ for ( int k = 0; k < 3; k++) {
+ G[i][j] += (g[i][k]*g[k][j] + g[j][k]*g[k][i]) / 2.; // The change
+ }
+ }
+ }
+
+ T trace = 0.;
+ for ( int i = 0; i < 3; i++) {
+ trace += (1./3.) * g[i][i] * g[i][i];
+ }
+
+ for ( int i = 0; i < 3; i++) {
+ G[i][i] -= trace;
+ }
+
+
+ // inner product of the traceless symmetric part of the square of the velocity gradient tensor
+ T G_ip = 0;
+ for ( int i = 0; i < 3; i++) {
+ for ( int j = 0; j < 3; j++) {
+ G_ip = G[i][j] * G[i][j];
+ }
+ }
+
+ // inner product of the strain rate
+ T s_ip = 0;
+ for ( int i = 0; i < 3; i++) {
+ for ( int j = 0; j < 3; j++) {
+ s_ip = s[i][j] * s[i][j];
+ }
+ }
+
+ // Turbulent relaxation time
+ T tau_turb = 3. * this->getPreFactor() * (pow(G_ip,1.5) / (pow(s_ip,2.5) + pow(G_ip,1.25)));
+ if ((pow(s_ip,2.5) + pow(G_ip,1.25)) == 0) {
+ tau_turb = 0.;
+ }
+
+ // Physical turbulent viscosity must be equal or higher that zero
+ if (tau_turb < 0.) {
+ tau_turb = 0.;
+ }
+
+ /// Molecular relaxation time
+ T tau_mol = 1. /this->getOmega();
+
+ /// Effective relaxation time
+ T tau_eff = tau_mol + tau_turb;
+ T omega_new = 1. / tau_eff;
+
+ return omega_new;
+
+}
+
+////////////////////// Class WALEForcedBGKdynamics //////////////////////////
+/** \param vs2_ speed of sound
+ * \param momenta_ a Momenta object to know how to compute velocity momenta
+ * \param momenta_ a Momenta object to know how to compute velocity momenta
+ */
+template
+WALEForcedBGKdynamics::WALEForcedBGKdynamics(T omega_,
+ Momenta& momenta_, T smagoConst_)
+ : SmagorinskyForcedBGKdynamics(omega_, momenta_, smagoConst_)
+{
+ this->preFactor = this->getSmagoConst()*this->getSmagoConst();
+}
+
+template
+T WALEForcedBGKdynamics::computePreFactor()
+{
+ return (T)this->getSmagoConst()*this->getSmagoConst();
+}
+
+template