Neko 1.99.9
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
dynamic_smagorinsky_nut_kernel.h
Go to the documentation of this file.
1#ifndef __COMMON_DYNAMIC_SMAGORINSKY_NUT_KERNEL_H__
2#define __COMMON_DYNAMIC_SMAGORINSKY_NUT_KERNEL_H__
3/*
4 Copyright (c) 2024, The Neko Authors
5 All rights reserved.
6
7 Redistribution and use in source and binary forms, with or without
8 modification, are permitted provided that the following conditions
9 are met:
10
11 * Redistributions of source code must retain the above copyright
12 notice, this list of conditions and the following disclaimer.
13
14 * Redistributions in binary form must reproduce the above
15 copyright notice, this list of conditions and the following
16 disclaimer in the documentation and/or other materials provided
17 with the distribution.
18
19 * Neither the name of the authors nor the names of its
20 contributors may be used to endorse or promote products derived
21 from this software without specific prior written permission.
22
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
26 FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
27 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
28 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
29 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
31 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
32 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
33 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
34 POSSIBILITY OF SUCH DAMAGE.
35*/
36
40#include <cmath>
41#include <algorithm>
42template< typename T>
43__global__ void s_abs_compute(T * __restrict__ s_abs,
44 const T * __restrict__ s11,
45 const T * __restrict__ s22,
46 const T * __restrict__ s33,
47 const T * __restrict__ s12,
48 const T * __restrict__ s13,
49 const T * __restrict__ s23,
50 const int n){
51
52 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
53 const int str = blockDim.x * gridDim.x;
54
55 for (int i = idx; i < n; i += str) {
56 const T s11_r = s11[i];
57 const T s22_r = s22[i];
58 const T s33_r = s33[i];
59 const T s12_r = s12[i];
60 const T s13_r = s13[i];
61 const T s23_r = s23[i];
62
63 s_abs[i] = sqrt(2.0 * (s11_r * s11_r +
64 s22_r * s22_r +
65 s33_r * s33_r) +
66 4.0 * (s12_r * s12_r +
67 s13_r * s13_r +
68 s23_r * s23_r));
69
70 }
71}
72
73template< typename T>
74__global__ void lij_compute_part1(T * __restrict__ l11,
75 T * __restrict__ l22,
76 T * __restrict__ l33,
77 T * __restrict__ l12,
78 T * __restrict__ l13,
79 T * __restrict__ l23,
80 const T * __restrict__ u,
81 const T * __restrict__ v,
82 const T * __restrict__ w,
83 const T * __restrict__ fu,
84 const T * __restrict__ fv,
85 const T * __restrict__ fw,
86 T * __restrict__ fuu,
87 T * __restrict__ fvv,
88 T * __restrict__ fww,
89 T * __restrict__ fuv,
90 T * __restrict__ fuw,
91 T * __restrict__ fvw,
92 const int n){
93
94 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
95 const int str = blockDim.x * gridDim.x;
96
97 for (int i = idx; i < n; i += str) {
98 const T u_r = u[i];
99 const T v_r = v[i];
100 const T w_r = w[i];
101 const T fu_r = fu[i];
102 const T fv_r = fv[i];
103 const T fw_r = fw[i];
104
105 l11[i] = fu_r * fu_r;
106 l22[i] = fv_r * fv_r;
107 l33[i] = fw_r * fw_r;
108 l12[i] = fu_r * fv_r;
109 l13[i] = fu_r * fw_r;
110 l23[i] = fv_r * fw_r;
111
112 fuu[i] = u_r * u_r;
113 fvv[i] = v_r * v_r;
114 fww[i] = w_r * w_r;
115 fuv[i] = u_r * v_r;
116 fuw[i] = u_r * w_r;
117 fvw[i] = v_r * w_r;
118 }
119}
120
121template< typename T>
122__global__ void lij_compute_part2(T * __restrict__ l11,
123 T * __restrict__ l22,
124 T * __restrict__ l33,
125 T * __restrict__ l12,
126 T * __restrict__ l13,
127 T * __restrict__ l23,
128 const T * __restrict__ fuu,
129 const T * __restrict__ fvv,
130 const T * __restrict__ fww,
131 const T * __restrict__ fuv,
132 const T * __restrict__ fuw,
133 const T * __restrict__ fvw,
134 const int n){
135
136 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
137 const int str = blockDim.x * gridDim.x;
138
139 for (int i = idx; i < n; i += str) {
140 l11[i] -= fuu[i];
141 l22[i] -= fvv[i];
142 l33[i] -= fww[i];
143 l12[i] -= fuv[i];
144 l13[i] -= fuw[i];
145 l23[i] -= fvw[i];
146 }
147}
148
149template< typename T>
150__global__ void mij_compute_part1(T * __restrict__ m11,
151 T * __restrict__ m22,
152 T * __restrict__ m33,
153 T * __restrict__ m12,
154 T * __restrict__ m13,
155 T * __restrict__ m23,
156 const T * __restrict__ s_abs,
157 const T * __restrict__ s11,
158 const T * __restrict__ s22,
159 const T * __restrict__ s33,
160 const T * __restrict__ s12,
161 const T * __restrict__ s13,
162 const T * __restrict__ s23,
163 const T * __restrict__ delta,
164 const T * __restrict__ fs11,
165 const T * __restrict__ fs22,
166 const T * __restrict__ fs33,
167 const T * __restrict__ fs12,
168 const T * __restrict__ fs13,
169 const T * __restrict__ fs23,
170 T * __restrict__ fsabss11,
171 T * __restrict__ fsabss22,
172 T * __restrict__ fsabss33,
173 T * __restrict__ fsabss12,
174 T * __restrict__ fsabss13,
175 T * __restrict__ fsabss23,
176 const T delta_ratio2,
177 const int n){
178
179 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
180 const int str = blockDim.x * gridDim.x;
181
182 for (int i = idx; i < n; i += str) {
183 const T s_abs_r = s_abs[i];
184 const T delta_r = delta[i];
185 const T delta2 = delta_r * delta_r;
186 const T fs_abs_r = sqrt(2.0 *
187 (fs11[i] * fs11[i] + fs22[i] * fs22[i] + fs33[i] * fs33[i]) +
188 4.0 * (fs12[i] * fs12[i] + fs13[i] * fs13[i] +
189 fs23[i] * fs23[i]));
190
191 m11[i] = delta_ratio2 * delta2 * fs_abs_r * fs11[i];
192 m22[i] = delta_ratio2 * delta2 * fs_abs_r * fs22[i];
193 m33[i] = delta_ratio2 * delta2 * fs_abs_r * fs33[i];
194 m12[i] = delta_ratio2 * delta2 * fs_abs_r * fs12[i];
195 m13[i] = delta_ratio2 * delta2 * fs_abs_r * fs13[i];
196 m23[i] = delta_ratio2 * delta2 * fs_abs_r * fs23[i];
197
198 fsabss11[i] = delta2 * s_abs_r * s11[i];
199 fsabss22[i] = delta2 * s_abs_r * s22[i];
200 fsabss33[i] = delta2 * s_abs_r * s33[i];
201 fsabss12[i] = delta2 * s_abs_r * s12[i];
202 fsabss13[i] = delta2 * s_abs_r * s13[i];
203 fsabss23[i] = delta2 * s_abs_r * s23[i];
204 }
205}
206
207template< typename T>
208__global__ void mij_nut_compute_part2(const T * __restrict__ m11,
209 const T * __restrict__ m22,
210 const T * __restrict__ m33,
211 const T * __restrict__ m12,
212 const T * __restrict__ m13,
213 const T * __restrict__ m23,
214 T * __restrict__ l11,
215 T * __restrict__ l22,
216 T * __restrict__ l33,
217 T * __restrict__ l12,
218 T * __restrict__ l13,
219 T * __restrict__ l23,
220 const T * __restrict__ fsabss11,
221 const T * __restrict__ fsabss22,
222 const T * __restrict__ fsabss33,
223 const T * __restrict__ fsabss12,
224 const T * __restrict__ fsabss13,
225 const T * __restrict__ fsabss23,
226 T * __restrict__ num,
227 T * __restrict__ den,
228 T * __restrict__ c_dyn,
229 const T * __restrict__ delta,
230 const T * __restrict__ s_abs,
231 T * __restrict__ nut,
232 const T alpha,
233 const int n){
234
235 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
236 const int str = blockDim.x * gridDim.x;
237
238 T num_curr, den_curr;
239
240 for (int i = idx; i < n; i += str) {
241 const T delta_r = delta[i];
242 const T delta2 = delta_r * delta_r;
243 T c_dyn_r;
244
245 const T m11_r = m11[i] - fsabss11[i];
246 const T m22_r = m22[i] - fsabss22[i];
247 const T m33_r = m33[i] - fsabss33[i];
248 const T m12_r = m12[i] - fsabss12[i];
249 const T m13_r = m13[i] - fsabss13[i];
250 const T m23_r = m23[i] - fsabss23[i];
251
252 num_curr = m11_r * l11[i]
253 + m22_r * l22[i]
254 + m33_r * l33[i]
255 + 2.0 * (m12_r * l12[i]
256 + m13_r * l13[i]
257 + m23_r * l23[i]);
258
259 den_curr = m11_r * m11_r
260 + m22_r * m22_r
261 + m33_r * m33_r
262 + 2.0 * (m12_r * m12_r
263 + m13_r * m13_r
264 + m23_r * m23_r);
265
266 const T num_r = alpha * num[i] + (1.0 - alpha) * num_curr;
267 const T den_r = alpha * den[i] + (1.0 - alpha) * den_curr;
268
269 num[i] = num_r;
270 den[i] = den_r;
271
272 if (den_r > 0.0) {
273 c_dyn_r = 0.5 * num_r / den_r;
274 } else {
275 c_dyn_r = 0.0;
276 }
277
278 c_dyn_r = max(c_dyn_r, 0.0);
279 c_dyn[i] = c_dyn_r;
280 nut[i] = c_dyn_r * delta2 * s_abs[i];
281
282 }
283}
284#endif // __COMMON_DYNAMIC_SMAGORINSKY_NUT_KERNEL_H__
ale_add_kinematics_kernel
__global__ void ale_add_kinematics_kernel(const int n, T *__restrict__ wx, T *__restrict__ wy, T *__restrict__ wz, const T *__restrict__ x_ref, const T *__restrict__ y_ref, const T *__restrict__ z_ref, const T *__restrict__ phi, const T *__restrict__ x, const T *__restrict__ y, const T *__restrict__ z, const kinematics_params_t kin_params)
Definition ale_kinematics_kernel.h:52
w
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ w
Definition ax_helm_compr_kernel.h:978
i
const int i
Definition ax_helm_compr_kernel.h:246
u
__global__ void const T *__restrict__ u
Definition ax_helm_compr_kernel.h:199
v
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ v
Definition ax_helm_compr_kernel.h:977
mij_compute_part1
__global__ void mij_compute_part1(T *__restrict__ m11, T *__restrict__ m22, T *__restrict__ m33, T *__restrict__ m12, T *__restrict__ m13, T *__restrict__ m23, const T *__restrict__ s_abs, const T *__restrict__ s11, const T *__restrict__ s22, const T *__restrict__ s33, const T *__restrict__ s12, const T *__restrict__ s13, const T *__restrict__ s23, const T *__restrict__ delta, const T *__restrict__ fs11, const T *__restrict__ fs22, const T *__restrict__ fs33, const T *__restrict__ fs12, const T *__restrict__ fs13, const T *__restrict__ fs23, T *__restrict__ fsabss11, T *__restrict__ fsabss22, T *__restrict__ fsabss33, T *__restrict__ fsabss12, T *__restrict__ fsabss13, T *__restrict__ fsabss23, const T delta_ratio2, const int n)
Definition dynamic_smagorinsky_nut_kernel.h:150
mij_nut_compute_part2
__global__ void mij_nut_compute_part2(const T *__restrict__ m11, const T *__restrict__ m22, const T *__restrict__ m33, const T *__restrict__ m12, const T *__restrict__ m13, const T *__restrict__ m23, T *__restrict__ l11, T *__restrict__ l22, T *__restrict__ l33, T *__restrict__ l12, T *__restrict__ l13, T *__restrict__ l23, const T *__restrict__ fsabss11, const T *__restrict__ fsabss22, const T *__restrict__ fsabss33, const T *__restrict__ fsabss12, const T *__restrict__ fsabss13, const T *__restrict__ fsabss23, T *__restrict__ num, T *__restrict__ den, T *__restrict__ c_dyn, const T *__restrict__ delta, const T *__restrict__ s_abs, T *__restrict__ nut, const T alpha, const int n)
Definition dynamic_smagorinsky_nut_kernel.h:208
lij_compute_part1
__global__ void lij_compute_part1(T *__restrict__ l11, T *__restrict__ l22, T *__restrict__ l33, T *__restrict__ l12, T *__restrict__ l13, T *__restrict__ l23, const T *__restrict__ u, const T *__restrict__ v, const T *__restrict__ w, const T *__restrict__ fu, const T *__restrict__ fv, const T *__restrict__ fw, T *__restrict__ fuu, T *__restrict__ fvv, T *__restrict__ fww, T *__restrict__ fuv, T *__restrict__ fuw, T *__restrict__ fvw, const int n)
Definition dynamic_smagorinsky_nut_kernel.h:74
lij_compute_part2
__global__ void lij_compute_part2(T *__restrict__ l11, T *__restrict__ l22, T *__restrict__ l33, T *__restrict__ l12, T *__restrict__ l13, T *__restrict__ l23, const T *__restrict__ fuu, const T *__restrict__ fvv, const T *__restrict__ fww, const T *__restrict__ fuv, const T *__restrict__ fuw, const T *__restrict__ fvw, const int n)
Definition dynamic_smagorinsky_nut_kernel.h:122
s_abs_compute
__global__ void s_abs_compute(T *__restrict__ s_abs, const T *__restrict__ s11, const T *__restrict__ s22, const T *__restrict__ s33, const T *__restrict__ s12, const T *__restrict__ s13, const T *__restrict__ s23, const int n)
Definition dynamic_smagorinsky_nut_kernel.h:43
max
#define max(a, b)
Definition tensor.cu:40