Neko 0.9.99
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 T * __restrict__ s11,
45 T * __restrict__ s22,
46 T * __restrict__ s33,
47 T * __restrict__ s12,
48 T * __restrict__ s13,
49 T * __restrict__ s23,
50 T * __restrict__ mult,
51 const int n){
52
53 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
54 const int str = blockDim.x * gridDim.x;
55
56 for (int i = idx; i < n; i += str) {
57 s11[i] = s11[i] * mult[i];
58 s22[i] = s22[i] * mult[i];
59 s33[i] = s33[i] * mult[i];
60 s12[i] = s12[i] * mult[i];
61 s13[i] = s13[i] * mult[i];
62 s23[i] = s23[i] * mult[i];
63
64 s_abs[i] = sqrt(2.0 * (s11[i] * s11[i] +
65 s22[i] * s22[i] +
66 s33[i] * s33[i]) +
67 4.0 * (s12[i] * s12[i] +
68 s13[i] * s13[i] +
69 s23[i] * s23[i]));
70
71 }
72}
73
74template< typename T>
75__global__ void lij_compute_part1(T * __restrict__ l11,
76 T * __restrict__ l22,
77 T * __restrict__ l33,
78 T * __restrict__ l12,
79 T * __restrict__ l13,
80 T * __restrict__ l23,
81 T * __restrict__ u,
82 T * __restrict__ v,
83 T * __restrict__ w,
84 T * __restrict__ fu,
85 T * __restrict__ fv,
86 T * __restrict__ fw,
87 T * __restrict__ fuu,
88 T * __restrict__ fvv,
89 T * __restrict__ fww,
90 T * __restrict__ fuv,
91 T * __restrict__ fuw,
92 T * __restrict__ fvw,
93 const int n){
94
95 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
96 const int str = blockDim.x * gridDim.x;
97
98 for (int i = idx; i < n; i += str) {
99 l11[i] = fu[i] * fu[i];
100 l22[i] = fv[i] * fv[i];
101 l33[i] = fw[i] * fw[i];
102 l12[i] = fu[i] * fv[i];
103 l13[i] = fu[i] * fw[i];
104 l23[i] = fv[i] * fw[i];
105
106 fuu[i] = u[i] * u[i];
107 fvv[i] = v[i] * v[i];
108 fww[i] = w[i] * w[i];
109 fuv[i] = u[i] * v[i];
110 fuw[i] = u[i] * w[i];
111 fvw[i] = v[i] * w[i];
112 }
113}
114
115template< typename T>
116__global__ void lij_compute_part2(T * __restrict__ l11,
117 T * __restrict__ l22,
118 T * __restrict__ l33,
119 T * __restrict__ l12,
120 T * __restrict__ l13,
121 T * __restrict__ l23,
122 T * __restrict__ fuu,
123 T * __restrict__ fvv,
124 T * __restrict__ fww,
125 T * __restrict__ fuv,
126 T * __restrict__ fuw,
127 T * __restrict__ fvw,
128 const int n){
129
130 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
131 const int str = blockDim.x * gridDim.x;
132
133 for (int i = idx; i < n; i += str) {
134 l11[i] -= fuu[i];
135 l22[i] -= fvv[i];
136 l33[i] -= fww[i];
137 l12[i] -= fuv[i];
138 l13[i] -= fuw[i];
139 l23[i] -= fvw[i];
140 }
141}
142
143template< typename T>
144__global__ void mij_compute_part1(T * __restrict__ m11,
145 T * __restrict__ m22,
146 T * __restrict__ m33,
147 T * __restrict__ m12,
148 T * __restrict__ m13,
149 T * __restrict__ m23,
150 T * __restrict__ s_abs,
151 T * __restrict__ s11,
152 T * __restrict__ s22,
153 T * __restrict__ s33,
154 T * __restrict__ s12,
155 T * __restrict__ s13,
156 T * __restrict__ s23,
157 T * __restrict__ fs_abs,
158 T * __restrict__ fs11,
159 T * __restrict__ fs22,
160 T * __restrict__ fs33,
161 T * __restrict__ fs12,
162 T * __restrict__ fs13,
163 T * __restrict__ fs23,
164 T * __restrict__ fsabss11,
165 T * __restrict__ fsabss22,
166 T * __restrict__ fsabss33,
167 T * __restrict__ fsabss12,
168 T * __restrict__ fsabss13,
169 T * __restrict__ fsabss23,
170 const T delta_ratio2,
171 const int n){
172
173 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
174 const int str = blockDim.x * gridDim.x;
175
176 for (int i = idx; i < n; i += str) {
177 m11[i] = delta_ratio2 * (fs_abs[i] * fs11[i]);
178 m22[i] = delta_ratio2 * (fs_abs[i] * fs22[i]);
179 m33[i] = delta_ratio2 * (fs_abs[i] * fs33[i]);
180 m12[i] = delta_ratio2 * (fs_abs[i] * fs12[i]);
181 m13[i] = delta_ratio2 * (fs_abs[i] * fs13[i]);
182 m23[i] = delta_ratio2 * (fs_abs[i] * fs23[i]);
183
184 fsabss11[i] = s_abs[i] * s11[i];
185 fsabss22[i] = s_abs[i] * s22[i];
186 fsabss33[i] = s_abs[i] * s33[i];
187 fsabss12[i] = s_abs[i] * s12[i];
188 fsabss13[i] = s_abs[i] * s13[i];
189 fsabss23[i] = s_abs[i] * s23[i];
190 }
191}
192
193template< typename T>
194__global__ void mij_nut_compute_part2(T * __restrict__ m11,
195 T * __restrict__ m22,
196 T * __restrict__ m33,
197 T * __restrict__ m12,
198 T * __restrict__ m13,
199 T * __restrict__ m23,
200 T * __restrict__ l11,
201 T * __restrict__ l22,
202 T * __restrict__ l33,
203 T * __restrict__ l12,
204 T * __restrict__ l13,
205 T * __restrict__ l23,
206 T * __restrict__ fsabss11,
207 T * __restrict__ fsabss22,
208 T * __restrict__ fsabss33,
209 T * __restrict__ fsabss12,
210 T * __restrict__ fsabss13,
211 T * __restrict__ fsabss23,
212 T * __restrict__ num,
213 T * __restrict__ den,
214 T * __restrict__ c_dyn,
215 T * __restrict__ delta,
216 T * __restrict__ s_abs,
217 T * __restrict__ nut,
218 const T alpha,
219 const int n){
220
221 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
222 const int str = blockDim.x * gridDim.x;
223
224 for (int i = idx; i < n; i += str) {
225 T num_curr, den_curr;
226 T delta2 = delta[i] * delta[i];
227
228 m11[i] = delta2 * (m11[i] - fsabss11[i]);
229 m22[i] = delta2 * (m22[i] - fsabss22[i]);
230 m33[i] = delta2 * (m33[i] - fsabss33[i]);
231 m12[i] = delta2 * (m12[i] - fsabss12[i]);
232 m13[i] = delta2 * (m13[i] - fsabss13[i]);
233 m23[i] = delta2 * (m23[i] - fsabss23[i]);
234
235 num_curr = m11[i] * l11[i]
236 + m22[i] * l22[i]
237 + m33[i] * l33[i]
238 + 2.0 * (m12[i] * l12[i]
239 + m13[i] * l13[i]
240 + m23[i] * l23[i]);
241
242 den_curr = m11[i] * m11[i]
243 + m22[i] * m22[i]
244 + m33[i] * m33[i]
245 + 2.0 * (m12[i] * m12[i]
246 + m13[i] * m13[i]
247 + m23[i] * m23[i]);
248
249 num[i] = alpha * num[i] + (1.0 - alpha) * num_curr;
250 den[i] = alpha * den[i] + (1.0 - alpha) * den_curr;
251
252 if (den[i] > 0.0) {
253 c_dyn[i] = 0.5 * num[i] / den[i];
254 } else {
255 c_dyn[i] = 0.0;
256 }
257
258 c_dyn[i] = max(c_dyn[i], 0.0);
259 nut[i] = c_dyn[i] * delta2 * s_abs[i];
260
261 }
262}
263#endif // __COMMON_DYNAMIC_SMAGORINSKY_NUT_KERNEL_H__
w
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ w
Definition ax_helm_full_kernel.h:48
i
const int i
Definition ax_helm_full_kernel.h:95
u
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ u
Definition ax_helm_full_kernel.h:46
v
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ v
Definition ax_helm_full_kernel.h:47
dirichlet_apply_scalar_kernel
__global__ void dirichlet_apply_scalar_kernel(const int *__restrict__ msk, T *__restrict__ x, const T g, const int m)
Definition dirichlet_kernel.h:42
mij_nut_compute_part2
__global__ void mij_nut_compute_part2(T *__restrict__ m11, T *__restrict__ m22, T *__restrict__ m33, T *__restrict__ m12, T *__restrict__ m13, T *__restrict__ m23, T *__restrict__ l11, T *__restrict__ l22, T *__restrict__ l33, T *__restrict__ l12, T *__restrict__ l13, T *__restrict__ l23, T *__restrict__ fsabss11, T *__restrict__ fsabss22, T *__restrict__ fsabss33, T *__restrict__ fsabss12, T *__restrict__ fsabss13, T *__restrict__ fsabss23, T *__restrict__ num, T *__restrict__ den, T *__restrict__ c_dyn, T *__restrict__ delta, T *__restrict__ s_abs, T *__restrict__ nut, const T alpha, const int n)
Definition dynamic_smagorinsky_nut_kernel.h:194
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, T *__restrict__ s_abs, T *__restrict__ s11, T *__restrict__ s22, T *__restrict__ s33, T *__restrict__ s12, T *__restrict__ s13, T *__restrict__ s23, T *__restrict__ fs_abs, T *__restrict__ fs11, T *__restrict__ fs22, T *__restrict__ fs33, T *__restrict__ fs12, T *__restrict__ fs13, 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:144
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, T *__restrict__ u, T *__restrict__ v, T *__restrict__ w, T *__restrict__ fu, T *__restrict__ fv, 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:75
s_abs_compute
__global__ void s_abs_compute(T *__restrict__ s_abs, T *__restrict__ s11, T *__restrict__ s22, T *__restrict__ s33, T *__restrict__ s12, T *__restrict__ s13, T *__restrict__ s23, T *__restrict__ mult, const int n)
Definition dynamic_smagorinsky_nut_kernel.h:43
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, 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:116
max
#define max(a, b)
Definition tensor.cu:40