Loading [MathJax]/extensions/tex2jax.js
Neko 0.9.99
A portable framework for high-order spectral element flow simulations
All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Pages
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>
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>
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,
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>
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>
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__ fs_abs,
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,
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 fs_abs_r = fs_abs[i];
185
186 m11[i] = delta_ratio2 * (fs_abs_r * fs11[i]);
187 m22[i] = delta_ratio2 * (fs_abs_r * fs22[i]);
188 m33[i] = delta_ratio2 * (fs_abs_r * fs33[i]);
189 m12[i] = delta_ratio2 * (fs_abs_r * fs12[i]);
190 m13[i] = delta_ratio2 * (fs_abs_r * fs13[i]);
191 m23[i] = delta_ratio2 * (fs_abs_r * fs23[i]);
192
193 fsabss11[i] = s_abs_r * s11[i];
194 fsabss22[i] = s_abs_r * s22[i];
195 fsabss33[i] = s_abs_r * s33[i];
196 fsabss12[i] = s_abs_r * s12[i];
197 fsabss13[i] = s_abs_r * s13[i];
198 fsabss23[i] = s_abs_r * s23[i];
199 }
200}
201
202template< typename T>
204 const T * __restrict__ m22,
205 const T * __restrict__ m33,
206 const T * __restrict__ m12,
207 const T * __restrict__ m13,
208 const T * __restrict__ m23,
215 const T * __restrict__ fsabss11,
216 const T * __restrict__ fsabss22,
217 const T * __restrict__ fsabss33,
218 const T * __restrict__ fsabss12,
219 const T * __restrict__ fsabss13,
220 const T * __restrict__ fsabss23,
221 T * __restrict__ num,
222 T * __restrict__ den,
223 T * __restrict__ c_dyn,
224 const T * __restrict__ delta,
225 const T * __restrict__ s_abs,
226 T * __restrict__ nut,
227 const T alpha,
228 const T * __restrict__ mult,
229 const int n){
230
231 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
232 const int str = blockDim.x * gridDim.x;
233
235
236 for (int i = idx; i < n; i += str) {
237 const T delta_r = delta[i];
238 const T delta2 = delta_r * delta_r;
239 T c_dyn_r;
240
241 const T m11_r = delta2 * (m11[i] - fsabss11[i]);
242 const T m22_r = delta2 * (m22[i] - fsabss22[i]);
243 const T m33_r = delta2 * (m33[i] - fsabss33[i]);
244 const T m12_r = delta2 * (m12[i] - fsabss12[i]);
245 const T m13_r = delta2 * (m13[i] - fsabss13[i]);
246 const T m23_r = delta2 * (m23[i] - fsabss23[i]);
247
248 num_curr = m11_r * l11[i]
249 + m22_r * l22[i]
250 + m33_r * l33[i]
251 + 2.0 * (m12_r * l12[i]
252 + m13_r * l13[i]
253 + m23_r * l23[i]);
254
256 + m22_r * m22_r
257 + m33_r * m33_r
258 + 2.0 * (m12_r * m12_r
259 + m13_r * m13_r
260 + m23_r * m23_r);
261
262 const T num_r = alpha * num[i] + (1.0 - alpha) * num_curr;
263 const T den_r = alpha * den[i] + (1.0 - alpha) * den_curr;
264
265 num[i] = num_r;
266 den[i] = den_r;
267
268 if (den_r > 0.0) {
269 c_dyn_r = 0.5 * num_r / den_r;
270 } else {
271 c_dyn_r = 0.0;
272 }
273
274 c_dyn_r = max(c_dyn_r, 0.0);
275 c_dyn[i] = c_dyn_r;
276 nut[i] = c_dyn_r * delta2 * s_abs[i] * mult[i];;
277
278 }
279}
280#endif // __COMMON_DYNAMIC_SMAGORINSKY_NUT_KERNEL_H__
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ w
const int i
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ u
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ v
__global__ void dirichlet_apply_scalar_kernel(const int *__restrict__ msk, T *__restrict__ x, const T g, const int m)
__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__ fs_abs, 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)
__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 T *__restrict__ mult, const int n)
__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)
__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)
__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)
#define max(a, b)
Definition tensor.cu:40