OpenLB 1.7
Loading...
Searching...
No Matches
setBouzidiVelocityBoundary2D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2020 Alexander Schulz
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24//This file contains the Bouzidi Velocity Boundary
25//This is an offLattice Boundary
26//This is a new version of the Boundary, which only contains free floating functions
27#ifndef SET_BOUZIDI_VELOCITY_BOUNDARY_2D_HH
28#define SET_BOUZIDI_VELOCITY_BOUNDARY_2D_HH
29
31
32
33namespace olb {
34
35namespace legacy {
36
38
40template<typename T, typename DESCRIPTOR, class MixinDynamics>
41void setBouzidiVelocityBoundary(SuperLattice<T, DESCRIPTOR>& sLattice, SuperGeometry<T,2>& superGeometry, int material, IndicatorF2D<T>& geometryIndicator,
42 std::vector<int> bulkMaterials)
43{
44
45 setBouzidiVelocityBoundary<T,DESCRIPTOR,MixinDynamics>(sLattice, superGeometry.getMaterialIndicator(material),
46 superGeometry.getMaterialIndicator(std::move(bulkMaterials)),
47 geometryIndicator);
48}
49
51template<typename T, typename DESCRIPTOR, class MixinDynamics>
53 FunctorPtr<SuperIndicatorF2D<T>>&& bulkIndicator,
54 IndicatorF2D<T>& geometryIndicator)
55{
56 T _epsFraction = 0.0001;
57 /* local boundaries: _overlap = 0;
58 * interp boundaries: _overlap = 1;
59 * bouzidi boundaries: _overlap = 1;
60 * extField boundaries: _overlap = 1;
61 * advectionDiffusion boundaries: _overlap = 1;
62 */
63 int _overlap = 1;
64
65 OstreamManager clout(std::cout, "setBouzidiVelocityBoundary");
66 clout << "epsFraction=" << _epsFraction << std::endl;
67 clout.setMultiOutput(true);
68 for (int iCloc = 0; iCloc < sLattice.getLoadBalancer().size(); ++iCloc) {
69 clout << "Cuboid globiC " << sLattice.getLoadBalancer().glob(iCloc)
70 << " starts to read distances for Velocity Boundary..." << std::endl;
71 setBouzidiVelocityBoundary<T,DESCRIPTOR,MixinDynamics>(sLattice.getBlock(iCloc),boundaryIndicator->getBlockIndicatorF(iCloc),
72 bulkIndicator->getBlockIndicatorF(iCloc),
73 geometryIndicator);
74 clout << "Cuboid globiC " << sLattice.getLoadBalancer().glob(iCloc)
75 << " finished reading distances for Velocity Boundary." << std::endl;
76 }
77 clout.setMultiOutput(false);
79 addPoints2CommBC<T,DESCRIPTOR>(sLattice, std::forward<decltype(boundaryIndicator)>(boundaryIndicator), _overlap);
80}
81
82
84
85
86template<typename T, typename DESCRIPTOR, class MixinDynamics>
88{
89 block.forSpatialLocations([&](auto iX, auto iY) {
90 if (boundaryIndicator(iX,iY)) {
91 setBouzidiVelocityBoundary<T,DESCRIPTOR,MixinDynamics>(block,
92 bulkIndicator.getBlockGeometry(),
93 iX, iY,
94 bulkIndicator,
95 geometryIndicator);
96 }
97 });
98}
99
100template<typename T, typename DESCRIPTOR, class MixinDynamics>
101void setBouzidiVelocityBoundary(BlockLattice<T,DESCRIPTOR>& block, BlockGeometry<T,2>& blockGeometryStructure, int iX, int iY, BlockIndicatorF2D<T>& bulkIndicator, IndicatorF2D<T>& geometryIndicator)
102{
103 T _epsFraction = 0.0001;
104 OstreamManager clout(std::cout, "setBouzidiVelocityBoundary");
105 T distances[DESCRIPTOR::q];
106 for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
107 distances[iPop] = -1;
108 }
109
110 for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
111 const int iXn = iX + descriptors::c<DESCRIPTOR >(iPop,0);
112 const int iYn = iY + descriptors::c<DESCRIPTOR >(iPop,1);
113 if (blockGeometryStructure.isInside(iXn,iYn) && bulkIndicator(iXn,iYn)) {
114 T dist = -1;
115 T physR[2];
116 blockGeometryStructure.getPhysR(physR,{iXn,iYn});
117 T voxelSize=blockGeometryStructure.getDeltaR();
118 Vector<T,2> physC(physR);
119
120 Vector<T,2> direction(-voxelSize*descriptors::c<DESCRIPTOR >(iPop,0),-voxelSize*descriptors::c<DESCRIPTOR >(iPop,1));
121 T cPhysNorm = voxelSize*util::sqrt(descriptors::c<DESCRIPTOR >(iPop,0)*descriptors::c<DESCRIPTOR >(iPop,0)+descriptors::c<DESCRIPTOR >(iPop,1)*descriptors::c<DESCRIPTOR >(iPop,1));
122
123 if (!geometryIndicator.distance(dist,physC,direction,blockGeometryStructure.getIcGlob() ) ) {
124 T epsX = voxelSize*descriptors::c<DESCRIPTOR >(iPop,0)*_epsFraction;
125 T epsY = voxelSize*descriptors::c<DESCRIPTOR >(iPop,1)*_epsFraction;
126
127 Vector<T,2> physC2(physC);
128 physC2[0] += epsX;
129 physC2[1] += epsY;
130 Vector<T,2> direction2(direction);
131 direction2[0] -= 2.*epsX;
132 direction2[1] -= 2.*epsY;
133
134 if ( !geometryIndicator.distance(dist,physC2,direction2,blockGeometryStructure.getIcGlob())) {
135 clout << "ERROR: no boundary found at (" << iXn << "," << iYn <<") ~ ("
136 << physR[0] << "," << physR[1] << "), "
137 << "in direction " << descriptors::opposite<DESCRIPTOR >(iPop)
138 << std::endl;
139 }
140 T distNew = (dist - util::sqrt(epsX*epsX+epsY*epsY))/cPhysNorm;
141 if (distNew < 0.5) {
142 dist = 0;
143 }
144 else {
145 dist = 0.5 * cPhysNorm;
146 clout << "WARNING: distance at (" << iXn << "," << iYn <<") ~ ("
147 << physR[0] << "," << physR[1] <<"), "
148 << "in direction " << descriptors::opposite<DESCRIPTOR >(iPop) << ": "
149 << distNew
150 << " rounded to "
151 << dist/cPhysNorm
152 << std::endl;
153 }
154 }
155 distances[descriptors::opposite<DESCRIPTOR >(iPop)] = dist/cPhysNorm;
156 } // bulk indicator
157 } // iPop loop
158 setBouzidiVelocityBoundary<T,DESCRIPTOR,MixinDynamics>(block, blockGeometryStructure, iX, iY, distances);
159
160}
161
162template<typename T, typename DESCRIPTOR, class MixinDynamics>
163void setBouzidiVelocityBoundary(BlockLattice<T,DESCRIPTOR>& block, BlockGeometry<T,2>& blockGeometryStructure, int x, int y, T distances[DESCRIPTOR::q])
164{
165 typedef DESCRIPTOR L;
166 T location[DESCRIPTOR::d];
167 blockGeometryStructure.getPhysR(location, {x,y});
168
169 T distancesCopy[L::q];
170 T spacing = blockGeometryStructure.getDeltaR();
171 for (int iPop = 1; iPop < L::q ; ++iPop) {
172 distancesCopy[iPop] = spacing*(1.-distances[iPop]);
173 if ( !util::nearZero(distances[iPop]+1) ) {
174 distancesCopy[iPop] = -1;
175 }
176 }
177 //can be universally used for Bouzidi and bounceBackVelocityBoundary
178 setOffDynamics<T,DESCRIPTOR,MixinDynamics>(block, x, y, location, distancesCopy);
179
180 for (int iPop = 1; iPop < L::q ; ++iPop) {
181 if ( !util::nearZero(distances[iPop]+1) ) {
182 setBouzidiVelocityBoundary<T,DESCRIPTOR,MixinDynamics>(block, blockGeometryStructure, x-descriptors::c<L>(iPop,0), y-descriptors::c<L>(iPop,1), iPop, distances[iPop]);
183 }
184 }
185
186}
187
188
189template<typename T, typename DESCRIPTOR, class MixinDynamics>
190void setBouzidiVelocityBoundary(BlockLattice<T,DESCRIPTOR>& block, BlockGeometry<T,2>& blockGeometryStructure, int x, int y, int iPop, T dist)
191{
193 if (blockGeometryStructure.getMaterial({x-descriptors::c<DESCRIPTOR >(iPop,0), y-descriptors::c<DESCRIPTOR >(iPop,1)}) != 1) {
195 postProcessor = new VelocityBounceBackPostProcessorGenerator2D<T,DESCRIPTOR>(x, y, iPop, dist);
196
197 }
198 else {
200 postProcessor =new VelocityBouzidiLinearPostProcessorGenerator2D<T,DESCRIPTOR>(x, y, iPop, dist);
201 }
202
203 if (postProcessor && !block.isPadding({x,y})) {
204 block.addPostProcessor(*postProcessor);
205 }
206}
207
208//set Dynamics for indicated Bouzidi and BounceBackVelocityBoundary cells
209template<typename T, typename DESCRIPTOR, class MixinDynamics>
210void setOffDynamics(BlockLattice<T,DESCRIPTOR>& block, int x, int y, T location[DESCRIPTOR::d], T distances[DESCRIPTOR::q])
211{
212 //Dynamics<T,DESCRIPTOR>* dynamics = BoundaryManager::getOffDynamics(location, distances);
213 Dynamics<T,DESCRIPTOR>* dynamics = new legacy::OffDynamics<T, DESCRIPTOR>(location, distances);
214 block.defineDynamics({x,y}, dynamics);
215}
216
217}
218
219}//namespace olb
220
221#endif
Representation of a block geometry.
Vector< T, D > getPhysR(LatticeR< D > latticeR)
int getMaterial(LatticeR< D > latticeR) const
returns the (iX,iY) entry in the 2D scalar field
int getIcGlob() const
Read only access to the global iC number which is given !=-1 if the block geometries are part of a su...
T getDeltaR() const
Read only access to the voxel size given in SI units (meter)
Base block indicator functor (discrete)
BlockGeometry< T, 2 > & getBlockGeometry()
Get underlying block geometry structure.
Platform-abstracted block lattice for external access and inter-block interaction.
virtual void addPostProcessor(std::type_index stage, LatticeR< DESCRIPTOR::d > latticeR, PostProcessorPromise< T, DESCRIPTOR > &&promise)=0
Schedule post processor for application to latticeR in stage.
void defineDynamics(LatticeR< DESCRIPTOR::d > latticeR, DynamicsPromise< T, DESCRIPTOR > &&promise)
Assign promised DYNAMICS to latticeR.
void forSpatialLocations(F f) const
bool isPadding(LatticeR< D > latticeR) const
Return whether location is valid.
bool isInside(LatticeR< D > latticeR) const
Return whether location is valid.
Smart pointer for managing the various ways of passing functors around.
Definition functorPtr.h:60
IndicatorF2D is an application from .
virtual bool distance(S &distance, const Vector< S, 2 > &origin, S precision, const Vector< S, 2 > &direction)
returns false or true and pos. distance if there was one found for an given origin and direction
class for marking output with some text
void setMultiOutput(bool b)
enable message output for all MPI processes, disabled by default
Representation of a statistic for a parallel 2D geometry.
std::unique_ptr< SuperIndicatorF< T, D > > getMaterialIndicator(std::vector< int > &&materials)
Returns a material indicator using the given vector of materials.
Super class maintaining block lattices for a cuboid decomposition.
BlockLattice< T, DESCRIPTOR > & getBlock(int locC)
Return BlockLattice with local index locC.
LoadBalancer< T > & getLoadBalancer()
Read and write access to the load balancer.
Plain old scalar vector.
Definition vector.h:47
Dynamics for offLattice boundary conditions OffDynamics are basically NoLatticeDynamics with the addi...
Definition dynamics.h:145
void setOffDynamics(BlockLattice< T, DESCRIPTOR > &block, int x, int y, T location[DESCRIPTOR::d], T distances[DESCRIPTOR::q])
void setBouzidiVelocityBoundary(SuperLattice< T, DESCRIPTOR > &sLattice, SuperGeometry< T, 2 > &superGeometry, int material, IndicatorF2D< T > &geometryIndicator, std::vector< int > bulkMaterials=std::vector< int >(1, 1))
Set offDynamics with boundary links and post processors using indicators.
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
Top level namespace for all of OpenLB.
Interface for per-cell dynamics.
Definition interface.h:54