VirtualFluids 0.2.0
Parallel CFD LBM Solver
Loading...
Searching...
No Matches
NoSlipInterpolatedIncompressible.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 NoSlipInterpolatedIncompressible_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 //////////////////////////////////////////////////////////////////////////
65 //! - Run for all indices in size of boundary condition (numberOfBCnodes)
66 //!
67 if(nodeIndex < numberOfBCnodes)
68 {
69
70 //////////////////////////////////////////////////////////////////////////
71 //! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on timestep is based on the esoteric twist algorithm \ref
72 //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
73 //!
74 Distributions27 dist;
75 getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
76
77 ////////////////////////////////////////////////////////////////////////////////
78 //! - Set local subgrid distances (q's)
79 //!
80 SubgridDistances27 subgridD;
81 getPointersToSubgridDistances(subgridD, subgridDistances, numberOfBCnodes);
82
83 ////////////////////////////////////////////////////////////////////////////////
84 //! - Set neighbor indices (necessary for indirect addressing)
85 //!
86 unsigned int indexOfBCnode = subgridDistanceIndices[nodeIndex];
87 unsigned int kzero= indexOfBCnode;
88 unsigned int ke = indexOfBCnode;
89 unsigned int kw = neighborX[indexOfBCnode];
90 unsigned int kn = indexOfBCnode;
91 unsigned int ks = neighborY[indexOfBCnode];
92 unsigned int kt = indexOfBCnode;
93 unsigned int kb = neighborZ[indexOfBCnode];
94 unsigned int ksw = neighborY[kw];
95 unsigned int kne = indexOfBCnode;
96 unsigned int kse = ks;
97 unsigned int knw = kw;
98 unsigned int kbw = neighborZ[kw];
99 unsigned int kte = indexOfBCnode;
100 unsigned int kbe = kb;
101 unsigned int ktw = kw;
102 unsigned int kbs = neighborZ[ks];
103 unsigned int ktn = indexOfBCnode;
104 unsigned int kbn = kb;
105 unsigned int kts = ks;
106 unsigned int ktse = ks;
107 unsigned int kbnw = kbw;
108 unsigned int ktnw = kw;
109 unsigned int kbse = kbs;
110 unsigned int ktsw = ksw;
111 unsigned int kbne = kb;
112 unsigned int ktne = indexOfBCnode;
113 unsigned int kbsw = neighborZ[ksw];
114
115 ////////////////////////////////////////////////////////////////////////////////
116 //! - Set local distributions
117 //!
118 real f_W = (dist.f[dP00])[ke ];
119 real f_E = (dist.f[dM00])[kw ];
120 real f_S = (dist.f[d0P0])[kn ];
121 real f_N = (dist.f[d0M0])[ks ];
122 real f_B = (dist.f[d00P])[kt ];
123 real f_T = (dist.f[d00M])[kb ];
124 real f_SW = (dist.f[dPP0])[kne ];
125 real f_NE = (dist.f[dMM0])[ksw ];
126 real f_NW = (dist.f[dPM0])[kse ];
127 real f_SE = (dist.f[dMP0])[knw ];
128 real f_BW = (dist.f[dP0P])[kte ];
129 real f_TE = (dist.f[dM0M])[kbw ];
130 real f_TW = (dist.f[dP0M])[kbe ];
131 real f_BE = (dist.f[dM0P])[ktw ];
132 real f_BS = (dist.f[d0PP])[ktn ];
133 real f_TN = (dist.f[d0MM])[kbs ];
134 real f_TS = (dist.f[d0PM])[kbn ];
135 real f_BN = (dist.f[d0MP])[kts ];
136 real f_BSW = (dist.f[dPPP])[ktne ];
137 real f_BNE = (dist.f[dMMP])[ktsw ];
138 real f_BNW = (dist.f[dPMP])[ktse ];
139 real f_BSE = (dist.f[dMPP])[ktnw ];
140 real f_TSW = (dist.f[dPPM])[kbne ];
141 real f_TNE = (dist.f[dMMM])[kbsw ];
142 real f_TNW = (dist.f[dPMM])[kbse ];
143 real f_TSE = (dist.f[dMPM])[kbnw ];
144
145 ////////////////////////////////////////////////////////////////////////////////
146 //! - Calculate macroscopic quantities
147 //!
148 real drho = f_TSE + f_TNW + f_TNE + f_TSW + f_BSE + f_BNW + f_BNE + f_BSW +
149 f_BN + f_TS + f_TN + f_BS + f_BE + f_TW + f_TE + f_BW + f_SE + f_NW + f_NE + f_SW +
150 f_T + f_B + f_N + f_S + f_E + f_W + ((dist.f[d000])[kzero]);
151
152 real vx1 = (((f_TSE - f_BNW) - (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
153 ((f_BE - f_TW) + (f_TE - f_BW)) + ((f_SE - f_NW) + (f_NE - f_SW)) +
154 (f_E - f_W));
155
156 real vx2 = ((-(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_SE - f_NW) + (f_NE - f_SW)) +
158 (f_N - f_S));
159
160 real vx3 = (((f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) + (f_TSW - f_BNE)) +
161 (-(f_BN - f_TS) + (f_TN - f_BS)) + ((f_TE - f_BW) - (f_BE - f_TW)) +
162 (f_T - f_B));
163
164 real cu_sq = c3o2 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
165
166 ////////////////////////////////////////////////////////////////////////////////
167 //! - change the pointer to write the results in the correct array
168 //!
169 getPointersToDistributions(dist, distributions, numberOfLBnodes, !isEvenTimestep);
170
171 ////////////////////////////////////////////////////////////////////////////////
172 //! - Update distributions with subgrid distance (q) between zero and one
173 //!
174 real feq, q, velocityLB;
175 q = (subgridD.q[dP00])[nodeIndex];
176 if (q>=c0o1 && q<=c1o1) // only update distribution for q between zero and one
177 {
178 velocityLB = vx1;
179 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
180 (dist.f[dM00])[kw] = getInterpolatedDistributionForNoSlipBC(q, f_E, f_W, feq, omega);
181 }
182
183 q = (subgridD.q[dM00])[nodeIndex];
184 if (q>=c0o1 && q<=c1o1)
185 {
186 velocityLB = -vx1;
187 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
188 (dist.f[dP00])[ke] = getInterpolatedDistributionForNoSlipBC(q, f_W, f_E, feq, omega);
189 }
190
191 q = (subgridD.q[d0P0])[nodeIndex];
192 if (q>=c0o1 && q<=c1o1)
193 {
194 velocityLB = vx2;
195 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
196 (dist.f[d0M0])[ks] = getInterpolatedDistributionForNoSlipBC(q, f_N, f_S, feq, omega);
197 }
198
199 q = (subgridD.q[d0M0])[nodeIndex];
200 if (q>=c0o1 && q<=c1o1)
201 {
202 velocityLB = -vx2;
203 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
204 (dist.f[d0P0])[kn] = getInterpolatedDistributionForNoSlipBC(q, f_S, f_N, feq, omega);
205 }
206
207 q = (subgridD.q[d00P])[nodeIndex];
208 if (q>=c0o1 && q<=c1o1)
209 {
210 velocityLB = vx3;
211 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
212 (dist.f[d00M])[kb] = getInterpolatedDistributionForNoSlipBC(q, f_T, f_B, feq, omega);
213 }
214
215 q = (subgridD.q[d00M])[nodeIndex];
216 if (q>=c0o1 && q<=c1o1)
217 {
218 velocityLB = -vx3;
219 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
220 (dist.f[d00P])[kt] = getInterpolatedDistributionForNoSlipBC(q, f_B, f_T, feq, omega);
221 }
222
223 q = (subgridD.q[dPP0])[nodeIndex];
224 if (q>=c0o1 && q<=c1o1)
225 {
226 velocityLB = vx1 + vx2;
227 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
228 (dist.f[dMM0])[ksw] = getInterpolatedDistributionForNoSlipBC(q, f_NE, f_SW, feq, omega);
229 }
230
231 q = (subgridD.q[dMM0])[nodeIndex];
232 if (q>=c0o1 && q<=c1o1)
233 {
234 velocityLB = -vx1 - vx2;
235 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
236 (dist.f[dPP0])[kne] = getInterpolatedDistributionForNoSlipBC(q, f_SW, f_NE, feq, omega);
237 }
238
239 q = (subgridD.q[dPM0])[nodeIndex];
240 if (q>=c0o1 && q<=c1o1)
241 {
242 velocityLB = vx1 - vx2;
243 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
244 (dist.f[dMP0])[knw] = getInterpolatedDistributionForNoSlipBC(q, f_SE, f_NW, feq, omega);
245 }
246
247 q = (subgridD.q[dMP0])[nodeIndex];
248 if (q>=c0o1 && q<=c1o1)
249 {
250 velocityLB = -vx1 + vx2;
251 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
252 (dist.f[dPM0])[kse] = getInterpolatedDistributionForNoSlipBC(q, f_NW, f_SE, feq, omega);
253 }
254
255 q = (subgridD.q[dP0P])[nodeIndex];
256 if (q>=c0o1 && q<=c1o1)
257 {
258 velocityLB = vx1 + vx3;
259 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
260 (dist.f[dM0M])[kbw] = getInterpolatedDistributionForNoSlipBC(q, f_TE, f_BW, feq, omega);
261 }
262
263 q = (subgridD.q[dM0M])[nodeIndex];
264 if (q>=c0o1 && q<=c1o1)
265 {
266 velocityLB = -vx1 - vx3;
267 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
268 (dist.f[dP0P])[kte] = getInterpolatedDistributionForNoSlipBC(q, f_BW, f_TE, feq, omega);
269 }
270
271 q = (subgridD.q[dP0M])[nodeIndex];
272 if (q>=c0o1 && q<=c1o1)
273 {
274 velocityLB = vx1 - vx3;
275 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
276 (dist.f[dM0P])[ktw] = getInterpolatedDistributionForNoSlipBC(q, f_BE, f_TW, feq, omega);
277 }
278
279 q = (subgridD.q[dM0P])[nodeIndex];
280 if (q>=c0o1 && q<=c1o1)
281 {
282 velocityLB = -vx1 + vx3;
283 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
284 (dist.f[dP0M])[kbe] = getInterpolatedDistributionForNoSlipBC(q, f_TW, f_BE, feq, omega);
285 }
286
287 q = (subgridD.q[d0PP])[nodeIndex];
288 if (q>=c0o1 && q<=c1o1)
289 {
290 velocityLB = vx2 + vx3;
291 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
292 (dist.f[d0MM])[kbs] = getInterpolatedDistributionForNoSlipBC(q, f_TN, f_BS, feq, omega);
293 }
294
295 q = (subgridD.q[d0MM])[nodeIndex];
296 if (q>=c0o1 && q<=c1o1)
297 {
298 velocityLB = -vx2 - vx3;
299 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
300 (dist.f[d0PP])[ktn] = getInterpolatedDistributionForNoSlipBC(q, f_BS, f_TN, feq, omega);
301 }
302
303 q = (subgridD.q[d0PM])[nodeIndex];
304 if (q>=c0o1 && q<=c1o1)
305 {
306 velocityLB = vx2 - vx3;
307 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
308 (dist.f[d0MP])[kts] = getInterpolatedDistributionForNoSlipBC(q, f_BN, f_TS, feq, omega);
309 }
310
311 q = (subgridD.q[d0MP])[nodeIndex];
312 if (q>=c0o1 && q<=c1o1)
313 {
314 velocityLB = -vx2 + vx3;
315 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
316 (dist.f[d0PM])[kbn] = getInterpolatedDistributionForNoSlipBC(q, f_TS, f_BN, feq, omega);
317 }
318
319 q = (subgridD.q[dPPP])[nodeIndex];
320 if (q>=c0o1 && q<=c1o1)
321 {
322 velocityLB = vx1 + vx2 + vx3;
323 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
324 (dist.f[dMMM])[kbsw] = getInterpolatedDistributionForNoSlipBC(q, f_TNE, f_BSW, feq, omega);
325 }
326
327 q = (subgridD.q[dMMM])[nodeIndex];
328 if (q>=c0o1 && q<=c1o1)
329 {
330 velocityLB = -vx1 - vx2 - vx3;
331 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
332 (dist.f[dPPP])[ktne] = getInterpolatedDistributionForNoSlipBC(q, f_BSW, f_TNE, feq, omega);
333 }
334
335 q = (subgridD.q[dPPM])[nodeIndex];
336 if (q>=c0o1 && q<=c1o1)
337 {
338 velocityLB = vx1 + vx2 - vx3;
339 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
340 (dist.f[dMMP])[ktsw] = getInterpolatedDistributionForNoSlipBC(q, f_BNE, f_TSW, feq, omega);
341 }
342
343 q = (subgridD.q[dMMP])[nodeIndex];
344 if (q>=c0o1 && q<=c1o1)
345 {
346 velocityLB = -vx1 - vx2 + vx3;
347 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
348 (dist.f[dPPM])[kbne] = getInterpolatedDistributionForNoSlipBC(q, f_TSW, f_BNE, feq, omega);
349 }
350
351 q = (subgridD.q[dPMP])[nodeIndex];
352 if (q>=c0o1 && q<=c1o1)
353 {
354 velocityLB = vx1 - vx2 + vx3;
355 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
356 (dist.f[dMPM])[kbnw] = getInterpolatedDistributionForNoSlipBC(q, f_TSE, f_BNW, feq, omega);
357 }
358
359 q = (subgridD.q[dMPM])[nodeIndex];
360 if (q>=c0o1 && q<=c1o1)
361 {
362 velocityLB = -vx1 + vx2 - vx3;
363 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
364 (dist.f[dPMP])[ktse] = getInterpolatedDistributionForNoSlipBC(q, f_BNW, f_TSE, feq, omega);
365 }
366
367 q = (subgridD.q[dPMM])[nodeIndex];
368 if (q>=c0o1 && q<=c1o1)
369 {
370 velocityLB = vx1 - vx2 - vx3;
371 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
372 (dist.f[dMPP])[ktnw] = getInterpolatedDistributionForNoSlipBC(q, f_BSE, f_TNW, feq, omega);
373 }
374
375 q = (subgridD.q[dMPP])[nodeIndex];
376 if (q>=c0o1 && q<=c1o1)
377 {
378 velocityLB = -vx1 + vx2 + vx3;
379 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
380 (dist.f[dPMM])[kbse] = getInterpolatedDistributionForNoSlipBC(q, f_TNW, f_BSE, feq, omega);
381 }
382 }
383}
384////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
385
386}
387
388//! \}