OpenLB 1.7
Loading...
Searching...
No Matches
powerLawBGKdynamics.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2012, 2015 Mathias J. Krause, Vojtech Cvrcek, Davide Dapelo
4 * 2022 Nando Suntoyo, 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
31#ifndef POWER_LAW_BGK_DYNAMICS_H
32#define POWER_LAW_BGK_DYNAMICS_H
33
34#include "dynamics/dynamics.h"
35#include "core/cell.h"
36
37#include "collisionLES.h"
38#include "porousBGKdynamics.h"
39
40namespace olb {
41
42// *INDENT-OFF*
43
44namespace powerlaw {
45
46struct OMEGA_MIN : public descriptors::FIELD_BASE<1> { };
47struct OMEGA_MAX : public descriptors::FIELD_BASE<1> { };
48struct M : public descriptors::FIELD_BASE<1> { };
49struct N : public descriptors::FIELD_BASE<1> { };
50// The following is used for Herschel-Bulkley only
53
55template <typename COLLISION, bool HERSCHELBULKLEY=false>
57 using parameters = typename COLLISION::parameters::template include<
59 >;
60
61 static std::string getName()
62 {
63 return "powerlaw::OmegaFromCell<" + COLLISION::getName() + ">";
64 }
65
66 template <typename CELL, typename PARAMETERS, typename OMEGA, typename RHO, typename PI, typename V=typename CELL::value_t>
67 static V computeOmega(CELL& cell, PARAMETERS& parameters, OMEGA& omega0, RHO& rho, PI& pi) any_platform
68 {
69 using DESCRIPTOR = typename CELL::descriptor_t;
70 V pre2 = V{0.5} * descriptors::invCs2<V,DESCRIPTOR>() * omega0 / rho; // strain rate tensor prefactor
71 pre2 *= pre2;
72 V gamma{};
73 if constexpr (DESCRIPTOR::template provides<descriptors::SHEAR_RATE_MAGNITUDE>()) {
74 gamma = cell.template getField<descriptors::SHEAR_RATE_MAGNITUDE>();
75 }
76 else {
77 if constexpr (DESCRIPTOR::template provides<descriptors::FORCE>()) {
78 // Cannot be done in just one line, it gives error - I don't know why. Davide Dapelo
79 const auto force = cell.template getField<descriptors::FORCE>();
80 gamma = util::sqrt(V{2}*pre2*lbm<DESCRIPTOR>::computePiNeqNormSqr(cell, force));
81 }
82 else {
83 gamma = util::sqrt(V{2}*pre2*lbm<DESCRIPTOR>::computePiNeqNormSqr(cell));
84 }
85 }
86 if constexpr(HERSCHELBULKLEY) {
87 gamma = util::max(parameters.template get<SHEAR_RATE_MIN>(), gamma);
88 }
89 V m = parameters.template get<M>();
90 V n = parameters.template get<N>();
91 V nuNew = m * util::pow(gamma, n-V{1}); // Ostwald-de Waele relation
92 if constexpr(HERSCHELBULKLEY) {
93 // Second term necessary for Herschel-Bulkley relation
94 nuNew += parameters.template get<YIELD_STRESS>() / gamma;
95 }
96 V newOmega = V{1} / (nuNew*descriptors::invCs2<V,DESCRIPTOR>() + V{0.5});
97 V omegaMax = parameters.template get<OMEGA_MAX>();
98 newOmega = util::min(newOmega, omegaMax);
99 V omegaMin = parameters.template get<OMEGA_MIN>();
100 newOmega = util::max(newOmega, omegaMin);
101 return newOmega;
102 }
103
104 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
105 struct type {
106 using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
107 using CollisionO = typename COLLISION::template type<DESCRIPTOR, MOMENTA, EQUILIBRIUM>;
108
109 template <typename CELL, typename PARAMETERS, typename V=typename CELL::value_t>
111 {
112 V rho, u[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR>::n] { };
113 MomentaF().computeAllMomenta(cell, rho, u, pi);
114 const V oldOmega = cell.template getField<descriptors::OMEGA>();
115 const V newOmega = computeOmega(cell, parameters, oldOmega, rho, pi);
116 cell.template setField<descriptors::OMEGA>(newOmega);
117 parameters.template set<descriptors::OMEGA>(newOmega);
118 return CollisionO().apply(cell, parameters);
119 }
120 };
121};
122
123
124template <int... NORMAL>
126
128template <int... NORMAL>
130 static std::string getName()
131 {
132 return "PeriodicPressureOffset";
133 }
134
135 template <typename DESCRIPTOR, typename MOMENTA>
136 using combined_momenta = typename MOMENTA::template type<DESCRIPTOR>;
137
138 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
139 using combined_equilibrium = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
140
141 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM, typename COLLISION>
143 using CollisionO = typename COLLISION::template type<DESCRIPTOR, MOMENTA, EQUILIBRIUM>;
144
145 template <typename CELL, typename PARAMETERS, typename V=typename CELL::value_t>
146 CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform
147 {
148 static constexpr auto populations = util::populationsContributingToDirection<DESCRIPTOR, NORMAL...>();
149
150 auto statistic = CollisionO().apply(cell, parameters);
151 const V densityOffset = parameters.template get<PRESSURE_OFFSET<NORMAL...>>();
152 for (unsigned iPop : populations) {
153 cell[iPop] += (cell[iPop] + descriptors::t<V,DESCRIPTOR>(iPop)) * densityOffset;
154 }
155 return statistic;
156 }
157 };
158
159 template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM, typename COLLISION>
160 using combined_parameters = typename COLLISION::parameters
161 ::template include<PRESSURE_OFFSET<NORMAL...>>;
162};
163
164}
165
167template <typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
169 T, DESCRIPTOR,
170 MOMENTA,
173>;
174
176template <typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
178 T, DESCRIPTOR,
179 MOMENTA,
183>;
184
186template <typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
188 T, DESCRIPTOR,
189 MOMENTA,
192>;
193
195template <typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
197 T, DESCRIPTOR,
198 MOMENTA,
202>;
203
205template<typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
207 T, DESCRIPTOR,
208 MOMENTA,
211>;
212
214template<typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
216 T, DESCRIPTOR,
217 MOMENTA,
221>;
222
224template<typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
226 T, DESCRIPTOR,
227 MOMENTA,
230>;
231
232// *INDENT-ON*
233
234}
235
236#endif
Definition of a LB cell – header file.
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
constexpr auto populationsContributingToDirection() any_platform
Return array of population indices where c[iPop][iD] == NORMAL[iD].
Definition util.h:226
Top level namespace for all of OpenLB.
#define any_platform
Define preprocessor macros for device-side functions, constant storage.
Definition platform.h:78
Return value of any collision.
Definition interface.h:43
Base of a field whose size is defined by [C,U_1,...,U_N]^T * [1,V_1,...V_N].
Dynamics constructed as a tuple of momenta, equilibrium and collision.
Definition interface.h:182
Dynamics combination rule implementing the forcing scheme by Guo et al.
Definition forcing.h:38
static V computePiNeqNormSqr(CELL &cell, const FORCE &force) any_platform
Computes squared norm of non-equilibrium part of 2nd momentum for forced dynamics.
Definition lbm.h:424
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
typename COLLISION::template type< DESCRIPTOR, MOMENTA, EQUILIBRIUM > CollisionO
typename MOMENTA::template type< DESCRIPTOR > MomentaF
Compute and update cell-wise OMEGA using Oswald-de-waele model.
static V computeOmega(CELL &cell, PARAMETERS &parameters, OMEGA &omega0, RHO &rho, PI &pi) any_platform
typename COLLISION::parameters::template include< descriptors::OMEGA, OMEGA_MIN, OMEGA_MAX, M, N, YIELD_STRESS, SHEAR_RATE_MIN > parameters
typename COLLISION::template type< DESCRIPTOR, MOMENTA, EQUILIBRIUM > CollisionO
CellStatistic< V > apply(CELL &cell, PARAMETERS &parameters) any_platform
Combination rule to realize a pressure drop at a periodic boundary.
typename MOMENTA::template type< DESCRIPTOR > combined_momenta
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > combined_equilibrium
typename COLLISION::parameters ::template include< PRESSURE_OFFSET< NORMAL... > > combined_parameters
Compute number of elements of a symmetric d-dimensional tensor.
Definition util.h:210