diff options
Diffstat (limited to 'src/boundary/advectionDiffusionBoundaries.hh')
| -rw-r--r-- | src/boundary/advectionDiffusionBoundaries.hh | 313 |
1 files changed, 313 insertions, 0 deletions
diff --git a/src/boundary/advectionDiffusionBoundaries.hh b/src/boundary/advectionDiffusionBoundaries.hh new file mode 100644 index 0000000..7c7f391 --- /dev/null +++ b/src/boundary/advectionDiffusionBoundaries.hh @@ -0,0 +1,313 @@ +/* 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 + * <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 ADVECTION_DIFFUSION_BOUNDARIES_HH +#define ADVECTION_DIFFUSION_BOUNDARIES_HH + +#include "advectionDiffusionBoundaries.h" +#include "dynamics/latticeDescriptors.h" +#include "core/util.h" +#include "dynamics/lbHelpers.h" + +namespace olb { + + + +//================================================================================================== +//==================== For regularized Advection Diffusion Boundary Condition ====================== +//============================================================================================ + + +// For flat Walls + +template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation> +AdvectionDiffusionBoundariesDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>:: +AdvectionDiffusionBoundariesDynamics( T omega_, Momenta<T,DESCRIPTOR>& momenta_) + : BasicDynamics<T,DESCRIPTOR>(momenta_), boundaryDynamics(omega_, momenta_) +{ +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation> +T AdvectionDiffusionBoundariesDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>:: +computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return lbHelpers<T,DESCRIPTOR>::equilibriumFirstOrder( iPop, rho, u ); +} + + +template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation> +void AdvectionDiffusionBoundariesDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>:: +collide(Cell<T,DESCRIPTOR>& cell,LatticeStatistics<T>& statistics) +{ + typedef DESCRIPTOR L; + typedef lbHelpers<T,DESCRIPTOR> lbH; + + T dirichletTemperature = this->_momenta.computeRho(cell); + T* u = cell.template getFieldPointer<descriptors::VELOCITY>(); + + std::vector<int> unknownIndexes = util::subIndexOutgoing<L, direction, + orientation>(); + std::vector<int> knownIndexes = util::remainingIndexes<L>(unknownIndexes); + + int missingNormal = 0; + + if ((L::d == 3 && L::q == 7)||(L::d == 2 && L::q == 5)) { + T sum = T(); + for (unsigned i = 0; i < knownIndexes.size(); ++i) { + sum += cell[knownIndexes[i]]; + } + + T difference = dirichletTemperature - (T) 1 - sum; // on cell there are non-shiftet values -> temperature has to be changed + + // here I know all missing and non missing f_i + for (unsigned i = 0; i < unknownIndexes.size(); ++i) { + int numOfNonNullComp = 0; + for (int iDim = 0; iDim < L::d; ++iDim) { + numOfNonNullComp += abs(descriptors::c<L>(unknownIndexes[i],iDim)); + } + if (numOfNonNullComp == 1) { + missingNormal = unknownIndexes[i]; + // here missing diagonal directions are erased + // just the normal direction stays (D3Q7) + unknownIndexes.erase(unknownIndexes.begin() + i); + break; + + } + } + cell[missingNormal] = difference; // on cell there are non-shiftet values -> temperature has to be changed + boundaryDynamics.collide(cell, statistics); // only for D3Q7 + } else { + // part for q=19 copied from AdvectionDiffusionEdgesDynamics.collide() + // but here just all missing directions, even at border of inlet area + // has to be checked! + for (unsigned iteratePop = 0; iteratePop < unknownIndexes.size(); + ++iteratePop) { + cell[unknownIndexes[iteratePop]] = + lbH::equilibriumFirstOrder(unknownIndexes[iteratePop], dirichletTemperature, u) + - (cell[util::opposite<L>(unknownIndexes[iteratePop])] + - lbH::equilibriumFirstOrder( + util::opposite<L>(unknownIndexes[iteratePop]), + dirichletTemperature, u)); + } + } +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation> +T AdvectionDiffusionBoundariesDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>:: +getOmega() const +{ + return boundaryDynamics.getOmega(); +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation> +void AdvectionDiffusionBoundariesDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>:: +setOmega(T omega_) +{ + boundaryDynamics.setOmega(omega_); +} + +//================================================================= +// For 2D Corners with regularized Dynamic ============================================== +//================================================================= +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal> +AdvectionDiffusionCornerDynamics2D<T,DESCRIPTOR,Dynamics,xNormal,yNormal>::AdvectionDiffusionCornerDynamics2D( + T omega_, Momenta<T,DESCRIPTOR>& momenta_) + : BasicDynamics<T,DESCRIPTOR>(momenta_), + boundaryDynamics(omega_, momenta_) +{ +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal> +T AdvectionDiffusionCornerDynamics2D<T,DESCRIPTOR,Dynamics,xNormal,yNormal>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return lbHelpers<T,DESCRIPTOR>::equilibriumFirstOrder( iPop, rho, u ); +} + + +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal> +void AdvectionDiffusionCornerDynamics2D<T,DESCRIPTOR,Dynamics,xNormal,yNormal>::collide(Cell<T,DESCRIPTOR>& cell,LatticeStatistics<T>& statistics) +{ + typedef DESCRIPTOR L; + typedef lbHelpers<T,DESCRIPTOR> lbH; + + T temperature = this->_momenta.computeRho(cell); + T* u = cell.template getFieldPointer<descriptors::VELOCITY>(); + // I need to get Missing information on the corners !!!! + std::vector<int> unknownIndexes = util::subIndexOutgoing2DonCorners<L,xNormal,yNormal>(); + // here I know all missing and non missing f_i + + + // The collision procedure for D2Q5 and D3Q7 lattice is the same ... + // Given the rule f_i_neq = -f_opposite(i)_neq + // I have the right number of equations for the number of unknowns using these lattices + + for (unsigned iPop = 0; iPop < unknownIndexes.size(); ++iPop) { + cell[unknownIndexes[iPop]] = lbH::equilibriumFirstOrder(unknownIndexes[iPop], temperature, u) + -(cell[util::opposite<L>(unknownIndexes[iPop])] + - lbH::equilibriumFirstOrder(util::opposite<L>(unknownIndexes[iPop]), temperature, u) ) ; + } + + // Once all the f_i are known, I can call the collision for the Regularized Model. + boundaryDynamics.collide(cell, statistics); + +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal> +T AdvectionDiffusionCornerDynamics2D<T,DESCRIPTOR,Dynamics,xNormal,yNormal>::getOmega() const +{ + return boundaryDynamics.getOmega(); +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal> +void AdvectionDiffusionCornerDynamics2D<T,DESCRIPTOR,Dynamics,xNormal,yNormal>::setOmega(T omega_) +{ + boundaryDynamics.setOmega(omega_); +} + + + +//================================================================= +// For 3D Corners with regularized Dynamic ============================================== +//================================================================= +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal, int zNormal> +AdvectionDiffusionCornerDynamics3D<T,DESCRIPTOR,Dynamics,xNormal,yNormal,zNormal>::AdvectionDiffusionCornerDynamics3D( + T omega_, Momenta<T,DESCRIPTOR>& momenta_) + : BasicDynamics<T,DESCRIPTOR>(momenta_), + boundaryDynamics(omega_, momenta_) +{ +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal, int zNormal> +T AdvectionDiffusionCornerDynamics3D<T,DESCRIPTOR,Dynamics,xNormal,yNormal,zNormal>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return lbHelpers<T, DESCRIPTOR>::equilibriumFirstOrder( iPop, rho, u ); +} + + +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal, int zNormal> +void AdvectionDiffusionCornerDynamics3D<T,DESCRIPTOR,Dynamics,xNormal,yNormal,zNormal>::collide(Cell<T,DESCRIPTOR>& cell,LatticeStatistics<T>& statistics) +{ + typedef DESCRIPTOR L; + typedef lbHelpers<T,DESCRIPTOR> lbH; + + T temperature = this->_momenta.computeRho(cell); + T* u = cell.template getFieldPointer<descriptors::VELOCITY>(); + // I need to get Missing information on the corners !!!! + std::vector<int> unknownIndexes = util::subIndexOutgoing3DonCorners<L,xNormal,yNormal,zNormal>(); + // here I know all missing and non missing f_i + + + // The collision procedure for D2Q5 and D3Q7 lattice is the same ... + // Given the rule f_i_neq = -f_opposite(i)_neq + // I have the right number of equations for the number of unknowns using these lattices + + for (unsigned iPop = 0; iPop < unknownIndexes.size(); ++iPop) { + cell[unknownIndexes[iPop]] = lbH::equilibriumFirstOrder(unknownIndexes[iPop], temperature, u) + -(cell[util::opposite<L>(unknownIndexes[iPop])] + - lbH::equilibriumFirstOrder(util::opposite<L>(unknownIndexes[iPop]), temperature, u) ) ; + } + + // Once all the f_i are known, I can call the collision for the Regularized Model. + boundaryDynamics.collide(cell, statistics); + +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal, int zNormal> +T AdvectionDiffusionCornerDynamics3D<T,DESCRIPTOR,Dynamics,xNormal,yNormal,zNormal>::getOmega() const +{ + return boundaryDynamics.getOmega(); +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int xNormal, int yNormal, int zNormal> +void AdvectionDiffusionCornerDynamics3D<T,DESCRIPTOR,Dynamics,xNormal,yNormal,zNormal>::setOmega(T omega_) +{ + boundaryDynamics.setOmega(omega_); +} + +//================================================================= +// For 3D Edges with regularized Dynamic ============================================== +//================================================================= +template<typename T, typename DESCRIPTOR, typename Dynamics, int plane, int normal1, int normal2> +AdvectionDiffusionEdgesDynamics<T,DESCRIPTOR,Dynamics,plane,normal1, normal2>::AdvectionDiffusionEdgesDynamics( + T omega_, Momenta<T,DESCRIPTOR>& momenta_) + : BasicDynamics<T,DESCRIPTOR>(momenta_), + boundaryDynamics(omega_, momenta_) +{ +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int plane, int normal1, int normal2> +T AdvectionDiffusionEdgesDynamics<T,DESCRIPTOR,Dynamics,plane,normal1, normal2>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return lbHelpers<T,DESCRIPTOR>::equilibriumFirstOrder( iPop, rho, u ); +} + + +template<typename T, typename DESCRIPTOR, typename Dynamics, int plane, int normal1, int normal2> +void AdvectionDiffusionEdgesDynamics<T,DESCRIPTOR,Dynamics,plane,normal1, normal2>::collide(Cell<T,DESCRIPTOR>& cell,LatticeStatistics<T>& statistics) +{ + typedef DESCRIPTOR L; + typedef lbHelpers<T,DESCRIPTOR> lbH; + + T temperature = this->_momenta.computeRho(cell); + T* u = cell.template getFieldPointer<descriptors::VELOCITY>(); + // I need to get Missing information on the corners !!!! + std::vector<int> unknownIndexes = util::subIndexOutgoing3DonEdges<L,plane,normal1, normal2>(); + // here I know all missing and non missing f_i + + + // The collision procedure for D2Q5 and D3Q7 lattice is the same ... + // Given the rule f_i_neq = -f_opposite(i)_neq + // I have the right number of equations for the number of unknowns using these lattices + + for (unsigned iPop = 0; iPop < unknownIndexes.size(); ++iPop) { + cell[unknownIndexes[iPop]] = lbH::equilibriumFirstOrder(unknownIndexes[iPop], temperature, u) + -(cell[util::opposite<L>(unknownIndexes[iPop])] + - lbH::equilibriumFirstOrder(util::opposite<L>(unknownIndexes[iPop]), temperature, u) ) ; + } + + // Once all the f_i are known, I can call the collision for the Regularized Model. + boundaryDynamics.collide(cell, statistics); + +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int plane, int normal1, int normal2> +T AdvectionDiffusionEdgesDynamics<T,DESCRIPTOR,Dynamics,plane,normal1, normal2>::getOmega() const +{ + return boundaryDynamics.getOmega(); +} + +template<typename T, typename DESCRIPTOR, typename Dynamics, int plane, int normal1, int normal2> +void AdvectionDiffusionEdgesDynamics<T,DESCRIPTOR,Dynamics,plane,normal1, normal2>::setOmega(T omega_) +{ + boundaryDynamics.setOmega(omega_); +} + + + + +} // namespace olb + + + + +#endif |