VirtualFluids 0.2.0
Parallel CFD LBM Solver
Loading...
Searching...
No Matches
NoSlipInterpolatedCompressible.cu
Go to the documentation of this file.
1//=======================================================================================
2// ____ ____ __ ______ __________ __ __ __ __
3// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
4// \ \ | | | | | |_) | | | | | | | / \ | |
5// \ \ | | | | | _ / | | | | | | / /\ \ | |
6// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
7// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
8// \ \ | | ________________________________________________________________
9// \ \ | | | ______________________________________________________________|
10// \ \| | | | __ __ __ __ ______ _______
11// \ | | |_____ | | | | | | | | | _ \ / _____)
12// \ | | _____| | | | | | | | | | | \ \ \_______
13// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
14// \ _____| |__| |________| \_______/ |__| |______/ (_______/
15//
16// This file is part of VirtualFluids. VirtualFluids is free software: you can
17// redistribute it and/or modify it under the terms of the GNU General Public
18// License as published by the Free Software Foundation, either version 3 of
19// the License, or (at your option) any later version.
20//
21// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
22// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
23// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
24// for more details.
25//
26// SPDX-License-Identifier: GPL-3.0-or-later
27// SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
28//
29//! \addtogroup gpu_BoundaryConditions BoundaryConditions
30//! \ingroup gpu_core core
31//! \{
32//! \author Martin Schoenherr, Anna Wellmann
33//======================================================================================
34#include "Calculation/Calculation.h"
35#include "lbm/constants/D3Q27.h"
36#include "basics/constants/NumericConstants.h"
37#include "Utilities/KernelUtilities.h"
38#include "cuda_helper/CudaIndexCalculation.h"
39
40using namespace vf::basics::constant;
41using namespace vf::lbm::dir;
42namespace vf::gpu {
43
44__global__ void NoSlipInterpolatedCompressible_Device(
45 real* distributions,
46 int* subgridDistanceIndices,
47 real* subgridDistances,
48 unsigned int numberOfBCnodes,
49 real omega,
50 unsigned int* neighborX,
51 unsigned int* neighborY,
52 unsigned int* neighborZ,
53 unsigned long long numberOfLBnodes,
54 bool isEvenTimestep)
55{
56 //////////////////////////////////////////////////////////////////////////
57 //! The no-slip boundary condition is executed in the following steps
58 //!
59 ////////////////////////////////////////////////////////////////////////////////
60 //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
61 //!
62 const unsigned nodeIndex = vf::cuda::get1DIndexFrom2DBlock();
63
64 if(nodeIndex < numberOfBCnodes)
65 {
66 //////////////////////////////////////////////////////////////////////////
67 //! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on timestep is based on the esoteric twist algorithm \ref
68 //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
69 //!
70 Distributions27 dist;
71 getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
72
73 ////////////////////////////////////////////////////////////////////////////////
74 //! - Set local subgrid distances (q's)
75 //!
76 SubgridDistances27 subgridD;
77 getPointersToSubgridDistances(subgridD, subgridDistances, numberOfBCnodes);
78
79 ////////////////////////////////////////////////////////////////////////////////
80 //! - Set neighbor indices (necessary for indirect addressing)
81 //!
82 unsigned int indexOfBCnode = subgridDistanceIndices[nodeIndex];
83 unsigned int kzero= indexOfBCnode;
84 unsigned int ke = indexOfBCnode;
85 unsigned int kw = neighborX[indexOfBCnode];
86 unsigned int kn = indexOfBCnode;
87 unsigned int ks = neighborY[indexOfBCnode];
88 unsigned int kt = indexOfBCnode;
89 unsigned int kb = neighborZ[indexOfBCnode];
90 unsigned int ksw = neighborY[kw];
91 unsigned int kne = indexOfBCnode;
92 unsigned int kse = ks;
93 unsigned int knw = kw;
94 unsigned int kbw = neighborZ[kw];
95 unsigned int kte = indexOfBCnode;
96 unsigned int kbe = kb;
97 unsigned int ktw = kw;
98 unsigned int kbs = neighborZ[ks];
99 unsigned int ktn = indexOfBCnode;
100 unsigned int kbn = kb;
101 unsigned int kts = ks;
102 unsigned int ktse = ks;
103 unsigned int kbnw = kbw;
104 unsigned int ktnw = kw;
105 unsigned int kbse = kbs;
106 unsigned int ktsw = ksw;
107 unsigned int kbne = kb;
108 unsigned int ktne = indexOfBCnode;
109 unsigned int kbsw = neighborZ[ksw];
110
111 ////////////////////////////////////////////////////////////////////////////////
112 //! - Set local distributions
113 //!
114 real f_W = (dist.f[dP00])[ke ];
115 real f_E = (dist.f[dM00])[kw ];
116 real f_S = (dist.f[d0P0])[kn ];
117 real f_N = (dist.f[d0M0])[ks ];
118 real f_B = (dist.f[d00P])[kt ];
119 real f_T = (dist.f[d00M])[kb ];
120 real f_SW = (dist.f[dPP0])[kne ];
121 real f_NE = (dist.f[dMM0])[ksw ];
122 real f_NW = (dist.f[dPM0])[kse ];
123 real f_SE = (dist.f[dMP0])[knw ];
124 real f_BW = (dist.f[dP0P])[kte ];
125 real f_TE = (dist.f[dM0M])[kbw ];
126 real f_TW = (dist.f[dP0M])[kbe ];
127 real f_BE = (dist.f[dM0P])[ktw ];
128 real f_BS = (dist.f[d0PP])[ktn ];
129 real f_TN = (dist.f[d0MM])[kbs ];
130 real f_TS = (dist.f[d0PM])[kbn ];
131 real f_BN = (dist.f[d0MP])[kts ];
132 real f_BSW = (dist.f[dPPP])[ktne ];
133 real f_BNE = (dist.f[dMMP])[ktsw ];
134 real f_BNW = (dist.f[dPMP])[ktse ];
135 real f_BSE = (dist.f[dMPP])[ktnw ];
136 real f_TSW = (dist.f[dPPM])[kbne ];
137 real f_TNE = (dist.f[dMMM])[kbsw ];
138 real f_TNW = (dist.f[dPMM])[kbse ];
139 real f_TSE = (dist.f[dMPM])[kbnw ];
140
141 ////////////////////////////////////////////////////////////////////////////////
142 //! - Calculate macroscopic quantities
143 //!
144 real drho = f_TSE + f_TNW + f_TNE + f_TSW + f_BSE + f_BNW + f_BNE + f_BSW +
145 f_BN + f_TS + f_TN + f_BS + f_BE + f_TW + f_TE + f_BW + f_SE + f_NW + f_NE + f_SW +
146 f_T + f_B + f_N + f_S + f_E + f_W + ((dist.f[d000])[kzero]);
147
148 real vx1 = (((f_TSE - f_BNW) - (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
149 ((f_BE - f_TW) + (f_TE - f_BW)) + ((f_SE - f_NW) + (f_NE - f_SW)) +
150 (f_E - f_W)) / (c1o1 + drho);
151
152 real vx2 = ((-(f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
153 ((f_BN - f_TS) + (f_TN - f_BS)) + (-(f_SE - f_NW) + (f_NE - f_SW)) +
154 (f_N - f_S)) / (c1o1 + drho);
155
156 real vx3 = (((f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) + (f_TSW - f_BNE)) +
157 (-(f_BN - f_TS) + (f_TN - f_BS)) + ((f_TE - f_BW) - (f_BE - f_TW)) +
158 (f_T - f_B)) / (c1o1 + drho);
159
160 real cu_sq = c3o2 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3) * (c1o1 + drho);
161
162 ////////////////////////////////////////////////////////////////////////////////
163 //! - change the pointer to write the results in the correct array
164 //!
165 getPointersToDistributions(dist, distributions, numberOfLBnodes, !isEvenTimestep);
166
167 ////////////////////////////////////////////////////////////////////////////////
168 //! - Update distributions with subgrid distance (q) between zero and one
169 real feq, q, velocityLB;
170 q = (subgridD.q[dP00])[nodeIndex];
171 if (q>=c0o1 && q<=c1o1) // only update distribution for q between zero and one
172 {
173 velocityLB = vx1;
174 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
175 (dist.f[dM00])[kw] = getInterpolatedDistributionForNoSlipBC(q, f_E, f_W, feq, omega);
176 }
177
178 q = (subgridD.q[dM00])[nodeIndex];
179 if (q>=c0o1 && q<=c1o1)
180 {
181 velocityLB = -vx1;
182 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
183 (dist.f[dP00])[ke] = getInterpolatedDistributionForNoSlipBC(q, f_W, f_E, feq, omega);
184 }
185
186 q = (subgridD.q[d0P0])[nodeIndex];
187 if (q>=c0o1 && q<=c1o1)
188 {
189 velocityLB = vx2;
190 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
191 (dist.f[d0M0])[ks] = getInterpolatedDistributionForNoSlipBC(q, f_N, f_S, feq, omega);
192 }
193
194 q = (subgridD.q[d0M0])[nodeIndex];
195 if (q>=c0o1 && q<=c1o1)
196 {
197 velocityLB = -vx2;
198 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
199 (dist.f[d0P0])[kn] = getInterpolatedDistributionForNoSlipBC(q, f_S, f_N, feq, omega);
200 }
201
202 q = (subgridD.q[d00P])[nodeIndex];
203 if (q>=c0o1 && q<=c1o1)
204 {
205 velocityLB = vx3;
206 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
207 (dist.f[d00M])[kb] = getInterpolatedDistributionForNoSlipBC(q, f_T, f_B, feq, omega);
208 }
209
210 q = (subgridD.q[d00M])[nodeIndex];
211 if (q>=c0o1 && q<=c1o1)
212 {
213 velocityLB = -vx3;
214 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
215 (dist.f[d00P])[kt] = getInterpolatedDistributionForNoSlipBC(q, f_B, f_T, feq, omega);
216 }
217
218 q = (subgridD.q[dPP0])[nodeIndex];
219 if (q>=c0o1 && q<=c1o1)
220 {
221 velocityLB = vx1 + vx2;
222 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
223 (dist.f[dMM0])[ksw] = getInterpolatedDistributionForNoSlipBC(q, f_NE, f_SW, feq, omega);
224 }
225
226 q = (subgridD.q[dMM0])[nodeIndex];
227 if (q>=c0o1 && q<=c1o1)
228 {
229 velocityLB = -vx1 - vx2;
230 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
231 (dist.f[dPP0])[kne] = getInterpolatedDistributionForNoSlipBC(q, f_SW, f_NE, feq, omega);
232 }
233
234 q = (subgridD.q[dPM0])[nodeIndex];
235 if (q>=c0o1 && q<=c1o1)
236 {
237 velocityLB = vx1 - vx2;
238 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
239 (dist.f[dMP0])[knw] = getInterpolatedDistributionForNoSlipBC(q, f_SE, f_NW, feq, omega);
240 }
241
242 q = (subgridD.q[dMP0])[nodeIndex];
243 if (q>=c0o1 && q<=c1o1)
244 {
245 velocityLB = -vx1 + vx2;
246 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
247 (dist.f[dPM0])[kse] = getInterpolatedDistributionForNoSlipBC(q, f_NW, f_SE, feq, omega);
248 }
249
250 q = (subgridD.q[dP0P])[nodeIndex];
251 if (q>=c0o1 && q<=c1o1)
252 {
253 velocityLB = vx1 + vx3;
254 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
255 (dist.f[dM0M])[kbw] = getInterpolatedDistributionForNoSlipBC(q, f_TE, f_BW, feq, omega);
256 }
257
258 q = (subgridD.q[dM0M])[nodeIndex];
259 if (q>=c0o1 && q<=c1o1)
260 {
261 velocityLB = -vx1 - vx3;
262 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
263 (dist.f[dP0P])[kte] = getInterpolatedDistributionForNoSlipBC(q, f_BW, f_TE, feq, omega);
264 }
265
266 q = (subgridD.q[dP0M])[nodeIndex];
267 if (q>=c0o1 && q<=c1o1)
268 {
269 velocityLB = vx1 - vx3;
270 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
271 (dist.f[dM0P])[ktw] = getInterpolatedDistributionForNoSlipBC(q, f_BE, f_TW, feq, omega);
272 }
273
274 q = (subgridD.q[dM0P])[nodeIndex];
275 if (q>=c0o1 && q<=c1o1)
276 {
277 velocityLB = -vx1 + vx3;
278 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
279 (dist.f[dP0M])[kbe] = getInterpolatedDistributionForNoSlipBC(q, f_TW, f_BE, feq, omega);
280 }
281
282 q = (subgridD.q[d0PP])[nodeIndex];
283 if (q>=c0o1 && q<=c1o1)
284 {
285 velocityLB = vx2 + vx3;
286 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
287 (dist.f[d0MM])[kbs] = getInterpolatedDistributionForNoSlipBC(q, f_TN, f_BS, feq, omega);
288 }
289
290 q = (subgridD.q[d0MM])[nodeIndex];
291 if (q>=c0o1 && q<=c1o1)
292 {
293 velocityLB = -vx2 - vx3;
294 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
295 (dist.f[d0PP])[ktn] = getInterpolatedDistributionForNoSlipBC(q, f_BS, f_TN, feq, omega);
296 }
297
298 q = (subgridD.q[d0PM])[nodeIndex];
299 if (q>=c0o1 && q<=c1o1)
300 {
301 velocityLB = vx2 - vx3;
302 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
303 (dist.f[d0MP])[kts] = getInterpolatedDistributionForNoSlipBC(q, f_BN, f_TS, feq, omega);
304 }
305
306 q = (subgridD.q[d0MP])[nodeIndex];
307 if (q>=c0o1 && q<=c1o1)
308 {
309 velocityLB = -vx2 + vx3;
310 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
311 (dist.f[d0PM])[kbn] = getInterpolatedDistributionForNoSlipBC(q, f_TS, f_BN, feq, omega);
312 }
313
314 q = (subgridD.q[dPPP])[nodeIndex];
315 if (q>=c0o1 && q<=c1o1)
316 {
317 velocityLB = vx1 + vx2 + vx3;
318 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
319 (dist.f[dMMM])[kbsw] = getInterpolatedDistributionForNoSlipBC(q, f_TNE, f_BSW, feq, omega);
320 }
321
322 q = (subgridD.q[dMMM])[nodeIndex];
323 if (q>=c0o1 && q<=c1o1)
324 {
325 velocityLB = -vx1 - vx2 - vx3;
326 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
327 (dist.f[dPPP])[ktne] = getInterpolatedDistributionForNoSlipBC(q, f_BSW, f_TNE, feq, omega);
328 }
329
330 q = (subgridD.q[dPPM])[nodeIndex];
331 if (q>=c0o1 && q<=c1o1)
332 {
333 velocityLB = vx1 + vx2 - vx3;
334 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
335 (dist.f[dMMP])[ktsw] = getInterpolatedDistributionForNoSlipBC(q, f_BNE, f_TSW, feq, omega);
336 }
337
338 q = (subgridD.q[dMMP])[nodeIndex];
339 if (q>=c0o1 && q<=c1o1)
340 {
341 velocityLB = -vx1 - vx2 + vx3;
342 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
343 (dist.f[dPPM])[kbne] = getInterpolatedDistributionForNoSlipBC(q, f_TSW, f_BNE, feq, omega);
344 }
345
346 q = (subgridD.q[dPMP])[nodeIndex];
347 if (q>=c0o1 && q<=c1o1)
348 {
349 velocityLB = vx1 - vx2 + vx3;
350 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
351 (dist.f[dMPM])[kbnw] = getInterpolatedDistributionForNoSlipBC(q, f_TSE, f_BNW, feq, omega);
352 }
353
354 q = (subgridD.q[dMPM])[nodeIndex];
355 if (q>=c0o1 && q<=c1o1)
356 {
357 velocityLB = -vx1 + vx2 - vx3;
358 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
359 (dist.f[dPMP])[ktse] = getInterpolatedDistributionForNoSlipBC(q, f_BNW, f_TSE, feq, omega);
360 }
361
362 q = (subgridD.q[dPMM])[nodeIndex];
363 if (q>=c0o1 && q<=c1o1)
364 {
365 velocityLB = vx1 - vx2 - vx3;
366 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
367 (dist.f[dMPP])[ktnw] = getInterpolatedDistributionForNoSlipBC(q, f_BSE, f_TNW, feq, omega);
368 }
369
370 q = (subgridD.q[dMPP])[nodeIndex];
371 if (q>=c0o1 && q<=c1o1)
372 {
373 velocityLB = -vx1 + vx2 + vx3;
374 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
375 (dist.f[dPMM])[kbse] = getInterpolatedDistributionForNoSlipBC(q, f_TNW, f_BSE, feq, omega);
376 }
377 }
378}
379
380}
381
382//! \}