Neko  0.9.99
A portable framework for high-order spectral element flow simulations
opgrad_kernel.h
Go to the documentation of this file.
1 #ifndef __MATH_OPGRAD_KERNEL_H__
2 #define __MATH_OPGRAD_KERNEL_H__
3 /*
4  Copyright (c) 2021-2023, 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 
41 template< typename T, const int LX, const int CHUNKS >
42 __global__ void opgrad_kernel_1d(T * __restrict__ ux,
43  T * __restrict__ uy,
44  T * __restrict__ uz,
45  const T * __restrict__ u,
46  const T * __restrict__ dx,
47  const T * __restrict__ dy,
48  const T * __restrict__ dz,
49  const T * __restrict__ drdx,
50  const T * __restrict__ dsdx,
51  const T * __restrict__ dtdx,
52  const T * __restrict__ drdy,
53  const T * __restrict__ dsdy,
54  const T * __restrict__ dtdy,
55  const T * __restrict__ drdz,
56  const T * __restrict__ dsdz,
57  const T * __restrict__ dtdz,
58  const T * __restrict__ w3) {
59 
60  __shared__ T shu[LX * LX * LX];
61 
62  __shared__ T shdx[LX * LX];
63  __shared__ T shdy[LX * LX];
64  __shared__ T shdz[LX * LX];
65 
66 
67  int i,j,k;
68 
69  const int e = blockIdx.x;
70  const int iii = threadIdx.x;
71  const int nchunks = (LX * LX * LX - 1) / CHUNKS + 1;
72 
73  if (iii < (LX * LX)) {
74  shdx[iii] = dx[iii];
75  shdy[iii] = dy[iii];
76  shdz[iii] = dz[iii];
77  }
78 
79  j = iii;
80  while(j < (LX * LX * LX)) {
81  shu[j] = u[j + e * LX * LX * LX];
82  j = j + CHUNKS;
83  }
84 
85  __syncthreads();
86 
87  for (int n = 0; n < nchunks; n++) {
88  const int ijk = iii + n * CHUNKS;
89  const int jk = ijk / LX;
90  i = ijk - jk * LX;
91  k = jk / LX;
92  j = jk - k * LX;
93  if ( i < LX && j < LX && k < LX ) {
94  T rtmp = 0.0;
95  T stmp = 0.0;
96  T ttmp = 0.0;
97  for (int l = 0; l < LX; l++) {
98  rtmp += shdx[i + l * LX] * shu[l + j * LX + k * LX * LX];
99  stmp += shdy[j + l * LX] * shu[i + l * LX + k * LX * LX];
100  ttmp += shdz[k + l * LX] * shu[i + j * LX + l * LX * LX];
101  }
102 
103  ux[ijk + e * LX * LX * LX] = w3[ijk]
104  * (drdx[ijk + e * LX * LX * LX] * rtmp
105  + dsdx[ijk + e * LX * LX * LX] * stmp
106  + dtdx[ijk + e * LX * LX * LX] * ttmp);
107 
108  uy[ijk + e * LX * LX * LX] = w3[ijk]
109  * (drdy[ijk + e * LX * LX * LX] * rtmp
110  + dsdy[ijk + e * LX * LX * LX] * stmp
111  + dtdy[ijk + e * LX * LX * LX] * ttmp);
112 
113  uz[ijk + e * LX * LX * LX] = w3[ijk]
114  * (drdz[ijk + e * LX * LX * LX] * rtmp
115  + dsdz[ijk + e * LX * LX * LX] * stmp
116  + dtdz[ijk + e * LX * LX * LX] * ttmp);
117 
118  }
119  }
120 
121 }
122 
123 template< typename T, const int LX >
124 __global__ void __launch_bounds__(LX*LX,3)
125 opgrad_kernel_kstep(T * __restrict__ ux,
126  T * __restrict__ uy,
127  T * __restrict__ uz,
128  const T * __restrict__ u,
129  const T * __restrict__ dx,
130  const T * __restrict__ dy,
131  const T * __restrict__ dz,
132  const T * __restrict__ drdx,
133  const T * __restrict__ dsdx,
134  const T * __restrict__ dtdx,
135  const T * __restrict__ drdy,
136  const T * __restrict__ dsdy,
137  const T * __restrict__ dtdy,
138  const T * __restrict__ drdz,
139  const T * __restrict__ dsdz,
140  const T * __restrict__ dtdz,
141  const T * __restrict__ w3) {
142 
143  __shared__ T shu[LX * LX];
144 
145  __shared__ T shdx[LX * LX];
146  __shared__ T shdy[LX * LX];
147  __shared__ T shdz[LX * LX];
148 
149  const int e = blockIdx.x;
150  const int j = threadIdx.y;
151  const int i = threadIdx.x;
152  const int ij = i + j * LX;
153  const int ele = e*LX*LX*LX;
154 
155  shdx[ij] = dx[ij];
156  shdy[ij] = dy[ij];
157  shdz[ij] = dz[ij];
158 
159  T ru[LX];
160 
161 #pragma unroll LX
162  for (int k = 0; k < LX; ++k) {
163  ru[k] = u[ij + k*LX*LX + ele];
164  }
165 
167 
168  #pragma unroll
169  for (int k = 0; k < LX; ++k) {
170  const int ijk = ij + k*LX*LX;
171  const T W3 = w3[ijk];
172  T ttmp = 0.0;
173  shu[ij] = ru[k];
174  for (int l = 0; l < LX; l++) {
175  ttmp += shdz[k+l*LX] * ru[l];
176  }
177  __syncthreads();
178 
179  T rtmp = 0.0;
180  T stmp = 0.0;
181 #pragma unroll
182  for (int l = 0; l < LX; l++) {
183  rtmp += shdx[i+l*LX] * shu[l+j*LX];
184  stmp += shdy[j+l*LX] * shu[i+l*LX];
185  }
186 
187  ux[ijk + ele] = W3 * (drdx[ijk + ele] * rtmp
188  + dsdx[ijk + ele] * stmp
189  + dtdx[ijk + ele] * ttmp);
190 
191  uy[ijk + ele] = W3 * (drdy[ijk + ele] * rtmp
192  + dsdy[ijk + ele] * stmp
193  + dtdy[ijk + ele] * ttmp);
194 
195  uz[ijk + ele] = W3 * (drdz[ijk + ele] * rtmp
196  + dsdz[ijk + ele] * stmp
197  + dtdz[ijk + ele] * ttmp);
198  __syncthreads();
199  }
200 }
201 
202 
203 #endif // __MATH_OPGRAD_KERNEL_H__
__shared__ T shu[LX *LX]
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ dsdx
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ dtdz
__global__ void __launch_bounds__(LX *LX, 3) opgrad_kernel_kstep(T *__restrict__ ux
__global__ void T *__restrict__ uy
__global__ void opgrad_kernel_1d(T *__restrict__ ux, T *__restrict__ uy, T *__restrict__ uz, const T *__restrict__ u, const T *__restrict__ dx, const T *__restrict__ dy, const T *__restrict__ dz, const T *__restrict__ drdx, const T *__restrict__ dsdx, const T *__restrict__ dtdx, const T *__restrict__ drdy, const T *__restrict__ dsdy, const T *__restrict__ dtdy, const T *__restrict__ drdz, const T *__restrict__ dsdz, const T *__restrict__ dtdz, const T *__restrict__ w3)
Definition: opgrad_kernel.h:42
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ dx
__global__ void T *__restrict__ T *__restrict__ uz
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ dsdz
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ dz
T ru[LX]
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ drdy
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ drdz
__shared__ T shdz[LX *LX]
const int i
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ dtdx
const int ij
__shared__ T shdy[LX *LX]
const int e
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ dtdy
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ drdx
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ dsdy
__shared__ T shdx[LX *LX]
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ dy
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ u
const int ele
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ const T *__restrict__ w3
const int j
__syncthreads()