diff options
Diffstat (limited to 'src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh')
-rw-r--r-- | src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh | 417 |
1 files changed, 417 insertions, 0 deletions
diff --git a/src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh b/src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh new file mode 100644 index 0000000..90867bc --- /dev/null +++ b/src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh @@ -0,0 +1,417 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani, Jonas Latt + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * <http://www.openlb.net/> + * + * 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. +*/ + +#ifndef NAVIER_STOKES_ADVECTION_DIFFUSION_COUPLING_POST_PROCESSOR_2D_HH +#define NAVIER_STOKES_ADVECTION_DIFFUSION_COUPLING_POST_PROCESSOR_2D_HH + +#include "latticeDescriptors.h" +#include "navierStokesAdvectionDiffusionCouplingPostProcessor2D.h" +#include "core/blockLattice2D.h" +#include "core/util.h" +#include "core/finiteDifference2D.h" + +using namespace std; + +namespace olb { + +//===================================================================================== +//============== NavierStokesAdvectionDiffusionCouplingPostProcessor2D =============== +//===================================================================================== + +template<typename T, typename DESCRIPTOR> +NavierStokesAdvectionDiffusionCouplingPostProcessor2D<T,DESCRIPTOR>:: +NavierStokesAdvectionDiffusionCouplingPostProcessor2D(int x0_, int x1_, int y0_, int y1_, + T gravity_, T T0_, T deltaTemp_, std::vector<T> dir_, + std::vector<SpatiallyExtendedObject2D* > partners_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), + gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), + dir(dir_), partners(partners_) +{ + // we normalize the direction of force vector + T normDir = T(); + for (unsigned iD = 0; iD < dir.size(); ++iD) { + normDir += dir[iD]*dir[iD]; + } + normDir = sqrt(normDir); + for (unsigned iD = 0; iD < dir.size(); ++iD) { + dir[iD] /= normDir; + } + + for (unsigned iD = 0; iD < dir.size(); ++iD) { + forcePrefactor[iD] = gravity * dir[iD]; + } + + tPartner = dynamic_cast<BlockLattice2D<T,descriptors::D2Q5<descriptors::VELOCITY>> *>(partners[0]); +} + +template<typename T, typename DESCRIPTOR> +void NavierStokesAdvectionDiffusionCouplingPostProcessor2D<T,DESCRIPTOR>:: +processSubDomain(BlockLattice2D<T,DESCRIPTOR>& blockLattice, + int x0_, int x1_, int y0_, int y1_) +{ + + int newX0, newX1, newY0, newY1; + if ( util::intersect ( + x0, x1, y0, y1, + x0_, x1_, y0_, y1_, + newX0, newX1, newY0, newY1 ) ) { + + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + // computation of the bousinessq force + T *force = blockLattice.get(iX,iY).template getFieldPointer<descriptors::FORCE>(); + T temperatureDifference = tPartner->get(iX,iY).computeRho() - T0; + for (unsigned iD = 0; iD < L::d; ++iD) { + force[iD] = forcePrefactor[iD] * temperatureDifference; + } + // Velocity coupling + T *u = tPartner->get(iX,iY).template getFieldPointer<descriptors::VELOCITY>(); + blockLattice.get(iX,iY).computeU(u); + } + } + } +} + +template<typename T, typename DESCRIPTOR> +void NavierStokesAdvectionDiffusionCouplingPostProcessor2D<T,DESCRIPTOR>:: +process(BlockLattice2D<T,DESCRIPTOR>& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1); +} + +/// LatticeCouplingGenerator for advectionDiffusion coupling + +template<typename T, typename DESCRIPTOR> +NavierStokesAdvectionDiffusionCouplingGenerator2D<T,DESCRIPTOR>:: +NavierStokesAdvectionDiffusionCouplingGenerator2D(int x0_, int x1_, int y0_, int y1_, + T gravity_, T T0_, T deltaTemp_, std::vector<T> dir_) + : LatticeCouplingGenerator2D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_), + gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), dir(dir_) +{ } + +template<typename T, typename DESCRIPTOR> +PostProcessor2D<T,DESCRIPTOR>* NavierStokesAdvectionDiffusionCouplingGenerator2D<T,DESCRIPTOR>::generate ( + std::vector<SpatiallyExtendedObject2D* > partners) const +{ + return new NavierStokesAdvectionDiffusionCouplingPostProcessor2D<T,DESCRIPTOR>( + this->x0,this->x1,this->y0,this->y1, gravity, T0, deltaTemp, dir,partners); +} + +template<typename T, typename DESCRIPTOR> +LatticeCouplingGenerator2D<T,DESCRIPTOR>* NavierStokesAdvectionDiffusionCouplingGenerator2D<T,DESCRIPTOR>::clone() const +{ + return new NavierStokesAdvectionDiffusionCouplingGenerator2D<T,DESCRIPTOR>(*this); +} + + +//===================================================================================== +//============== SmagorinskyBoussinesqCouplingPostProcessor2D =============== +//===================================================================================== + +template<typename T, typename DESCRIPTOR> +SmagorinskyBoussinesqCouplingPostProcessor2D<T,DESCRIPTOR>:: +SmagorinskyBoussinesqCouplingPostProcessor2D(int x0_, int x1_, int y0_, int y1_, + T gravity_, T T0_, T deltaTemp_, std::vector<T> dir_, T PrTurb_, + std::vector<SpatiallyExtendedObject2D* > partners_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), + gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), + dir(dir_), PrTurb(PrTurb_), partners(partners_) +{ + // we normalize the direction of force vector + T normDir = T(); + for (unsigned iD = 0; iD < dir.size(); ++iD) { + normDir += dir[iD]*dir[iD]; + } + normDir = sqrt(normDir); + for (unsigned iD = 0; iD < dir.size(); ++iD) { + dir[iD] /= normDir; + } + + for (unsigned iD = 0; iD < dir.size(); ++iD) { + forcePrefactor[iD] = gravity * dir[iD]; + } + + tauTurbADPrefactor = descriptors::invCs2<T,descriptors::D2Q5<descriptors::VELOCITY,descriptors::TAU_EFF>>() / descriptors::invCs2<T,DESCRIPTOR>() / PrTurb; + tPartner = dynamic_cast<BlockLattice2D<T,descriptors::D2Q5<descriptors::VELOCITY,descriptors::TAU_EFF>> *>(partners[0]); +} + +template<typename T, typename DESCRIPTOR> +void SmagorinskyBoussinesqCouplingPostProcessor2D<T,DESCRIPTOR>:: +processSubDomain(BlockLattice2D<T,DESCRIPTOR>& blockLattice, + int x0_, int x1_, int y0_, int y1_) +{ + + int newX0, newX1, newY0, newY1; + if ( util::intersect ( + x0, x1, y0, y1, + x0_, x1_, y0_, y1_, + newX0, newX1, newY0, newY1 ) ) { + + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + + // computation of the bousinessq force + T *force = blockLattice.get(iX,iY).template getFieldPointer<descriptors::FORCE>(); + T temperatureDifference = tPartner->get(iX,iY).computeRho() - T0; + for (unsigned iD = 0; iD < L::d; ++iD) { + force[iD] = forcePrefactor[iD] * temperatureDifference; + } + + // Velocity coupling + T *u = tPartner->get(iX,iY).template getFieldPointer<descriptors::VELOCITY>(); + + // tau coupling + T *tauNS = blockLattice.get(iX,iY).template getFieldPointer<descriptors::TAU_EFF>(); + T *tauAD = tPartner->get(iX,iY).template getFieldPointer<descriptors::TAU_EFF>(); + + T rho, pi[util::TensorVal<DESCRIPTOR >::n]; + blockLattice.get(iX,iY).computeAllMomenta(rho, u, pi); + T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2]; + if (util::TensorVal<DESCRIPTOR >::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_NS = 1. / blockLattice.get(iX,iY).getDynamics()->getOmega(); + T tau_mol_AD = 1. / tPartner->get(iX,iY).getDynamics()->getOmega(); + /// Turbulent realaxation time + T tau_turb_NS = 0.5*(sqrt(tau_mol_NS*tau_mol_NS + dynamic_cast<SmagorinskyDynamics<T,DESCRIPTOR>*>(blockLattice.get(iX,iY).getDynamics())->getPreFactor()/rho*PiNeqNorm) - tau_mol_NS); + /// Effective realaxation time + tauNS[0] = tau_mol_NS+tau_turb_NS; + + T tau_turb_AD = tau_turb_NS * tauTurbADPrefactor; + tauAD[0] = tau_mol_AD+tau_turb_AD; + } + } + } + +} + +template<typename T, typename DESCRIPTOR> +void SmagorinskyBoussinesqCouplingPostProcessor2D<T,DESCRIPTOR>:: +process(BlockLattice2D<T,DESCRIPTOR>& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1); +} + +/// LatticeCouplingGenerator for advectionDiffusion coupling + +template<typename T, typename DESCRIPTOR> +SmagorinskyBoussinesqCouplingGenerator2D<T,DESCRIPTOR>:: +SmagorinskyBoussinesqCouplingGenerator2D(int x0_, int x1_, int y0_, int y1_, + T gravity_, T T0_, T deltaTemp_, std::vector<T> dir_, T PrTurb_) + : LatticeCouplingGenerator2D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_), + gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), dir(dir_), PrTurb(PrTurb_) +{ } + +template<typename T, typename DESCRIPTOR> +PostProcessor2D<T,DESCRIPTOR>* SmagorinskyBoussinesqCouplingGenerator2D<T,DESCRIPTOR>::generate ( + std::vector<SpatiallyExtendedObject2D* > partners) const +{ + return new SmagorinskyBoussinesqCouplingPostProcessor2D<T,DESCRIPTOR>( + this->x0,this->x1,this->y0,this->y1, gravity, T0, deltaTemp, dir, PrTurb, partners); +} + +template<typename T, typename DESCRIPTOR> +LatticeCouplingGenerator2D<T,DESCRIPTOR>* SmagorinskyBoussinesqCouplingGenerator2D<T,DESCRIPTOR>::clone() const +{ + return new SmagorinskyBoussinesqCouplingGenerator2D<T,DESCRIPTOR>(*this); +} + + +//===================================================================================== +//============== MixedScaleBoussinesqCouplingPostProcessor2D =============== +//===================================================================================== + +template<typename T, typename DESCRIPTOR> +MixedScaleBoussinesqCouplingPostProcessor2D<T,DESCRIPTOR>:: +MixedScaleBoussinesqCouplingPostProcessor2D(int x0_, int x1_, int y0_, int y1_, + T gravity_, T T0_, T deltaTemp_, std::vector<T> dir_, T PrTurb_, + std::vector<SpatiallyExtendedObject2D* > partners_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), + gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), + dir(dir_), PrTurb(PrTurb_), partners(partners_) +{ + // we normalize the direction of force vector + T normDir = T(); + for (unsigned iD = 0; iD < dir.size(); ++iD) { + normDir += dir[iD]*dir[iD]; + } + normDir = sqrt(normDir); + for (unsigned iD = 0; iD < dir.size(); ++iD) { + dir[iD] /= normDir; + } + + for (unsigned iD = 0; iD < dir.size(); ++iD) { + forcePrefactor[iD] = gravity * dir[iD]; + } + + tauTurbADPrefactor = descriptors::invCs2<T,descriptors::D2Q5<descriptors::VELOCITY,descriptors::TAU_EFF,descriptors::CUTOFF_HEAT_FLUX>>() / descriptors::invCs2<T,DESCRIPTOR>() / PrTurb; + tPartner = dynamic_cast<BlockLattice2D<T,descriptors::D2Q5<descriptors::VELOCITY,descriptors::TAU_EFF,descriptors::CUTOFF_HEAT_FLUX>> *>(partners[0]); +} + +template<typename T, typename DESCRIPTOR> +void MixedScaleBoussinesqCouplingPostProcessor2D<T,DESCRIPTOR>:: +processSubDomain(BlockLattice2D<T,DESCRIPTOR>& blockLattice, + int x0_, int x1_, int y0_, int y1_) +{ + + const T C_nu = 0.04; + const T C_alpha = 0.5; + const T deltaT = 1.0; + + const T invCs2_g = descriptors::invCs2<T,descriptors::D2Q5<descriptors::VELOCITY,descriptors::TAU_EFF,descriptors::CUTOFF_HEAT_FLUX>>(); + + int newX0, newX1, newY0, newY1; + if ( util::intersect ( + x0, x1, y0, y1, + x0_, x1_, y0_, y1_, + newX0, newX1, newY0, newY1 ) ) { + + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + T *cutoffKinEnergy_14 = blockLattice.get(iX, iY).template getFieldPointer<descriptors::CUTOFF_KIN_ENERGY>(); + T *cutoffHeatFlux_14 = tPartner->get(iX, iY).template getFieldPointer<descriptors::CUTOFF_HEAT_FLUX>(); + + // Velocity coupling + T *u = tPartner->get(iX,iY).template getFieldPointer<descriptors::VELOCITY>(); + + // tau coupling + T *tauNS = blockLattice.get(iX,iY).template getFieldPointer<descriptors::TAU_EFF>(); + T *tauAD = tPartner->get(iX,iY).template getFieldPointer<descriptors::TAU_EFF>(); + + /// Molecular realaxation time + T tau_mol_NS = 1. / blockLattice.get(iX,iY).getDynamics()->getOmega(); + T tau_mol_AD = 1. / tPartner->get(iX,iY).getDynamics()->getOmega(); + + const T temperature = tPartner->get(iX,iY).computeRho(); + + // computation of the bousinessq force + T *force = blockLattice.get(iX,iY).template getFieldPointer<descriptors::FORCE>(); + T temperatureDifference = temperature - T0; + for (unsigned iD = 0; iD < L::d; ++iD) { + force[iD] = forcePrefactor[iD] * temperatureDifference; + } + + T rho, pi[util::TensorVal<DESCRIPTOR>::n], j[DESCRIPTOR::d]; + blockLattice.get(iX,iY).computeAllMomenta(rho, u, pi); + + int iPi = 0; + for (int Alpha=0; Alpha<DESCRIPTOR::d; ++Alpha) { + for (int Beta=Alpha; Beta<DESCRIPTOR::d; ++Beta) { + pi[iPi] += rho/2.*(force[Alpha]*u[Beta] + u[Alpha]*force[Beta]); + ++iPi; + } + } + const T piSqr[3] = {pi[0]*pi[0], pi[1]*pi[1], pi[2]*pi[2]}; + const T PiNeqNormSqr = piSqr[0] + 2.0*piSqr[1] + piSqr[2]; + const T PiNeqNorm = sqrt(PiNeqNormSqr); + + tPartner->get(iX,iY).computeJ(j); + const T tmp_preFactor = invCs2_g / rho / tauAD[0]; + const T jNeq[2] = {(j[0] - temperature * u[0]), (j[1] - temperature * u[1])}; + const T jNeqSqr[2] = {jNeq[0]*jNeq[0], jNeq[1]*jNeq[1]}; + const T jNeqSqr_prefacor = 2. * 0.25 * (jNeq[0] + jNeq[1]) * (jNeq[0] + jNeq[1]); + + const T TnormSqr = jNeqSqr_prefacor*PiNeqNormSqr; + const T Tnorm = sqrt(TnormSqr); + + /// Turbulent realaxation time + // T tau_turb_NS = 0.5*(sqrt(tau_mol_NS*tau_mol_NS + dynamic_cast<SmagorinskyDynamics<T,DESCRIPTOR>*>(blockLattice.get(iX,iY).getDynamics())->getPreFactor()/rho*PiNeqNorm) - tau_mol_NS); + + // const T tmp_A = C_nu * sqrt(sqrt(2.)/2.) * descriptors::invCs2<T,DESCRIPTOR>() * descriptors::invCs2<T,DESCRIPTOR>() * sqrt(PiNeqNorm / rho) * cutoffKinEnergy_14[0]; + // const T tmp_A_2 = tmp_A * tmp_A; + // const T tmp_A_4 = tmp_A_2 * tmp_A_2; + + // const T tau_mol_NS_2 = tau_mol_NS * tau_mol_NS; + // const T tau_mol_NS_3 = tau_mol_NS_2 * tau_mol_NS; + + // const T tmp_1_3 = 1./3.; + // const T tmp_2_13 = pow(2., tmp_1_3); + // const T tmp_3_3_12 = 3. * sqrt(3.); + + // const T tmp_sqrtA = sqrt(27.*tmp_A_4-4.*tmp_A_2*tau_mol_NS_3); + + // // T tau_turb_NS = 1/3 ((27 A^2 + 3 sqrt(3) sqrt(27 A^4 - 4 A^2 b^3) - 2 b^3)^(1/3)/2^(1/3) + (2^(1/3) b^2)/(27 A^2 + 3 sqrt(3) sqrt(27 A^4 - 4 A^2 b^3) - 2 b^3)^(1/3) - b) + // T tau_turb_NS = ( pow(27.*tmp_A_2 + tmp_3_3_12*sqrt(27.*tmp_A_4-4.*tmp_A_2*tau_mol_NS_3)-2.*tau_mol_NS_3, tmp_1_3) / tmp_2_13 + // + (tmp_2_13*tau_mol_NS_2) / pow(27.*tmp_A_2+tmp_3_3_12*sqrt(27.*tmp_A_4-4.*tmp_A_2*tau_mol_NS_3) - 2.*tau_mol_NS_3, tmp_1_3) + // - tau_mol_NS + // ) * tmp_1_3; + + // if ( tau_turb_NS != tau_turb_NS ) + // tau_turb_NS = 0.; + + //cout << tau_turb_NS << " " << 27. * tmp_A_2 << " " << 4. * tau_mol_NS_3 << " " << PiNeqNorm << " " << " " << rho << endl; + + const T tmp_A = C_nu * sqrt(sqrt(2.)/2.) * descriptors::invCs2<T,DESCRIPTOR>() * descriptors::invCs2<T,DESCRIPTOR>() * sqrt(PiNeqNorm / rho / tauNS[0]) * cutoffKinEnergy_14[0]; + const T tau_turb_NS = tmp_A; + + // T tau_turb_AD = tau_turb_NS * tauTurbADPrefactor; + const T tmp_B = C_alpha * descriptors::invCs2<T,DESCRIPTOR>() / rho * sqrt(2.0 * Tnorm * invCs2_g / tauNS[0] / tauAD[0]) * cutoffHeatFlux_14[0]; + const T tau_turb_AD = tmp_B; + // cout << jNeq[0] << " " << jNeq[1] << " " << sqrt(Tnorm * invCs2_g / tauNS[0] / tauAD[0]) << " " << TnormSqr << endl; + + /// Effective realaxation time + tauNS[0] = tau_mol_NS+tau_turb_NS; + tauAD[0] = tau_mol_AD+tau_turb_AD; + + } + } + } + +} + +template<typename T, typename DESCRIPTOR> +void MixedScaleBoussinesqCouplingPostProcessor2D<T,DESCRIPTOR>:: +process(BlockLattice2D<T,DESCRIPTOR>& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1); +} + +/// LatticeCouplingGenerator for advectionDiffusion coupling + +template<typename T, typename DESCRIPTOR> +MixedScaleBoussinesqCouplingGenerator2D<T,DESCRIPTOR>:: +MixedScaleBoussinesqCouplingGenerator2D(int x0_, int x1_, int y0_, int y1_, + T gravity_, T T0_, T deltaTemp_, std::vector<T> dir_, T PrTurb_) + : LatticeCouplingGenerator2D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_), + gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), dir(dir_), PrTurb(PrTurb_) +{ } + +template<typename T, typename DESCRIPTOR> +PostProcessor2D<T,DESCRIPTOR>* MixedScaleBoussinesqCouplingGenerator2D<T,DESCRIPTOR>::generate ( + std::vector<SpatiallyExtendedObject2D* > partners) const +{ + return new MixedScaleBoussinesqCouplingPostProcessor2D<T,DESCRIPTOR>( + this->x0,this->x1,this->y0,this->y1, gravity, T0, deltaTemp, dir, PrTurb, partners); +} + +template<typename T, typename DESCRIPTOR> +LatticeCouplingGenerator2D<T,DESCRIPTOR>* MixedScaleBoussinesqCouplingGenerator2D<T,DESCRIPTOR>::clone() const +{ + return new MixedScaleBoussinesqCouplingGenerator2D<T,DESCRIPTOR>(*this); +} + +} // namespace olb + +#endif |