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/ADMlatticeDescriptors.h | 63 + src/dynamics/MakeHeader | 33 + .../SmagorinskyPorousParticleBGKdynamics.h | 71 + .../SmagorinskyPorousParticleBGKdynamics.hh | 119 ++ src/dynamics/SmagorinskyPowerLawBGKdynamics.h | 73 + src/dynamics/SmagorinskyPowerLawBGKdynamics.hh | 154 ++ .../SmagorinskyPowerLawPorousBGKdynamics.h | 66 + .../SmagorinskyPowerLawPorousBGKdynamics.hh | 133 ++ src/dynamics/WALELatticeDescriptors.h | 54 + src/dynamics/advectionDiffusionDynamics.h | 201 ++ src/dynamics/advectionDiffusionDynamics.hh | 442 +++++ src/dynamics/advectionDiffusionForces.h | 108 ++ src/dynamics/advectionDiffusionForces.hh | 130 ++ src/dynamics/advectionDiffusionMomenta.h | 88 + src/dynamics/advectionDiffusionMomenta.hh | 159 ++ src/dynamics/chopardDynamics.cpp | 39 + src/dynamics/chopardDynamics.h | 67 + src/dynamics/chopardDynamics.hh | 139 ++ src/dynamics/d3q13Helpers.h | 141 ++ src/dynamics/descriptorBase.h | 110 ++ src/dynamics/descriptorField.h | 201 ++ src/dynamics/descriptorFunction.h | 190 ++ src/dynamics/descriptorTag.h | 75 + src/dynamics/dynOmegaLatticeDescriptors.h | 73 + src/dynamics/dynSmagorinskyLatticeDescriptors.h | 47 + src/dynamics/dynamics.cpp | 105 + src/dynamics/dynamics.h | 995 ++++++++++ src/dynamics/dynamics.hh | 2026 ++++++++++++++++++++ src/dynamics/dynamics2D.h | 64 + src/dynamics/dynamics2D.hh | 56 + src/dynamics/dynamics3D.h | 68 + src/dynamics/dynamics3D.hh | 59 + src/dynamics/entropicDynamics.cpp | 41 + src/dynamics/entropicDynamics.h | 143 ++ src/dynamics/entropicDynamics.hh | 433 +++++ src/dynamics/entropicLbHelpers.h | 79 + src/dynamics/entropicLbHelpers2D.h | 102 + src/dynamics/entropicLbHelpers3D.h | 106 + src/dynamics/firstOrderLbHelpers.h | 132 ++ src/dynamics/firstOrderLbHelpers2D.h | 105 + src/dynamics/firstOrderLbHelpers3D.h | 285 +++ src/dynamics/freeEnergyDynamics.cpp | 48 + src/dynamics/freeEnergyDynamics.h | 120 ++ src/dynamics/freeEnergyDynamics.hh | 243 +++ src/dynamics/freeEnergyPostProcessor2D.h | 242 +++ src/dynamics/freeEnergyPostProcessor2D.hh | 511 +++++ src/dynamics/freeEnergyPostProcessor3D.h | 231 +++ src/dynamics/freeEnergyPostProcessor3D.hh | 631 ++++++ src/dynamics/guoZhaoDynamics.h | 74 + src/dynamics/guoZhaoDynamics.hh | 143 ++ src/dynamics/guoZhaoLatticeDescriptors.h | 67 + src/dynamics/guoZhaoLbHelpers.h | 166 ++ src/dynamics/interactionPotential.cpp | 62 + src/dynamics/interactionPotential.h | 140 ++ src/dynamics/interactionPotential.hh | 194 ++ src/dynamics/latticeDescriptors.cpp | 51 + src/dynamics/latticeDescriptors.h | 404 ++++ src/dynamics/lbHelpers.h | 678 +++++++ src/dynamics/lbHelpersD2Q5.h | 222 +++ src/dynamics/lbHelpersD2Q9.h | 343 ++++ src/dynamics/lbHelpersD3Q15.h | 305 +++ src/dynamics/lbHelpersD3Q19.h | 567 ++++++ src/dynamics/lbHelpersD3Q27.h | 548 ++++++ src/dynamics/lbHelpersD3Q7.h | 305 +++ src/dynamics/module.mk | 27 + src/dynamics/mrtDynamics.cpp | 75 + src/dynamics/mrtDynamics.h | 97 + src/dynamics/mrtDynamics.hh | 203 ++ src/dynamics/mrtHelpers.h | 182 ++ src/dynamics/mrtHelpers2D.h | 229 +++ src/dynamics/mrtHelpers3D.h | 586 ++++++ src/dynamics/mrtLatticeDescriptors.h | 477 +++++ ...okesAdvectionDiffusionCouplingPostProcessor2D.h | 184 ++ ...kesAdvectionDiffusionCouplingPostProcessor2D.hh | 417 ++++ ...okesAdvectionDiffusionCouplingPostProcessor3D.h | 244 +++ ...kesAdvectionDiffusionCouplingPostProcessor3D.hh | 485 +++++ ...sAdvectionDiffusionMRTCouplingPostProcessor2D.h | 91 + ...AdvectionDiffusionMRTCouplingPostProcessor2D.hh | 134 ++ ...sAdvectionDiffusionMRTCouplingPostProcessor3D.h | 149 ++ ...AdvectionDiffusionMRTCouplingPostProcessor3D.hh | 281 +++ src/dynamics/porousAdvectionDiffusionDescriptors.h | 43 + src/dynamics/porousAdvectionDiffusionDynamics.h | 63 + src/dynamics/porousAdvectionDiffusionDynamics.hh | 117 ++ src/dynamics/porousBGKdynamics.h | 273 +++ src/dynamics/porousBGKdynamics.hh | 671 +++++++ src/dynamics/porousForcedBGKDynamics.h | 65 + src/dynamics/porousForcedBGKDynamics.hh | 104 + src/dynamics/porousLatticeDescriptors.h | 87 + src/dynamics/powerLawBGKdynamics.h | 93 + src/dynamics/powerLawBGKdynamics.hh | 157 ++ src/dynamics/rtlbmDescriptors.h | 104 + src/dynamics/rtlbmDynamics.h | 90 + src/dynamics/rtlbmDynamics.hh | 187 ++ src/dynamics/shanChenDynGForcedPostProcessor2D.h | 87 + src/dynamics/shanChenDynGForcedPostProcessor2D.hh | 205 ++ .../shanChenDynOmegaForcedPostProcessor2D.h | 85 + .../shanChenDynOmegaForcedPostProcessor2D.hh | 197 ++ .../shanChenDynOmegaForcedPostProcessor3D.h | 85 + .../shanChenDynOmegaForcedPostProcessor3D.hh | 206 ++ src/dynamics/shanChenForcedLatticeDescriptors.h | 70 + src/dynamics/shanChenForcedPostProcessor2D.h | 85 + src/dynamics/shanChenForcedPostProcessor2D.hh | 197 ++ src/dynamics/shanChenForcedPostProcessor3D.h | 85 + src/dynamics/shanChenForcedPostProcessor3D.hh | 206 ++ .../shanChenForcedSingleComponentPostProcessor2D.h | 87 + ...shanChenForcedSingleComponentPostProcessor2D.hh | 171 ++ .../shanChenForcedSingleComponentPostProcessor3D.h | 85 + ...shanChenForcedSingleComponentPostProcessor3D.hh | 176 ++ src/dynamics/shearSmagorinskyLatticeDescriptors.h | 77 + src/dynamics/smagorinskyBGKdynamics.cpp | 44 + src/dynamics/smagorinskyBGKdynamics.h | 374 ++++ src/dynamics/smagorinskyBGKdynamics.hh | 1361 +++++++++++++ src/dynamics/smagorinskyGuoZhaoDynamics.h | 75 + src/dynamics/smagorinskyGuoZhaoDynamics.hh | 118 ++ src/dynamics/smagorinskyMRTdynamics.h | 95 + src/dynamics/smagorinskyMRTdynamics.hh | 164 ++ src/dynamics/stochasticSGSdynamics.h | 102 + src/dynamics/stochasticSGSdynamics.hh | 293 +++ src/dynamics/superGuoZhaoPostProcessor2D.h | 54 + src/dynamics/superGuoZhaoPostProcessor2D.hh | 92 + src/dynamics/wallFunctionLatticeDescriptors.h | 42 + 121 files changed, 24903 insertions(+) create mode 100644 src/dynamics/ADMlatticeDescriptors.h create mode 100644 src/dynamics/MakeHeader create mode 100644 src/dynamics/SmagorinskyPorousParticleBGKdynamics.h create mode 100644 src/dynamics/SmagorinskyPorousParticleBGKdynamics.hh create mode 100644 src/dynamics/SmagorinskyPowerLawBGKdynamics.h create mode 100644 src/dynamics/SmagorinskyPowerLawBGKdynamics.hh create mode 100644 src/dynamics/SmagorinskyPowerLawPorousBGKdynamics.h create mode 100644 src/dynamics/SmagorinskyPowerLawPorousBGKdynamics.hh create mode 100644 src/dynamics/WALELatticeDescriptors.h create mode 100644 src/dynamics/advectionDiffusionDynamics.h create mode 100644 src/dynamics/advectionDiffusionDynamics.hh create mode 100644 src/dynamics/advectionDiffusionForces.h create mode 100644 src/dynamics/advectionDiffusionForces.hh create mode 100644 src/dynamics/advectionDiffusionMomenta.h create mode 100644 src/dynamics/advectionDiffusionMomenta.hh create mode 100644 src/dynamics/chopardDynamics.cpp create mode 100644 src/dynamics/chopardDynamics.h create mode 100644 src/dynamics/chopardDynamics.hh create mode 100644 src/dynamics/d3q13Helpers.h create mode 100644 src/dynamics/descriptorBase.h create mode 100644 src/dynamics/descriptorField.h create mode 100644 src/dynamics/descriptorFunction.h create mode 100644 src/dynamics/descriptorTag.h create mode 100644 src/dynamics/dynOmegaLatticeDescriptors.h create mode 100644 src/dynamics/dynSmagorinskyLatticeDescriptors.h create mode 100644 src/dynamics/dynamics.cpp create mode 100644 src/dynamics/dynamics.h create mode 100644 src/dynamics/dynamics.hh create mode 100644 src/dynamics/dynamics2D.h create mode 100644 src/dynamics/dynamics2D.hh create mode 100644 src/dynamics/dynamics3D.h create mode 100644 src/dynamics/dynamics3D.hh create mode 100644 src/dynamics/entropicDynamics.cpp create mode 100644 src/dynamics/entropicDynamics.h create mode 100644 src/dynamics/entropicDynamics.hh create mode 100644 src/dynamics/entropicLbHelpers.h create mode 100644 src/dynamics/entropicLbHelpers2D.h create mode 100644 src/dynamics/entropicLbHelpers3D.h create mode 100644 src/dynamics/firstOrderLbHelpers.h create mode 100644 src/dynamics/firstOrderLbHelpers2D.h create mode 100644 src/dynamics/firstOrderLbHelpers3D.h create mode 100644 src/dynamics/freeEnergyDynamics.cpp create mode 100644 src/dynamics/freeEnergyDynamics.h create mode 100644 src/dynamics/freeEnergyDynamics.hh create mode 100644 src/dynamics/freeEnergyPostProcessor2D.h create mode 100644 src/dynamics/freeEnergyPostProcessor2D.hh create mode 100644 src/dynamics/freeEnergyPostProcessor3D.h create mode 100644 src/dynamics/freeEnergyPostProcessor3D.hh create mode 100644 src/dynamics/guoZhaoDynamics.h create mode 100644 src/dynamics/guoZhaoDynamics.hh create mode 100644 src/dynamics/guoZhaoLatticeDescriptors.h create mode 100644 src/dynamics/guoZhaoLbHelpers.h create mode 100644 src/dynamics/interactionPotential.cpp create mode 100644 src/dynamics/interactionPotential.h create mode 100644 src/dynamics/interactionPotential.hh create mode 100644 src/dynamics/latticeDescriptors.cpp create mode 100644 src/dynamics/latticeDescriptors.h create mode 100644 src/dynamics/lbHelpers.h create mode 100644 src/dynamics/lbHelpersD2Q5.h create mode 100644 src/dynamics/lbHelpersD2Q9.h create mode 100644 src/dynamics/lbHelpersD3Q15.h create mode 100644 src/dynamics/lbHelpersD3Q19.h create mode 100644 src/dynamics/lbHelpersD3Q27.h create mode 100644 src/dynamics/lbHelpersD3Q7.h create mode 100644 src/dynamics/module.mk create mode 100644 src/dynamics/mrtDynamics.cpp create mode 100644 src/dynamics/mrtDynamics.h create mode 100644 src/dynamics/mrtDynamics.hh create mode 100644 src/dynamics/mrtHelpers.h create mode 100644 src/dynamics/mrtHelpers2D.h create mode 100644 src/dynamics/mrtHelpers3D.h create mode 100644 src/dynamics/mrtLatticeDescriptors.h create mode 100644 src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.h create mode 100644 src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh create mode 100644 src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor3D.h create mode 100644 src/dynamics/navierStokesAdvectionDiffusionCouplingPostProcessor3D.hh create mode 100644 src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor2D.h create mode 100644 src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor2D.hh create mode 100644 src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.h create mode 100644 src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh create mode 100644 src/dynamics/porousAdvectionDiffusionDescriptors.h create mode 100644 src/dynamics/porousAdvectionDiffusionDynamics.h create mode 100644 src/dynamics/porousAdvectionDiffusionDynamics.hh create mode 100644 src/dynamics/porousBGKdynamics.h create mode 100644 src/dynamics/porousBGKdynamics.hh create mode 100644 src/dynamics/porousForcedBGKDynamics.h create mode 100644 src/dynamics/porousForcedBGKDynamics.hh create mode 100644 src/dynamics/porousLatticeDescriptors.h create mode 100644 src/dynamics/powerLawBGKdynamics.h create mode 100644 src/dynamics/powerLawBGKdynamics.hh create mode 100644 src/dynamics/rtlbmDescriptors.h create mode 100644 src/dynamics/rtlbmDynamics.h create mode 100644 src/dynamics/rtlbmDynamics.hh create mode 100644 src/dynamics/shanChenDynGForcedPostProcessor2D.h create mode 100644 src/dynamics/shanChenDynGForcedPostProcessor2D.hh create mode 100644 src/dynamics/shanChenDynOmegaForcedPostProcessor2D.h create mode 100644 src/dynamics/shanChenDynOmegaForcedPostProcessor2D.hh create mode 100644 src/dynamics/shanChenDynOmegaForcedPostProcessor3D.h create mode 100644 src/dynamics/shanChenDynOmegaForcedPostProcessor3D.hh create mode 100644 src/dynamics/shanChenForcedLatticeDescriptors.h create mode 100644 src/dynamics/shanChenForcedPostProcessor2D.h create mode 100644 src/dynamics/shanChenForcedPostProcessor2D.hh create mode 100644 src/dynamics/shanChenForcedPostProcessor3D.h create mode 100644 src/dynamics/shanChenForcedPostProcessor3D.hh create mode 100644 src/dynamics/shanChenForcedSingleComponentPostProcessor2D.h create mode 100644 src/dynamics/shanChenForcedSingleComponentPostProcessor2D.hh create mode 100644 src/dynamics/shanChenForcedSingleComponentPostProcessor3D.h create mode 100644 src/dynamics/shanChenForcedSingleComponentPostProcessor3D.hh create mode 100644 src/dynamics/shearSmagorinskyLatticeDescriptors.h create mode 100644 src/dynamics/smagorinskyBGKdynamics.cpp create mode 100644 src/dynamics/smagorinskyBGKdynamics.h create mode 100644 src/dynamics/smagorinskyBGKdynamics.hh create mode 100644 src/dynamics/smagorinskyGuoZhaoDynamics.h create mode 100644 src/dynamics/smagorinskyGuoZhaoDynamics.hh create mode 100644 src/dynamics/smagorinskyMRTdynamics.h create mode 100644 src/dynamics/smagorinskyMRTdynamics.hh create mode 100644 src/dynamics/stochasticSGSdynamics.h create mode 100644 src/dynamics/stochasticSGSdynamics.hh create mode 100644 src/dynamics/superGuoZhaoPostProcessor2D.h create mode 100644 src/dynamics/superGuoZhaoPostProcessor2D.hh create mode 100644 src/dynamics/wallFunctionLatticeDescriptors.h (limited to 'src/dynamics') diff --git a/src/dynamics/ADMlatticeDescriptors.h b/src/dynamics/ADMlatticeDescriptors.h new file mode 100644 index 0000000..d6bb59d --- /dev/null +++ b/src/dynamics/ADMlatticeDescriptors.h @@ -0,0 +1,63 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012 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. +*/ + +#ifndef ADM_BGK_DYNAMICS_DESCRIPTOR_H +#define ADM_BGK_DYNAMICS_DESCRIPTOR_H + +#include "dynamics/latticeDescriptors.h" +//#include + + +namespace olb { + +namespace descriptors { + +// 2D Descriptors for ADM + +using ADMD2Q9Descriptor = D2Q9; + +//////////////////////////////////////////////////////////////////////////////// +// extended descriptor for ADM + +using ADMD3Q19Descriptor = D3Q19; + +//////////////////////////////////////// +/// ADM Descriptors for forced fields + +//// Forced 2D ADM scheme +using ForcedADMD2Q9Descriptor = D2Q9; + +//// Forced 3D ADM scheme + +using ForcedADMD3Q19Descriptor = D3Q19; + +//// Forced adapted 3D ADM scheme + +using ForcedAdaptiveADMD3Q19Descriptor = D3Q19; + + +} // namespace descriptors + +} // namespace olb + +#endif diff --git a/src/dynamics/MakeHeader b/src/dynamics/MakeHeader new file mode 100644 index 0000000..b1acde0 --- /dev/null +++ b/src/dynamics/MakeHeader @@ -0,0 +1,33 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias 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. + + +generic := + +precompiled := chopardDynamics \ + dynamics \ + entropicDynamics \ + latticeDescriptors \ + smagorinskyBGKdynamics \ + freeEnergyDynamics + + diff --git a/src/dynamics/SmagorinskyPorousParticleBGKdynamics.h b/src/dynamics/SmagorinskyPorousParticleBGKdynamics.h new file mode 100644 index 0000000..414408f --- /dev/null +++ b/src/dynamics/SmagorinskyPorousParticleBGKdynamics.h @@ -0,0 +1,71 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2017 Davide Dapelo + * 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 + * Smagorinsky BGK Dynamics for porous media -- header file. + */ +#ifndef SMAGORINSKY_POROUS_PARTICLE_BGK_DYNAMICS_H +#define SMAGORINSKY_POROUS_PARTICLE_BGK_DYNAMICS_H + +#include "dynamics/porousBGKdynamics.h" + +namespace olb { + +/// Implementation of the BGK collision step for a porosity model +template +class SmagorinskyPorousParticleBGKdynamics : public PorousParticleBGKdynamics { +public: + /// Constructor + SmagorinskyPorousParticleBGKdynamics(T omega_, Momenta& momenta_, T smagoConst_, + T dx_ = 1, T dt_ = 1); + + /// Collision step + virtual void collide(Cell& cell, + LatticeStatistics& statistics_); + + /// set relaxation parameter + void setOmega(T omega_); + /// Get local smagorinsky relaxation parameter of the dynamics + virtual T getSmagorinskyOmega(Cell& cell_); + + +protected: + /// Computes a constant prefactor in order to speed up the computation + T computePreFactor(T omega_, T smagoConst_); + /// Computes the local smagorinsky relaxation parameter + T computeOmega(T omega0, T preFactor_, T rho, + T pi[util::TensorVal::n] ); + + /// effective collision time based upon Smagorisnky approach + T tau_eff; + /// Smagorinsky constant + T smagoConst; + /// Precomputed constant which speeeds up the computation + T preFactor; + T dx; + T dt; +}; + +} // olb + +#endif diff --git a/src/dynamics/SmagorinskyPorousParticleBGKdynamics.hh b/src/dynamics/SmagorinskyPorousParticleBGKdynamics.hh new file mode 100644 index 0000000..482c7db --- /dev/null +++ b/src/dynamics/SmagorinskyPorousParticleBGKdynamics.hh @@ -0,0 +1,119 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2017 Davide Dapelo + * 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 + * Smagorinsky BGK Dynamics for porous -- generic implementation. + */ +#ifndef SMAGORINSKY_POROUS_PARTICLE_BGK_DYNAMICS_HH +#define SMAGORINSKY_POROUS_PARTICLE_BGK_DYNAMICS_HH + +#include "dynamics/porousBGKdynamics.hh" + +namespace olb { + +////////////////////// Class PorousParticleBGKdynamics ////////////////////////// + +template +SmagorinskyPorousParticleBGKdynamics::SmagorinskyPorousParticleBGKdynamics(T omega_, + Momenta& momenta_, T smagoConst_, T dx_, T dt_ ) + : PorousParticleBGKdynamics(omega_,momenta_), smagoConst(smagoConst_), + preFactor(computePreFactor(omega_,smagoConst_) ) +{ } + +template +void SmagorinskyPorousParticleBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell, rho, u, pi); +/* + T* porosity = cell.template getFieldPointer(); + for (int i=0; i(); + T* velNumerator = cell.template getFieldPointer(); + T* porosity = cell.template getFieldPointer(); + if (*velDenominator > std::numeric_limits::epsilon()) { + *porosity = 1.-*porosity; // 1-prod(1-smoothInd) + for (int i=0; i < DESCRIPTOR::d; i++) { + u[i] += *porosity * (*(velNumerator+i) / *velDenominator - u[i]); + } + } + T newOmega = computeOmega(this->getOmega(), preFactor, rho, pi); + T uSqr = lbHelpers::bgkCollision(cell, rho, u, newOmega); + statistics.incrementStats(rho, uSqr); + + cell.template defineField(1.0); + cell.template defineField(0.0); + cell.template defineField(0.0); +} + +template +void SmagorinskyPorousParticleBGKdynamics::setOmega(T omega_) +{ + PorousParticleBGKdynamics::setOmega(omega_); + preFactor = computePreFactor(omega_, smagoConst); +} + +template +T SmagorinskyPorousParticleBGKdynamics::getSmagorinskyOmega(Cell& cell ) +{ + T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell, rho, uTemp, pi); + T newOmega = computeOmega(this->getOmega(), preFactor, rho, pi); + return newOmega; +} + +template +T SmagorinskyPorousParticleBGKdynamics::computePreFactor(T omega_, T smagoConst_) +{ + return (T)smagoConst_*smagoConst_*descriptors::invCs2()*descriptors::invCs2()*2*sqrt(2); +} + + + +template +T SmagorinskyPorousParticleBGKdynamics::computeOmega(T omega0, T preFactor_, T rho, + T pi[util::TensorVal::n] ) +{ + 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. /omega0; + /// Turbulent realaxation time + T tau_turb = 0.5*(sqrt(tau_mol*tau_mol + preFactor_/rho*PiNeqNorm) - tau_mol); + /// Effective realaxation time + tau_eff = tau_mol+tau_turb; + T omega_new= 1./tau_eff; + return omega_new; +} + +} // olb + +#endif diff --git a/src/dynamics/SmagorinskyPowerLawBGKdynamics.h b/src/dynamics/SmagorinskyPowerLawBGKdynamics.h new file mode 100644 index 0000000..e0dad17 --- /dev/null +++ b/src/dynamics/SmagorinskyPowerLawBGKdynamics.h @@ -0,0 +1,73 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012, 2015 Mathias J. Krause, Vojtech Cvrcek, Davide Dapelo + * 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 + * Porous-particle BGK Dynamics with adjusted omega + * and Smagorinsky turbulence model -- header file. + * Strain rate similar to "J.Boyd, J. Buick and S.Green: A second-order accurate lattice Boltzmann non-Newtonian flow model" + * Power Law similar to "Huidan Yu, Sharath S. Girimaji, Li-Shi Luo - DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method" + * + */ +#ifndef SMAGORINSKY_POWER_LAW_BGK_DYNAMICS_H +#define SMAGORINSKY_POWER_LAW_BGK_DYNAMICS_H + +#include "core/cell.h" + +namespace olb { + +/// Implementation of the BGK collision step +template +class SmagorinskyPowerLawBGKdynamics : public SmagorinskyBGKdynamics, PowerLawDynamics { +public: + /// Constructor + /// m,n...parameter in the power law model + SmagorinskyPowerLawBGKdynamics(T omega, Momenta& momenta, T m=0.1, T n=.5, T nuMin=T(2.9686e-3), T nuMax=T(3.1667), T smagoConst=T(0.14)); + /// Collision step + virtual void collide(Cell& cell, LatticeStatistics& statistics) override; + /// Computes the local smagorinsky relaxation parameter with omega0 as an input + virtual T computeEffectiveOmega(Cell& cell, T omega0); +protected: + /// Computes squared norm of non-equilibrium part of 2nd momentum (for power-law) + virtual T PiNeqNormSqr(Cell& cell ) override; +}; + +/// Implementation of the ForcedBGK collision step +template +class SmagorinskyPowerLawForcedBGKdynamics : public SmagorinskyForcedBGKdynamics, public PowerLawDynamics { +public: + /// Constructor + /// m,n...parameter in the power law model + SmagorinskyPowerLawForcedBGKdynamics(T omega, Momenta& momenta, T m=0.1, T n=.5, T nuMin=T(2.9686e-3), T nuMax=T(3.1667), T smagoConst=T(0.14)); + /// Collision step + virtual void collide(Cell& cell, LatticeStatistics& statistics) override; + +protected: + /// Computes the local smagorinsky relaxation parameter with omega0 as an input + virtual T computeEffectiveOmega(Cell& cell, T omega0); + /// Computes squared norm of non-equilibrium part of 2nd momentum (for power-law) + virtual T PiNeqNormSqr(Cell& cell ) override; +}; + +} + +#endif diff --git a/src/dynamics/SmagorinskyPowerLawBGKdynamics.hh b/src/dynamics/SmagorinskyPowerLawBGKdynamics.hh new file mode 100644 index 0000000..181f4d6 --- /dev/null +++ b/src/dynamics/SmagorinskyPowerLawBGKdynamics.hh @@ -0,0 +1,154 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012, 2015 Mathias J. Krause, Vojtech Cvrcekt, Davide Dapelo + * 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 + * Porous-particle BGK Dynamics with adjusted omega + * and Smagorinsky turbulence model -- generic implementation. + * Strain rate similar to "J.Boyd, J. Buick and S.Green: A second-order accurate lattice Boltzmann non-Newtonian flow model" + * Power Law similar to "Huidan Yu, Sharath S. Girimaji, Li-Shi Luo - DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method" + */ +#ifndef SMAGORINSKY_POWER_LAW_BGK_DYNAMICS_HH +#define SMAGORINSKY_POWER_LAW_BGK_DYNAMICS_HH + +#include "../dynamics/powerLawBGKdynamics.h" +#include "SmagorinskyPowerLawBGKdynamics.h" +#include "math.h" + +namespace olb { + +////////////////////// Class SmagorinskyPowerLawBGKdynamics ////////////////////////// + +/** \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 +SmagorinskyPowerLawBGKdynamics::SmagorinskyPowerLawBGKdynamics ( + T omega, Momenta& momenta, T m, T n , T nuMin, T nuMax, T smagoConst) + : SmagorinskyBGKdynamics(omega, momenta, smagoConst), + PowerLawDynamics(m, n, nuMin, nuMax) +{ } + +template +void SmagorinskyPowerLawBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell, rho, u, pi); + + // Computation of the power-law omega. + // An external is used in place of BGKdynamics::_omega to keep generality and flexibility. + T oldOmega = cell.template getFieldPointer()[0]; + T intOmega = this->computeOmegaPL(cell, oldOmega, rho, pi); + T newOmega = computeEffectiveOmega(cell, intOmega); // turbulent omega + + T uSqr = lbHelpers::bgkCollision(cell, rho, u, newOmega); + cell.template getFieldPointer()[0] = intOmega; // updating omega + statistics.incrementStats(rho, uSqr); +} + +template +T SmagorinskyPowerLawBGKdynamics::computeEffectiveOmega(Cell& cell, T omega0) +{ + T rho = this->_momenta.computeRho(cell); + T PiNeqNorm = sqrt(PiNeqNormSqr(cell)); + /// Molecular realaxation time + T tau_mol = 1. /omega0; + /// 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; +} + +template +T SmagorinskyPowerLawBGKdynamics::PiNeqNormSqr(Cell& cell ) +{ + return lbHelpers::computePiNeqNormSqr(cell); +} + + +////////////////////// Class SmagorinskyPowerLawForcedBGKdynamics ////////////////////////// + +/** \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 +SmagorinskyPowerLawForcedBGKdynamics::SmagorinskyPowerLawForcedBGKdynamics ( + T omega, Momenta& momenta, T m, T n , T nuMin, T nuMax, T smagoConst) + : SmagorinskyForcedBGKdynamics(omega, momenta, smagoConst), + PowerLawDynamics(m, n, nuMin, nuMax) +{ } + +template +void SmagorinskyPowerLawForcedBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell, rho, u, pi); + + // Computation of the power-law omega. + // An external is used in place of BGKdynamics::_omega to keep generality and flexibility. + T oldOmega = cell.template getFieldPointer()[0]; + T intOmega = this->computeOmegaPL(cell, oldOmega, rho, pi); + T newOmega = computeEffectiveOmega(cell, intOmega); // turbulent omega + + T* force = cell.template getFieldPointer(); + for (int iVel=0; iVel::bgkCollision(cell, rho, u, newOmega); + cell.template getFieldPointer()[0] = intOmega; // updating omega + lbHelpers::addExternalForce(cell, u, newOmega, rho); + statistics.incrementStats(rho, uSqr); +} + +template +T SmagorinskyPowerLawForcedBGKdynamics::computeEffectiveOmega(Cell& cell, T omega0) +{ + T rho = this->_momenta.computeRho(cell); + T PiNeqNorm = sqrt(PiNeqNormSqr(cell)); + /// Molecular realaxation time + T tau_mol = 1. /omega0; + /// 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; +} + +template +T SmagorinskyPowerLawForcedBGKdynamics::PiNeqNormSqr(Cell& cell ) +{ + return lbHelpers::computeForcedPiNeqNormSqr(cell); +} + +} + +#endif diff --git a/src/dynamics/SmagorinskyPowerLawPorousBGKdynamics.h b/src/dynamics/SmagorinskyPowerLawPorousBGKdynamics.h new file mode 100644 index 0000000..d3a43c7 --- /dev/null +++ b/src/dynamics/SmagorinskyPowerLawPorousBGKdynamics.h @@ -0,0 +1,66 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012, 2015 Mathias J. Krause, Vojtech Cvrcek, Davide Dapelo + * 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 + * Porous-particle BGK Dynamics with adjusted omega + * and Smagorinsky turbulence model -- header file. + * Strain rate similar to "J.Boyd, J. Buick and S.Green: A second-order accurate lattice Boltzmann non-Newtonian flow model" + * Power Law similar to "Huidan Yu, Sharath S. Girimaji, Li-Shi Luo - DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method" + * + */ +#ifndef SMAGORINSKY_POWER_LAW_POROUS_BGK_DYNAMICS_H +#define SMAGORINSKY_POWER_LAW_POROUS_BGK_DYNAMICS_H + +#include "SmagorinskyPorousParticleBGKdynamics.h" +#include "core/cell.h" + +namespace olb { + +/// Implementation of the BGK collision step +template +class SmagorinskyPowerLawPorousParticleBGKdynamics : public SmagorinskyPorousParticleBGKdynamics { +public: + /// Constructor + /// m,n...parameter in the power law model, dt...the explizit typed time step from LBM model + SmagorinskyPowerLawPorousParticleBGKdynamics(T omega_, Momenta& momenta_, T m_=0.1, T n_=.5, T dtPL_=T(0.0016), T nuMin=T(2.9686e-3), T nuMax=T(3.1667), T smagoConst_=T(0.14)); + + /// Collision step + virtual void collide(Cell& cell, + LatticeStatistics& statistics_); + +private: + + /// Computes the local powerLaw relaxation parameter + T computeOmegaPL(T omega0_, T rho_, T pi_[util::TensorVal::n] ); + +private: + T m; + T n; + T dtPL; + T omegaMin; + T omegaMax; +}; + +} + +#endif diff --git a/src/dynamics/SmagorinskyPowerLawPorousBGKdynamics.hh b/src/dynamics/SmagorinskyPowerLawPorousBGKdynamics.hh new file mode 100644 index 0000000..7568070 --- /dev/null +++ b/src/dynamics/SmagorinskyPowerLawPorousBGKdynamics.hh @@ -0,0 +1,133 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012, 2015 Mathias J. Krause, Vojtech Cvrcekt, Davide Dapelo + * 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 + * Porous-particle BGK Dynamics with adjusted omega + * and Smagorinsky turbulence model -- generic implementation. + * Strain rate similar to "J.Boyd, J. Buick and S.Green: A second-order accurate lattice Boltzmann non-Newtonian flow model" + * Power Law similar to "Huidan Yu, Sharath S. Girimaji, Li-Shi Luo - DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method" + */ +#ifndef SMAGORINSKY_POWER_LAW_POROUS_BGK_DYNAMICS_HH +#define SMAGORINSKY_POWER_LAW_POROUS_BGK_DYNAMICS_HH + +#include "SmagorinskyPowerLawPorousBGKdynamics.h" +#include "SmagorinskyPorousParticleBGKdynamics.hh" +#include "math.h" + +namespace olb { + +////////////////////// Class SmagorinskyPowerLawPorousParticleBGKdynamics ////////////////////////// + +/** \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 +SmagorinskyPowerLawPorousParticleBGKdynamics::SmagorinskyPowerLawPorousParticleBGKdynamics ( + T omega_, Momenta& momenta_, T m_, T n_ , T dtPL_, T nuMin, T nuMax, T smagoConst_) + : SmagorinskyPorousParticleBGKdynamics(omega_,momenta_,smagoConst_), + m(m_), + n(n_), + dtPL(dtPL_) + //preFactor(computePreFactor(omega_,smagoConstPL_) ) +{ + omegaMin = 2./(nuMax*2.*descriptors::invCs2() + 1.); + omegaMax = 2./(nuMin*2.*descriptors::invCs2() + 1.); +} + +template +void SmagorinskyPowerLawPorousParticleBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell, rho, u, pi); + // load old omega from dyn. omega descriptor +// T oldOmega = this->getOmega(); //compute with constant omega + T oldOmega = cell.template getFieldPointer()[0]; //compute with dynamic omega + T OmegaPL = computeOmegaPL(oldOmega, rho, pi); + T* velDenominator = cell.template getFieldPointer(); + T* velNumerator = cell.template getFieldPointer(); + T* porosity = cell.template getFieldPointer(); + if (*velDenominator > std::numeric_limits::epsilon()) { + *porosity = 1.-*porosity; // 1-prod(1-smoothInd) + for (int i=0; i < DESCRIPTOR::d; i++) { + u[i] += *porosity * (*(velNumerator+i) / *velDenominator - u[i]); + } + } + + T newOmega = this->computeOmega(OmegaPL, this->preFactor, rho, pi); + + T uSqr = lbHelpers::bgkCollision(cell, rho, u, newOmega); + // save new omega to dyn. omega descriptor + cell.template getFieldPointer()[0] = newOmega; //compute with dynamic omega + statistics.incrementStats(rho, uSqr); + + cell.template defineField(1.0); + cell.template defineField(0.0); + cell.template defineField(0.0); +} + +template +T SmagorinskyPowerLawPorousParticleBGKdynamics::computeOmegaPL(T omega0, T rho, T pi[util::TensorVal::n] ) +{ + + // strain rate tensor without prefactor + T PiNeqNormSqr = pi[0]*pi[0] + 2.*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 pre2 = pow(descriptors::invCs2()/2./dtPL* omega0/rho,2.); // prefactor to the strain rate tensor + T D = pre2*PiNeqNormSqr; // Strain rate tensor + T gamma = sqrt(2.*D); // shear rate + + T nuNew = m*pow(gamma,n-1.); //nu for non-Newtonian fluid + //T newOmega = 2./(nuNew*6.+1.); + T newOmega = 2./(nuNew*2.*descriptors::invCs2() + 1.); + + /* + * problem if newOmega too small or too big is see e.g. "Daniel Conrad , Andreas Schneider, Martin Böhle: + * A viscosity adaption method for Lattice Boltzmann simulations" + */ + //if (newOmega>1.965) { + // newOmega = 1.965; //std::cout << newOmega << std::endl; + //} + //if (newOmega<0.1) { + // newOmega = 0.1; //std::cout << newOmega << std::endl; + //} + if (newOmega>omegaMax) { + newOmega = omegaMax; //std::cout << newOmega << std::endl; + } + if (newOmega + * + * 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 WALE_LATTICE_DESCRIPTOR_H +#define WALE_LATTICE_DESCRIPTOR_H + +#include "dynamics/latticeDescriptors.h" + + +namespace olb { + +namespace descriptors { + + + +///////////////////////////////////////////////////////////////////////////////// +// 3D Descriptors for flow with Wall Adaptive Local Eddy Viscosity (WALE) + +using WALED3Q19Descriptor = D3Q19; + +using WALED3Q27Descriptor = D3Q27; + +/// WALE 3D Forced + +using WALEForcedD3Q19Descriptor = D3Q19; + +using WALEForcedD3Q27Descriptor = D3Q27; + + +} // namespace descriptors + +} // namespace olb + +#endif diff --git a/src/dynamics/advectionDiffusionDynamics.h b/src/dynamics/advectionDiffusionDynamics.h new file mode 100644 index 0000000..0ea7142 --- /dev/null +++ b/src/dynamics/advectionDiffusionDynamics.h @@ -0,0 +1,201 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani + * 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 -- header file. + */ +#ifndef ADVECTION_DIFFUSION_DYNAMICS_H +#define ADVECTION_DIFFUSION_DYNAMICS_H + +#include "latticeDescriptors.h" +#include "dynamics/dynamics.h" +#include "core/unitConverter.h" + +namespace olb { + +// ========= the RLB advection diffusion dynamics ========// +/// it uses the regularized approximation that can be found in the thesis of J. Latt (2007). +template +class AdvectionDiffusionRLBdynamics : public BasicDynamics { +public: + /// Constructor + AdvectionDiffusionRLBdynamics( T omega_, Momenta& momenta_ ); + /// Compute equilibrium distribution function + T computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const override; + /// Collision step + void collide( Cell& cell, LatticeStatistics& statistics ) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega( T omega_ ) override; +private: + T omega; ///< relaxation parameter +}; + +/// Implementation of Regularized BGK collision, followed by any Dynamics +template +class CombinedAdvectionDiffusionRLBdynamics : public BasicDynamics { +public: + /// Constructor + CombinedAdvectionDiffusionRLBdynamics(T omega, Momenta& momenta); + /// Compute equilibrium distribution function + T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const override; + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +private: + Dynamics _boundaryDynamics; +}; + +// ========= the BGK advection diffusion dynamics ========// +/// This approach contains a slight error in the diffusion term. +template +class AdvectionDiffusionBGKdynamics : public BasicDynamics { +public: + /// Constructor + AdvectionDiffusionBGKdynamics( T omega, Momenta& momenta ); + AdvectionDiffusionBGKdynamics( const UnitConverter& converter, Momenta& momenta ); + /// Compute equilibrium distribution function + T computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const override; + /// Collision step + void collide( Cell& cell, LatticeStatistics& statistics ) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega( T omega ) override; +protected: + T _omega; ///< relaxation parameter +}; + + + +// ========= the TRT advection diffusion dynamics ========// +/// This approach contains a slight error in the diffusion term. +template +class AdvectionDiffusionTRTdynamics : public AdvectionDiffusionBGKdynamics { +public: + /// Constructor + AdvectionDiffusionTRTdynamics( T omega, Momenta& momenta, T magicParameter ); + /// Collision step + void collide( Cell& cell, LatticeStatistics& statistics ) override; +protected: + T _omega2; /// relaxation parameter for odd moments + T _magicParameter; +}; + + +// ======= BGK advection diffusion dynamics with source term ======// +// following Seta, T. (2013). Implicit temperature-correction-based +// immersed-boundary thermal lattice Boltzmann method for the simulation +// of natural convection. Physical Review E, 87(6), 063304. +template +class SourcedAdvectionDiffusionBGKdynamics : public AdvectionDiffusionBGKdynamics { +public: + /// Constructor + SourcedAdvectionDiffusionBGKdynamics(T omega_, Momenta& momenta_ ); + /// Collision step + virtual void collide(Cell& cell, LatticeStatistics& statistics ) override; + /// Compute Density + T computeRho(Cell const& cell) const override; + /// Compute fluid velocity and particle density on the cell. + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; +private: + const T _omegaMod; +}; + +// ========= the BGK advection diffusion Stokes drag dynamics with a Smagorinsky turbulence model ========// +/// This approach contains a slight error in the diffusion term. +template +class SmagorinskyParticleAdvectionDiffusionBGKdynamics : public olb::AdvectionDiffusionBGKdynamics { +public: + /// Constructor + SmagorinskyParticleAdvectionDiffusionBGKdynamics(T omega_, Momenta& momenta_, T smagoConst_, T dx_, T dt_); + /// Collision step + virtual void collide(Cell& cell, LatticeStatistics& statistics ); + /// Get local smagorinsky relaxation parameter of the dynamics + virtual T getSmagorinskyOmega(Cell& cell); + /// Set local relaxation parameter of the dynamics + virtual void setOmega(T omega_); + +private: + /// Computes a constant prefactor in order to speed up the computation + T computePreFactor(T omega_, T smagoConst_, T dx_, T dt_); + /// Computes the local smagorinsky relaxation parameter + T computeOmega(T omega0, T preFacto_r, T rho, T pi[util::TensorVal::n] ); + + /// effective collision time based upon Smagorisnky approach + T tau_eff; + /// Smagorinsky constant + T smagoConst; + /// Precomputed constant which speeeds up the computation + T preFactor; + T dx; + T dt; +}; + +// ========= the BGK advection diffusion Stokes drag dynamics ========// +/// This approach contains a slight error in the diffusion term. +template +class ParticleAdvectionDiffusionBGKdynamics : public olb::AdvectionDiffusionBGKdynamics { +public: + /// Constructor + ParticleAdvectionDiffusionBGKdynamics(T omega_, Momenta& momenta_); + /// Collision step + void collide(Cell& cell, LatticeStatistics& statistics ) override; +private: + T omega; ///< relaxation parameter +}; + + +// ========= the MRT advection diffusion dynamics ========// +/// This approach is based on the multi-distribution LBM model. +/// The couplintg is done using the Boussinesq approximation +template +class AdvectionDiffusionMRTdynamics : public BasicDynamics { +public: + /// Constructor + AdvectionDiffusionMRTdynamics( T omega, Momenta& momenta ); + /// Compute equilibrium distribution function + T computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const override; + /// Collision step + void collide( Cell& cell, LatticeStatistics& statistics ) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega( T omega ) override; +private: + T _omega; ///< relaxation parameter +protected: + T invM_S[DESCRIPTOR::q][DESCRIPTOR::q]; ///< inverse relaxation times matrix +}; + +} // namespace olb + +#endif diff --git a/src/dynamics/advectionDiffusionDynamics.hh b/src/dynamics/advectionDiffusionDynamics.hh new file mode 100644 index 0000000..54a3ab7 --- /dev/null +++ b/src/dynamics/advectionDiffusionDynamics.hh @@ -0,0 +1,442 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani + * 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 ADVECTION_DIFFUSION_DYNAMICS_HH +#define ADVECTION_DIFFUSION_DYNAMICS_HH + +#include +#include +#include "advectionDiffusionDynamics.h" +#include "dynamics/lbHelpers.h" + +namespace olb { + + +////////////////////// Class AdvectionDiffusionRLBdynamics ////////////////////////// + +/** \param omega_ relaxation parameter, related to the dynamic viscosity + * \param momenta_ a Momenta object to know how to compute velocity momenta + */ +//==================================================================// +//============= Regularized Model for Advection diffusion===========// +//==================================================================// + +template +AdvectionDiffusionRLBdynamics::AdvectionDiffusionRLBdynamics ( + T omega_, Momenta& momenta_ ) + : BasicDynamics( momenta_ ), + omega( omega_ ) +{ } + +template +T AdvectionDiffusionRLBdynamics::computeEquilibrium( int iPop, T rho, + const T u[DESCRIPTOR::d], T uSqr ) const +{ + return lbHelpers::equilibriumFirstOrder( iPop, rho, u ); +} + + +template +void AdvectionDiffusionRLBdynamics::collide( Cell& cell, + LatticeStatistics& statistics ) +{ + T temperature = this->_momenta.computeRho( cell ); + + const T* u = cell.template getFieldPointer(); + + T uSqr = lbHelpers:: + rlbCollision( cell, temperature, u, omega ); + + statistics.incrementStats( temperature, uSqr ); +} + +template +T AdvectionDiffusionRLBdynamics::getOmega() const +{ + return omega; +} + +template +void AdvectionDiffusionRLBdynamics::setOmega( T omega_ ) +{ + omega = omega_; +} + + +////////////////////// Class CombinedAdvectionDiffusionRLBdynamics ///////////////////////// + +template +CombinedAdvectionDiffusionRLBdynamics::CombinedAdvectionDiffusionRLBdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics(momenta), + _boundaryDynamics(omega, momenta) +{ } + +template +T CombinedAdvectionDiffusionRLBdynamics:: +computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return lbHelpers::equilibriumFirstOrder( iPop, rho, u ); +} + +template +void CombinedAdvectionDiffusionRLBdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + typedef DESCRIPTOR L; + + T temperature = this->_momenta.computeRho( cell ); + const T* u = cell.template getFieldPointer(); + + T jNeq[DESCRIPTOR::d]; + // this->_momenta.computeJ( cell, jNeq ); + dynamic_cast&>(this->_momenta).computeJneq( cell, jNeq ); + // cout << jNeq[0] << " " << jNeq[1] << " " << u[0] << " " << u[1] << endl; + // stripe of equilibrium part u * T + // for ( int iD = 0; iD < DESCRIPTOR::d; ++iD ) { + // jNeq[iD] -= u[iD] * temperature; + // } + // cout << jNeq[0] << " " << jNeq[1] << " " << u[0] << " " << u[1] << endl; + + for (int iPop = 0; iPop < L::q; ++iPop) { + cell[iPop] = lbHelpers::equilibriumFirstOrder( iPop, temperature, u ) + + firstOrderLbHelpers::fromJneqToFneq(iPop, jNeq); + // cout << firstOrderLbHelpers::fromJneqToFneq(iPop,jNeq) << " "; + } + // cout << endl; + // lbHelpers::computeJ(cell, jNeq); + // cout << jNeq[0] << " " << jNeq[1] << " " << u[0] << " " << u[1] << endl; + + _boundaryDynamics.collide(cell, statistics); + // lbHelpers::computeJ(cell, jNeq); + // cout << jNeq[0] << " " << jNeq[1] << " " << u[0] << " " << u[1] << endl; +} + +template +T CombinedAdvectionDiffusionRLBdynamics::getOmega() const +{ + return _boundaryDynamics.getOmega(); +} + +template +void CombinedAdvectionDiffusionRLBdynamics::setOmega(T omega) +{ + _boundaryDynamics.setOmega(omega); +} + +//==================================================================// +//============= BGK Model for Advection diffusion===========// +//==================================================================// + +template +AdvectionDiffusionBGKdynamics::AdvectionDiffusionBGKdynamics ( + T omega, Momenta& momenta ) + : BasicDynamics( momenta ), + _omega(omega) +{ } + +template +AdvectionDiffusionBGKdynamics::AdvectionDiffusionBGKdynamics ( + const UnitConverter& converter, Momenta& momenta ) + : BasicDynamics( momenta ), + _omega(converter.getLatticeRelaxationFrequency()) +{ } + +template +T AdvectionDiffusionBGKdynamics::computeEquilibrium( int iPop, T rho, + const T u[DESCRIPTOR::d], T uSqr ) const +{ + return lbHelpers::equilibriumFirstOrder( iPop, rho, u ); +} + + +template +void AdvectionDiffusionBGKdynamics::collide( Cell& cell, + LatticeStatistics& statistics ) +{ + T temperature = this->_momenta.computeRho( cell ); + const T* u = cell.template getFieldPointer(); + + T uSqr = lbHelpers:: + bgkCollision( cell, temperature, u, _omega ); + + statistics.incrementStats( temperature, uSqr ); +} + +template +T AdvectionDiffusionBGKdynamics::getOmega() const +{ + return _omega; +} + +template +void AdvectionDiffusionBGKdynamics::setOmega( T omega ) +{ + _omega = omega; +} + + +//==================================================================// +//============= TRT Model for Advection diffusion===========// +//==================================================================// + +template +AdvectionDiffusionTRTdynamics::AdvectionDiffusionTRTdynamics ( + T omega, Momenta& momenta, T magicParameter ) + : AdvectionDiffusionBGKdynamics( omega, momenta ), + _omega2(1/(magicParameter/(1/omega-0.5)+0.5)), _magicParameter(magicParameter) +{ } + +template +void AdvectionDiffusionTRTdynamics::collide( Cell& cell, + LatticeStatistics& statistics ) +{ + T temperature = this->_momenta.computeRho( cell ); + const T* u = cell.template getFieldPointer(); + + 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::equilibriumFirstOrder(iPop, temperature, u); + } + 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] -= _omega2 * (fPlus[