diff options
Diffstat (limited to 'src/functors/lattice/blockLatticeLocalF3D.hh')
-rw-r--r-- | src/functors/lattice/blockLatticeLocalF3D.hh | 1556 |
1 files changed, 1556 insertions, 0 deletions
diff --git a/src/functors/lattice/blockLatticeLocalF3D.hh b/src/functors/lattice/blockLatticeLocalF3D.hh new file mode 100644 index 0000000..6359b55 --- /dev/null +++ b/src/functors/lattice/blockLatticeLocalF3D.hh @@ -0,0 +1,1556 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012 Lukas Baron, Mathias J. Krause, Albert Mink + * 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 BLOCK_LATTICE_LOCAL_F_3D_HH +#define BLOCK_LATTICE_LOCAL_F_3D_HH + + +#include<cmath> +#include<math.h> + +#include "blockLatticeLocalF3D.h" +#include "blockBaseF3D.h" +#include "core/blockLatticeStructure3D.h" +#include "dynamics/lbHelpers.h" // for computation of lattice rho and velocity +#include "communication/mpiManager.h" +#include "utilities/vectorHelpers.h" + +namespace olb { + +template<typename T, typename DESCRIPTOR> +BlockLatticeFpop3D<T, DESCRIPTOR>::BlockLatticeFpop3D(BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, DESCRIPTOR::q) +{ + this->getName() = "fPop"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeFpop3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + output[iPop] = + this->_blockLattice.get(input[0], input[1], input[2])[iPop] + + descriptors::t<T,DESCRIPTOR>(iPop); + } + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeDissipation3D<T, DESCRIPTOR>::BlockLatticeDissipation3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1), _converter(converter) +{ + this->getName() = "dissipation"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeDissipation3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n]; + this->_blockLattice.get(input[0], input[1], input[2]).computeAllMomenta(rho, + uTemp, + pi); + + T PiNeqNormSqr = pi[0] * pi[0] + 2. * 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 nuLattice = _converter.getLatticeViscosity(); + T omega = 1. / _converter.getLatticeRelaxationTime(); + output[0] = PiNeqNormSqr * nuLattice + * pow(omega * descriptors::invCs2<T,DESCRIPTOR>(), 2) / rho / 2.; + + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticePhysDissipation3D<T, DESCRIPTOR>::BlockLatticePhysDissipation3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, + int overlap, + const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1), + _overlap(overlap), + _converter(converter) +{ + this->getName() = "physDissipation"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysDissipation3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n]; + this->_blockLattice.get( + input[0]+_overlap, input[1]+_overlap, input[2]+_overlap + ).computeAllMomenta(rho, uTemp, pi); + + T PiNeqNormSqr = pi[0] * pi[0] + 2. * 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 nuLattice = this->_converter.getLatticeViscosity(); + T omega = 1. / this->_converter.getLatticeRelaxationTime(); + T dt = this->_converter.getConversionFactorTime(); + output[0] = PiNeqNormSqr * nuLattice + * pow(omega * descriptors::invCs2<T,DESCRIPTOR>() / rho, 2) / 2. + * this->_converter.getPhysViscosity() / nuLattice / dt / dt; + + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeEffevtiveDissipation3D<T, DESCRIPTOR>::BlockLatticeEffevtiveDissipation3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter, T smagoConst, + LESDynamics<T, DESCRIPTOR>& LESdynamics) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1), + _converter(converter), _smagoConst(smagoConst), _LESdynamics(LESdynamics) +{ + this->getName() = "EffevtiveDissipation"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeEffevtiveDissipation3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n]; + this->_blockLattice.get(input[0], input[1], input[2]).computeAllMomenta(rho, + uTemp, + pi); + + T PiNeqNormSqr = pi[0] * pi[0] + 2. * 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 omegaEff = _LESdynamics.getEffectiveOmega(this->_blockLattice.get(input[0], input[1], input[2])); + T nuEff = ((1./omegaEff)-0.5)/descriptors::invCs2<T,DESCRIPTOR>(); + + output[0] = PiNeqNormSqr * (nuEff) + * pow(omegaEff * descriptors::invCs2<T,DESCRIPTOR>() / rho, 2) / 2.; + + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticePhysEffevtiveDissipation3D<T, DESCRIPTOR>::BlockLatticePhysEffevtiveDissipation3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter, T smagoConst, + LESDynamics<T, DESCRIPTOR>& LESdynamics) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1), + _converter(converter), _smagoConst(smagoConst), _LESdynamics(LESdynamics) +{ + this->getName() = "physEffevtiveDissipation"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysEffevtiveDissipation3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n]; + this->_blockLattice.get(input[0], input[1], input[2]).computeAllMomenta(rho, + uTemp, + pi); + + T PiNeqNormSqr = pi[0] * pi[0] + 2. * 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 dt = 1./ _converter.getConversionFactorTime(); + T omegaEff = _LESdynamics.getEffectiveOmega(this->_blockLattice.get(input[0], input[1], input[2])); + T nuEff = ((1./omegaEff)-0.5)/descriptors::invCs2<T,DESCRIPTOR>(); // BGK shear viscosity definition + + output[0] = PiNeqNormSqr * nuEff + * pow(omegaEff * descriptors::invCs2<T,DESCRIPTOR>() / rho, 2) / 2. + * _converter.getPhysViscosity() / _converter.getLatticeViscosity() / dt / dt; + + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeDensity3D<T, DESCRIPTOR>::BlockLatticeDensity3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1) +{ + this->getName() = "density"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeDensity3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + output[0] = this->_blockLattice.get(input[0], input[1], input[2]).computeRho(); + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeVelocity3D<T, DESCRIPTOR>::BlockLatticeVelocity3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3) +{ + this->getName() = "velocity"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeVelocity3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T rho; + this->_blockLattice.get(input[0], input[1], input[2]).computeRhoU(rho, output); + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeExternalVelocity3D<T, DESCRIPTOR>::BlockLatticeExternalVelocity3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3) +{ + this->getName() = "externalVelocity"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeExternalVelocity3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T* ExtVel = this->_blockLattice.get(input[0], input[1], input[2]).template getFieldPointer<descriptors::VELOCITY>(); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + output[iVel] = ExtVel[iVel]; + } + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeFlux3D<T, DESCRIPTOR>::BlockLatticeFlux3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3) +{ + this->getName() = "flux"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeFlux3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + this->_blockLattice.get(input[0], input[1], input[2]).computeJ(output); + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeStrainRate3D<T, DESCRIPTOR>::BlockLatticeStrainRate3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 9) +{ + this->getName() = "strainRate"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeStrainRate3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n]; + this->_blockLattice.get(input[0], input[1], input[2]).computeAllMomenta(rho, + uTemp, + pi); + + T omega = 1. / this->_converter.getLatticeRelaxationTime(); + + output[0] = -pi[0] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + output[1] = -pi[1] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + output[2] = -pi[2] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + output[3] = -pi[1] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + output[4] = -pi[3] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + output[5] = -pi[4] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + output[6] = -pi[2] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + output[7] = -pi[4] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + output[8] = -pi[5] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.; + + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticePhysStrainRate3D<T, DESCRIPTOR>::BlockLatticePhysStrainRate3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, + int overlap, + const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 9), + _overlap(overlap) +{ + this->getName() = "strainRate"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysStrainRate3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n]; + this->_blockLattice.get( + input[0]+_overlap, input[1]+_overlap, input[2]+_overlap + ).computeAllMomenta(rho, uTemp, pi); + + T omega = 1. / this->_converter.getLatticeRelaxationTime(); + T dt = this->_converter.getConversionFactorTime(); + + output[0] = -pi[0] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + output[1] = -pi[1] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + output[2] = -pi[2] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + output[3] = -pi[1] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + output[4] = -pi[3] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + output[5] = -pi[4] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + output[6] = -pi[2] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + output[7] = -pi[4] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + output[8] = -pi[5] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt; + + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeGeometry3D<T, DESCRIPTOR>::BlockLatticeGeometry3D(BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, + BlockGeometryStructure3D<T>& blockGeometry, int material) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1), + _blockGeometry(blockGeometry), + _material(material) +{ + this->getName() = "geometry"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeGeometry3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + output[0] = _blockGeometry.getMaterial(input[0], input[1], input[2]); + + if (_material != -1) { + if ( util::nearZero(_material-output[0]) ) { + output[0] = 1.; + return true; + } + else { + output[0] = 0.; + return true; + } + } + return true; +} + + +template<typename T, typename DESCRIPTOR> +BlockLatticeRank3D<T,DESCRIPTOR>::BlockLatticeRank3D( + BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice) + : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 1) +{ + this->getName() = "rank"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeRank3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + output[0] = singleton::mpi().getRank() + 1; + return true; +} + + +template<typename T, typename DESCRIPTOR> +BlockLatticeCuboid3D<T,DESCRIPTOR>::BlockLatticeCuboid3D( + BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, int iC) + : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 1), _iC(iC) +{ + this->getName() = "cuboid"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeCuboid3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + output[0] = _iC + 1; + return true; +} + + + +template<typename T, typename DESCRIPTOR> +BlockLatticePhysPressure3D<T, DESCRIPTOR>::BlockLatticePhysPressure3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, + int overlap, + const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 1), + _overlap(overlap) +{ + this->getName() = "physPressure"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysPressure3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + // lattice pressure = c_s^2 ( rho -1 ) + T latticePressure = ( this->_blockLattice.get(input[0]+_overlap, input[1]+_overlap, input[2]+_overlap).computeRho() - 1.0) / descriptors::invCs2<T,DESCRIPTOR>(); + output[0] = this->_converter.getPhysPressure(latticePressure); + + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticePhysVelocity3D<T, DESCRIPTOR>::BlockLatticePhysVelocity3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, + int overlap, + const UnitConverter<T,DESCRIPTOR>& converter, + bool print) + : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 3), + _overlap(overlap), + _print(print) +{ + this->getName() = "physVelocity"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysVelocity3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + if (_print) { + std::cout << input[0] << " " << input[1] << " " << input[2] << " | " + << singleton::mpi().getRank() << std::endl; + } + + T rho; + this->_blockLattice.get( + input[0]+_overlap, input[1]+_overlap, input[2]+_overlap).computeRhoU(rho, output); + output[0] = this->_converter.getPhysVelocity(output[0]); + output[1] = this->_converter.getPhysVelocity(output[1]); + output[2] = this->_converter.getPhysVelocity(output[2]); + + return true; +} + +template <typename T, typename DESCRIPTOR> +BlockLatticePhysExternalPorosity3D<T,DESCRIPTOR>::BlockLatticePhysExternalPorosity3D( + BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, + int overlap, + const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice,converter,2), + _overlap(overlap) +{ + this->getName() = "ExtPorosityField"; +} + +template <typename T, typename DESCRIPTOR> +bool BlockLatticePhysExternalPorosity3D<T,DESCRIPTOR>::operator() (T output[], const int input[]) +{ + this->_blockLattice.get( + input[0]+_overlap, input[1]+_overlap, input[2]+_overlap + ).template computeField<descriptors::POROSITY>(output); + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticePhysExternalVelocity3D<T, DESCRIPTOR>::BlockLatticePhysExternalVelocity3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 3) +{ + this->getName() = "physVelExtField"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysExternalVelocity3D<T, DESCRIPTOR>::operator()( + T output[], const int input[]) +{ + this->_blockLattice.get(input[0], input[1], input[2]).template computeField<descriptors::VELOCITY>(output); + output[0] = this->_converter.getPhysVelocity(output[0]); + output[1] = this->_converter.getPhysVelocity(output[1]); + output[2] = this->_converter.getPhysVelocity(output[2]); + return true; +} + +template <typename T, typename DESCRIPTOR> +BlockLatticePhysExternalParticleVelocity3D<T,DESCRIPTOR>::BlockLatticePhysExternalParticleVelocity3D +(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice,converter,2) +{ + this->getName() = "ExtParticleVelocityField"; +} + +template <typename T, typename DESCRIPTOR> +bool BlockLatticePhysExternalParticleVelocity3D<T,DESCRIPTOR>::operator() (T output[], const int input[]) +{ + const T* velocity_numerator = this->blockLattice.get(input).template getFieldPointer<descriptors::VELOCITY_NUMERATOR>(); + const T* velocity_denominator = this->blockLattice.get(input).template getFieldPointer<descriptors::VELOCITY_DENOMINATOR>(); + + if (velocity_denominator[0] > std::numeric_limits<T>::epsilon()) { + output[0]=this->_converter.getPhysVelocity(velocity_numerator[0]/velocity_denominator[0]); + output[1]=this->_converter.getPhysVelocity(velocity_numerator[1]/velocity_denominator[0]); + output[2]=this->_converter.getPhysVelocity(velocity_numerator[2]/velocity_denominator[0]); + return true; + } + output[0]=this->_converter.getPhysVelocity(velocity_numerator[0]); + output[1]=this->_converter.getPhysVelocity(velocity_numerator[1]); + output[2]=this->_converter.getPhysVelocity(velocity_numerator[2]); + return true; +} + + +template<typename T, typename DESCRIPTOR> +BlockLatticePhysExternal3D<T, DESCRIPTOR>::BlockLatticePhysExternal3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, + T convFactorToPhysUnits, int offset, int size) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3), + _convFactorToPhysUnits(convFactorToPhysUnits), + _offset(offset), _size(size) +{ + this->getName() = "physExtField"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysExternal3D<T, DESCRIPTOR>::operator()( + T output[], const int input[]) +{ + const auto& cell = this->_blockLattice.get(input[0], input[1], input[2]); + output[0] = _convFactorToPhysUnits*cell[_offset+0]; + output[1] = _convFactorToPhysUnits*cell[_offset+1]; + output[2] = _convFactorToPhysUnits*cell[_offset+2]; + return true; +} + +template <typename T, typename DESCRIPTOR> +BlockLatticePhysBoundaryForce3D<T,DESCRIPTOR>::BlockLatticePhysBoundaryForce3D( + BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, + BlockIndicatorF3D<T>& indicatorF, + const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice, converter, 3), + _indicatorF(indicatorF), + _blockGeometry(indicatorF.getBlockGeometryStructure()) +{ + this->getName() = "physBoundaryForce"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysBoundaryForce3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] = T(); + } + + if (_indicatorF(input)) { + for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) { + // Get direction + const Vector<int,3> c = descriptors::c<DESCRIPTOR>(iPop); + // Get next cell located in the current direction + // Check if the next cell is a fluid node + if (_blockGeometry.get(input[0] + c[0], input[1] + c[1], input[2] + c[2]) == 1) { + // Get f_q of next fluid cell where l = opposite(q) + T f = this->_blockLattice.get(input[0] + c[0], input[1] + c[1], input[2] + c[2])[iPop]; + // Get f_l of the boundary cell + // Add f_q and f_opp + f += this->_blockLattice.get(input)[util::opposite<DESCRIPTOR>(iPop)]; + // Update force + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] -= c[i] * f; + } + } + } + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] = this->_converter.getPhysForce(output[i]); + } + } + return true; +} + +template <typename T, typename DESCRIPTOR> +BlockLatticePSMPhysForce3D<T,DESCRIPTOR>::BlockLatticePSMPhysForce3D( + BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, + const UnitConverter<T,DESCRIPTOR>& converter, + int mode_) + : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice, converter, 3) +{ + this->getName() = "physPSMForce"; + mode = (Mode) mode_; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePSMPhysForce3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] = T(); + } + + T epsilon = 1. - *(this->_blockLattice.get(input).template getFieldPointer<descriptors::POROSITY>()); + + //if ((epsilon > 1e-5 && epsilon < 1 - 1e-5)) { + if ((epsilon > 1e-5)) { + T rho, u[DESCRIPTOR::d], u_s[DESCRIPTOR::d]; + + for (int i = 0; i < DESCRIPTOR::d; i++) { + u_s[i] = (this->_blockLattice.get(input).template getFieldPointer<descriptors::VELOCITY_SOLID>())[i]; + } + T paramA = this->_converter.getLatticeRelaxationTime() - 0.5; + // speed up paramB + T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA); + + T omega_s; + T omega = 1. / this->_converter.getLatticeRelaxationTime(); + + this->_blockLattice.get(input).computeRhoU(rho, u); + + const T uSqr_s = util::normSqr<T,DESCRIPTOR::d>(u_s); + T uSqr = util::normSqr<T,DESCRIPTOR::d>(u); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + switch (mode) { + case M2: + omega_s = (lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u_s, uSqr_s) + - this->_blockLattice.get(input[0], input[1], input[2])[iPop]) + + (1 - omega) + * (this->_blockLattice.get(input[0], input[1], input[2])[iPop] + - lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr)); + break; + case M3: + omega_s = + (this->_blockLattice.get(input[0], input[1], input[2])[descriptors::opposite<DESCRIPTOR>(iPop)] + - lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium( + descriptors::opposite<DESCRIPTOR>(iPop), rho, u_s, uSqr_s)) + - (this->_blockLattice.get(input[0], input[1], input[2])[iPop] + - lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u_s, uSqr_s)); + } + + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] -= descriptors::c<DESCRIPTOR>(iPop,i) * omega_s; + } + } + + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] = this->_converter.getPhysForce(output[i] * paramB); + } + } + return true; +} + +template <typename T, typename DESCRIPTOR> +BlockLatticePhysWallShearStress3D<T,DESCRIPTOR>::BlockLatticePhysWallShearStress3D( + BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, + BlockGeometryStructure3D<T>& blockGeometry, + int overlap, + int material, + const UnitConverter<T,DESCRIPTOR>& converter, + IndicatorF3D<T>& indicator) + : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice,converter,1), + _blockGeometry(blockGeometry), + _overlap(overlap), + _material(material) + { + this->getName() = "physWallShearStress"; + const T scaling = this->_converter.getConversionFactorLength() * 0.1; + const T omega = 1. / this->_converter.getLatticeRelaxationTime(); + const T dt = this->_converter.getConversionFactorTime(); + _physFactor = -omega * descriptors::invCs2<T,DESCRIPTOR>() / dt * this->_converter.getPhysDensity() * this->_converter.getPhysViscosity(); + std::vector<int> discreteNormalOutwards(4, 0); + + for (int iX = 1 ; iX < _blockGeometry.getNx() - 1; iX++) { + _discreteNormal.resize(_blockGeometry.getNx() - 2); + _normal.resize(_blockGeometry.getNx() - 2); + + for (int iY = 1; iY < _blockGeometry.getNy() - 1; iY++) { + _discreteNormal[iX-1].resize(_blockGeometry.getNy() - 2); + _normal[iX-1].resize(_blockGeometry.getNy() - 2); + + for (int iZ = 1; iZ < _blockGeometry.getNz() - 1; iZ++) { + _discreteNormal[iX-1][iY-1].resize(_blockGeometry.getNz() - 2); + _normal[iX-1][iY-1].resize(_blockGeometry.getNz() - 2); + + if (_blockGeometry.get(iX, iY, iZ) == _material) { + discreteNormalOutwards = _blockGeometry.getStatistics().getType(iX, iY, iZ); + _discreteNormal[iX - 1][iY - 1][iZ - 1].resize(3); + _normal[iX - 1][iY - 1][iZ - 1].resize(3); + + _discreteNormal[iX - 1][iY- 1][iZ- 1][0] = -discreteNormalOutwards[1]; + _discreteNormal[iX- 1][iY- 1][iZ- 1][1] = -discreteNormalOutwards[2]; + _discreteNormal[iX- 1][iY- 1][iZ- 1][2] = -discreteNormalOutwards[3]; + + T physR[3]; + _blockGeometry.getPhysR(physR,iX, iY, iZ); + Vector<T,3> origin(physR[0],physR[1],physR[2]); + Vector<T,3> direction(-_discreteNormal[iX- 1][iY- 1][iZ- 1][0] * scaling, + -_discreteNormal[iX- 1][iY- 1][iZ- 1][1] * scaling, + -_discreteNormal[iX- 1][iY- 1][iZ- 1][2] * scaling); + Vector<T,3> normal(0.,0.,0.); + origin[0] = physR[0]; + origin[1] = physR[1]; + origin[2] = physR[2]; + + indicator.normal(normal, origin, direction); + normal.normalize(); + + _normal[iX- 1][iY- 1][iZ- 1][0] = normal[0]; + _normal[iX- 1][iY- 1][iZ- 1][1] = normal[1]; + _normal[iX- 1][iY- 1][iZ- 1][2] = normal[2]; + } + } + } + } +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysWallShearStress3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + output[0] = T(); + + if (input[0] + _overlap < 1 || + input[1] + _overlap < 1 || + input[2] + _overlap < 1 || + input[0] + _overlap >= _blockGeometry.getNx()-1 || + input[1] + _overlap >= _blockGeometry.getNy()-1 || + input[2] + _overlap >= _blockGeometry.getNz()-1 ){ + +#ifdef OLB_DEBUG + std::cout << "Input address not mapped by _discreteNormal, overlap too small" << std::endl; +#endif + return true; + } + + if (_blockGeometry.get(input[0]+_overlap,input[1]+_overlap,input[2]+_overlap) == _material) { + + T traction[3]; + T stress[6]; + T rho = this->_blockLattice.get(input[0] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0], + input[1] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1], + input[2] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]).computeRho(); + this->_blockLattice.get(input[0] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0], + input[1] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1], + input[2] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]).computeStress(stress); + + traction[0] = stress[0]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0] + + stress[1]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1] + + stress[2]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]; + traction[1] = stress[1]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0] + + stress[3]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1] + + stress[4]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]; + traction[2] = stress[2]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0] + + stress[4]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1] + + stress[5]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]; + + T traction_normal_SP; + T tractionNormalComponent[3]; + // scalar product of traction and normal vector + traction_normal_SP = traction[0] * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0] + + traction[1] * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1] + + traction[2] * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]; + tractionNormalComponent[0] = traction_normal_SP * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0]; + tractionNormalComponent[1] = traction_normal_SP * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1]; + tractionNormalComponent[2] = traction_normal_SP * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]; + + T WSS[3]; + WSS[0] = traction[0] - tractionNormalComponent[0]; + WSS[1] = traction[1] - tractionNormalComponent[1]; + WSS[2] = traction[2] - tractionNormalComponent[2]; + // magnitude of the wall shear stress vector + output[0] = sqrt(WSS[0]*WSS[0] + WSS[1]*WSS[1] + WSS[2]*WSS[2]); + + return true; + } + else { + return true; + } +} + + +template <typename T, typename DESCRIPTOR> +BlockLatticePhysCorrBoundaryForce3D<T,DESCRIPTOR>::BlockLatticePhysCorrBoundaryForce3D( + BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, + BlockIndicatorF3D<T>& indicatorF, + const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice, converter, 3), + _indicatorF(indicatorF), + _blockGeometry(indicatorF.getBlockGeometryStructure()) +{ + this->getName() = "physCorrBoundaryForce"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysCorrBoundaryForce3D<T, DESCRIPTOR>::operator()( + T output[], const int input[]) +{ + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] = T(); + } + + if (_indicatorF(input)) { + for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) { + // Get direction + const Vector<int,3> c = descriptors::c<DESCRIPTOR>(iPop); + // Get next cell located in the current direction + // Check if the next cell is a fluid node + if (_blockGeometry.get(input[0] + c[0], input[1] + c[1], input[2] + c[2]) == 1) { + // Get f_q of next fluid cell where l = opposite(q) + T f = this->_blockLattice.get(input[0] + c[0], input[1] + c[1], input[2] + c[2])[iPop]; + // Get f_l of the boundary cell + // Add f_q and f_opp + f += this->_blockLattice.get(input)[util::opposite<DESCRIPTOR>(iPop)]; + // Update force + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] -= c[i] * (f - 2. * descriptors::t<T,DESCRIPTOR>(iPop)); + } + } + } + for (int i = 0; i < this->getTargetDim(); ++i) { + output[i] = this->_converter.getPhysForce(output[i]); + } + } + return true; +} + +template<typename T, typename DESCRIPTOR, typename FIELD> +BlockLatticeField3D<T,DESCRIPTOR,FIELD>::BlockLatticeField3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, DESCRIPTOR::template size<FIELD>()) +{ + this->getName() = "externalField"; +} + +template<typename T, typename DESCRIPTOR, typename FIELD> +bool BlockLatticeField3D<T,DESCRIPTOR,FIELD>::operator()(T output[], const int input[]) +{ + this->_blockLattice.get(input[0], input[1], input[2]).template computeField<FIELD>(output); + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticePorosity3D<T, DESCRIPTOR>::BlockLatticePorosity3D(BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1) +{ + this->getName() = "porosity"; +} + +// under construction +template<typename T, typename DESCRIPTOR> +bool BlockLatticePorosity3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ +#ifndef excludeDualDynamics + this->_blockLattice.get(input[0], input[1], input[2]).template computeField<descriptors::POROSITY>( + output); +#else + this->_blockLattice.get(input[0], input[1], input[2]).template computeField<descriptors::POROSITY>( + output); + +#endif + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticeVolumeFractionApproximation3D<T, DESCRIPTOR>::BlockLatticeVolumeFractionApproximation3D(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, + BlockGeometryStructure3D<T>& blockGeometry, + IndicatorF3D<T>& indicator, + int refinementLevel, + const UnitConverter<T,DESCRIPTOR>& converter, + bool insideOut) + : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1), + _blockGeometry(blockGeometry), _indicator(indicator), _refinementLevel(refinementLevel), _converter(converter), _insideOut(insideOut), + _physSubGridMinPhysRshift((1./ T(_refinementLevel * 2.) - 0.5) * _converter.getPhysDeltaX()), + _physSubGridDeltaX(1. / T(_refinementLevel) * _converter.getPhysDeltaX()), + _latticeSubGridVolume(1. / (_refinementLevel * _refinementLevel * _refinementLevel)) +{ + this->getName() = "volumeFractionApproximation"; +} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticeVolumeFractionApproximation3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + output[0] = 0.; + T physR[3]; + bool inside[1]; + _blockGeometry.getPhysR(physR, input[0], input[1], input[2]); + + T subGridMinPhysR[3]; + subGridMinPhysR[0] = physR[0] + _physSubGridMinPhysRshift; + subGridMinPhysR[1] = physR[1] + _physSubGridMinPhysRshift; + subGridMinPhysR[2] = physR[2] + _physSubGridMinPhysRshift; + + for (int x = 0; x < _refinementLevel; x++) { + for (int y = 0; y < _refinementLevel; y++) { + for (int z = 0; z < _refinementLevel; z++) { + physR[0] = subGridMinPhysR[0] + x * _physSubGridDeltaX; + physR[1] = subGridMinPhysR[1] + y * _physSubGridDeltaX; + physR[2] = subGridMinPhysR[2] + z * _physSubGridDeltaX; + _indicator(inside, physR); + if (!_insideOut) { + if (inside[0]) { + output[0] += _latticeSubGridVolume; + } + } + else { + if (!inside[0]) { + output[0] += _latticeSubGridVolume; + } + } + } + } + } + return true; +} + +template<typename T, typename DESCRIPTOR> +BlockLatticePhysPermeability3D<T, DESCRIPTOR>::BlockLatticePhysPermeability3D( + BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter) + : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 1) +{ + this->getName() = "permeability"; +} + +//TODO: consistency with 2D (181219) +//template<typename T, typename DESCRIPTOR> +//BlockLatticePhysPermeability3D<T, DESCRIPTOR>::BlockLatticePhysPermeability3D( +// BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, +// BlockGeometry3D<T>& blockGeometry, int material, +// const UnitConverter<T>& converter) +// : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 1), +// _blockGeometry(blockGeometry), +// _material(material) +//{ +// this->getName() = "permeability"; +//} + +template<typename T, typename DESCRIPTOR> +bool BlockLatticePhysPermeability3D<T, DESCRIPTOR>::operator()(T output[], const int input[]) +{ + T value; |