VirtualFluids 0.2.0
Parallel CFD LBM Solver
Loading...
Searching...
No Matches
SurfaceLayer_Device.cuh
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//! \author Henry Korb
30//=======================================================================================
31#ifndef SurfaceLayer_Device_H
32#define SurfaceLayer_Device_H
33
34#include <cmath>
35
36#include <basics/DataTypes.h>
37#include <basics/constants/NumericConstants.h>
38
39#include <cuda_helper/CudaIndexCalculation.h>
40
41#include <lbm/MacroscopicQuantities.h>
42#include <lbm/advectionDiffusion/BoundaryConditions.h>
43#include <lbm/collision/TurbulentViscosity.h>
44#include <lbm/constants/D3Q27.h>
45
46#include "BoundaryConditions/BoundaryConditionFactory.h"
47#include "Calculation/Calculation.h"
48#include "Stress.h"
49#include "SurfaceLayer.h"
50#include "Utilities/KernelUtilities.h"
51#include "inverseMomentumExchange.cuh"
52#include "wallModelMoninObukhov.h"
53
54namespace vf::gpu {
55
56template <BoundaryConditionFactory::StressBC stressBCType, BoundaryConditionFactory::SurfaceLayerBC heatFluxBCtype,
57 bool useDelayedBounceBack>
58__global__ void
59SurfaceLayerDevice27(GridParameter gridParams, QforBoundaryConditions boundaryParams, WallModelParameters wallModelParams,
60 TemperatureWallModelParameters temperatureWallModelParams, TemperatureParameters temperatureParams)
61{
62 using namespace vf::basics::constant;
63 using namespace vf::lbm::dir;
64 using namespace vf::lbm::advection_diffusion;
65
66 constexpr real filterFrequency = 1e-3F;
67 constexpr uint maxIter = 100;
68 constexpr real convergenceCriteria = 1e-3F;
69 constexpr real zero = c0o1; // need for some std functions
70
71 using StressBC = BoundaryConditionFactory::StressBC;
72 using SurfaceLayerBC = BoundaryConditionFactory::SurfaceLayerBC;
73
74 const uint nodeIndex = vf::cuda::get1DIndexFrom2DBlock();
75
76 if (nodeIndex >= boundaryParams.numberOfBCnodes)
77 return;
78
79 ///////////////////////////////////////////////////////////
80 // Load and compute momentum inputs
81 /////////////////////////////////////////////////////////
82
83 Distributions27 populationReferences =
84 getDistributionReferences27(gridParams.distributions, gridParams.numberOfNodes, gridParams.isEvenTimestep);
85 SubgridDistances27 subgridDistances;
86 getPointersToSubgridDistances(subgridDistances, boundaryParams.q27[0], boundaryParams.numberOfBCnodes);
87
88 const uint k_000 = boundaryParams.k[nodeIndex];
89 const ListIndices listIndices(k_000, gridParams.neighborX, gridParams.neighborY, gridParams.neighborZ);
90
91 real populations[NUMBER_Of_DIRECTIONS];
92 getPostCollisionDistribution(populations, populationReferences, listIndices);
93
94 const real drho = vf::lbm::getDensity(populations);
95 const real3 velocityNode = { vf::lbm::getCompressibleVelocityX1(populations, drho),
96 vf::lbm::getCompressibleVelocityX2(populations, drho),
97 vf::lbm::getCompressibleVelocityX3(populations, drho) };
98
99 const real density = c1o1;
100
101 const real3 wallNormal { boundaryParams.normalX[nodeIndex], boundaryParams.normalY[nodeIndex],
102 boundaryParams.normalZ[nodeIndex] };
103 const real3 velocityNodeTangential = computeTangentialVector(velocityNode, wallNormal);
104
105 const real velocityNodeMeanTangentialMagnitude = smoothAndSaveMean(
106 computeMagnitude(velocityNodeTangential), filterFrequency, wallModelParams.velocityMagnitudeNode[nodeIndex]);
107
108 const uint samplingIndex = wallModelParams.samplingIndices[nodeIndex];
109
110 const real3 velocitySample = { gridParams.velocityX[samplingIndex], gridParams.velocityY[samplingIndex],
111 gridParams.velocityZ[samplingIndex] };
112
113 const real velocitySampleTangentialMagnitude = computeMagnitude(computeTangentialVector(velocitySample, wallNormal));
114
115 const real velocitySampleMeanTangentialMagnitude = smoothAndSaveMean(velocitySampleTangentialMagnitude, filterFrequency,
116 wallModelParams.velocityMagnitudeSample[nodeIndex]);
117
118 ///////////////////////////////////////////////////////////
119 // Load and compute temperature inputs
120 ///////////////////////////////////////////////////////////
121
122 auto populationReferencesTemperature = getDistributionReferences27(
123 temperatureParams.distributionsTemperature, gridParams.numberOfNodes, gridParams.isEvenTimestep);
124 real populationsTemperature[NUMBER_Of_DIRECTIONS];
125 getPostCollisionDistribution(populationsTemperature, populationReferencesTemperature, listIndices);
126
127 const real temperatureRelativeNode = vf::lbm::getDensity(populationsTemperature);
128 temperatureWallModelParams.temperatureNode[nodeIndex] = temperatureRelativeNode;
129
130 const real temperatureRelativeSample = temperatureParams.temperature[samplingIndex];
131 const real temperatureRelativeSampleMean = smoothAndSaveMean(temperatureRelativeSample, filterFrequency,
132 temperatureWallModelParams.temperatureSample[nodeIndex]);
133
134 ///////////////////////////////////////////////////////////
135 // Load wall model parameters
136 ///////////////////////////////////////////////////////////
137
138 const real samplingDistance = wallModelParams.samplingDistance[nodeIndex];
139 const real roughnessLength = wallModelParams.roughnessLength[nodeIndex];
140 const real vonKarmanConstant = wallModelParams.vonKarmanConstant[nodeIndex];
141 const real roughnessLengthTemperature = temperatureWallModelParams.roughnessLength[nodeIndex];
142
143 const real heatingRate = temperatureWallModelParams.heatingRate[nodeIndex];
144 const real surfaceTemperature = temperatureWallModelParams.surfaceTemperature[nodeIndex] + heatingRate;
145 temperatureWallModelParams.surfaceTemperature[nodeIndex] = surfaceTemperature;
146 const real temperatureDifference = temperatureRelativeSampleMean - surfaceTemperature;
147
148 ///////////////////////////////////////////////////////////
149 // Compute stress and heat flux from wall model
150 ///////////////////////////////////////////////////////////
151 real frictionVelocity = wallModelParams.frictionVelocity[nodeIndex];
152 real surfaceHeatFlux = temperatureWallModelParams.surfaceHeatFlux[nodeIndex];
153
154 uint iteration = 0;
155 bool converged = false;
156 do {
157 const real stabilityParameter =
158 computeStabilityParameter(samplingDistance, temperatureParams.gravity, surfaceHeatFlux, frictionVelocity,
159 temperatureParams.referenceTemperature, vonKarmanConstant);
160 const real stabilityCorrection = computeStabilityCorrectionMomentum(stabilityParameter);
161 const real frictionVelocityOld = frictionVelocity;
162 frictionVelocity = std::max(zero, computeFrictionVelocity(velocitySampleMeanTangentialMagnitude, vonKarmanConstant,
163 samplingDistance, roughnessLength, stabilityCorrection));
164 if (heatFluxBCtype == SurfaceLayerBC::SurfaceTemperature)
165 surfaceHeatFlux = computeSurfaceHeatFlux(temperatureDifference, frictionVelocity, vonKarmanConstant,
166 samplingDistance, roughnessLengthTemperature,
167 computeStabilityCorrectionTemperature(stabilityParameter));
168 iteration++;
169 converged = std::abs(frictionVelocity - frictionVelocityOld) < convergenceCriteria * std::abs(frictionVelocityOld);
170 } while ((iteration < maxIter) && !converged);
171
172 wallModelParams.frictionVelocity[nodeIndex] = frictionVelocity;
173 if (heatFluxBCtype == SurfaceLayerBC::SurfaceTemperature)
174 temperatureWallModelParams.surfaceHeatFlux[nodeIndex] = surfaceHeatFlux;
175
176 const real3 wallShearStress =
177 computeWallShearStress(frictionVelocity, velocityNodeTangential, velocityNodeMeanTangentialMagnitude, density);
178
179 ///////////////////////////////////////////////////////////
180 // Apply inverse Momentum Exchange
181 ///////////////////////////////////////////////////////////
182
183 real populationsBouncedBack[NUMBER_Of_DIRECTIONS];
184 bool linkIsCut[NUMBER_Of_DIRECTIONS];
185
186 switch (stressBCType) {
187 case StressBC::StressBounceBackCompressible:
188 computeBouncedBackDistributionsBB(subgridDistances, populations, linkIsCut, populationsBouncedBack, nodeIndex);
189 break;
190 case StressBC::StressBounceBackWithPressureCompressible:
191 computeBouncedBackDistributionsBBPressure(subgridDistances, drho, populations, linkIsCut, populationsBouncedBack,
192 nodeIndex);
193 break;
194 case StressBC::StressInterpolatedCompressible: {
195 const real relaxationFrequency = vf::lbm::calculateOmegaWithTurbulentViscosity(
196 gridParams.relaxationFrequency, gridParams.turbulentViscosity[k_000]);
197 computeBouncedBackDistributionsInterpolated(subgridDistances, velocityNode, drho, relaxationFrequency,
198 populations, linkIsCut, populationsBouncedBack, nodeIndex);
199 } break;
200 }
201
202 real3 wallMomentum = computeWallMomentumBounceBack(linkIsCut, populationsBouncedBack, populations);
203
204 const real wallArea = c1o1;
205 const real subgridDistance = (subgridDistances.q[d00M])[nodeIndex];
206 const real interpolationFactor =
207 stressBCType == StressBC::StressInterpolatedCompressible ? c1o1 + subgridDistance : c1o1;
208
209 const real3 fakeWallVelocity = computeFakeWallVelocity(wallNormal, velocitySample, wallShearStress, density,
210 interpolationFactor, wallArea, wallMomentum);
211 if (!useDelayedBounceBack)
212 populationReferences = getDistributionReferences27(gridParams.distributions, gridParams.numberOfNodes,
213 !gridParams.isEvenTimestep);
214 wallMomentum += writeDistributionsBB(populationReferences, linkIsCut, populationsBouncedBack, fakeWallVelocity, density,
215 listIndices);
216
217 wallModelParams.forceX[nodeIndex] = wallMomentum.x;
218 wallModelParams.forceY[nodeIndex] = wallMomentum.y;
219 wallModelParams.forceZ[nodeIndex] = wallMomentum.z;
220
221 ///////////////////////////////////////////////////////////
222 // Apply Heat Flux Boundary Condition
223 ///////////////////////////////////////////////////////////
224
225 const real3 diffusiveFlux = real3 { vf::lbm::getIncompressibleVelocityX1(populationsTemperature),
226 vf::lbm::getIncompressibleVelocityX2(populationsTemperature),
227 vf::lbm::getIncompressibleVelocityX3(populationsTemperature) } -
228 velocityNode * temperatureRelativeNode;
229 const real normalDiffusiveFlux = dot(diffusiveFlux, wallNormal);
230 const real3 wallFlux = diffusiveFlux + wallNormal * (surfaceHeatFlux - normalDiffusiveFlux);
231
232 if (!useDelayedBounceBack)
233 populationReferencesTemperature = getDistributionReferences27(
234 temperatureParams.distributionsTemperature, gridParams.numberOfNodes, !gridParams.isEvenTimestep);
235
236 switch (stressBCType) {
237 case StressBC::StressBounceBackCompressible:
238 case StressBC::StressBounceBackWithPressureCompressible:
239 forEachNonRestDirection([&](auto direction) {
240 if (!linkIsCut[direction])
241 return;
242 const real population = computePopulationSimpleBounceBackWithFlux<direction>(
243 populationsTemperature, wallFlux.x, wallFlux.y, wallFlux.z);
244 writeInInverseDirection<direction>(population, listIndices, populationReferencesTemperature);
245 });
246 break;
247 case StressBC::StressInterpolatedCompressible:
248 const real diffusivityNode = temperatureParams.diffusivity + temperatureParams.turbulentDiffusivity[k_000];
249 const real relaxationFrequency = vf::lbm::computeRelaxationFrequency(diffusivityNode);
250 forEachNonRestDirection([&](auto direction) {
251 const real subgridDistance = (subgridDistances.q[direction])[nodeIndex];
252 if (subgridDistance < c0o1 || subgridDistance > c1o1)
253 return;
254 const real population = computePopulationInterpolatedBounceBackWithFlux<direction>(
255 subgridDistance, populationsTemperature, velocityNode.x, velocityNode.y, velocityNode.z,
256 relaxationFrequency, temperatureRelativeNode, wallFlux.x, wallFlux.y, wallFlux.z);
257 writeInInverseDirection<direction>(population, listIndices, populationReferencesTemperature);
258 });
259 break;
260 }
261}
262
263}
264
265#endif