OpenLB 1.7
Loading...
Searching...
No Matches
collision.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2021 Adrian Kummerlaender
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#ifndef DYNAMICS_COLLISION_H
25#define DYNAMICS_COLLISION_H
26
27#include "lbm.h"
28#include "descriptorField.h"
29#include "rtlbmDescriptors.h"
30
32
33namespace olb {
34
35namespace collision {
36
37struct None {
38 using parameters = typename meta::list<>;
39
40 static std::string getName() {
41 return "None";
42 }
43
44 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
45 struct type {
46 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
47 CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform {
48 return {-1, -1};
49 };
50 };
51};
52
53struct BGK {
55
56 static std::string getName() {
57 return "BGK";
58 }
59
60 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
61 struct type {
62 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
63
64 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
65 CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform {
66 V fEq[DESCRIPTOR::q] { };
67 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
68 const V omega = parameters.template get<descriptors::OMEGA>();
69 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
70 cell[iPop] *= V{1} - omega;
71 cell[iPop] += omega * fEq[iPop];
72 }
73 return statistic;
74 };
75 };
76};
77
78struct RLB {
80
81 static std::string getName() {
82 return "RLB";
83 }
84
85 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
86 struct type {
87 using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
88 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
89
90 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
91 CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform {
93 MomentaF().computeStress(cell, pi);
94 V fEq[DESCRIPTOR::q] { };
95 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
96 const V omega = parameters.template get<descriptors::OMEGA>();
97 cell[0] = fEq[0] + (V{1} - omega) * equilibrium<DESCRIPTOR>::template fromPiToFneq<V>(0, pi);
98 for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) {
99 cell[iPop] = fEq[iPop];
100 cell[iPop+DESCRIPTOR::q/2] = fEq[iPop+DESCRIPTOR::q/2];
101 V fNeq = (V{1} - omega) * equilibrium<DESCRIPTOR>::template fromPiToFneq<V>(iPop, pi);
102 cell[iPop] += fNeq;
103 cell[iPop+DESCRIPTOR::q/2] += fNeq;
104 }
105 return statistic;
106 };
107 };
108};
109
112
113 static std::string getName() {
114 return "RLB";
115 }
116
117 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
118 struct type {
119 using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
120 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
121
122 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
124 V fEq[DESCRIPTOR::q] { };
125 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
126 V j1[DESCRIPTOR::d] { };
127 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop ) {
128 for (int iD=0; iD < DESCRIPTOR::d; ++iD ) {
129 j1[iD] += descriptors::c<DESCRIPTOR>(iPop,iD) * (cell[iPop] - fEq[iPop]);
130 }
131 }
132 const V omega = parameters.template get<descriptors::OMEGA>();
133 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop ) {
134 V fNeq{};
135 for (int iD=0; iD < DESCRIPTOR::d; ++iD) {
136 fNeq += descriptors::c<DESCRIPTOR>(iPop,iD) * j1[iD];
137 }
138 fNeq *= descriptors::t<V,DESCRIPTOR>(iPop) * descriptors::invCs2<V,DESCRIPTOR>();
139 cell[iPop] = fEq[iPop] + (V{1} - omega) * fNeq;
140 }
141 return statistic;
142 };
143 };
144};
145
148
149 static std::string getName() {
150 return "ConstRhoBGK";
151 }
152
153 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
154 struct type {
155 using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
156 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
157
158 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
160 V fEq[DESCRIPTOR::q];
161 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
162
163 V deltaRho = V{1} - parameters.template get<statistics::AVERAGE_RHO>();
164 V ratioRho = V{1} + deltaRho / statistic.rho;
165
166 const V omega = parameters.template get<descriptors::OMEGA>();
167
168 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
169 cell[iPop] = ratioRho * (fEq[iPop] + descriptors::t<V,DESCRIPTOR>(iPop))
170 - descriptors::t<V,DESCRIPTOR>(iPop)
171 + (V{1} - omega) * (cell[iPop] - fEq[iPop]);
172 }
173 return {statistic.rho+deltaRho, statistic.uSqr};
174 };
175 };
176};
177
179 struct OMEGA : public descriptors::FIELD_BASE<0,0,1> { };
180
182
183 static std::string getName() {
184 return "PerPopulationBGK";
185 }
186
187 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
188 struct type {
189 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
190
191 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
193 V fEq[DESCRIPTOR::q] { };
194 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
195 const auto omega = parameters.template get<OMEGA>();
196 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
197 cell[iPop] *= V{1} - omega[iPop];
198 cell[iPop] += omega[iPop] * fEq[iPop];
199 }
200 return statistic;
201 };
202 };
203};
204
205struct TRT {
206 struct MAGIC : public descriptors::FIELD_BASE<1> { };
207
209
210 static std::string getName() {
211 return "TRT";
212 }
213
214 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
215 struct type {
216 using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
217 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
218
219 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
221 V fPlus[DESCRIPTOR::q], fMinus[DESCRIPTOR::q];
222 V fEq[DESCRIPTOR::q], fEqPlus[DESCRIPTOR::q], fEqMinus[DESCRIPTOR::q];
223 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
224 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
225 fPlus[iPop] = 0.5 * (cell[iPop] + cell[descriptors::opposite<DESCRIPTOR>(iPop)]);
226 fMinus[iPop] = 0.5 * (cell[iPop] - cell[descriptors::opposite<DESCRIPTOR>(iPop)]);
227 }
228 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
229 fEqPlus[iPop] = 0.5 * (fEq[iPop] + fEq[descriptors::opposite<DESCRIPTOR>(iPop)]);
230 fEqMinus[iPop] = 0.5 * (fEq[iPop] - fEq[descriptors::opposite<DESCRIPTOR>(iPop)]);
231 }
232
233 const V omega = parameters.template get<descriptors::OMEGA>();
234 const V magic = parameters.template get<MAGIC>();
235 const V omega2 = V{1} / (magic / (V{1} / omega - 0.5) + V{0.5});
236
237 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
238 cell[iPop] -= omega * (fPlus[iPop] - fEqPlus[iPop]) + omega2 * (fMinus[iPop] - fEqMinus[iPop]);
239 }
240 return statistic;
241 };
242 };
243};
244
245struct Revert {
246 using parameters = typename meta::list<>;
247
248 static std::string getName() {
249 return "Revert";
250 }
251
252 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
253 struct type {
254 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
256 for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) {
257 V cell_iPop = cell[iPop];
258 cell[iPop] = cell[descriptors::opposite<DESCRIPTOR>(iPop)];
259 cell[descriptors::opposite<DESCRIPTOR>(iPop)] = cell_iPop;
260 }
261 return {V{-1}, V{-1}};
262 };
263 };
264};
265
267
274template <typename COLLISION>
276 using parameters = typename COLLISION::parameters;
277
278 static std::string getName() {
279 return "NguyenLaddCorrection<" + COLLISION::getName() + ">";
280 }
281
282 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
283 struct type {
284 using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
285 using CollisionO = typename COLLISION::template type<DESCRIPTOR,MOMENTA,EQUILIBRIUM>;
286
287 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
289 CollisionO().apply(cell, parameters);
290 V rho, u[DESCRIPTOR::d];
291 MomentaF().computeRhoU(cell, rho, u);
292 for (int iPop=1; iPop < DESCRIPTOR::q; ++iPop) {
293 for (int iD=0; iD < DESCRIPTOR::d; ++iD) {
294 cell[iPop] += 2 * rho * u[iD]*descriptors::c<DESCRIPTOR>(iPop,iD) * descriptors::t<V,DESCRIPTOR>(iPop) * descriptors::invCs2<V,DESCRIPTOR>();
295 }
296 }
297 return {rho, util::normSqr<V,DESCRIPTOR::d>(u)};
298 };
299 };
300};
301
303 struct RF : public descriptors::FIELD_BASE<1> { };
304 using parameters = typename meta::list<RF>;
305
306 static std::string getName() {
307 return "PartialBounceBack";
308 }
309
310 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
311 struct type {
312 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
314 for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) {
315 std::swap(cell[iPop], cell[iPop+DESCRIPTOR::q/2]);
316 }
317 const V rf = parameters.template get<RF>();
318 for (int iPop=1; iPop < DESCRIPTOR::q; ++iPop) {
319 cell[iPop] = (rf - V{1}) * (cell[iPop] + descriptors::t<V,DESCRIPTOR>(iPop))
320 - descriptors::t<V,DESCRIPTOR>(iPop);
321 }
322 return {-1, -1};
323 };
324 };
325};
326
327struct Poisson {
328 struct SINK : public descriptors::FIELD_BASE<1> { };
330
331 static std::string getName() {
332 return "Poisson";
333 }
334
335 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
336 struct type {
337 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
338
339 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
341 V fEq[DESCRIPTOR::q] { };
342 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
343 const V omega = parameters.template get<descriptors::OMEGA>();
344 const V sink = parameters.template get<SINK>();
345 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
346 cell[iPop] = (cell[iPop] + descriptors::t<V,DESCRIPTOR>(iPop))
347 - omega * ((cell[iPop] + descriptors::t<V,DESCRIPTOR>(iPop)) - descriptors::t<V,DESCRIPTOR>(iPop) * statistic.rho)
348 - sink * (cell[iPop] + descriptors::t<V,DESCRIPTOR>(iPop))
349 - descriptors::t<V,DESCRIPTOR>(iPop);
350 }
351 return statistic;
352 };
353 };
354};
355
356struct P1 {
357 struct SCATTERING : public descriptors::FIELD_BASE<1> { };
358 struct ABSORPTION : public descriptors::FIELD_BASE<1> { };
360
361 static std::string getName() {
362 return "P1";
363 }
364
365 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
366 struct type {
367 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
368
369 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
371 V fEq[DESCRIPTOR::q] { };
372 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
373 const V scattering = parameters.template get<SCATTERING>();
374 const V absorption = parameters.template get<ABSORPTION>();
375 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
376 cell[iPop] -= scattering * descriptors::norm_c<V,DESCRIPTOR>(iPop)*(cell[iPop] - fEq[iPop])
377 + absorption * descriptors::norm_c<V,DESCRIPTOR>(iPop)*(cell[iPop] + descriptors::t<V,DESCRIPTOR>(iPop));
378 }
379 return statistic;
380 };
381 };
382};
383
385 using parameters = typename meta::list<>;
386
387 static std::string getName() {
388 return "FixedEquilibrium";
389 }
390
391 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
392 struct type {
393 using EquilibriumF = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
394
395 template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
397 V fEq[DESCRIPTOR::q] { };
398 const auto statistic = EquilibriumF().compute(cell, parameters, fEq);
399 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
400 cell[iPop] = fEq[iPop];
401 }
402 return statistic;
403 };
404 };
405};
406
407}
408
409}
410
411#endif
412
413#include "collision.cse.h"
Interface for post-processing steps – header file.
Top level namespace for all of OpenLB.
#define any_platform
Define preprocessor macros for device-side functions, constant storage.
Definition platform.h:78
– header file
Return value of any collision.
Definition interface.h:43
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:123
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:120
typename MOMENTA::template type< DESCRIPTOR > MomentaF
Definition collision.h:119
typename meta::list< descriptors::OMEGA > parameters
Definition collision.h:111
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:65
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:62
typename meta::list< descriptors::OMEGA > parameters
Definition collision.h:54
static std::string getName()
Definition collision.h:56
typename MOMENTA::template type< DESCRIPTOR > MomentaF
Definition collision.h:155
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:156
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:159
static std::string getName()
Definition collision.h:149
typename meta::list< descriptors::OMEGA, statistics::AVERAGE_RHO > parameters
Definition collision.h:147
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:393
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:396
typename meta::list<> parameters
Definition collision.h:385
static std::string getName()
Definition collision.h:387
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:288
typename MOMENTA::template type< DESCRIPTOR > MomentaF
Definition collision.h:284
typename COLLISION::template type< DESCRIPTOR, MOMENTA, EQUILIBRIUM > CollisionO
Definition collision.h:285
Nguyen-Ladd Velocity Correction using momenta-defined velocity.
Definition collision.h:275
typename COLLISION::parameters parameters
Definition collision.h:276
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:47
static std::string getName()
Definition collision.h:40
typename meta::list<> parameters
Definition collision.h:38
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:367
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:370
static std::string getName()
Definition collision.h:361
typename meta::list< SCATTERING, ABSORPTION > parameters
Definition collision.h:359
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:313
static std::string getName()
Definition collision.h:306
typename meta::list< RF > parameters
Definition collision.h:304
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:192
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:189
typename meta::list< OMEGA > parameters
Definition collision.h:181
static std::string getName()
Definition collision.h:183
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:337
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:340
static std::string getName()
Definition collision.h:331
typename meta::list< descriptors::OMEGA, SINK > parameters
Definition collision.h:329
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:91
typename MOMENTA::template type< DESCRIPTOR > MomentaF
Definition collision.h:87
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:88
typename meta::list< descriptors::OMEGA > parameters
Definition collision.h:79
static std::string getName()
Definition collision.h:81
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:255
static std::string getName()
Definition collision.h:248
typename meta::list<> parameters
Definition collision.h:246
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition collision.h:220
typename MOMENTA::template type< DESCRIPTOR > MomentaF
Definition collision.h:216
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
Definition collision.h:217
typename meta::list< descriptors::OMEGA, MAGIC > parameters
Definition collision.h:208
static std::string getName()
Definition collision.h:210
Base of a field whose size is defined by [C,U_1,...,U_N]^T * [1,V_1,...V_N].
Plain wrapper for list of types.
Definition meta.h:276
Compute number of elements of a symmetric d-dimensional tensor.
Definition util.h:210