OpenLB 1.7
Loading...
Searching...
No Matches
entropicDynamics.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2006, 2007, 2017 Orestis Malaspinas, Jonas Latt, Mathias J. Krause
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
28#ifndef ENTROPIC_LB_DYNAMICS_HH
29#define ENTROPIC_LB_DYNAMICS_HH
30
31#include <algorithm>
32#include <limits>
33#include "lbm.h"
34#include "entropicDynamics.h"
35#include "entropicLbHelpers.h"
36
37namespace olb {
38
39//==============================================================================//
41//==============================================================================//
45template<typename T, typename DESCRIPTOR, typename MOMENTA>
46EntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::EntropicEqDynamics (T omega_)
47 : legacy::BasicDynamics<T,DESCRIPTOR,MOMENTA>(),
48 omega(omega_)
49{
50 this->getName() = "EntropicEqDynamics";
51}
52
53template<typename T, typename DESCRIPTOR, typename MOMENTA>
54T EntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
55{
56 return entropicLbHelpers<T,DESCRIPTOR>::equilibrium(iPop,rho,u);
57}
58
59template<typename T, typename DESCRIPTOR, typename MOMENTA>
60CellStatistic<T> EntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::collide (
61 Cell<T,DESCRIPTOR>& cell)
62{
63 typedef DESCRIPTOR L;
64 typedef entropicLbHelpers<T,DESCRIPTOR> eLbH;
65
66 T rho, u[DESCRIPTOR::d];
67 MOMENTA().computeRhoU(cell, rho, u);
68 T uSqr = util::normSqr<T,L::d>(u);
69
70 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
71 cell[iPop] += omega * (eLbH::equilibrium(iPop,rho,u) - cell[iPop]);
72 }
73
74 //statistics.incrementStats(rho, uSqr);
75}
76
77template<typename T, typename DESCRIPTOR, typename MOMENTA>
78T EntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::getOmega() const
79{
80 return omega;
81}
82
83template<typename T, typename DESCRIPTOR, typename MOMENTA>
84void EntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::setOmega(T omega_)
85{
86 omega = omega_;
87}
88
89
90//====================================================================//
92//====================================================================//
93
96template<typename T, typename DESCRIPTOR, typename MOMENTA>
97ForcedEntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::ForcedEntropicEqDynamics (T omega_)
98 : legacy::BasicDynamics<T,DESCRIPTOR,MOMENTA>(),
99 omega(omega_)
100{
101 this->getName() = "ForcedEntropicEqDynamics";
102}
103
104template<typename T, typename DESCRIPTOR, typename MOMENTA>
105T ForcedEntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
106{
107 return entropicLbHelpers<T,DESCRIPTOR>::equilibrium(iPop,rho,u);
108}
109
110
111template<typename T, typename DESCRIPTOR, typename MOMENTA>
112CellStatistic<T> ForcedEntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::collide (
113 Cell<T,DESCRIPTOR>& cell)
114{
115 typedef DESCRIPTOR L;
116 typedef entropicLbHelpers<T,DESCRIPTOR> eLbH;
117
118 T rho, u[DESCRIPTOR::d];
119 MOMENTA().computeRhoU(cell, rho, u);
120
121 T* force = cell.template getFieldPointer<descriptors::FORCE>();
122 for (int iDim=0; iDim<DESCRIPTOR::d; ++iDim) {
123 u[iDim] += force[iDim] / (T)2.;
124 }
125 T uSqr = util::normSqr<T,L::d>(u);
126
127 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
128 cell[iPop] += omega * (eLbH::equilibrium(iPop,rho,u) - cell[iPop]);
129 }
130
131 lbm<DESCRIPTOR>::addExternalForce(cell, u, omega);
132
133 //statistics.incrementStats(rho, uSqr);
134}
135
136template<typename T, typename DESCRIPTOR, typename MOMENTA>
141
142template<typename T, typename DESCRIPTOR, typename MOMENTA>
143void ForcedEntropicEqDynamics<T,DESCRIPTOR,MOMENTA>::setOmega(T omega_)
144{
145 omega = omega_;
146}
147
148
149//==============================================================================//
151//==============================================================================//
155template<typename T, typename DESCRIPTOR, typename MOMENTA>
156EntropicDynamics<T,DESCRIPTOR,MOMENTA>::EntropicDynamics (T omega_)
157 : legacy::BasicDynamics<T,DESCRIPTOR,MOMENTA>(),
158 omega(omega_)
159{
160 this->getName() = "EntropicDynamics";
161}
162
163template<typename T, typename DESCRIPTOR, typename MOMENTA>
164T EntropicDynamics<T,DESCRIPTOR,MOMENTA>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
165{
166 return entropicLbHelpers<T,DESCRIPTOR>::equilibrium(iPop,rho,u);
167}
168
169template<typename T, typename DESCRIPTOR, typename MOMENTA>
170CellStatistic<T> EntropicDynamics<T,DESCRIPTOR,MOMENTA>::collide (
171 Cell<T,DESCRIPTOR>& cell)
172{
173 typedef DESCRIPTOR L;
174 typedef entropicLbHelpers<T,DESCRIPTOR> eLbH;
175
176 T rho, u[DESCRIPTOR::d];
177 MOMENTA().computeRhoU(cell, rho, u);
178 T uSqr = util::normSqr<T,L::d>(u);
179
180 T f[L::q], fEq[L::q], fNeq[L::q];
181 for (int iPop = 0; iPop < L::q; ++iPop) {
182 fEq[iPop] = eLbH::equilibrium(iPop,rho,u);
183 fNeq[iPop] = cell[iPop] - fEq[iPop];
184 f[iPop] = cell[iPop] + descriptors::t<T,L>(iPop);
185 fEq[iPop] += descriptors::t<T,L>(iPop);
186 }
187 //==============================================================================//
188 //============= Evaluation of alpha using a Newton Raphson algorithm ===========//
189 //==============================================================================//
190
191 T alpha = 2.0;
192 bool converged = getAlpha(alpha,f,fNeq);
193 if (!converged) {
194 std::cout << "Newton-Raphson failed to converge.\n";
195 exit(1);
196 }
197
198 OLB_ASSERT(converged,"Entropy growth failed to converge!");
199
200 T omegaTot = omega / 2.0 * alpha;
201 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
202 cell[iPop] *= (T)1-omegaTot;
203 cell[iPop] += omegaTot * (fEq[iPop]-descriptors::t<T,L>(iPop));
204 }
205
206 //statistics.incrementStats(rho, uSqr);
207}
208
209template<typename T, typename DESCRIPTOR, typename MOMENTA>
210T EntropicDynamics<T,DESCRIPTOR,MOMENTA>::getOmega() const
211{
212 return omega;
213}
214
215template<typename T, typename DESCRIPTOR, typename MOMENTA>
216void EntropicDynamics<T,DESCRIPTOR,MOMENTA>::setOmega(T omega_)
217{
218 omega = omega_;
219}
220
221template<typename T, typename DESCRIPTOR, typename MOMENTA>
222T EntropicDynamics<T,DESCRIPTOR,MOMENTA>::computeEntropy(const T f[])
223{
224 typedef DESCRIPTOR L;
225 T entropy = T();
226 for (int iPop = 0; iPop < L::q; ++iPop) {
227 OLB_ASSERT(f[iPop] > T(), "f[iPop] <= 0");
228 entropy += f[iPop]*util::log(f[iPop]/descriptors::t<T,L>(iPop));
229 }
230
231 return entropy;
232}
233
234template<typename T, typename DESCRIPTOR, typename MOMENTA>
235T EntropicDynamics<T,DESCRIPTOR,MOMENTA>::computeEntropyGrowth(const T f[], const T fNeq[], const T &alpha)
236{
237 typedef DESCRIPTOR L;
238
239 T fAlphaFneq[L::q];
240 for (int iPop = 0; iPop < L::q; ++iPop) {
241 fAlphaFneq[iPop] = f[iPop] - alpha*fNeq[iPop];
242 }
243
244 return computeEntropy(f) - computeEntropy(fAlphaFneq);
245}
246
247template<typename T, typename DESCRIPTOR, typename MOMENTA>
248T EntropicDynamics<T,DESCRIPTOR,MOMENTA>::computeEntropyGrowthDerivative(const T f[], const T fNeq[], const T &alpha)
249{
250 typedef DESCRIPTOR L;
251
252 T entropyGrowthDerivative = T();
253 for (int iPop = 0; iPop < L::q; ++iPop) {
254 T tmp = f[iPop] - alpha*fNeq[iPop];
255 OLB_ASSERT(tmp > T(), "f[iPop] - alpha*fNeq[iPop] <= 0");
256 entropyGrowthDerivative += fNeq[iPop]*(util::log(tmp/descriptors::t<T,L>(iPop)));
257 }
258
259 return entropyGrowthDerivative;
260}
261
262template<typename T, typename DESCRIPTOR, typename MOMENTA>
263bool EntropicDynamics<T,DESCRIPTOR,MOMENTA>::getAlpha(T &alpha, const T f[], const T fNeq[])
264{
265 const T epsilon = std::numeric_limits<T>::epsilon();
266
267 T alphaGuess = T();
268 const T var = 100.0;
269 const T errorMax = epsilon*var;
270 T error = 1.0;
271 int count = 0;
272 for (count = 0; count < 10000; ++count) {
273 T entGrowth = computeEntropyGrowth(f,fNeq,alpha);
274 T entGrowthDerivative = computeEntropyGrowthDerivative(f,fNeq,alpha);
275 if ((error < errorMax) || (util::fabs(entGrowth) < var*epsilon)) {
276 return true;
277 }
278 alphaGuess = alpha - entGrowth /
279 entGrowthDerivative;
280 error = util::fabs(alpha-alphaGuess);
281 alpha = alphaGuess;
282 }
283 return false;
284}
285
286//====================================================================//
288//====================================================================//
289
292template<typename T, typename DESCRIPTOR, typename MOMENTA>
293ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::ForcedEntropicDynamics (T omega_)
294 : legacy::BasicDynamics<T,DESCRIPTOR,MOMENTA>(),
295 omega(omega_)
296{
297 this->getName() = "ForcedEntropicDynamics";
298}
299
300template<typename T, typename DESCRIPTOR, typename MOMENTA>
301T ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
302{
303 return entropicLbHelpers<T,DESCRIPTOR>::equilibrium(iPop,rho,u);
304}
305
306
307template<typename T, typename DESCRIPTOR, typename MOMENTA>
308CellStatistic<T> ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::collide (
309 Cell<T,DESCRIPTOR>& cell)
310{
311 typedef DESCRIPTOR L;
312 typedef entropicLbHelpers<T,DESCRIPTOR> eLbH;
313
314 T rho, u[DESCRIPTOR::d];
315 MOMENTA().computeRhoU(cell, rho, u);
316 T uSqr = util::normSqr<T,L::d>(u);
317
318 T f[L::q], fEq[L::q], fNeq[L::q];
319 for (int iPop = 0; iPop < L::q; ++iPop) {
320 fEq[iPop] = eLbH::equilibrium(iPop,rho,u);
321 fNeq[iPop] = cell[iPop] - fEq[iPop];
322 f[iPop] = cell[iPop] + descriptors::t<T,L>(iPop);
323 fEq[iPop] += descriptors::t<T,L>(iPop);
324 }
325 //==============================================================================//
326 //============= Evaluation of alpha using a Newton Raphson algorithm ===========//
327 //==============================================================================//
328
329 T alpha = 2.0;
330 bool converged = getAlpha(alpha,f,fNeq);
331 if (!converged) {
332 std::cout << "Newton-Raphson failed to converge.\n";
333 exit(1);
334 }
335
336 OLB_ASSERT(converged,"Entropy growth failed to converge!");
337
338 T* force = cell.template getFieldPointer<descriptors::FORCE>();
339 for (int iDim=0; iDim<DESCRIPTOR::d; ++iDim) {
340 u[iDim] += force[iDim] / (T)2.;
341 }
342 uSqr = util::normSqr<T,L::d>(u);
343 T omegaTot = omega / 2.0 * alpha;
344 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
345 cell[iPop] *= (T)1-omegaTot;
346 cell[iPop] += omegaTot * eLbH::equilibrium(iPop,rho,u);
347 }
348 lbm<DESCRIPTOR>::addExternalForce(cell, u, omegaTot);
349
350 //statistics.incrementStats(rho, uSqr);
351}
352
353template<typename T, typename DESCRIPTOR, typename MOMENTA>
354T ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::getOmega() const
355{
356 return omega;
357}
358
359template<typename T, typename DESCRIPTOR, typename MOMENTA>
360void ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::setOmega(T omega_)
361{
362 omega = omega_;
363}
364
365template<typename T, typename DESCRIPTOR, typename MOMENTA>
366T ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::computeEntropy(const T f[])
367{
368 typedef DESCRIPTOR L;
369 T entropy = T();
370 for (int iPop = 0; iPop < L::q; ++iPop) {
371 OLB_ASSERT(f[iPop] > T(), "f[iPop] <= 0");
372 entropy += f[iPop]*util::log(f[iPop]/descriptors::t<T,L>(iPop));
373 }
374
375 return entropy;
376}
377
378template<typename T, typename DESCRIPTOR, typename MOMENTA>
379T ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::computeEntropyGrowth(const T f[], const T fNeq[], const T &alpha)
380{
381 typedef DESCRIPTOR L;
382
383 T fAlphaFneq[L::q];
384 for (int iPop = 0; iPop < L::q; ++iPop) {
385 fAlphaFneq[iPop] = f[iPop] - alpha*fNeq[iPop];
386 }
387
388 return computeEntropy(f) - computeEntropy(fAlphaFneq);
389}
390
391template<typename T, typename DESCRIPTOR, typename MOMENTA>
392T ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::computeEntropyGrowthDerivative(const T f[], const T fNeq[], const T &alpha)
393{
394 typedef DESCRIPTOR L;
395
396 T entropyGrowthDerivative = T();
397 for (int iPop = 0; iPop < L::q; ++iPop) {
398 T tmp = f[iPop] - alpha*fNeq[iPop];
399 OLB_ASSERT(tmp > T(), "f[iPop] - alpha*fNeq[iPop] <= 0");
400 entropyGrowthDerivative += fNeq[iPop]*util::log(tmp/descriptors::t<T,L>(iPop));
401 }
402
403 return entropyGrowthDerivative;
404}
405
406template<typename T, typename DESCRIPTOR, typename MOMENTA>
407bool ForcedEntropicDynamics<T,DESCRIPTOR,MOMENTA>::getAlpha(T &alpha, const T f[], const T fNeq[])
408{
409 const T epsilon = std::numeric_limits<T>::epsilon();
410
411 T alphaGuess = T();
412 const T var = 100.0;
413 const T errorMax = epsilon*var;
414 T error = 1.0;
415 int count = 0;
416 for (count = 0; count < 10000; ++count) {
417 T entGrowth = computeEntropyGrowth(f,fNeq,alpha);
418 T entGrowthDerivative = computeEntropyGrowthDerivative(f,fNeq,alpha);
419 if ((error < errorMax) || (util::fabs(entGrowth) < var*epsilon)) {
420 return true;
421 }
422 alphaGuess = alpha - entGrowth /
423 entGrowthDerivative;
424 error = util::fabs(alpha-alphaGuess);
425 alpha = alphaGuess;
426 }
427 return false;
428}
429
430}
431
432#endif
433
olb::Cell
Highest-level interface to Cell data.
Definition cell.h:148
olb::EntropicDynamics
Implementation of the entropic collision step.
Definition entropicDynamics.h:88
olb::EntropicDynamics::computeEquilibrium
T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const override
Compute equilibrium distribution function.
Definition entropicDynamics.hh:164
olb::EntropicDynamics::collide
CellStatistic< T > collide(Cell< T, DESCRIPTOR > &cell) override
Collision step.
Definition entropicDynamics.hh:170
olb::EntropicDynamics::setOmega
void setOmega(T omega_)
Set local relaxation parameter of the dynamics.
Definition entropicDynamics.hh:216
olb::EntropicDynamics::getOmega
T getOmega() const
Get local relaxation parameter of the dynamics.
Definition entropicDynamics.hh:210
olb::EntropicDynamics::EntropicDynamics
EntropicDynamics(T omega_)
Constructor.
Definition entropicDynamics.hh:156
olb::EntropicEqDynamics::EntropicEqDynamics
EntropicEqDynamics(T omega_)
Constructor.
Definition entropicDynamics.hh:46
olb::EntropicEqDynamics::computeEquilibrium
T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const override
Compute equilibrium distribution function.
Definition entropicDynamics.hh:54
olb::EntropicEqDynamics::collide
CellStatistic< T > collide(Cell< T, DESCRIPTOR > &cell) override
Collision step.
Definition entropicDynamics.hh:60
olb::EntropicEqDynamics::setOmega
void setOmega(T omega_)
Set local relaxation parameter of the dynamics.
Definition entropicDynamics.hh:84
olb::EntropicEqDynamics::getOmega
T getOmega() const
Get local relaxation parameter of the dynamics.
Definition entropicDynamics.hh:78
olb::ForcedEntropicDynamics
Implementation of the forced entropic collision step.
Definition entropicDynamics.h:119
olb::ForcedEntropicDynamics::setOmega
virtual void setOmega(T omega_)
Set local relaxation parameter of the dynamics.
Definition entropicDynamics.hh:360
olb::ForcedEntropicDynamics::ForcedEntropicDynamics
ForcedEntropicDynamics(T omega_)
Constructor.
Definition entropicDynamics.hh:293
olb::ForcedEntropicDynamics::getOmega
virtual T getOmega() const
Get local relaxation parameter of the dynamics.
Definition entropicDynamics.hh:354
olb::ForcedEntropicDynamics::collide
virtual CellStatistic< T > collide(Cell< T, DESCRIPTOR > &cell)
Collision step.
Definition entropicDynamics.hh:308
olb::ForcedEntropicDynamics::computeEquilibrium
virtual T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
Compute equilibrium distribution function.
Definition entropicDynamics.hh:301
olb::ForcedEntropicEqDynamics::computeEquilibrium
virtual T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
Compute equilibrium distribution function.
Definition entropicDynamics.hh:105
olb::ForcedEntropicEqDynamics::setOmega
virtual void setOmega(T omega_)
Set local relaxation parameter of the dynamics.
Definition entropicDynamics.hh:143
olb::ForcedEntropicEqDynamics::collide
virtual CellStatistic< T > collide(Cell< T, DESCRIPTOR > &cell)
Collision step.
Definition entropicDynamics.hh:112
olb::ForcedEntropicEqDynamics::getOmega
virtual T getOmega() const
Get local relaxation parameter of the dynamics.
Definition entropicDynamics.hh:137
olb::ForcedEntropicEqDynamics::ForcedEntropicEqDynamics
ForcedEntropicEqDynamics(T omega_)
Constructor.
Definition entropicDynamics.hh:97
entropicDynamics.h
A collection of entropic dynamics classes (e.g.
entropicLbHelpers.h
A collection of dynamics classes (e.g.
olb::util::log
ADf< T, DIM > log(const ADf< T, DIM > &a)
Definition aDiff.h:475
olb::util::fabs
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Definition pack.h:106
olb
Top level namespace for all of OpenLB.
Definition boundaryPostProcessors2D.h:34
OLB_ASSERT
#define OLB_ASSERT(COND, MESSAGE)
Definition olbDebug.h:45
olb::CellStatistic
Return value of any collision.
Definition interface.h:43
olb::Dynamics::getName
virtual std::string getName() const
Return human-readable name.
Definition interface.h:63
olb::entropicLbHelpers::equilibrium
static T equilibrium(int iPop, T rho, const T u[DESCRIPTOR::d])
Computation of equilibrium distribution.
Definition entropicLbHelpers.h:38
olb::lbm::addExternalForce
static void addExternalForce(CELL &cell, const RHO &rho, const U &u, const OMEGA &omega, const FORCE &force) any_platform
Add a force term after BGK collision.
Definition lbm.h:463
Generated on Tue Mar 5 2024 09:57:46 for OpenLB by doxygen 1.10.0