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 --- .../twoWayCouplings/twoWayHelperFunctionals.hh | 431 +++++++++++++++++++++ 1 file changed, 431 insertions(+) create mode 100644 src/particles/twoWayCouplings/twoWayHelperFunctionals.hh (limited to 'src/particles/twoWayCouplings/twoWayHelperFunctionals.hh') diff --git a/src/particles/twoWayCouplings/twoWayHelperFunctionals.hh b/src/particles/twoWayCouplings/twoWayHelperFunctionals.hh new file mode 100644 index 0000000..b7db1b3 --- /dev/null +++ b/src/particles/twoWayCouplings/twoWayHelperFunctionals.hh @@ -0,0 +1,431 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2019 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. + */ + +/* Helper functionals for Lagrangian two-way coupling methods -- generic implementation. + */ + +#ifndef LB_TWO_WAY_HELPER_FUNCTIONALS_HH +#define LB_TWO_WAY_HELPER_FUNCTIONALS_HH + +namespace olb { + +////////////////////// Class ParticleReynoldsNumberBase //////////////////////// + +template class Particle> +ParticleReynoldsNumberBase::ParticleReynoldsNumberBase(UnitConverter& converter) + : _converter(converter) +{} + + +////////////////////// Class NewtonianParticleReynoldsNumber //////////////////////// + +template class Particle> +NewtonianParticleReynoldsNumber::NewtonianParticleReynoldsNumber(UnitConverter& converter) + : ParticleReynoldsNumberBase(converter) +{} + +template class Particle> +T NewtonianParticleReynoldsNumber::operator() ( Particle* p, T magU, int globicFull[]) +{ + T ReP = 2. * p->getRad() * magU / this->_converter.getPhysViscosity(); + return ReP > this->_RePmin ? ReP : this->_RePmin; +} + + +////////////////////// Class PowerLawParticleReynoldsNumber //////////////////////// + +template class Particle> +PowerLawParticleReynoldsNumber::PowerLawParticleReynoldsNumber ( + UnitConverter& converter, SuperLattice3D& sLattice ) + : ParticleReynoldsNumberBase(converter), + _sLattice(sLattice) +{} + +template class Particle> +T PowerLawParticleReynoldsNumber::operator() ( Particle* p, T magU, int globicFull[]) +{ + // loc() indicates the local cuboid number locIC of the actual processing thread, + // for given global cuboid number iC + // this is to get appropriate particle system on locIC + int locIC = _sLattice.getLoadBalancer().loc(globicFull[0]); + + T physPosP[3] = {T(), T(), T()}; // particle's physical position + physPosP[0] = (p->getPos()[0]); + physPosP[1] = (p->getPos()[1]); + physPosP[2] = (p->getPos()[2]); + + // particle's dimensionless position, rounded at neighbouring voxel + int latticeRoundedPosP[3] = { globicFull[1], globicFull[2], globicFull[3] }; + + // getting the power-law relaxation frequency form the dynamics's external field + T omega = _sLattice.getBlockLattice(locIC).get ( + latticeRoundedPosP[0], + latticeRoundedPosP[1], + latticeRoundedPosP[2] ).template getFieldPointer()[0]; + + // physical viscosity from relaxation time + T nu = this->_converter.getPhysViscosity ( + (1./omega - 0.5) / descriptors::invCs2() ); + + T ReP = 2. * p->getRad() * magU / nu; + return ReP > this->_RePmin ? ReP : this->_RePmin; +} + + +////////////////////// Class TwoWayHelperFunctional //////////////////////// + +template +TwoWayHelperFunctional::TwoWayHelperFunctional ( + UnitConverter& converter, + SuperLattice3D& sLattice ) + : _converter(converter), + _sLattice(sLattice) +{} + +template +TwoWayHelperFunctional::~TwoWayHelperFunctional() +{ + _interpLatticeDensity.reset(); + _interpLatticeVelocity.reset(); +} + + +////////////////////// Class NaiveMomentumExchange //////////////////////// + +template +NaiveMomentumExchange::NaiveMomentumExchange ( + UnitConverter& converter, + SuperLattice3D& sLattice, + std::shared_ptr > interpLatticeDensity ) + : TwoWayHelperFunctional(converter, sLattice) +{ + this->_interpLatticeDensity = interpLatticeDensity; +} + +template +bool NaiveMomentumExchange::operator() ( T gF[], T latticeVelF[], T latticeVelP[], + T physPosP[], int latticeRoundedP[], + int globic ) +{ + T magU = sqrt( pow(latticeVelF[0] - latticeVelP[0],2) + + pow(latticeVelF[1] - latticeVelP[1],2) + + pow(latticeVelF[2] - latticeVelP[2],2) ); + + // Interpolated/exact pdf (depending on the functional) + T f[Lattice::q] = { T() }; + this->_interpLatticeDensity->operator()(f, physPosP, globic); + + // Mock cell with the same dynamics of the cell containing the particle + int locIC = this->_sLattice.getLoadBalancer().loc(globic); + Cell cell ( this->_sLattice.getBlockLattice(locIC).get ( + latticeRoundedP[0], + latticeRoundedP[1], + latticeRoundedP[2] ).getDynamics() ); + + // Filling the mock cell with the interpolated pdf + for (unsigned iPop = 0; iPop < Lattice::q; ++iPop) { + cell[iPop] = f[iPop]; + } + + T rhoL = cell.computeRho(); + + gF[0] = rhoL * magU; + gF[1] = rhoL * magU; + gF[2] = rhoL * magU; + + return true; +} + + +////////////////////// Class LaddMomentumExchange //////////////////////// + +template +LaddMomentumExchange::LaddMomentumExchange ( + UnitConverter& converter, + SuperLattice3D& sLattice, + std::shared_ptr > interpLatticeDensity, + std::shared_ptr > interpLatticeVelocity ) + : TwoWayHelperFunctional(converter, sLattice) +{ + this->_interpLatticeDensity = interpLatticeDensity; + this->_interpLatticeVelocity = interpLatticeVelocity; +} + +template +bool LaddMomentumExchange::operator() ( T gF[], T latticeVelF[], T latticeVelP[], + T physPosP[], int latticeRoundedP[], + int globic ) +{ + T physLatticeL = this->_converter.getConversionFactorLength(); + + // force density gF + gF[0] = T(); + gF[1] = T(); + gF[2] = T(); + + T fiPop = T(); + T sp = T(); // dot product for particle velocity + T faPos[3] = {T(), T(), T()}; // fAlphaPosition = particle position + T fbPos[3] = {T(), T(), T()}; // fBetaPosition = neighbor position to particle position in direction iPop + + T fa[Lattice::q] = { T() }; // fAlpha = interpolated density to fAlphaPosition + T fb[Lattice::q] = { T() }; // fBeta = interpolated density to fBetaPosition + T lFU[3] = {T(), T(), T()}; + + // runs through all q discrete velocity directions + for (unsigned iPop = 0; iPop < Lattice::q; ++iPop) { + // physical position on neighbor to get pre-streaming collision part + faPos[0] = physPosP[0] + physLatticeL * descriptors::c(iPop,0); + faPos[1] = physPosP[1] + physLatticeL * descriptors::c(iPop,1); + faPos[2] = physPosP[2] + physLatticeL * descriptors::c(iPop,2); + // Lagrange interpolated polynomial to get density on particle position + this->_interpLatticeDensity->operator() (fa, faPos, globic); + + // physical position on neighbor to get pre-streaming collision part + fbPos[0] = physPosP[0] - physLatticeL * descriptors::c(iPop,0); + fbPos[1] = physPosP[1] - physLatticeL * descriptors::c(iPop,1); + fbPos[2] = physPosP[2] - physLatticeL * descriptors::c(iPop,2); + // Lagrange interpolated polynomial to get density on particle position + this->_interpLatticeDensity->operator() (fb, fbPos, globic); + + // fiPop = density on fBetaPosition in direction iPop + fiPop = fb[util::opposite(iPop)]; + // Get f_l of the boundary cell + // add density fAlphaL of opposite direction to iPop + fiPop -= fa[iPop]; + + // physical velocity + lFU[0] = -descriptors::c(iPop,0) * fiPop; + lFU[1] = -descriptors::c(iPop,1) * fiPop; + lFU[2] = -descriptors::c(iPop,2) * fiPop; + + // point product + sp = descriptors::c(iPop,0) * latticeVelP[0] + descriptors::c(iPop,1) * latticeVelP[1] + + descriptors::c(iPop,2) * latticeVelP[2]; + sp *= 2. * descriptors::invCs2() + * descriptors::t(iPop); + + // external force density that acts on particles + gF[0] += (lFU[0] - descriptors::c(iPop,0) * (sp)); + gF[1] += (lFU[1] - descriptors::c(iPop,1) * (sp)); + gF[2] += (lFU[2] - descriptors::c(iPop,2) * (sp)); + } + gF[0] = fabs(gF[0]); + gF[1] = fabs(gF[1]); + gF[2] = fabs(gF[2]); + + return true; +} + + +////////////////////// Class SmoothingFunctional //////////////////////// + +template +SmoothingFunctional::SmoothingFunctional ( + T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice ) + : _kernelLength(kernelLength), + _converter(converter), + _sLattice(sLattice) +{} + +template +int SmoothingFunctional::getSize() +{ + return _latticePosAndWeight.size(); +} + +template +void SmoothingFunctional::getLatticePos(int latticePos[], int i) +{ + latticePos[0] = _latticePosAndWeight[i].latticePos[0]; + latticePos[1] = _latticePosAndWeight[i].latticePos[1]; + latticePos[2] = _latticePosAndWeight[i].latticePos[2]; +} + +template +int SmoothingFunctional::getGlobic(int i) +{ + return _latticePosAndWeight[i].globic; +} + +template +T SmoothingFunctional::getWeight(int i) +{ + return _latticePosAndWeight[i].weight; +} + +template +bool SmoothingFunctional::update(T physPosP[], int globic) +{ + // Bottom-left corner of a cube centered at the particle, with side 2*_kernelLength + T physPosMin[3] = {T(), T(), T()}; + physPosMin[0] = physPosP[0] - _kernelLength; + physPosMin[1] = physPosP[1] - _kernelLength; + physPosMin[2] = physPosP[2] - _kernelLength; + int latticePosMin[3] = {0, 0, 0}; + this->_sLattice.getCuboidGeometry().get(globic).getLatticeR ( + latticePosMin, physPosMin ); + + // Top-right corner of a cube centered at the particle, with side 2*_kernelLength + T physPosMax[3] = {T(), T(), T()}; + physPosMax[0] = physPosP[0] + _kernelLength; + physPosMax[1] = physPosP[1] + _kernelLength; + physPosMax[2] = physPosP[2] + _kernelLength; + int latticePosMax[3] = {0, 0, 0}; + this->_sLattice.getCuboidGeometry().get(globic).getLatticeR ( + latticePosMax, physPosMax ); + + // Clearing the _latticePosAndWeight list + _latticePosAndWeight.clear(); + + T normalizer = T(); + int iLatticePos[3] = {0, 0, 0}; + // Cycling all the cells on a cube containing a sphee centered in bubble's position and with kernel smoothing length as radius + for (iLatticePos[0]=latticePosMin[0]; iLatticePos[0]<=latticePosMax[0]; iLatticePos[0]++) { + for (iLatticePos[1]=latticePosMin[1]; iLatticePos[1]<=latticePosMax[1]; iLatticePos[1]++) { + for (iLatticePos[2]=latticePosMin[2]; iLatticePos[2]<=latticePosMax[2]; iLatticePos[2]++) { + + T iPhysPos[3] = {T(), T(), T()}; + this->_sLattice.getCuboidGeometry().get(globic).getPhysR ( + iPhysPos, iLatticePos ); + + // Is the voxel within a smooting kernel length from the bubble's position? + if ( pow(physPosP[0] - iPhysPos[0], 2) + + pow(physPosP[1] - iPhysPos[1], 2) + + pow(physPosP[2] - iPhysPos[2], 2) < pow(_kernelLength, 2) ) { + + // Adding the voxel's position (and relative weight) to the _latticePosAndWeight list + LatticePosAndWeight item; + item.latticePos[0] = iLatticePos[0]; + item.latticePos[1] = iLatticePos[1]; + item.latticePos[2] = iLatticePos[2]; + item.weight = this->compute(physPosP, iPhysPos); + + normalizer += item.weight; + _latticePosAndWeight.push_back(item); + } + } + } + } + + // If normalizer is zero, then no voxels are within a kernel smoothing length from the bubble's location. + // And it is a problem. + if (normalizer == T()) { + std::cout << "ERROR: SmoothingFunctional::update(...):" << std::endl + << "[smoothingFunctional] physPosP: " + << physPosP[0] << " " + << physPosP[1] << " " + << physPosP[2] << std::endl + << "[smoothingFunctional] physPosMin: " + << physPosMin[0] << " " + << physPosMin[1] << " " + << physPosMin[2] << std::endl + << "[smoothingFunctional] physPosMax: " + << latticePosMax[0] << " " + << latticePosMax[1] << " " + << latticePosMax[2] << std::endl + << "[smoothingFunctional] normalizer: " + << normalizer << std::endl; + return false; + } + + // Normalizing to one + for (int i=0; i +LinearAveragingSmoothingFunctional::LinearAveragingSmoothingFunctional ( + T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice ) + : SmoothingFunctional(kernelLength, converter, sLattice) +{} + +template +T LinearAveragingSmoothingFunctional::compute(T physPosP[], T physPosL[]) +{ + return this->smoothingFunction(physPosP[0] - physPosL[0]) + * this->smoothingFunction(physPosP[1] - physPosL[1]) + * this->smoothingFunction(physPosP[2] - physPosL[2]); +} + + +////////////////////// Class VolumeAveragingSmoothingFunctional //////////////////////// + +template +VolumeAveragingSmoothingFunctional::VolumeAveragingSmoothingFunctional ( + T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice ) + : SmoothingFunctional(kernelLength, converter, sLattice) +{} + +template +T VolumeAveragingSmoothingFunctional::compute(T physPosP[], T physPosL[]) +{ + return this->smoothingFunction ( sqrt ( + pow(physPosP[0] - physPosL[0], 2) + + pow(physPosP[1] - physPosL[1], 2) + + pow(physPosP[2] - physPosL[2], 2) ) ); +} + + +////////////////////// Class DeenSmoothingFunctional //////////////////////// + +template +DeenSmoothingFunctional::DeenSmoothingFunctional ( + T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice ) + : LinearAveragingSmoothingFunctional(kernelLength, converter, sLattice) +{} + +template +T DeenSmoothingFunctional::smoothingFunction(T delta) +{ + return ( pow(delta, 4)/pow(this->_kernelLength, 5) + - 2.*pow(delta, 2)/pow(this->_kernelLength, 3) + + 1./this->_kernelLength + ); +} + + +////////////////////// Class StepSmoothingFunctional //////////////////////// + +template +StepSmoothingFunctional::StepSmoothingFunctional ( + T kernelLength, UnitConverter& converter, SuperLattice3D& sLattice ) + : VolumeAveragingSmoothingFunctional(kernelLength, converter, sLattice) +{} + +template +T StepSmoothingFunctional::smoothingFunction(T delta) +{ + return 1.; +} + + +} + +#endif -- cgit v1.2.3