OpenLB 1.7
Loading...
Searching...
No Matches
forwardCouplingModels3D.hh
Go to the documentation of this file.
1/* Lattice Boltzmann sample, written in C++, using the OpenLB
2 * library
3 *
4 * Copyright (C) 2019 Davide Dapelo
5 * E-mail contact: info@openlb.net
6 * The most recent release of OpenLB can be downloaded at
7 * <http://www.openlb.net/>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public
20 * License along with this program; if not, write to the Free
21 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 * Boston, MA 02110-1301, USA.
23 */
24
25/* Models for Lagrangian forward-coupling methods -- generic implementation.
26 */
27
28#ifndef LB_FORWARD_COUPLING_MODELS_HH
29#define LB_FORWARD_COUPLING_MODELS_HH
30
31namespace olb {
32
33
35
36template<typename T, template<typename V> class Particle>
38 SuperGeometry<T,3>& sGeometry,
39 std::shared_ptr<DragModel<T, Particle>> dragModel )
40 : _sGeometry(sGeometry),
41 _dragModel(dragModel)
42{}
43
44
46
47template<typename T, typename Lattice, template<typename V> class Particle>
51 SuperGeometry<T,3>& sGeometry,
52 std::shared_ptr<DragModel<T, Particle>> dragModel )
53 : ForwardCouplingModel<T,Particle>(sGeometry, dragModel),
54 _converter(converter),
55 _sLattice(sLattice)
56{
57 this->_interpLatticeDensity = std::make_shared<SuperLatticeInterpDensity3Degree3D<T, Lattice> > (
58 this->_sLattice, this->_sGeometry, this->_converter );
59 this->_interpLatticeVelocity = std::make_shared<SuperLatticeInterpPhysVelocity3D<T, Lattice> > (
60 this->_sLattice, this->_converter );
61}
62
63
65
66template<typename T, typename Lattice, template<typename V> class Particle>
70 SuperGeometry<T,3>& sGeometry,
71 std::shared_ptr<DragModel<T, Particle>> dragModel )
72 : BaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel)
73{}
74
75template<typename T, typename Lattice, template<typename V> class Particle>
77{
79 T physPosP[3] = { p->getPos()[0],
80 p->getPos()[1],
81 p->getPos()[2]
82 }; // particle's physical position
83 // particle's dimensionless position, rounded at neighbouring voxel
84 int latticeRoundedPosP[3] = {0, 0, 0};
85 this->_sLattice.getCuboidGeometry().get(globic).getLatticeR (
86 latticeRoundedPosP, physPosP );
87 // { globic, latticeRoundedP[0, ..., 2] }
88 int globicFull[4] = { globic,
89 latticeRoundedPosP[0],
90 latticeRoundedPosP[1],
91 latticeRoundedPosP[2]
92 };
93
94 // Particle's velocity
95 T physVelP[3] = { p->getVel()[0],
96 p->getVel()[1],
97 p->getVel()[2]
98 }; // Physical
99 // Lattice
100 T latticeVelP[3] = { this->_converter.getLatticeVelocity(physVelP[0]),
101 this->_converter.getLatticeVelocity(physVelP[1]),
102 this->_converter.getLatticeVelocity(physVelP[2])
103 }; // particle's dimensionless velocity
104
105 // Lattice's velocity at particle's location
106 T physVelF[3] = {T(), T(), T()}; // Physical
107 this->_interpLatticeVelocity->operator() (physVelF, physPosP, globic);
108 // Lattice
109 T latticeVelF[3] = { this->_converter.getLatticeVelocity(physVelF[0]),
110 this->_converter.getLatticeVelocity(physVelF[1]),
111 this->_converter.getLatticeVelocity(physVelF[2])
112 }; // Lattice's dimensionless velocity at particle's location
113
114 // Computing fluid-particle momentum transfer
115 T gF[3] = {T(), T(), T()}; // force density gF
116 this->_momentumExchange->operator() (gF, latticeVelF, latticeVelP, physPosP, latticeRoundedPosP, globic);
117
118 // Computing drag coefficient
119 T Cd = this->_dragModel->operator()(p, latticeVelF, latticeVelP, globicFull);
120#ifdef VERBOSE
121 std:: cout << Cd
122 << " physPosP=(" << physPosP[0] << ", " << physPosP[1] << ", " << physPosP[2] << ") "
123 << " physVelP=(" << physVelP[0] << ", " << physVelP[1] << ", " << physVelP[2] << ") "
124 << " physVelF=(" << physVelF[0] << ", " << physVelF[1] << ", " << physVelF[2] << ") "
125 << std::endl;
126#endif
127 /*
128 if (Cd > 100.)
129 throw std::range_error ( "LocalBaseForwardCouplingModel::operator(). Cd="+std::to_string(Cd)
130 + "\nphysVelP=(" + std::to_string(physVelP[0]) + ", " + std::to_string(physVelP[1]) + ", " + std::to_string(physVelP[2]) + ")"
131 + "\nphysVelF=(" + std::to_string(physVelF[0]) + ", " + std::to_string(physVelF[1]) + ", " + std::to_string(physVelF[2]) + ")"
132 + "\n");
133 */
134
136 T latticePRad = p->getRad() / this->_converter.getConversionFactorLength();
137 T latticeForceP[3] = { .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[0] * (latticeVelF[0] - latticeVelP[0]),
138 .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[1] * (latticeVelF[1] - latticeVelP[1]),
139 .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[2] * (latticeVelF[2] - latticeVelP[2])
140 }; // dimensionless force acting on the particle
141
143 std::vector<T> physForceP(3, T()); // physical force acting on the particle
144 physForceP[0] = latticeForceP[0] * this->_converter.getConversionFactorForce();
145 physForceP[1] = latticeForceP[1] * this->_converter.getConversionFactorForce();
146 physForceP[2] = latticeForceP[2] * this->_converter.getConversionFactorForce();
147
149 p->setStoreForce(physForceP);
150 return true;
151}
152
153/*
154template<typename T, typename Lattice, template<typename V> class Particle>
155bool LocalBaseForwardCouplingModel<T,Lattice,Particle>::operator() (SuperParticleSystem3D<T,Particle>& spSys)
156{
157 return true;
158}
159*/
160
161
163
164template<typename T, typename Lattice, template<typename V> class Particle>
166 UnitConverter<T, Lattice>& converter,
167 SuperLattice<T, Lattice>& sLattice,
168 SuperGeometry<T,3>& sGeometry,
169 std::shared_ptr<DragModel<T, Particle>> dragModel )
170 : LocalBaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel)
171{
172 this->_momentumExchange = std::make_shared<NaiveMomentumExchange<T, Lattice> > (
173 this->_converter, this->_sLattice, this->_interpLatticeDensity );
174}
175
176
178
179template<typename T, typename Lattice, template<typename V> class Particle>
181 UnitConverter<T, Lattice>& converter,
182 SuperLattice<T, Lattice>& sLattice,
183 SuperGeometry<T,3>& sGeometry,
184 std::shared_ptr<DragModel<T, Particle>> dragModel )
185 : LocalBaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel)
186{
187 this->_momentumExchange = std::make_shared<LaddMomentumExchange<T, Lattice> > (
188 this->_converter, this->_sLattice,
190}
191
192
194
195template<typename T, typename Lattice, template<typename V> class Particle>
197 UnitConverter<T, Lattice>& converter,
198 SuperLattice<T, Lattice>& sLattice,
199 SuperGeometry<T,3>& sGeometry,
200 std::shared_ptr<DragModel<T, Particle>> dragModel,
201 std::shared_ptr<SmoothingFunctional<T, Lattice>> smoothingFunctional )
202 : BaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel),
203 _smoothingFunctional(smoothingFunctional)
204{}
205
206template<typename T, typename Lattice, template<typename V> class Particle>
208{
210 T physPosP[3] = { p->getPos()[0],
211 p->getPos()[1],
212 p->getPos()[2]
213 }; // particle's physical position
214 // particle's dimensionless position, rounded at neighbouring voxel
215 int latticeRoundedPosP[3] = {0, 0, 0};
216 this->_sLattice.getCuboidGeometry().get(globic).getLatticeR (
217 latticeRoundedPosP, physPosP );
218 // { globic, latticeRoundedP[0, ..., 2] }
219 int globicFull[4] = { globic,
220 latticeRoundedPosP[0],
221 latticeRoundedPosP[1],
222 latticeRoundedPosP[2]
223 };
224
225 // Particle's velocity
226 T physVelP[3] = { p->getVel()[0],
227 p->getVel()[1],
228 p->getVel()[2]
229 }; // Physical
230 // Lattice
231 T latticeVelP[3] = { this->_converter.getLatticeVelocity(physVelP[0]),
232 this->_converter.getLatticeVelocity(physVelP[1]),
233 this->_converter.getLatticeVelocity(physVelP[2])
234 }; // particle's dimensionless velocity
235
236 // Update the smoothing functional
237 if ( ! this->_smoothingFunctional->update(physPosP, globic)) {
238 std::cout << "ERROR: no lattice point enclosed in particle's kernel length!" << std::endl;
239 return false;
240 }
241
243 T latticeForceP[3] = {T(), T(), T()}; // dimensionless force acting on the particle
244 for (auto&& i : this->_smoothingFunctional->getData()) {
245
246 // Position of the iterated voxel
247 int iLatticePosF[3] = {i.latticePos[0], i.latticePos[1], i.latticePos[2]};
248 // Physical
249 T iPhysPosF[3] = {T(), T(), T()};
250 this->_sLattice.getCuboidGeometry().get(globic).getPhysR (
251 iPhysPosF, iLatticePosF );
252
253 // Fluid velocity at the iterated voxel
254 T iPhysVelF[3] = {T(), T(), T()}; // Physical
255 this->_interpLatticeVelocity->operator() (iPhysVelF, iPhysPosF, globic);
256 // Lattice
257 T iLatticeVelF[3] = { this->_converter.getLatticeVelocity(iPhysVelF[0]),
258 this->_converter.getLatticeVelocity(iPhysVelF[1]),
259 this->_converter.getLatticeVelocity(iPhysVelF[2])
260 }; // Lattice's dimensionless velocity at particle's location
261
262 // Computing fluid-particle momentum transfer
263 T gF[3] = {T(), T(), T()}; // force density gF
264 this->_momentumExchange->operator() ( gF, iLatticeVelF, latticeVelP, physPosP, iLatticePosF, globic);
265
266 // Computing drag coefficient
267 T Cd = this->_dragModel->operator()(p, iLatticeVelF, latticeVelP, globicFull);
268
270 T latticePRad = p->getRad() / this->_converter.getConversionFactorLength();
271 latticeForceP[0] += .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[0] * (iLatticeVelF[0] - latticeVelP[0]) * i.weight;
272 latticeForceP[1] += .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[1] * (iLatticeVelF[1] - latticeVelP[1]) * i.weight;
273 latticeForceP[2] += .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[2] * (iLatticeVelF[2] - latticeVelP[2]) * i.weight;
274 }
275
277 std::vector<T> physForceP(3, T()); // physical force acting on the particle
278 physForceP[0] = latticeForceP[0] * this->_converter.getConversionFactorForce();
279 physForceP[1] = latticeForceP[1] * this->_converter.getConversionFactorForce();
280 physForceP[2] = latticeForceP[2] * this->_converter.getConversionFactorForce();
281
283 p->setStoreForce(physForceP);
284 return true;
285}
286
287
289
290template<typename T, typename Lattice, template<typename V> class Particle>
292 UnitConverter<T, Lattice>& converter,
293 SuperLattice<T, Lattice>& sLattice,
294 SuperGeometry<T,3>& sGeometry,
295 std::shared_ptr<DragModel<T, Particle>> dragModel,
296 std::shared_ptr<SmoothingFunctional<T, Lattice>> smoothingFunctional )
297 : NonLocalBaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel, smoothingFunctional)
298{
299 this->_momentumExchange = std::make_shared<NaiveMomentumExchange<T, Lattice> > (
300 this->_converter, this->_sLattice, this->_interpLatticeDensity );
301}
302
303
305
306template<typename T, typename Lattice, template<typename V> class Particle>
308 UnitConverter<T, Lattice>& converter,
309 SuperLattice<T, Lattice>& sLattice,
310 SuperGeometry<T,3>& sGeometry,
311 std::shared_ptr<DragModel<T, Particle>> dragModel,
312 std::shared_ptr<SmoothingFunctional<T, Lattice>> smoothingFunctional, int nVoxelInterpPoints )
313 : BaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel),
314 _smoothingFunctional(smoothingFunctional)
315{
316 this->_momentumExchange = std::make_shared<NaiveMomentumExchange<T, Lattice> > (
317 this->_converter, this->_sLattice, this->_interpLatticeDensity );
318}
319
320template<typename T, typename Lattice, template<typename V> class Particle>
322{
324 T physPosP[3] = { p->getPos()[0],
325 p->getPos()[1],
326 p->getPos()[2]
327 }; // particle's physical position
328 // particle's dimensionless position, rounded at neighbouring voxel
329 int latticeRoundedPosP[3] = {0, 0, 0};
330 this->_sLattice.getCuboidGeometry().get(globic).getLatticeR (
331 latticeRoundedPosP, physPosP );
332 // { globic, latticeRoundedP[0, ..., 2] }
333 int globicFull[4] = { globic,
334 latticeRoundedPosP[0],
335 latticeRoundedPosP[1],
336 latticeRoundedPosP[2]
337 };
338
339 // Particle's velocity
340 T physVelP[3] = { p->getVel()[0],
341 p->getVel()[1],
342 p->getVel()[2]
343 }; // Physical
344 // Lattice
345 T latticeVelP[3] = { this->_converter.getLatticeVelocity(physVelP[0]),
346 this->_converter.getLatticeVelocity(physVelP[1]),
347 this->_converter.getLatticeVelocity(physVelP[2])
348 }; // particle's dimensionless velocity
349
350 // Lattice's velocity at particle's location
351 T physVelF[3] = {T(), T(), T()}; // Physical
352 this->_interpLatticeVelocity->operator() (physVelF, physPosP, globic);
353 // Lattice
354 T latticeVelF[3] = { this->_converter.getLatticeVelocity(physVelF[0]),
355 this->_converter.getLatticeVelocity(physVelF[1]),
356 this->_converter.getLatticeVelocity(physVelF[2])
357 }; // Lattice's dimensionless velocity at particle's location
358
359 // Computing fluid-particle momentum transfer
360 T gF[3] = {T(), T(), T()}; // force density gF
361 this->_momentumExchange->operator() (gF, latticeVelF, latticeVelP, physPosP, latticeRoundedPosP, globic);
362
363 // Update the smoothing functional
364 if ( ! this->_smoothingFunctional->update(physPosP, globic)) {
365 std::cout << "ERROR: no lattice point enclosed in particle's kernel length!" << std::endl;
366 return false;
367 }
368
369 // Computing continuous phase fraction across kernel
370 T continuousPhaseFraction = T();
371 for (auto&& i : this->_smoothingFunctional->getData()) {
372 continuousPhaseFraction += i.continuousPhaseFraction * i.weight;
373 }
374
375 // Computing drag coefficient
376 T Cd = this->_dragModel->operator()(p, latticeVelF, latticeVelP, globicFull, continuousPhaseFraction);
377#ifdef VERBOSE
378 std:: cout << Cd
379 << " physPosP=(" << physPosP[0] << ", " << physPosP[1] << ", " << physPosP[2] << ") "
380 << " physVelP=(" << physVelP[0] << ", " << physVelP[1] << ", " << physVelP[2] << ") "
381 << " physVelF=(" << physVelF[0] << ", " << physVelF[1] << ", " << physVelF[2] << ") "
382 << std::endl;
383#endif
384 /*
385 if (Cd > 100.)
386 throw std::range_error ( "LocalBaseForwardCouplingModel::operator(). Cd="+std::to_string(Cd)
387 + "\nphysVelP=(" + std::to_string(physVelP[0]) + ", " + std::to_string(physVelP[1]) + ", " + std::to_string(physVelP[2]) + ")"
388 + "\nphysVelF=(" + std::to_string(physVelF[0]) + ", " + std::to_string(physVelF[1]) + ", " + std::to_string(physVelF[2]) + ")"
389 + "\n");
390 */
391
393 T latticePRad = p->getRad() / this->_converter.getConversionFactorLength();
394 T latticeForceP[3] = { .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[0] * (latticeVelF[0] - latticeVelP[0]) / util::pow(continuousPhaseFraction, 1.65),
395 .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[1] * (latticeVelF[1] - latticeVelP[1]) / util::pow(continuousPhaseFraction, 1.65),
396 .5 * Cd * M_PI*util::pow(latticePRad,2) * gF[2] * (latticeVelF[2] - latticeVelP[2]) / util::pow(continuousPhaseFraction, 1.65)
397 };
398
400 std::vector<T> physForceP(3, T()); // physical force acting on the particle
401 physForceP[0] = latticeForceP[0] * this->_converter.getConversionFactorForce();
402 physForceP[1] = latticeForceP[1] * this->_converter.getConversionFactorForce();
403 physForceP[2] = latticeForceP[2] * this->_converter.getConversionFactorForce();
404
406 p->setStoreForce(physForceP);
407 return true;
408}
409
410}
411
412#endif
#define M_PI
Abstact base class for all the local/non-local forward-coupling models.
UnitConverter< T, Lattice > & _converter
std::shared_ptr< TwoWayHelperFunctional< T, Lattice > > _momentumExchange
std::shared_ptr< SuperLatticeInterpPhysVelocity3D< T, Lattice > > _interpLatticeVelocity
BaseForwardCouplingModel(UnitConverter< T, Lattice > &converter, SuperLattice< T, Lattice > &sLattice, SuperGeometry< T, 3 > &sGeometry, std::shared_ptr< DragModel< T, Particle > > dragModel)
Constructor.
std::shared_ptr< SuperLatticeInterpDensity3Degree3D< T, Lattice > > _interpLatticeDensity
SuperLattice< T, Lattice > & _sLattice
Abstact base class for DragModelBase.
Abstact base class for all the forward-coupling models Its raison d'etre consists of not being temple...
ForwardCouplingModel(SuperGeometry< T, 3 > &sGeometry, std::shared_ptr< DragModel< T, Particle > > dragModel)
Constructor.
SuperGeometry< T, 3 > & _sGeometry
LaddForwardCouplingModel(UnitConverter< T, Lattice > &converter, SuperLattice< T, Lattice > &sLattice, SuperGeometry< T, 3 > &sGeometry, std::shared_ptr< DragModel< T, Particle > > dragModel)
Constructor.
Abstact class for all the local forward-coupling models, viz., momentum coupling from fluid to partic...
LocalBaseForwardCouplingModel(UnitConverter< T, Lattice > &converter, SuperLattice< T, Lattice > &sLattice, SuperGeometry< T, 3 > &sGeometry, std::shared_ptr< DragModel< T, Particle > > dragModel)
Constructor.
virtual bool operator()(Particle< T > *p, int globic) override
Class operator to apply the coupling, for overload.
NaiveForwardCouplingModel(UnitConverter< T, Lattice > &converter, SuperLattice< T, Lattice > &sLattice, SuperGeometry< T, 3 > &sGeometry, std::shared_ptr< DragModel< T, Particle > > dragModel)
Constructor.
NaiveNonLocalForwardCouplingModel(UnitConverter< T, Lattice > &converter, SuperLattice< T, Lattice > &sLattice, SuperGeometry< T, 3 > &sGeometry, std::shared_ptr< DragModel< T, Particle > > dragModel, std::shared_ptr< SmoothingFunctional< T, Lattice > > smoothingFunctional)
Constructor.
Abstact class for all the non-local forward-coupling models, viz., momentum coupling from fluid to pa...
NonLocalBaseForwardCouplingModel(UnitConverter< T, Lattice > &converter, SuperLattice< T, Lattice > &sLattice, SuperGeometry< T, 3 > &sGeometry, std::shared_ptr< DragModel< T, Particle > > dragModel, std::shared_ptr< SmoothingFunctional< T, Lattice > > smoothingFunctional)
Constructor.
virtual bool operator()(Particle< T > *p, int globic) override
Class operator to apply the coupling, for overload.
Abstact class for all the smoothing functionals.
Representation of a statistic for a parallel 2D geometry.
Super class maintaining block lattices for a cuboid decomposition.
Conversion between physical and lattice units, as well as discretization.
virtual bool operator()(Particle< T > *p, int globic) override
Class operator to apply the coupling, for overload.
vanWachemForwardCouplingModel(UnitConverter< T, Lattice > &converter, SuperLattice< T, Lattice > &sLattice, SuperGeometry< T, 3 > &sGeometry, std::shared_ptr< DragModel< T, Particle > > dragModel, std::shared_ptr< SmoothingFunctional< T, Lattice > > smoothingFunctional, int nVoxelInterpPoints)
Constructor.
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
Top level namespace for all of OpenLB.