blob: bb55d00420ac3de05e8c4369c661bf02f33ef527 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
|
#pragma once
#include "descriptor.hpp"
namespace kel {
namespace lbm {
template<typename T, typename Descriptor>
saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw::data<sch::FixedArray<T,Descriptor::D>> vel){
using dfi = df_info<T, Descriptor>;
saw::data<sch::FixedArray<T,Descriptor::Q>> eq;
// ^
// 0.0
// / \
// | |
//
// Velocity * Velocity meaning || vel ||_2^2 or <vel,vel>_2
saw::data<T> vel_vel{0.0};
for(uint64_t j = 0u; j < Descriptor::D; ++j){
vel_vel = vel_vel + vel.at(j) * vel.at(j);
}
/**
* Calculate equilibrium
*/
for(uint64_t i = 0u; i < eq.template get_dim_size<0u>(); ++i){
saw::data<T> vel_c{};
for(uint64_t j = 0u; j < Descriptor::D; ++j){
// <vel,c_i>_2
vel_c = vel_c + (vel.at(j) * saw::data<T>{static_cast<saw::native_data_type<T>::type>(dfi::directions[i][j])});
}
auto vel_c_cs2 = vel_c * saw::data<T>{dfi::inv_cs2};
eq.at(i).set(
dfi::weights[i] * rho.get() *
(
1.0
+ vel_c_cs2.get()
- dfi::inv_cs2 * 0.5 * vel_vel.get()
+ vel_c_cs2.get() * vel_c_cs2.get() * 0.5
)
);
}
return eq;
}
}
}
|
