-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmodel_host_device.h
executable file
·132 lines (127 loc) · 5.08 KB
/
model_host_device.h
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#ifndef MODEL_HOST_DEVICE_H
#define MODEL_HOST_DEVICE_H
#ifndef OPS_FUN_PREFIX
#define OPS_FUN_PREFIX
#endif
enum VariableTypes {
Variable_Rho = 0,
Variable_U = 1,
Variable_V = 2,
Variable_W = 3,
Variable_T = 4,
Variable_Qx = 5,
Variable_Qy = 6,
Variable_Qz = 7,
// This is for the force correction needed by using He1998 scheme.
Variable_U_Force = 8,
Variable_V_Force = 9,
Variable_W_Force = 10,
};
/*
* Calculate the first-order force term
* Author: Jianping Meng 22-Feb-2019
*/
static inline OPS_FUN_PREFIX Real CalcBodyForce(const int xiIndex, const Real rho,
const Real* acceleration) {
Real cf{0};
for (int i = 0; i < LATTDIM; i++) {
cf += CS * XI[xiIndex * LATTDIM + i] * acceleration[i];
}
return WEIGHTS[xiIndex] * rho * cf;
}
static inline OPS_FUN_PREFIX Real CalcBGKFeq(const int l, const Real rho, const Real u, const Real v,
const Real T, const int polyOrder) {
Real cu{(CS * XI[l * LATTDIM] * u + CS * XI[l * LATTDIM + 1] * v)};
Real c2{(CS * XI[l * LATTDIM] * CS * XI[l * LATTDIM] +
CS * XI[l * LATTDIM + 1] * CS * XI[l * LATTDIM + 1])};
Real cu2{cu * cu};
Real u2{u * u + v * v};
Real res = 1.0 + cu + 0.5 * (cu2 - u2 + (T - 1.0) * (c2 - LATTDIM));
if ((polyOrder) >= 3) {
res = res +
cu * (cu2 - 3.0 * u2 + 3.0 * (T - 1.0) * (c2 - LATTDIM - 2.0)) /
6.0;
}
if ((polyOrder) >= 4) {
res =
res + (cu2 * cu2 - 6.0 * cu2 * u2 + 3.0 * u2 * u2) / 24.0 +
(T - 1.0) * ((c2 - (LATTDIM + 2)) * (cu2 - u2) - 2.0 * cu2) / 4.0 +
(T - 1.0) * (T - 1.0) *
(c2 * c2 - 2.0 * (LATTDIM + 2) * c2 + LATTDIM * (LATTDIM + 2)) /
8.0;
}
return WEIGHTS[l] * rho * res;
}
static inline OPS_FUN_PREFIX Real CalcBGKFeqpseudo3D(const int l, const Real rho, const Real u, const Real v,
const Real w, const Real* g, const Real tau) {
Real cu{(XI[l * LATTDIM] * u + XI[l * LATTDIM + 1] * v +
XI[l * LATTDIM + 2] * w)};
Real cF{(XI[l * LATTDIM] * g[0] + XI[l * LATTDIM + 1] * g[1] +
XI[l * LATTDIM + 2] * g[2])};
Real uF{(u * g[0] + v * g[1] +
w * g[2])};
Real B = 1.0-0.5/tau;
Real C = 1.0-0.5/tau;
Real res = B*cF*CS*CS + C* (cF*cu*CS*CS*CS*CS - uF*CS*CS);
return WEIGHTS[l] * res;
}
static inline OPS_FUN_PREFIX Real CalcBGKBodyForce(const int l, const Real rho, const Real u, const Real v,
const Real w, const Real T, const int polyOrder) {
Real cu{(CS * XI[l * LATTDIM] * u * CS + CS * XI[l * LATTDIM + 1] * v * CS +
CS * XI[l * LATTDIM + 2] * w * CS)};
Real cu2{cu * cu};
Real u2{CS * u * u * CS + CS * v * v * CS + CS * w * w * CS};
Real res = 1.0 + cu + 0.5 * (cu2 - u2);
return WEIGHTS[l] * rho * res;
}
static inline OPS_FUN_PREFIX Real CalcBGKFeq(const int l, const Real rho, const Real u, const Real v,
const Real w, const Real T, const int polyOrder) {
Real cu{(CS * XI[l * LATTDIM] * u * CS + CS * XI[l * LATTDIM + 1] * v * CS +
CS * XI[l * LATTDIM + 2] * w * CS)};
Real c2{(CS * XI[l * LATTDIM] * CS * XI[l * LATTDIM] +
CS * XI[l * LATTDIM + 1] * CS * XI[l * LATTDIM + 1] +
CS * XI[l * LATTDIM + 2] * CS * XI[l * LATTDIM + 2])};
Real cu2{cu * cu};
Real u2{CS * u * u * CS + CS * v * v * CS + CS * w * w * CS};
Real res = 1.0 + cu + 0.5 * (cu2 - u2);
if ((polyOrder) >= 3) {
res = res +
cu * (cu2 - 3.0 * u2 + 3.0 * (T - 1.0) * (c2 - LATTDIM - 2.0)) /
6.0;
}
if ((polyOrder) >= 4) {
res =
res + (cu2 * cu2 - 6.0 * cu2 * u2 + 3.0 * u2 * u2) / 24.0 +
(T - 1.0) * ((c2 - (LATTDIM + 2)) * (cu2 - u2) - 2.0 * cu2) / 4.0 +
(T - 1.0) * (T - 1.0) *
(c2 * c2 - 2.0 * (LATTDIM + 2) * c2 + LATTDIM * (LATTDIM + 2)) /
8.0;
}
return WEIGHTS[l] * rho * res;
}
static inline OPS_FUN_PREFIX Real CalcSWEFeq(const int l, const Real h, const Real u, const Real v,
const int polyOrder) {
// Implementing the model derived in Please refer to Meng, Gu Emerson, Peng
// and Zhang, IJMPC 2018(29):1850080
Real cu{(CS * XI[l * LATTDIM] * u + CS * XI[l * LATTDIM + 1] * v)};
Real c2{(CS * XI[l * LATTDIM] * CS * XI[l * LATTDIM] +
CS * XI[l * LATTDIM + 1] * CS * XI[l * LATTDIM + 1])};
Real cu2{cu * cu};
Real u2{u * u + v * v};
Real res = 1.0 + cu + 0.5 * (cu2 - u2 + (h - 1.0) * (c2 - LATTDIM));
if (polyOrder >= 3) {
res = res +
cu * (cu2 - 3.0 * u2 + 3.0 * (h - 1.0) * (c2 - LATTDIM - 2.0)) /
6.0;
}
if (polyOrder >= 4) {
res =
res + (cu2 * cu2 - 6.0 * cu2 * u2 + 3.0 * u2 * u2) / 24.0 +
(h - 1.0) * ((c2 - (LATTDIM + 2)) * (cu2 - u2) - 2.0 * cu2) / 4.0 +
(h - 1.0) * (h - 1.0) *
(c2 * c2 - 2.0 * (LATTDIM + 2) * c2 + LATTDIM * (LATTDIM + 2)) /
8.0;
}
return WEIGHTS[l] * h * res;
}
#endif //MODEL_HOST_DEVICE_H