OpenLB 1.7
Loading...
Searching...
No Matches
Classes | Public Types | Public Member Functions | Static Public Attributes | List of all members
olb::MCMPForcedPostProcessor< POTENTIAL, N_COMPONENTS > Struct Template Reference

Multi-Component-Multi-Phase Shan-Chen force with thermodynamic equation of state based on. More...

#include <shanChenForcedPostProcessor.h>

+ Collaboration diagram for olb::MCMPForcedPostProcessor< POTENTIAL, N_COMPONENTS >:

Classes

struct  A
 
struct  ALPHA
 
struct  B
 
struct  CHI
 
struct  DEVI
 
struct  EPSILON
 
struct  G
 
struct  GI
 
struct  GII
 
struct  K
 
struct  MM
 
struct  SIGMA
 
struct  TCRIT
 
struct  TEMPERATURE
 

Public Types

using parameters = meta::list<CHI,TEMPERATURE,G,SIGMA,EPSILON,K,A,B,MM,TCRIT,DEVI,ALPHA,GI,GII>
 

Public Member Functions

int getPriority () const
 
template<typename CELLS , typename PARAMETERS >
void apply (CELLS &cells, PARAMETERS &parameters) any_platform
 

Static Public Attributes

static constexpr OperatorScope scope = OperatorScope::PerCellWithParameters
 

Detailed Description

template<typename POTENTIAL, unsigned N_COMPONENTS>
struct olb::MCMPForcedPostProcessor< POTENTIAL, N_COMPONENTS >

Multi-Component-Multi-Phase Shan-Chen force with thermodynamic equation of state based on.

  1. Czelusniak L E, Mapelli V P, Guzella M S, Cabezas-Gómez L, Wagner A J (2020) Force approach for the pseudopotential lattice Boltzmann method. Physical Review E, 102(3), 033307. DOI: 10.1103/PhysRevE.102.033307
  2. Peng C, Ayala L F, Ayala O M (2021) A thermodynamically consistent pseudo- potential lattice Boltzmann model for multi-component, multiphase, partially miscible mixtures. Journal of Computational Physics, 429, 110018. DOI: 10.1016/j.jcp.2020.110018

Definition at line 152 of file shanChenForcedPostProcessor.h.

Member Typedef Documentation

◆ parameters

template<typename POTENTIAL , unsigned N_COMPONENTS>
using olb::MCMPForcedPostProcessor< POTENTIAL, N_COMPONENTS >::parameters = meta::list<CHI,TEMPERATURE,G,SIGMA,EPSILON,K,A,B,MM,TCRIT,DEVI,ALPHA,GI,GII>

Definition at line 170 of file shanChenForcedPostProcessor.h.

Member Function Documentation

◆ apply()

template<typename POTENTIAL , unsigned N_COMPONENTS>
template<typename CELLS , typename PARAMETERS >
void olb::MCMPForcedPostProcessor< POTENTIAL, N_COMPONENTS >::apply ( CELLS & cells,
PARAMETERS & parameters )
inline

Definition at line 177 of file shanChenForcedPostProcessor.h.

