1//=======================================================================================
2// ____ ____ __ ______ __________ __ __ __ __
3// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
4// \ \ | | | | | |_) | | | | | | | / \ | |
5// \ \ | | | | | _ / | | | | | | / /\ \ | |
6// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
7// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
8// \ \ | | ________________________________________________________________
9// \ \ | | | ______________________________________________________________|
10// \ \| | | | __ __ __ __ ______ _______
11// \ | | |_____ | | | | | | | | | _ \ / _____)
12// \ | | _____| | | | | | | | | | | \ \ \_______
13// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
14// \ _____| |__| |________| \_______/ |__| |______/ (_______/
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.
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
26// SPDX-License-Identifier: GPL-3.0-or-later
27// SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
29//! \addtogroup gpu_BoundaryConditions BoundaryConditions
30//! \ingroup gpu_core core
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"
40using namespace vf::basics::constant;
41using namespace vf::lbm::dir;
45//////////////////////////////////////////////////////////////////////////////
46__global__ void SlipCompressible_Device(
48 int* subgridDistanceIndices,
49 real* subgridDistances,
50 unsigned int numberOfBCnodes,
52 unsigned int* neighborX,
53 unsigned int* neighborY,
54 unsigned int* neighborZ,
55 unsigned long long numberOfLBnodes,
58 //! The slip boundary condition is executed in the following steps
61 ////////////////////////////////////////////////////////////////////////////////
62 //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
64 const unsigned nodeIndex = vf::cuda::get1DIndexFrom2DBlock();
66 if(nodeIndex < numberOfBCnodes)
68 //////////////////////////////////////////////////////////////////////////
69 //! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on timestep is based on the esoteric twist algorithm \ref
70 //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
73 getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
75 ////////////////////////////////////////////////////////////////////////////////
76 //! - Set local subgrid distances (q's)
78 SubgridDistances27 subgridD;
79 getPointersToSubgridDistances(subgridD, subgridDistances, numberOfBCnodes);
81 ////////////////////////////////////////////////////////////////////////////////
82 //! - Set neighbor indices (necessary for indirect addressing)
84 unsigned int indexOfBCnode = subgridDistanceIndices[nodeIndex];
85 unsigned int kzero= indexOfBCnode;
86 unsigned int ke = indexOfBCnode;
87 unsigned int kw = neighborX[indexOfBCnode];
88 unsigned int kn = indexOfBCnode;
89 unsigned int ks = neighborY[indexOfBCnode];
90 unsigned int kt = indexOfBCnode;
91 unsigned int kb = neighborZ[indexOfBCnode];
92 unsigned int ksw = neighborY[kw];
93 unsigned int kne = indexOfBCnode;
94 unsigned int kse = ks;
95 unsigned int knw = kw;
96 unsigned int kbw = neighborZ[kw];
97 unsigned int kte = indexOfBCnode;
98 unsigned int kbe = kb;
99 unsigned int ktw = kw;
100 unsigned int kbs = neighborZ[ks];
101 unsigned int ktn = indexOfBCnode;
102 unsigned int kbn = kb;
103 unsigned int kts = ks;
104 unsigned int ktse = ks;
105 unsigned int kbnw = kbw;
106 unsigned int ktnw = kw;
107 unsigned int kbse = kbs;
108 unsigned int ktsw = ksw;
109 unsigned int kbne = kb;
110 unsigned int ktne = indexOfBCnode;
111 unsigned int kbsw = neighborZ[ksw];
113 ////////////////////////////////////////////////////////////////////////////////
114 //! - Set local distributions
116 real f_W = (dist.f[dP00])[ke ];
117 real f_E = (dist.f[dM00])[kw ];
118 real f_S = (dist.f[d0P0])[kn ];
119 real f_N = (dist.f[d0M0])[ks ];
120 real f_B = (dist.f[d00P])[kt ];
121 real f_T = (dist.f[d00M])[kb ];
122 real f_SW = (dist.f[dPP0])[kne ];
123 real f_NE = (dist.f[dMM0])[ksw ];
124 real f_NW = (dist.f[dPM0])[kse ];
125 real f_SE = (dist.f[dMP0])[knw ];
126 real f_BW = (dist.f[dP0P])[kte ];
127 real f_TE = (dist.f[dM0M])[kbw ];
128 real f_TW = (dist.f[dP0M])[kbe ];
129 real f_BE = (dist.f[dM0P])[ktw ];
130 real f_BS = (dist.f[d0PP])[ktn ];
131 real f_TN = (dist.f[d0MM])[kbs ];
132 real f_TS = (dist.f[d0PM])[kbn ];
133 real f_BN = (dist.f[d0MP])[kts ];
134 real f_BSW = (dist.f[dPPP])[ktne ];
135 real f_BNE = (dist.f[dMMP])[ktsw ];
136 real f_BNW = (dist.f[dPMP])[ktse ];
137 real f_BSE = (dist.f[dMPP])[ktnw ];
138 real f_TSW = (dist.f[dPPM])[kbne ];
139 real f_TNE = (dist.f[dMMM])[kbsw ];
140 real f_TNW = (dist.f[dPMM])[kbse ];
141 real f_TSE = (dist.f[dMPM])[kbnw ];
143 ////////////////////////////////////////////////////////////////////////////////
144 //! - Calculate macroscopic quantities
146 real drho = f_TSE + f_TNW + f_TNE + f_TSW + f_BSE + f_BNW + f_BNE + f_BSW +
147 f_BN + f_TS + f_TN + f_BS + f_BE + f_TW + f_TE + f_BW + f_SE + f_NW + f_NE + f_SW +
148 f_T + f_B + f_N + f_S + f_E + f_W + ((dist.f[d000])[kzero]);
150 real vx1 = (((f_TSE - f_BNW) - (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
151 ((f_BE - f_TW) + (f_TE - f_BW)) + ((f_SE - f_NW) + (f_NE - f_SW)) +
152 (f_E - f_W)) / (c1o1 + drho);
154 real vx2 = ((-(f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
155 ((f_BN - f_TS) + (f_TN - f_BS)) + (-(f_SE - f_NW) + (f_NE - f_SW)) +
156 (f_N - f_S)) / (c1o1 + drho);
158 real vx3 = (((f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) + (f_TSW - f_BNE)) +
159 (-(f_BN - f_TS) + (f_TN - f_BS)) + ((f_TE - f_BW) - (f_BE - f_TW)) +
160 (f_T - f_B)) / (c1o1 + drho);
162 real cu_sq = c3o2 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3) * (c1o1 + drho);
164 ////////////////////////////////////////////////////////////////////////////////
165 //! - change the pointer to write the results in the correct array
167 getPointersToDistributions(dist, distributions, numberOfLBnodes, !isEvenTimestep);
169 ////////////////////////////////////////////////////////////////////////////////
170 //! - Multiply the local velocities by the slipLength
172 real slipLength = c1o1;
173 real VeloX = slipLength*vx1;
174 real VeloY = slipLength*vx2;
175 real VeloZ = slipLength*vx3;
177 ////////////////////////////////////////////////////////////////////////////////
178 //! - Update distributions with subgrid distance (q) between zero and one
180 real feq, q, velocityLB, velocityBC;
186 q = (subgridD.q[dP00])[nodeIndex];
187 if (q>=c0o1 && q<=c1o1) // only update distribution for q between zero and one
193 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
195 (dist.f[dM00])[kw] = getInterpolatedDistributionForVeloBC(q, f_E, f_W, feq, omega, velocityBC, c2o27);
198 q = (subgridD.q[dM00])[nodeIndex];
199 if (q>=c0o1 && q<=c1o1)
205 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
207 (dist.f[dP00])[ke] = getInterpolatedDistributionForVeloBC(q, f_W, f_E, feq, omega, velocityBC, c2o27);
210 q = (subgridD.q[d0P0])[nodeIndex];
211 if (q>=c0o1 && q<=c1o1)
217 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
219 (dist.f[d0M0])[ks] = getInterpolatedDistributionForVeloBC(q, f_N, f_S, feq, omega, velocityBC, c2o27);
222 q = (subgridD.q[d0M0])[nodeIndex];
223 if (q>=c0o1 && q<=c1o1)
229 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
231 (dist.f[d0P0])[kn] = getInterpolatedDistributionForVeloBC(q, f_S, f_N, feq, omega, velocityBC, c2o27);
234 q = (subgridD.q[d00P])[nodeIndex];
235 if (q>=c0o1 && q<=c1o1)
241 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
243 (dist.f[d00M])[kb] = getInterpolatedDistributionForVeloBC(q, f_T, f_B, feq, omega, velocityBC, c2o27);
246 q = (subgridD.q[d00M])[nodeIndex];
247 if (q>=c0o1 && q<=c1o1)
253 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
255 (dist.f[d00P])[kt] = getInterpolatedDistributionForVeloBC(q, f_B, f_T, feq, omega, velocityBC, c2o27);
258 q = (subgridD.q[dPP0])[nodeIndex];
259 if (q>=c0o1 && q<=c1o1)
261 VeloX = slipLength*vx1;
262 VeloY = slipLength*vx2;
263 if (x == true) VeloX = c0o1;
264 if (y == true) VeloY = c0o1;
266 velocityLB = vx1 + vx2;
267 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
268 velocityBC = VeloX + VeloY;
269 (dist.f[dMM0])[ksw] = getInterpolatedDistributionForVeloBC(q, f_NE, f_SW, feq, omega, velocityBC, c1o54);
272 q = (subgridD.q[dMM0])[nodeIndex];
273 if (q>=c0o1 && q<=c1o1)
275 VeloX = slipLength*vx1;
276 VeloY = slipLength*vx2;
277 if (x == true) VeloX = c0o1;
278 if (y == true) VeloY = c0o1;
280 velocityLB = -vx1 - vx2;
281 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
282 velocityBC = -VeloX - VeloY;
283 (dist.f[dPP0])[kne] = getInterpolatedDistributionForVeloBC(q, f_SW, f_NE, feq, omega, velocityBC, c1o54);
286 q = (subgridD.q[dPM0])[nodeIndex];
287 if (q>=c0o1 && q<=c1o1)
289 VeloX = slipLength*vx1;
290 VeloY = slipLength*vx2;
291 if (x == true) VeloX = c0o1;
292 if (y == true) VeloY = c0o1;
294 velocityLB = vx1 - vx2;
295 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
296 velocityBC = VeloX - VeloY;
297 (dist.f[dMP0])[knw] = getInterpolatedDistributionForVeloBC(q, f_SE, f_NW, feq, omega, velocityBC, c1o54);
300 q = (subgridD.q[dMP0])[nodeIndex];
301 if (q>=c0o1 && q<=c1o1)
303 VeloX = slipLength*vx1;
304 VeloY = slipLength*vx2;
305 if (x == true) VeloX = c0o1;
306 if (y == true) VeloY = c0o1;
308 velocityLB = -vx1 + vx2;
309 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
310 velocityBC = -VeloX + VeloY;
311 (dist.f[dPM0])[kse] = getInterpolatedDistributionForVeloBC(q, f_NW, f_SE, feq, omega, velocityBC, c1o54);
314 q = (subgridD.q[dP0P])[nodeIndex];
315 if (q>=c0o1 && q<=c1o1)
317 VeloX = slipLength*vx1;
318 VeloZ = slipLength*vx3;
319 if (x == true) VeloX = c0o1;
320 if (z == true) VeloZ = c0o1;
322 velocityLB = vx1 + vx3;
323 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
324 velocityBC = VeloX + VeloZ;
325 (dist.f[dM0M])[kbw] = getInterpolatedDistributionForVeloBC(q, f_TE, f_BW, feq, omega, velocityBC, c1o54);
328 q = (subgridD.q[dM0M])[nodeIndex];
329 if (q>=c0o1 && q<=c1o1)
331 VeloX = slipLength*vx1;
332 VeloZ = slipLength*vx3;
333 if (x == true) VeloX = c0o1;
334 if (z == true) VeloZ = c0o1;
336 velocityLB = -vx1 - vx3;
337 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
338 velocityBC = -VeloX - VeloZ;
339 (dist.f[dP0P])[kte] = getInterpolatedDistributionForVeloBC(q, f_BW, f_TE, feq, omega, velocityBC, c1o54);
342 q = (subgridD.q[dP0M])[nodeIndex];
343 if (q>=c0o1 && q<=c1o1)
345 VeloX = slipLength*vx1;
346 VeloZ = slipLength*vx3;
347 if (x == true) VeloX = c0o1;
348 if (z == true) VeloZ = c0o1;
350 velocityLB = vx1 - vx3;
351 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
352 velocityBC = VeloX - VeloZ;
353 (dist.f[dM0P])[ktw] = getInterpolatedDistributionForVeloBC(q, f_BE, f_TW, feq, omega, velocityBC, c1o54);
356 q = (subgridD.q[dM0P])[nodeIndex];
357 if (q>=c0o1 && q<=c1o1)
359 VeloX = slipLength*vx1;
360 VeloZ = slipLength*vx3;
361 if (x == true) VeloX = c0o1;
362 if (z == true) VeloZ = c0o1;
364 velocityLB = -vx1 + vx3;
365 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
366 velocityBC = -VeloX + VeloZ;
367 (dist.f[dP0M])[kbe] = getInterpolatedDistributionForVeloBC(q, f_TW, f_BE, feq, omega, velocityBC, c1o54);
370 q = (subgridD.q[d0PP])[nodeIndex];
371 if (q>=c0o1 && q<=c1o1)
373 VeloY = slipLength*vx2;
374 VeloZ = slipLength*vx3;
375 if (y == true) VeloY = c0o1;
376 if (z == true) VeloZ = c0o1;
378 velocityLB = vx2 + vx3;
379 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
380 velocityBC = VeloY + VeloZ;
381 (dist.f[d0MM])[kbs] = getInterpolatedDistributionForVeloBC(q, f_TN, f_BS, feq, omega, velocityBC, c1o54);
384 q = (subgridD.q[d0MM])[nodeIndex];
385 if (q>=c0o1 && q<=c1o1)
387 VeloY = slipLength*vx2;
388 VeloZ = slipLength*vx3;
389 if (y == true) VeloY = c0o1;
390 if (z == true) VeloZ = c0o1;
392 velocityLB = -vx2 - vx3;
393 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
394 velocityBC = -VeloY - VeloZ;
395 (dist.f[d0PP])[ktn] = getInterpolatedDistributionForVeloBC(q, f_BS, f_TN, feq, omega, velocityBC, c1o54);
399 q = (subgridD.q[d0PM])[nodeIndex];
400 if (q>=c0o1 && q<=c1o1)
402 VeloY = slipLength*vx2;
403 VeloZ = slipLength*vx3;
404 if (y == true) VeloY = c0o1;
405 if (z == true) VeloZ = c0o1;
407 velocityLB = vx2 - vx3;
408 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
409 velocityBC = VeloY - VeloZ;
410 (dist.f[d0MP])[kts] = getInterpolatedDistributionForVeloBC(q, f_BN, f_TS, feq, omega, velocityBC, c1o54);
413 q = (subgridD.q[d0MP])[nodeIndex];
414 if (q>=c0o1 && q<=c1o1)
416 VeloY = slipLength*vx2;
417 VeloZ = slipLength*vx3;
418 if (y == true) VeloY = c0o1;
419 if (z == true) VeloZ = c0o1;
421 velocityLB = -vx2 + vx3;
422 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
423 velocityBC = -VeloY + VeloZ;
424 (dist.f[d0PM])[kbn] = getInterpolatedDistributionForVeloBC(q, f_TS, f_BN, feq, omega, velocityBC, c1o54);
427 q = (subgridD.q[dPPP])[nodeIndex];
428 if (q>=c0o1 && q<=c1o1)
430 VeloX = slipLength*vx1;
431 VeloY = slipLength*vx2;
432 VeloZ = slipLength*vx3;
433 if (x == true) VeloX = c0o1;
434 if (y == true) VeloY = c0o1;
435 if (z == true) VeloZ = c0o1;
436 velocityLB = vx1 + vx2 + vx3;
437 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
438 velocityBC = VeloX + VeloY + VeloZ;
439 (dist.f[dMMM])[kbsw] = getInterpolatedDistributionForVeloBC(q, f_TNE, f_BSW, feq, omega, velocityBC, c1o216);
442 q = (subgridD.q[dMMM])[nodeIndex];
443 if (q>=c0o1 && q<=c1o1)
445 VeloX = slipLength*vx1;
446 VeloY = slipLength*vx2;
447 VeloZ = slipLength*vx3;
448 if (x == true) VeloX = c0o1;
449 if (y == true) VeloY = c0o1;
450 if (z == true) VeloZ = c0o1;
451 velocityLB = -vx1 - vx2 - vx3;
452 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
453 velocityBC = -VeloX - VeloY - VeloZ;
454 (dist.f[dPPP])[ktne] = getInterpolatedDistributionForVeloBC(q, f_BSW, f_TNE, feq, omega, velocityBC, c1o216);
458 q = (subgridD.q[dPPM])[nodeIndex];
459 if (q>=c0o1 && q<=c1o1)
461 VeloX = slipLength*vx1;
462 VeloY = slipLength*vx2;
463 VeloZ = slipLength*vx3;
464 if (x == true) VeloX = c0o1;
465 if (y == true) VeloY = c0o1;
466 if (z == true) VeloZ = c0o1;
467 velocityLB = vx1 + vx2 - vx3;
468 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
469 velocityBC = VeloX + VeloY - VeloZ;
470 (dist.f[dMMP])[ktsw] = getInterpolatedDistributionForVeloBC(q, f_BNE, f_TSW, feq, omega, velocityBC, c1o216);
473 q = (subgridD.q[dMMP])[nodeIndex];
474 if (q>=c0o1 && q<=c1o1)
476 VeloX = slipLength*vx1;
477 VeloY = slipLength*vx2;
478 VeloZ = slipLength*vx3;
479 if (x == true) VeloX = c0o1;
480 if (y == true) VeloY = c0o1;
481 if (z == true) VeloZ = c0o1;
482 velocityLB = -vx1 - vx2 + vx3;
483 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
484 velocityBC = -VeloX - VeloY + VeloZ;
485 (dist.f[dPPM])[kbne] = getInterpolatedDistributionForVeloBC(q, f_TSW, f_BNE, feq, omega, velocityBC, c1o216);
488 q = (subgridD.q[dPMP])[nodeIndex];
489 if (q>=c0o1 && q<=c1o1)
491 VeloX = slipLength*vx1;
492 VeloY = slipLength*vx2;
493 VeloZ = slipLength*vx3;
494 if (x == true) VeloX = c0o1;
495 if (y == true) VeloY = c0o1;
496 if (z == true) VeloZ = c0o1;
497 velocityLB = vx1 - vx2 + vx3;
498 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
499 velocityBC = VeloX - VeloY + VeloZ;
500 (dist.f[dMPM])[kbnw] = getInterpolatedDistributionForVeloBC(q, f_TSE, f_BNW, feq, omega, velocityBC, c1o216);
503 q = (subgridD.q[dMPM])[nodeIndex];
504 if (q>=c0o1 && q<=c1o1)
506 VeloX = slipLength*vx1;
507 VeloY = slipLength*vx2;
508 VeloZ = slipLength*vx3;
509 if (x == true) VeloX = c0o1;
510 if (y == true) VeloY = c0o1;
511 if (z == true) VeloZ = c0o1;
512 velocityLB = -vx1 + vx2 - vx3;
513 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
514 velocityBC = -VeloX + VeloY - VeloZ;
515 (dist.f[dPMP])[ktse] = getInterpolatedDistributionForVeloBC(q, f_BNW, f_TSE, feq, omega, velocityBC, c1o216);
518 q = (subgridD.q[dPMM])[nodeIndex];
519 if (q>=c0o1 && q<=c1o1)
521 VeloX = slipLength*vx1;
522 VeloY = slipLength*vx2;
523 VeloZ = slipLength*vx3;
524 if (x == true) VeloX = c0o1;
525 if (y == true) VeloY = c0o1;
526 if (z == true) VeloZ = c0o1;
527 velocityLB = vx1 - vx2 - vx3;
528 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
529 velocityBC = VeloX - VeloY - VeloZ;
530 (dist.f[dMPP])[ktnw] = getInterpolatedDistributionForVeloBC(q, f_BSE, f_TNW, feq, omega, velocityBC, c1o216);
533 q = (subgridD.q[dMPP])[nodeIndex];
534 if (q>=c0o1 && q<=c1o1)
536 VeloX = slipLength*vx1;
537 VeloY = slipLength*vx2;
538 VeloZ = slipLength*vx3;
539 if (x == true) VeloX = c0o1;
540 if (y == true) VeloY = c0o1;
541 if (z == true) VeloZ = c0o1;
542 velocityLB = -vx1 + vx2 + vx3;
543 feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
544 velocityBC = -VeloX + VeloY + VeloZ;
545 (dist.f[dPMM])[kbse] = getInterpolatedDistributionForVeloBC(q, f_TNW, f_BSE, feq, omega, velocityBC, c1o216);
549////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////