Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
spalding_kernel.h
Go to the documentation of this file.
1#ifndef __COMMON_SPALDING_KERNEL_H__
2#define __COMMON_SPALDING_KERNEL_H__
3/*
4 Copyright (c) 2025, 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__device__ T solve(const T u, const T y, const T guess, const T nu,
44 const T kappa, const T B);
45
49template<typename T>
51 const T * __restrict__ v_d,
52 const T * __restrict__ w_d,
53 const int * __restrict__ ind_r_d,
54 const int * __restrict__ ind_s_d,
55 const int * __restrict__ ind_t_d,
56 const int * __restrict__ ind_e_d,
57 const T * __restrict__ n_x_d,
58 const T * __restrict__ n_y_d,
59 const T * __restrict__ n_z_d,
60 const T * __restrict__ nu_d,
61 const T * __restrict__ h_d,
65 const int n_nodes,
66 const int lx,
67 const T kappa,
68 const T B,
69 const int tstep) {
70
71 const int idx = blockIdx.x * blockDim.x + threadIdx.x;
72 const int str = blockDim.x * gridDim.x;
73 for (int i = idx; i < n_nodes; i += str) {
74 // Sample the velocity
75 const int index = (ind_e_d[i] - 1) * lx * lx * lx +
76 (ind_t_d[i] - 1) * lx * lx +
77 (ind_s_d[i] - 1) * lx +
78 (ind_r_d[i] - 1);
79
80 T ui = u_d[index];
81 T vi = v_d[index];
82 T wi = w_d[index];
83
84 // Load normal vectors and wall shear stress values once
85 T nx = n_x_d[i];
86 T ny = n_y_d[i];
87 T nz = n_z_d[i];
88 T h = h_d[i];
89
90 // Project on tangential direction
91 T normu = ui * nx + vi * ny + wi * nz;
92
93 ui -= normu * nx;
94 vi -= normu * ny;
95 wi -= normu * nz;
96
97 T magu = sqrt(ui * ui + vi * vi + wi * wi);
98
99 // Get initial guess for Newton solver
100 T guess;
101 if (tstep == 1) {
102 guess = sqrt(magu * nu_d[i] / h);
103 } else {
104 guess = tau_x_d[i] * tau_x_d[i] +
105 tau_y_d[i] * tau_y_d[i] +
106 tau_z_d[i] * tau_z_d[i];
107 guess = sqrt(sqrt(guess));
108 }
109
110 // Solve for utau using Newton's method
111 T utau = solve(magu, h, guess, nu_d[i], kappa, B);
112
113 // Distribute according to the velocity vector
114 tau_x_d[i] = -utau * utau * ui / magu;
115 tau_y_d[i] = -utau * utau * vi / magu;
116 tau_z_d[i] = -utau * utau * wi / magu;
117 }
118}
119
123template<typename T>
124__device__ T solve(const T u, const T y, const T guess, const T nu,
125 const T kappa, const T B) {
126 T utau = guess;
127 T yp, up, f, df, old, error;
128 const int maxiter = 100;
129
130 for (int k = 0; k < maxiter; ++k) {
131 up = u / utau;
132 yp = y * utau / nu;
133 old = utau;
134
135 // Evaluate function and its derivative
136 f = (up + exp(-kappa * B) *
137 (exp(kappa * up) - 1.0 - kappa * up -
138 0.5 * (kappa * up) * (kappa * up) -
139 (1.0 / 6.0) * (kappa * up) * (kappa * up) * (kappa * up)) -
140 yp);
141
142 df = (-y / nu - u / (utau * utau) -
143 kappa * up / utau * exp(-kappa * B) *
144 (exp(kappa * up) - 1.0 - kappa * up -
145 0.5 * (kappa * up) * (kappa * up)));
146
147 // Update solution
148 utau -= f / df;
149
150 error = fabs((old - utau) / old);
151
152 if (error < 1e-3) {
153 break;
154 }
155 }
156
157 return utau;
158}
159
160#endif // __COMMON_SPALDING_KERNEL_H__
const int i
__global__ void T *__restrict__ T *__restrict__ const T *__restrict__ u
const int e
__global__ void dirichlet_apply_scalar_kernel(const int *__restrict__ msk, T *__restrict__ x, const T g, const int m)
__global__ void spalding_compute(const T *__restrict__ u_d, const T *__restrict__ v_d, const T *__restrict__ w_d, const int *__restrict__ ind_r_d, const int *__restrict__ ind_s_d, const int *__restrict__ ind_t_d, const int *__restrict__ ind_e_d, const T *__restrict__ n_x_d, const T *__restrict__ n_y_d, const T *__restrict__ n_z_d, const T *__restrict__ nu_d, const T *__restrict__ h_d, T *__restrict__ tau_x_d, T *__restrict__ tau_y_d, T *__restrict__ tau_z_d, const int n_nodes, const int lx, const T kappa, const T B, const int tstep)
solve(case_descr)
Definition intf.py:300