178 {
179 using V = typename CELLS::template value_t<names::Component1>::value_t;
180 using DESCRIPTOR = typename CELLS::template value_t<names::Component1>::descriptor_t;
181
182 auto chi = parameters.template get<CHI>();
183 auto t = parameters.template get<TEMPERATURE>();
184 auto g = parameters.template get<G>();
185 auto sig = parameters.template get<SIGMA>();
186 auto eps = parameters.template get<EPSILON>();
187 auto k = parameters.template get<K>();
188 auto a = parameters.template get<A>();
189 auto b = parameters.template get<B>();
190 auto M = parameters.template get<MM>();
191 auto T_c = parameters.template get<TCRIT>();
192 auto m = parameters.template get<DEVI>();
193 auto alpha = parameters.template get<ALPHA>();
194 auto g_I = parameters.template get<GI>();
195 auto g_II = parameters.template get<GII>();
196
197 Vector<V,10> rhoField{};
198 V rhoTot = 0;
199 Vector<V,DESCRIPTOR::d> u_sum{};
200 auto& cell1 = cells.template get<names::Component1>();
201 rhoField[0] = cell1.template getFieldComponent<descriptors::STATISTIC>(0);
202 Vector<V,DESCRIPTOR::d> u1{};
203 lbm<DESCRIPTOR>::computeJ(cell1, u1);
204 u_sum += u1;
205
206 auto& cell2 = cells.template get<names::Component2>();
207 rhoField[1] = cell2.template getFieldComponent<descriptors::STATISTIC>(0);
208 Vector<V,DESCRIPTOR::d> u2{};
209 lbm<DESCRIPTOR>::computeJ(cell2, u2);
210 u_sum += u2;
211
212#ifdef THIRD_COMPONENT
213 auto& cell3 = cells.template get<names::Component3>();
214 rhoField[2] = cell3.template getFieldComponent<descriptors::STATISTIC>(0);
215 Vector<V,DESCRIPTOR::d> u3{};
216 lbm<DESCRIPTOR>::computeJ(cell3, u3);
217 u_sum += u3;
218#endif
219#ifdef FOURTH_COMPONENT
220 auto& cell4 = cells.template get<names::Component4>();
221 rhoField[3] = cell4.template getFieldComponent<descriptors::STATISTIC>(0);
222 Vector<V,DESCRIPTOR::d> u4{};
223 lbm<DESCRIPTOR>::computeJ(cell4, u4);
224 u_sum += u4;
225#endif
226 for (unsigned n = 0; n < N_COMPONENTS; ++n) {
227 rhoTot += rhoField[n];
228 }
229
230 // Computation of the interaction potential
231 Vector<V,DESCRIPTOR::d> totalForce{};
232 Vector<V,DESCRIPTOR::d> A_ij{};
233 //V psi = POTENTIAL().compute(rhoField, t, k, a, b, T_c, m, alpha, g_I, g_II, M);
234 V p = POTENTIAL().computeP(rhoField, t, a, b, T_c, m, alpha, g_I, g_II, M);
235 V psi = util::sqrt(util::abs(6.*(k*p - rhoTot/3.)));
236 //V psi = util::sqrt(6.*p);
237
238 for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
239 Vector<V,10> rhoField_i{};
240 auto nextCell1 = cell1.neighbor(descriptors::c<DESCRIPTOR>(iPop));
241 rhoField_i[0] = nextCell1.template getFieldComponent<descriptors::STATISTIC>(0);
242 auto nextCell2 = cell2.neighbor(descriptors::c<DESCRIPTOR>(iPop));
243 rhoField_i[1] = nextCell2.template getFieldComponent<descriptors::STATISTIC>(0);
244#ifdef THIRD_COMPONENT
245 auto nextCell3 = cell3.neighbor(descriptors::c<DESCRIPTOR>(iPop));
246 rhoField_i[2] = nextCell3.template getFieldComponent<descriptors::STATISTIC>(0);
247#endif
248#ifdef FOURTH_COMPONENT
249 auto nextCell4 = cell4.neighbor(descriptors::c<DESCRIPTOR>(iPop));
250 rhoField_i[3] = nextCell4.template getFieldComponent<descriptors::STATISTIC>(0);
251#endif
252 V rhoTot_i = 0;
253 for (unsigned n = 0; n < N_COMPONENTS; ++n) {
254 rhoTot_i += rhoField_i[n];
255 }
256 V p_i = POTENTIAL().computeP(rhoField_i, t, a, b, T_c, m, alpha, g_I, g_II, M);
257 //V psi_i = POTENTIAL().compute(rhoField_i, t, k, a, b, T_c, m, alpha, g_I, g_II, M);
258 V phi_i = 6.*(k*p_i - rhoTot_i/3.)/g;
259 //V phi_i = 6.*p_i;
260 V psi_i = 0;
261 if (phi_i >= 0){
262 psi_i = util::sqrt(util::abs(phi_i));
263 }
264 else {
265 psi_i = -1*util::sqrt(util::abs(phi_i));
266 }
267 totalForce += -g * psi * psi_i * descriptors::c<DESCRIPTOR>(iPop) * descriptors::t<V,DESCRIPTOR>(iPop);
268 for (int jPop = 0; jPop < DESCRIPTOR::q; ++jPop) {
269 Vector<V,10> rhoField_j{};
270 auto nextCell1 = cell1.neighbor(descriptors::c<DESCRIPTOR>(jPop));
271 rhoField_j[0] = nextCell1.template getFieldComponent<descriptors::STATISTIC>(0);
272 auto nextCell2 = cell2.neighbor(descriptors::c<DESCRIPTOR>(jPop));
273 rhoField_j[1] = nextCell2.template getFieldComponent<descriptors::STATISTIC>(0);
274#ifdef THIRD_COMPONENT
275 auto nextCell3 = cell3.neighbor(descriptors::c<DESCRIPTOR>(jPop));
276 rhoField_j[2] = nextCell3.template getFieldComponent<descriptors::STATISTIC>(0);
277#endif
278#ifdef FOURTH_COMPONENT
279 auto nextCell4 = cell4.neighbor(descriptors::c<DESCRIPTOR>(jPop));
280 rhoField_j[3] = nextCell4.template getFieldComponent<descriptors::STATISTIC>(0);
281#endif
282 V rhoTot_j = 0;
283 for (unsigned n = 0; n < N_COMPONENTS; ++n) {
284 rhoTot_j += rhoField_j[n];
285 }
286 V p_j = POTENTIAL().computeP(rhoField_j, t, a, b, T_c, m, alpha, g_I, g_II, M);
287 //V psi_j = POTENTIAL().compute(rhoField_j, t, k, a, b, T_c, m, alpha, g_I, g_II, M);
288 V phi_j = 6.*(k*p_j - (rhoTot_j)/3.)/g;
289 //V phi_j = 6.*p_j;
290 V psi_j = 0;
291 if (phi_j >= 0){
292 psi_j = util::sqrt(util::abs(phi_j));
293 }
294 else {
295 psi_j = -1*util::sqrt(util::abs(phi_j));
296 }
297 Vector<V,DESCRIPTOR::d> c1{};
298 Vector<V,DESCRIPTOR::d> c2{};
299 for (int m = 0; m < DESCRIPTOR::d; ++m) {
300 Vector<V,DESCRIPTOR::d> C1{};
301 V C2 = 0.;
302 for (int n = 0; n < DESCRIPTOR::d; ++n) {
303 C1[n] += descriptors::c<DESCRIPTOR>(jPop)[m] * descriptors::c<DESCRIPTOR>(jPop)[n];
304 C2 += descriptors::c<DESCRIPTOR>(jPop)[n] * descriptors::c<DESCRIPTOR>(jPop)[n] - 1./descriptors::invCs2<V,DESCRIPTOR>();
305 if (n==m){
306 C1[n] -= 1./descriptors::invCs2<V,DESCRIPTOR>();
307 }
308 }
309 for (int n = 0; n < DESCRIPTOR::d; ++n) {
310 c1[m] += descriptors::c<DESCRIPTOR>(iPop)[n] * C1[n];
311 }
312 c2[m] += descriptors::c<DESCRIPTOR>(iPop)[m] * C2;
313 }
314 A_ij = 0.5*g*descriptors::t<V,DESCRIPTOR>(iPop)*descriptors::t<V,DESCRIPTOR>(jPop)*descriptors::invCs2<V,DESCRIPTOR>()*
315 ((3./2.*eps-(sig-1))*c1 + (sig-1)*c2);
316 totalForce += A_ij * psi_i * psi_j;
317 }
318
319 }
320 // Computation of the common velocity, shared among the two populations, consisting
321 // of u and the momentum difference due to interaction potential plus external force
322 auto externalForce1 = cell1.template getField<descriptors::EXTERNAL_FORCE>();
323 Vector<V,DESCRIPTOR::d> uTot = (u_sum + totalForce/2.) / rhoTot + externalForce1/2.;
324 cell1.template setField<descriptors::FORCE>(externalForce1 + chi[0]*totalForce/rhoField[0]);
325 cell1.template setField<descriptors::VELOCITY>(uTot);
326
327 auto externalForce2 = cell2.template getField<descriptors::EXTERNAL_FORCE>();
328 cell2.template setField<descriptors::FORCE>(externalForce2 + chi[1]*totalForce/rhoField[1]);
329 cell2.template setField<descriptors::VELOCITY>(uTot);
330
331#ifdef THIRD_COMPONENT
332 auto externalForce3 = cell3.template getField<descriptors::EXTERNAL_FORCE>();
333 cell3.template setField<descriptors::FORCE>(externalForce3 + chi[2]*totalForce/rhoField[2]);
334 cell3.template setField<descriptors::VELOCITY>(uTot);
335#endif
336#ifdef FOURTH_COMPONENT
337 auto externalForce4 = cell4.template getField<descriptors::EXTERNAL_FORCE>();
338 cell4.template setField<descriptors::FORCE>(externalForce4 + chi[3]*totalForce/rhoField[3]);
339 cell4.template setField<descriptors::VELOCITY>(uTot);
340#endif
341
342 cell1.template setField<descriptors::SCALAR>(p);
343 cell2.template setField<descriptors::SCALAR>(psi);
344 }
olb::descriptors::data::t
platform_constant Fraction t[Q]
Definition descriptorFunction.h:63
olb::descriptors::mrt_data::M
platform_constant Fraction M[Q][Q]
Definition mrtLatticeDescriptors.h:65
olb::descriptors::m
constexpr T m(unsigned iPop, unsigned jPop, tag::MRT)
Definition mrtLatticeDescriptors.h:394
olb::util::sqrt
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
olb::util::abs
ADf< T, DIM > abs(const ADf< T, DIM > &a)
Definition aDiff.h:1019
olb::MCMPForcedPostProcessor::parameters
meta::list< CHI, TEMPERATURE, G, SIGMA, EPSILON, K, A, B, MM, TCRIT, DEVI, ALPHA, GI, GII > parameters
Definition shanChenForcedPostProcessor.h:170
olb::lbm::computeJ
static void computeJ(CELL &cell, J &j) any_platform
Computation of momentum.
Definition lbm.h:197

References olb::util::abs(), olb::lbm< DESCRIPTOR >::computeJ(), and olb::util::sqrt().

+ Here is the call graph for this function:

◆ getPriority()

template<typename POTENTIAL , unsigned N_COMPONENTS>
int olb::MCMPForcedPostProcessor< POTENTIAL, N_COMPONENTS >::getPriority ( ) const
inline

Definition at line 172 of file shanChenForcedPostProcessor.h.

172 {
173 return 0;
174 }

Member Data Documentation

◆ scope

template<typename POTENTIAL , unsigned N_COMPONENTS>
constexpr OperatorScope olb::MCMPForcedPostProcessor< POTENTIAL, N_COMPONENTS >::scope = OperatorScope::PerCellWithParameters
staticconstexpr

Definition at line 153 of file shanChenForcedPostProcessor.h.

The documentation for this struct was generated from the following file:
Generated on Tue Mar 5 2024 09:58:02 for OpenLB by doxygen 1.10.0