OpenLB 1.7
Loading...
Searching...
No Matches
dynamics.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2006, 2007 Jonas Latt, Mathias J. Krause
4 * 2021 Adrian Kummerlaender
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#ifndef LB_LEGACY_DYNAMICS_H
26#define LB_LEGACY_DYNAMICS_H
27
28#include "dynamics/interface.h"
29
30namespace olb {
31
32namespace legacy {
33
34template <typename T, typename DESCRIPTOR, typename MOMENTA>
35struct BasicDynamics : public dynamics::CustomCollision<T,DESCRIPTOR,MOMENTA> {
36 T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d]) const override {
37 return equilibrium<DESCRIPTOR>::secondOrder(iPop, rho, u);
38 };
39
40 std::type_index id() override {
42 }
43
45 throw std::bad_function_call();
46 }
47
48};
49
50template<typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
51class BGKdynamics : public BasicDynamics<T,DESCRIPTOR,MOMENTA> {
52public:
53 BGKdynamics(T omega);
56
57 virtual T getOmega() const {
58 return _omega;
59 }
60
61private:
62 using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
63
64 T _omega;
65};
66
67template<typename T, typename DESCRIPTOR, typename MOMENTA>
69 : BasicDynamics<T,DESCRIPTOR,MOMENTA>(),
70 _omega(omega)
71{ }
72
73template<typename T, typename DESCRIPTOR, typename MOMENTA>
75{
76 T rho, u[DESCRIPTOR::d];
77 MomentaF().computeRhoU(cell, rho, u);
78 T uSqr = lbm<DESCRIPTOR>::bgkCollision(cell, rho, u, _omega);
79 return {rho, uSqr};
80}
81
82template<typename T, typename DESCRIPTOR>
83class NoLatticeDynamics : public Dynamics<T,DESCRIPTOR> {
84public:
86 NoLatticeDynamics(T rho = T(1) );
88 T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d]) const override any_platform;
92 T computeRho(ConstCell<T,DESCRIPTOR>& cell) const override;
94 void computeU (
96 T u[DESCRIPTOR::d] ) const override;
98 void computeJ (
100 T j[DESCRIPTOR::d] ) const override;
102 void computeStress (
104 T rho, const T u[DESCRIPTOR::d],
105 T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
106 void computeRhoU (
108 T& rho, T u[DESCRIPTOR::d]) const override;
109 void computeAllMomenta (
111 T& rho, T u[DESCRIPTOR::d],
112 T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
114 void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
116 void defineU(Cell<T,DESCRIPTOR>& cell,
117 const T u[DESCRIPTOR::d]) override;
119 void defineRhoU (
120 Cell<T,DESCRIPTOR>& cell,
121 T rho, const T u[DESCRIPTOR::d]) override;
123 void defineAllMomenta (
124 Cell<T,DESCRIPTOR>& cell,
125 T rho, const T u[DESCRIPTOR::d],
126 const T pi[util::TensorVal<DESCRIPTOR >::n] ) override;
127
128 std::type_index id() override {
129 return typeid(NoLatticeDynamics<T,DESCRIPTOR>);
130 }
131
133 throw std::bad_function_call();
134 }
135
136private:
138 T _rho;
139};
140
144template<typename T, typename DESCRIPTOR>
145class OffDynamics : public NoLatticeDynamics<T,DESCRIPTOR> {
146public:
148 OffDynamics(const T _location[DESCRIPTOR::d]);
150 OffDynamics(const T _location[DESCRIPTOR::d], T _distances[DESCRIPTOR::q]);
152 T computeRho(ConstCell<T,DESCRIPTOR>& cell) const override;
154 void computeU(ConstCell<T,DESCRIPTOR>& cell, T u[DESCRIPTOR::d] ) const override;
156 void setBoundaryIntersection(int iPop, T distance);
158 bool getBoundaryIntersection(int iPop, T intersection[DESCRIPTOR::d]);
160 void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
162 void defineRho(int iPop, T rho);
164 void defineU(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d]) override;
166 void defineU(const T u[DESCRIPTOR::d]);
168 void defineU(int iPop, const T u[DESCRIPTOR::d]);
170 T getVelocityCoefficient(int iPop);
171
172 std::type_index id() override {
173 return typeid(OffDynamics<T,DESCRIPTOR>);
174 }
175
176
177private:
178 T _rho;
179 T _u[DESCRIPTOR::q][DESCRIPTOR::d];
180 T location[DESCRIPTOR::d];
181 T distances[DESCRIPTOR::q];
182 T boundaryIntersection[DESCRIPTOR::q][DESCRIPTOR::d];
183 T velocityCoefficient[DESCRIPTOR::q];
184};
185
187
188template<typename T, typename DESCRIPTOR>
190{
191 this->getName() = "NoLatticeDynamics";
192}
193
194template<typename T, typename DESCRIPTOR>
195T NoLatticeDynamics<T,DESCRIPTOR>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d]) const
196{
197 return T();
198}
199
200template<typename T, typename DESCRIPTOR>
204
205template<typename T, typename DESCRIPTOR>
210
211template<typename T, typename DESCRIPTOR>
214 T u[DESCRIPTOR::d]) const
215{
216 for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
217 u[iD] = T();
218 }
219}
220
221template<typename T, typename DESCRIPTOR>
224 T j[DESCRIPTOR::d]) const
225{
226 for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
227 j[iD] = T();
228 }
229}
230
231template<typename T, typename DESCRIPTOR>
234 T rho, const T u[DESCRIPTOR::d],
236{
237 for (int iPi=0; iPi<util::TensorVal<DESCRIPTOR >::n; ++iPi) {
238 pi[iPi] = T();//std::numeric_limits<T>::signaling_NaN();
239 }
240}
241
242template<typename T, typename DESCRIPTOR>
245 T& rho, T u[DESCRIPTOR::d]) const
246{
247 rho = computeRho(cell);
248 computeU(cell, u);
249}
250
251template<typename T, typename DESCRIPTOR>
254 T& rho, T u[DESCRIPTOR::d],
256{
257 computeRhoU(cell, rho, u);
258 computeStress(cell, rho, u, pi);
259}
260
261template<typename T, typename DESCRIPTOR>
264
265template<typename T, typename DESCRIPTOR>
267 Cell<T,DESCRIPTOR>& cell,
268 const T u[DESCRIPTOR::d])
269{ }
270
271template<typename T, typename DESCRIPTOR>
273 Cell<T,DESCRIPTOR>& cell,
274 T rho, const T u[DESCRIPTOR::d])
275{ }
276
277template<typename T, typename DESCRIPTOR>
279 Cell<T,DESCRIPTOR>& cell,
280 T rho, const T u[DESCRIPTOR::d],
282{ }
283
285
286template<typename T, typename DESCRIPTOR>
287OffDynamics<T,DESCRIPTOR>::OffDynamics(const T _location[DESCRIPTOR::d])
288{
289 this->getName() = "OffDynamics";
290 typedef DESCRIPTOR L;
291 for (int iD = 0; iD < L::d; iD++) {
292 location[iD] = _location[iD];
293 }
294 for (int iPop = 0; iPop < L::q; iPop++) {
295 distances[iPop] = -1;
296 velocityCoefficient[iPop] = 0;
297 for (int iD = 0; iD < L::d; iD++) {
298 boundaryIntersection[iPop][iD] = _location[iD];
299 _u[iPop][iD] = T();
300 }
301 }
302 _rho=T(1);
303}
304
305template<typename T, typename DESCRIPTOR>
306OffDynamics<T,DESCRIPTOR>::OffDynamics(const T _location[DESCRIPTOR::d], T _distances[DESCRIPTOR::q])
307{
308 this->getName() = "OffDynamics";
309 typedef DESCRIPTOR L;
310 for (int iD = 0; iD < L::d; iD++) {
311 location[iD] = _location[iD];
312 }
313 for (int iPop = 0; iPop < L::q; iPop++) {
314 distances[iPop] = _distances[iPop];
315 velocityCoefficient[iPop] = 0;
316 for (int iD = 0; iD < L::d; iD++) {
317 boundaryIntersection[iPop][iD] = _location[iD] - _distances[iPop]*descriptors::c<L>(iPop,iD);
318 _u[iPop][iD] = T();
319 }
320 }
321 _rho=T(1);
322}
323
324template<typename T, typename DESCRIPTOR>
326{
327 /*typedef DESCRIPTOR L;
328 T rhoTmp = T();
329 T counter = T();
330 int counter2 = int();
331 for (int iPop = 0; iPop < L::q; iPop++) {
332 if (distances[iPop] != -1) {
333 rhoTmp += (cell[iPop] + descriptors::t<T,L>(iPop))*descriptors::t<T,L>(iPop);
334 counter += descriptors::t<T,L>(iPop);
335 counter2++;
336 }
337 }
338 //if (rhoTmp/counter + 1<0.1999) std::cout << rhoTmp/counter2 + 1 <<std::endl;
339 //if (rhoTmp/counter + 1>1.001) std::cout << rhoTmp/counter2 + 1 <<std::endl;
340 return rhoTmp/counter/counter;*/
341 return _rho;
342}
343
344template<typename T, typename DESCRIPTOR>
346{
347 typedef DESCRIPTOR L;
348 for (int iD = 0; iD < L::d; iD++) {
349 u[iD] = T();
350 }
351 int counter = 0;
352 for (int iPop = 0; iPop < L::q; iPop++) {
353 if ( !util::nearZero(distances[iPop]+1) ) {
354 for (int iD = 0; iD < L::d; iD++) {
355 u[iD] += _u[iPop][iD];
356 }
357 counter++;
358 }
359 }
360 if (counter!=0) {
361 for (int iD = 0; iD < L::d; iD++) {
362 u[iD] /= counter;
363 }
364 }
365 return;
366}
367
368template<typename T, typename DESCRIPTOR>
370{
373 typedef DESCRIPTOR L;
374 distances[iPop] = distance;
375 for (int iD = 0; iD < L::d; iD++) {
376 boundaryIntersection[iPop][iD] = location[iD] - distance*descriptors::c<L>(iPop,iD);
377 }
378}
379
380template<typename T, typename DESCRIPTOR>
381bool OffDynamics<T,DESCRIPTOR>::getBoundaryIntersection(int iPop, T intersection[DESCRIPTOR::d])
382{
383 typedef DESCRIPTOR L;
384 if ( !util::nearZero(distances[iPop]+1) ) {
385 for (int iD = 0; iD < L::d; iD++) {
386 intersection[iD] = boundaryIntersection[iPop][iD];
387 }
388 return true;
389 }
390 return false;
391}
392
393template<typename T, typename DESCRIPTOR>
395{
396 _rho=rho;
397}
398
399template<typename T, typename DESCRIPTOR>
401{
402 _rho=rho;
403}
404
405template<typename T, typename DESCRIPTOR>
407 Cell<T,DESCRIPTOR>& cell,
408 const T u[DESCRIPTOR::d])
409{
410 defineU(u);
411}
412
413template<typename T, typename DESCRIPTOR>
414void OffDynamics<T,DESCRIPTOR>::defineU(const T u[DESCRIPTOR::d])
415{
416 typedef DESCRIPTOR L;
417 for (int iPop = 0; iPop < L::q; iPop++) {
418 if ( !util::nearZero(distances[iPop]+1) ) {
419 defineU(iPop, u);
420 }
421 }
422}
423
425
429template<typename T, typename DESCRIPTOR>
431 int iPop, const T u[DESCRIPTOR::d])
432{
433 OLB_PRECONDITION(distances[iPop] != -1)
434 typedef DESCRIPTOR L;
435 velocityCoefficient[iPop] = 0;
436 // scalar product of c(iPop) and u
437 for (int sum = 0; sum < L::d; sum++) { // +/- problem because of first stream than postprocess
438 velocityCoefficient[iPop] -= descriptors::c<L>(iPop,sum)*u[sum];
439 }
440 // compute summand for boundary condition
441 velocityCoefficient[iPop] *= 2*descriptors::invCs2<T,L>() * descriptors::t<T,L>(iPop);
442
443 for (int iD = 0; iD < L::d; iD++) {
444 _u[iPop][iD] = u[iD];
445 }
446}
447
448template<typename T, typename DESCRIPTOR>
450{
451 return velocityCoefficient[iPop];
452}
453
454}
455
456}
457
458#endif
Platform-abstracted block lattice for external access and inter-block interaction.
Highest-level interface to Cell data.
Definition cell.h:148
Highest-level interface to read-only Cell data.
Definition cell.h:53
virtual T getOmega() const
Definition dynamics.h:57
CellStatistic< T > collide(Cell< T, DESCRIPTOR > &cell) override
Collision step.
Definition dynamics.h:74
T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d]) const override any_platform
Yields 0;.
Definition dynamics.h:195
CellStatistic< T > collide(Cell< T, DESCRIPTOR > &cell) override
Collision step.
Definition dynamics.h:201
std::type_index id() override
Expose unique type-identifier for RTTI.
Definition dynamics.h:128
void defineAllMomenta(Cell< T, DESCRIPTOR > &cell, T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal< DESCRIPTOR >::n]) override
Does nothing.
Definition dynamics.h:278
void computeAllMomenta(ConstCell< T, DESCRIPTOR > &cell, T &rho, T u[DESCRIPTOR::d], T pi[util::TensorVal< DESCRIPTOR >::n]) const override
Compute all momenta up to second order.
Definition dynamics.h:252
void computeJ(ConstCell< T, DESCRIPTOR > &cell, T j[DESCRIPTOR::d]) const override
Yields 0;.
Definition dynamics.h:222
void defineRho(Cell< T, DESCRIPTOR > &cell, T rho) override
Does nothing.
Definition dynamics.h:262
void computeRhoU(ConstCell< T, DESCRIPTOR > &cell, T &rho, T u[DESCRIPTOR::d]) const override
Compute fluid velocity and particle density.
Definition dynamics.h:243
NoLatticeDynamics(T rho=T(1))
You may fix a fictitious density value on no dynamics node via this constructor.
Definition dynamics.h:189
void computeStress(ConstCell< T, DESCRIPTOR > &cell, T rho, const T u[DESCRIPTOR::d], T pi[util::TensorVal< DESCRIPTOR >::n]) const override
Yields NaN.
Definition dynamics.h:232
T computeRho(ConstCell< T, DESCRIPTOR > &cell) const override
Yields 1;.
Definition dynamics.h:206
void defineRhoU(Cell< T, DESCRIPTOR > &cell, T rho, const T u[DESCRIPTOR::d]) override
Does nothing.
Definition dynamics.h:272
void computeU(ConstCell< T, DESCRIPTOR > &cell, T u[DESCRIPTOR::d]) const override
Yields 0;.
Definition dynamics.h:212
AbstractParameters< T, DESCRIPTOR > & getParameters(BlockLattice< T, DESCRIPTOR > &block) override
Parameters access for legacy post processors.
Definition dynamics.h:132
void defineU(Cell< T, DESCRIPTOR > &cell, const T u[DESCRIPTOR::d]) override
Does nothing.
Definition dynamics.h:266
Dynamics for offLattice boundary conditions OffDynamics are basically NoLatticeDynamics with the addi...
Definition dynamics.h:145
OffDynamics(const T _location[DESCRIPTOR::d])
Constructor.
Definition dynamics.h:287
T getVelocityCoefficient(int iPop)
Get VelocitySummand for Bouzidi-Boundary Condition.
Definition dynamics.h:449
void defineRho(Cell< T, DESCRIPTOR > &cell, T rho) override
Set particle density on the cell.
Definition dynamics.h:394
void computeU(ConstCell< T, DESCRIPTOR > &cell, T u[DESCRIPTOR::d]) const override
Returns an average of the locally stored u.
Definition dynamics.h:345
std::type_index id() override
Expose unique type-identifier for RTTI.
Definition dynamics.h:172
T computeRho(ConstCell< T, DESCRIPTOR > &cell) const override
Returns local stored rho which is updated if the bc is used as velocity!=0 condition.
Definition dynamics.h:325
void defineU(Cell< T, DESCRIPTOR > &cell, const T u[DESCRIPTOR::d]) override
Set fluid velocity on the cell.
Definition dynamics.h:406
void setBoundaryIntersection(int iPop, T distance)
Set Intersection of the link and the boundary.
Definition dynamics.h:369
bool getBoundaryIntersection(int iPop, T intersection[DESCRIPTOR::d])
Get Intersection of the link and the boundary.
Definition dynamics.h:381
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
Top level namespace for all of OpenLB.
dynamics::Tuple< T, DESCRIPTOR, MOMENTA, std::conditional_t<(DESCRIPTOR::d==3 &&DESCRIPTOR::q==7)||(DESCRIPTOR::d==2 &&DESCRIPTOR::q==5), equilibria::FirstOrder, equilibria::SecondOrder >, collision::BGK > BGKdynamics
Common BGK collision step.
Definition dynamics.h:72
#define OLB_PRECONDITION(COND)
Definition olbDebug.h:46
#define any_platform
Define preprocessor macros for device-side functions, constant storage.
Definition platform.h:78
Dynamic access interface for FIELD-valued parameters.
Return value of any collision.
Definition interface.h:43
Interface for per-cell dynamics.
Definition interface.h:54
virtual std::string getName() const
Return human-readable name.
Definition interface.h:63
virtual CellStatistic< T > collide(Cell< T, DESCRIPTOR > &cell)
Perform purely-local collision step on Cell interface (legacy, to be deprecated)
Definition interface.h:74
static V secondOrder(int iPop, const RHO &rho, const U &u, const USQR &uSqr) any_platform
Computation of equilibrium distribution, second order in u.
Definition lbm.h:51
static V bgkCollision(CELL &cell, const RHO &rho, const VELOCITY &u, const OMEGA &omega) any_platform
BGK collision step.
Definition lbm.h:290
T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d]) const override
Return iPop equilibrium for given first and second momenta.
Definition dynamics.h:36
AbstractParameters< T, DESCRIPTOR > & getParameters(BlockLattice< T, DESCRIPTOR > &block) override
Parameters access for legacy post processors.
Definition dynamics.h:44
std::type_index id() override
Expose unique type-identifier for RTTI.
Definition dynamics.h:40
Compute number of elements of a symmetric d-dimensional tensor.
Definition util.h:210