OpenLB 1.7
Loading...
Searching...
No Matches
blockLatticeInteraction.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2006-2008 Jonas Latt
4 * 2008-2020 Mathias Krause
5 * 2020 Adrian Kummerlaender
6 * 2021 Nicolas Hafen
7 * E-mail contact: info@openlb.net
8 * The most recent release of OpenLB can be downloaded at
9 * <http://www.openlb.net/>
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public
22 * License along with this program; if not, write to the Free
23 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
24 * Boston, MA 02110-1301, USA.
25*/
26
27#ifndef BLOCK_LATTICE_INTERACTION_HH
28#define BLOCK_LATTICE_INTERACTION_HH
29
31// TODO: Use with c++20
32//#include <numbers>
33
34namespace olb {
35
36namespace particles {
37
38template <typename T, unsigned D>
40 T invDeltaX, LatticeR<D>& start,
41 LatticeR<D>& end, PhysR<T, D> position,
42 T circumRadius)
43{
45 LatticeR<D> blockMin(0);
46 LatticeR<D> blockMax(blockGeometry.getExtent() - 1);
47
49 PhysR<T, D> particleMin(position - circumRadius);
50 PhysR<T, D> relStart(particleMin - blockGeometry.getOrigin());
51 PhysR<T, D> particleMax(position + circumRadius);
52 PhysR<T, D> relEnd(particleMax - blockGeometry.getOrigin());
53
55 bool validRange = true;
56 for (unsigned iDim = 0; iDim < D; ++iDim) {
57 start[iDim] =
58 util::max(util::floor(invDeltaX * relStart[iDim]), blockMin[iDim]);
59 end[iDim] = util::min(util::ceil(invDeltaX * relEnd[iDim]), blockMax[iDim]);
60 int intersectionRange = end[iDim] - start[iDim];
61 validRange = validRange && (intersectionRange >= 0.);
62 }
63
64 return validRange;
65}
66
67template <typename T, unsigned D>
69 bool& surfaceOutOfGeometry, PhysR<T, D>& ghostPos,
70 const PhysR<T, D>& cellMin, const PhysR<T, D>& cellMax,
71 const PhysR<T, D>& position, T circumRadius,
72 const Vector<bool, D>& periodic)
73{
74 ghostPos = position;
75 surfaceOutOfGeometry = false;
76 for (unsigned i = 0; i < D; ++i) {
77 T const particleMin = position[i] - circumRadius;
78 T const particleMax = position[i] + circumRadius;
79 if (particleMin < cellMin[i] && periodic[i]) {
80 surfaceOutOfGeometry = true;
82 position[i], cellMax[i], cellMin[i]);
83 }
84 else if (particleMax > cellMax[i] && periodic[i]) {
85 surfaceOutOfGeometry = true;
87 position[i], cellMax[i], cellMin[i]);
88 }
89 }
90}
91
94template <typename T, typename DESCRIPTOR, typename PARTICLETYPE>
96 BlockLattice<T, DESCRIPTOR>& blockLattice,
97 UnitConverter<T, DESCRIPTOR> const& converter,
98 const PhysR<T, DESCRIPTOR::d>& cellMin,
99 const PhysR<T, DESCRIPTOR::d>& cellMax,
101 const Vector<bool, DESCRIPTOR::d>& periodic)
102{
103 constexpr unsigned D = DESCRIPTOR::d;
104 using namespace descriptors;
105 const T circumRadius = particles::access::getRadius(particle);
106 const PhysR<T, D> originalPosition(particles::access::getPosition(particle));
107 // Checking if the particle's surface leaves the domain
108 bool surfaceOutOfGeometry;
109 PhysR<T, D> ghostPos;
110 checkSmoothIndicatorOutOfGeometry(surfaceOutOfGeometry, ghostPos, cellMin,
111 cellMax, originalPosition, circumRadius,
112 periodic);
113
114 //Do the normal routine if the particle is in the geometry
115 if (!surfaceOutOfGeometry) {
116 setBlockParticleField(blockGeometry, blockLattice, converter, particle);
117 }
118 else {
119 //sets the Particle to ghost position on the other side of the domain and sets the field
120 particle.template setField<GENERAL, POSITION>(ghostPos);
121 setBlockParticleField(blockGeometry, blockLattice, converter, particle);
122 //Reverting Particle to its Previous position and setting the field
123 particle.template setField<GENERAL, POSITION>(originalPosition);
124 setBlockParticleField(blockGeometry, blockLattice, converter, particle);
125
126 // For parallized particles, the update of the position is processed during relocation
127 if constexpr(!particles::access::providesParallelization<PARTICLETYPE>()) {
128 // Check if center of mass left the simulation domain
129 bool centerOfMassOutOfDomain = false;
130 for (unsigned iDim = 0; iDim < D; ++iDim) {
131 centerOfMassOutOfDomain = centerOfMassOutOfDomain ||
132 originalPosition[iDim] < cellMin[iDim] ||
133 originalPosition[iDim] > cellMax[iDim];
134 }
135
136 // Replace the position of the center of mass with the ghost position
137 // (If the center of mass is outside the domain)
138 if (centerOfMassOutOfDomain) {
139 particle.template setField<GENERAL, POSITION>(ghostPos);
140 }
141 }
142 }
143}
144
145template <typename T, typename DESCRIPTOR, typename PARTICLETYPE,
146 typename PARTICLECONTACTTYPE, typename WALLCONTACTTYPE, typename F>
148 const BlockGeometry<T, DESCRIPTOR::d>& blockGeometry,
149 BlockLattice<T, DESCRIPTOR>& blockLattice,
150 UnitConverter<T, DESCRIPTOR> const& converter,
151 const PhysR<T, DESCRIPTOR::d>& cellMin,
152 const PhysR<T, DESCRIPTOR::d>& cellMax,
153 ParticleSystem<T, PARTICLETYPE>& particleSystem,
155 contactContainer,
156 size_t iP, Particle<T, PARTICLETYPE>& particle,
157 std::vector<SolidBoundary<T, DESCRIPTOR::d>>& solidBoundaries,
158 F getSetupPeriodicity)
159{
160 constexpr unsigned D = DESCRIPTOR::d;
161 using namespace descriptors;
162 const T circumRadius = particles::access::getRadius(particle);
163 // Checking if the Particle leaves the Domain
164 bool surfaceOutOfGeometry;
165 PhysR<T, D> originalPosition(particles::access::getPosition(particle));
166 PhysR<T, D> ghostPos;
168 surfaceOutOfGeometry, ghostPos, cellMin, cellMax, originalPosition,
169 circumRadius, getSetupPeriodicity());
170
171 //Do the normal routine if the particle is in the geometry
172 if (!surfaceOutOfGeometry) {
173 setBlockParticleField(blockGeometry, blockLattice, converter,
174 particleSystem, contactContainer, iP, particle,
175 solidBoundaries);
176 }
177 else {
178 //sets the Particle to ghost position on the other side of the domain and sets the field
179 particle.template setField<GENERAL, POSITION>(ghostPos);
180 setBlockParticleField(blockGeometry, blockLattice, converter,
181 particleSystem, contactContainer, iP, particle,
182 solidBoundaries, cellMin, cellMax,
183 getSetupPeriodicity);
184 //Reverting Particle to its Previous position and setting the field
185 particle.template setField<GENERAL, POSITION>(originalPosition);
186 setBlockParticleField(blockGeometry, blockLattice, converter,
187 particleSystem, contactContainer, iP, particle,
188 solidBoundaries, cellMin, cellMax,
189 getSetupPeriodicity);
190 }
191}
192
193//Interation for block particle intersection for provided lambda function F
194template <typename T, typename DESCRIPTOR, typename F>
196 const BlockGeometry<T, DESCRIPTOR::d>& blockGeometry,
197 BlockLattice<T, DESCRIPTOR>& blockLattice, int padding,
198 Vector<T, DESCRIPTOR::d> position, T circumRadius, F f)
199{
200 constexpr unsigned D = DESCRIPTOR::d;
201 LatticeR<D> start;
202 LatticeR<D> end;
203 if (getBlockParticleIntersection(blockGeometry,
204 T {1} / blockGeometry.getDeltaR(), start,
205 end, position, circumRadius)) {
206 //Extend range by padding if desired
207 start -= LatticeR<D>(padding);
208 end += LatticeR<D>(padding);
209 //For all cells in subdomain defined by start/end
210 blockLattice.forSpatialLocations(start, end,
211 [&](const LatticeR<D>& latticeR) {
212 f(latticeR);
213 });
214 }
215}
216
217//TODO: remove material 1 from hardcoded version
218//- Can be done with unordered_set, analogous to other version below
219template <typename T, typename DESCRIPTOR, typename PARTICLETYPE>
221 BlockLattice<T, DESCRIPTOR>& blockLattice,
222 UnitConverter<T, DESCRIPTOR> const& converter,
224{
225 constexpr unsigned D = DESCRIPTOR::d;
226
227 using namespace descriptors;
228 auto circumRadius =
229 particle.template getField<SURFACE, SINDICATOR>()->getCircumRadius();
230 auto position = particle.template getField<GENERAL, POSITION>();
231
232 //For all cells in block particle intersection
233 int padding = blockGeometry.getPadding();
235 blockGeometry, blockLattice, padding, position, circumRadius,
236 [&](const LatticeR<D>& latticeR) {
237 //Get phys position
238 T physR[D] = {};
239 blockGeometry.getPhysR(physR, latticeR);
240 //Retrieve porosity
242 particles::resolved::evalSolidVolumeFraction(porosity.data(), physR, particle);
243
244 //For cells containing particle bits
245 if (!util::nearZero(porosity[0])) {
246 auto material = blockGeometry.get(latticeR);
247 //TODO: remove material 1 from hardcoded version
248 if (material == 1) {
249 //Retrieve local velocity
252 particle, PhysR<T, D>(physR));
253 for (unsigned iDim = 0; iDim != D; ++iDim) {
254 velocity[iDim] = converter.getLatticeVelocity(velocity[iDim]);
255 }
256 //Calculate solid volume fraction
257 T solidVolumeFraction = 1. - porosity[0];
258 //Write to fields
259 {
260 auto cell = blockLattice.get(latticeR);
261 if constexpr (!DESCRIPTOR::template provides<descriptors::VELOCITY_SOLID>()) {
262 cell.template getFieldPointer<
263 descriptors::VELOCITY_NUMERATOR>() += porosity[0] * velocity;
264 cell.template getFieldPointer<
266 }
267 else {
268 cell.template getFieldPointer<
269 descriptors::VELOCITY_SOLID>() += velocity;
270 }
271 cell.template getFieldPointer<descriptors::POROSITY>() *=
272 solidVolumeFraction;
273 }
274 } //if (material==1)
275 } //if (!util::nearZero(porosity[0]))
276 });
277}
278
279template <typename T, typename DESCRIPTOR, typename PARTICLETYPE,
280 typename PARTICLECONTACTTYPE, typename WALLCONTACTTYPE, typename F>
282 const BlockGeometry<T, DESCRIPTOR::d>& blockGeometry,
283 BlockLattice<T, DESCRIPTOR>& blockLattice,
284 UnitConverter<T, DESCRIPTOR> const& converter,
285 ParticleSystem<T, PARTICLETYPE>& particleSystem,
287 contactContainer,
288 size_t iP, Particle<T, PARTICLETYPE>& particle,
289 std::vector<SolidBoundary<T, DESCRIPTOR::d>>& solidBoundaries,
290 const PhysR<T, DESCRIPTOR::d>& cellMin,
291 const PhysR<T, DESCRIPTOR::d>& cellMax, F getSetupPeriodicity)
292{
293 using namespace descriptors;
294 constexpr unsigned D = DESCRIPTOR::d;
295 const T deltaX = blockGeometry.getDeltaR();
296 auto sIndicator = particle.template getField<SURFACE, SINDICATOR>();
297 const T detectionDistance = particles::contact::evalContactDetectionDistance(particle, deltaX);
298 const T circumRadius = particles::contact::evalCircumRadius(
299 detectionDistance, sIndicator->getCircumRadius(), sIndicator->getEpsilon());
300 auto position = particle.template getField<GENERAL, POSITION>();
301 //Retrieve wallMaterials
302 std::unordered_set<int> wallMaterials = getLatticeMaterials( solidBoundaries );
303
304 //For all cells in block particle intersection
305 int padding = blockGeometry.getPadding();
307 blockGeometry, blockLattice, padding, position, circumRadius,
308 [&](const LatticeR<D>& latticeR) {
309 //Get phys position
310 T physR[D] = {};
311 blockGeometry.getPhysR(physR, latticeR);
313 particle, PhysR<T, D>(physR));
314 //Retrieve porosity
316 sdf::evalSolidVolumeFraction(porosity.data(), signedDist,
317 sIndicator->getEpsilon());
318
319 //For cells containing particle bits
320 if (!util::nearZero(porosity[0])) {
321 //Retrieve material at cell and ensure location outside of boundary
322 auto material = blockGeometry.getMaterial(latticeR);
323 if (material != 0 && wallMaterials.find(material)==wallMaterials.end()) {
324 //Retrieve local velocity
327 particle, PhysR<T, D>(physR));
328 for (unsigned iDim = 0; iDim != D; ++iDim) {
329 velocity[iDim] = converter.getLatticeVelocity(velocity[iDim]);
330 }
331 //Calculate solid volume fraction
332 T solidVolumeFraction = 1. - porosity[0];
333 //Write to fields
334 {
335 auto cell = blockLattice.get(latticeR);
336 if constexpr (!DESCRIPTOR::template provides<descriptors::VELOCITY_SOLID>()) {
337 cell.template getFieldPointer<descriptors::VELOCITY_NUMERATOR>() +=
338 porosity[0] * velocity;
339 cell.template getFieldPointer<
341 }
342 else {
343 cell.template getFieldPointer<descriptors::VELOCITY_SOLID>() += velocity;
344 }
345 cell.template getFieldPointer<descriptors::POROSITY>() *=
346 solidVolumeFraction;
347 }
348 } //if (wallMaterials.find(material)==wallMaterials.end())
349 } //if (!util::nearZero(porosity[0]))
350
351 // For contacts we have another distance
352 // Particle-particle detection (0.5 * deltaX layer for each equals deltaX which is then the same as for the wall detection)
353 if (signedDist < detectionDistance) {
354 // We must ignore cells on padding
355 bool ignoreCell = false;
356 for (unsigned i = 0; i < D; ++i) {
357 if (latticeR[i] < 0 ||
358 latticeR[i] > blockGeometry.getExtent()[i] - 1) {
359 ignoreCell = true;
360 }
361 }
362
363 // processing contact detection
364 if (!ignoreCell) {
365 auto contactHelper = blockLattice.get(latticeR)
366 .template getFieldPointer<
368 // this checks if another particle is already on this cell, if yes, the contacts are updated. If no, the ID (shifted by 1) is stored.
369 std::size_t iP2 = contactHelper[0] - 1;
370 if (contactHelper[0] > 0 && iP2 != iP) {
371 std::size_t localiP2 = iP2;
373 PARTICLETYPE>()) {
374 for (std::size_t localiP = 0; localiP < particleSystem.size();
375 ++localiP) {
376 const std::size_t globalParticleID =
377 particleSystem.get(localiP)
378 .template getField<PARALLELIZATION, ID>();
379 if (globalParticleID == iP2) {
380 localiP2 = localiP;
381 break;
382 }
383 }
384 }
385 auto particle2 = particleSystem.get(localiP2);
386
387 // Check that at least one particle should not ignore contacts
388 bool detectParticleContacts = true;
389 if constexpr (access::providesComputeContact<PARTICLETYPE>()) {
390 detectParticleContacts =
391 access::isContactComputationEnabled(particle, particle2);
392 }
393
394 if (detectParticleContacts) {
396 contactContainer, std::array<std::size_t, 2>({iP2, iP}),
397 PhysR<T, D>(physR), particle2, particle, cellMin, cellMax,
398 getSetupPeriodicity, deltaX);
399 }
400 }
401 else {
402 contactHelper[0] = iP + 1;
403 }
404
405 // Check if contacts should be ignored
406 bool detectWallContacts = true;
407 if constexpr (access::providesComputeContact<PARTICLETYPE>()) {
408 detectWallContacts =
410 }
411
412 if (detectWallContacts) {
413 // TODO: Use the following line with c++20
414 //for (unsigned short solidBoundaryID = 0; SolidBoundary<T, DESCRIPTOR::d>& solidBoundary : solidBoundaries)
415 unsigned solidBoundaryID = 0;
416 for (SolidBoundary<T, DESCRIPTOR::d>& solidBoundary :
417 solidBoundaries) {
418 const T solidBoundaryDetectionDistance =
419 0.5 * util::sqrt(D) * deltaX +
420 solidBoundary.getEnlargementForContact();
421 if (solidBoundary.getIndicator()->signedDistance(
422 PhysR<T, D>(physR)) < solidBoundaryDetectionDistance) {
424 contactContainer, iP, solidBoundaryID, PhysR<T, D>(physR),
425 particle, cellMin, cellMax, getSetupPeriodicity, deltaX);
426 }
427 ++solidBoundaryID;
428 }
429 }
430 }
431 }
432 });
433}
434
435template <typename T, typename DESCRIPTOR>
437 const BlockGeometry<T, DESCRIPTOR::d>& blockGeometry,
438 BlockLattice<T, DESCRIPTOR>& blockLattice)
439{
440 constexpr unsigned D = DESCRIPTOR::d;
441 int padding = blockGeometry.getPadding();
442 LatticeR<D> blockMin(0. - padding);
443 LatticeR<D> blockMax(blockGeometry.getExtent() - 1 + padding);
444
445 using namespace descriptors;
446 blockLattice.forSpatialLocations(
447 blockMin, blockMax, [&](const LatticeR<D>& latticeR) {
448 auto cell = blockLattice.get(latticeR);
449 particles::resetAllParticleRelatedFields<DESCRIPTOR, decltype(cell), typename decltype(cell)::value_t>(cell);
450 });
451}
452
453//TODO: HOTFIX ONLY
454template <typename T, typename DESCRIPTOR>
456 const BlockGeometry<T, DESCRIPTOR::d>& blockGeometry,
457 BlockLattice<T, DESCRIPTOR>& blockLattice)
458{
459 constexpr unsigned D = DESCRIPTOR::d;
460 int padding = blockGeometry.getPadding();
461 LatticeR<D> blockMin(0. - padding);
462 LatticeR<D> blockMax(blockGeometry.getExtent() - 1 + padding);
463
464 using namespace descriptors;
465 blockLattice.forSpatialLocations(
466 blockMin, blockMax, [&](LatticeR<D> latticeR) {
467 auto cell = blockLattice.get(latticeR);
468 particles::resetParticleContactRelatedFields<DESCRIPTOR, decltype(cell), typename decltype(cell)::value_t>(cell);
469 });
470}
471
472} //namespace particles
473
474} // namespace olb
475
476#endif
Representation of a block geometry.
Vector< T, D > getPhysR(LatticeR< D > latticeR)
int getMaterial(LatticeR< D > latticeR) const
returns the (iX,iY) entry in the 2D scalar field
T getDeltaR() const
Read only access to the voxel size given in SI units (meter)
Vector< T, D > getOrigin() const
Read only access to the origin position given in SI units (meter)
Vector< int, D > getExtent() const
Returns the extend of the block in lattice units.
std::enable_if_t< sizeof...(L)==D, int > get(L... latticeR) const
Read-only access to a material number.
Platform-abstracted block lattice for external access and inter-block interaction.
Cell< T, DESCRIPTOR > get(CellID iCell)
Get Cell interface for index iCell.
void forSpatialLocations(F f) const
int getPadding() const
Read only access to padding.
Conversion between physical and lattice units, as well as discretization.
constexpr T getLatticeVelocity(T physVelocity) const
conversion from physical to lattice velocity
Plain old scalar vector.
Definition vector.h:47
auto & get()
Expose container.
constexpr std::size_t size()
Size of ParticleSystem.
constexpr bool providesParallelization()
Provides group PARALLELIZATION.
T getRadius(Particle< T, PARTICLETYPE > &particle)
bool isContactComputationEnabled(Particle< T, PARTICLETYPE > &particle)
Check if contact should be regarded (specification for a single particle)
Vector< T, PARTICLETYPE::d > getPosition(Particle< T, PARTICLETYPE > particle)
T movePositionToEnd(const T position, const T max, const T min)
Definition utilities.h:260
T movePositionToStart(const T position, const T max, const T min)
Definition utilities.h:254
T evalCircumRadius(T const contactDetectionDistance, T const circumRadius, T const epsilon)
void updateContacts(ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer, std::array< std::size_t, 2 > &&ids, const PhysR< T, PARTICLETYPE::d > &pos, Particle< T, PARTICLETYPE > &particle1, Particle< T, PARTICLETYPE > &particle2, const PhysR< T, PARTICLETYPE::d > &cellMin, const PhysR< T, PARTICLETYPE::d > &cellMax, F getSetupPeriodicity, T deltaX)
T evalContactDetectionDistance(Particle< T, PARTICLETYPE > &particle, T const physDeltaX)
constexpr Vector< T, PARTICLETYPE::d > calculateLocalVelocity(Particle< T, PARTICLETYPE > &particle, const PhysR< T, PARTICLETYPE::d > &input)
const S signedDistanceToParticle(Particle< T, PARTICLETYPE > &particle, const PhysR< S, PARTICLETYPE::d > &input)
bool evalSolidVolumeFraction(T output[], const S input[], Particle< T, PARTICLETYPE > &particle)
void resetBlockParticleField(const BlockGeometry< T, DESCRIPTOR::d > &blockGeometry, BlockLattice< T, DESCRIPTOR > &blockLattice)
bool getBlockParticleIntersection(const BlockGeometry< T, D > &blockGeometry, T invDeltaX, LatticeR< D > &start, LatticeR< D > &end, Vector< T, D > position, T circumRadius)
void setBlockParticleField(const BlockGeometry< T, DESCRIPTOR::d > &blockGeometry, BlockLattice< T, DESCRIPTOR > &blockLattice, UnitConverter< T, DESCRIPTOR > const &converter, Particle< T, PARTICLETYPE > &particle)
void resetBlockContactField(const BlockGeometry< T, DESCRIPTOR::d > &blockGeometry, BlockLattice< T, DESCRIPTOR > &blockLattice)
void checkSmoothIndicatorOutOfGeometry(bool &outOfGeometry, Vector< T, D > &ghostPos, const PhysR< T, D > &cellMin, const PhysR< T, D > &cellMax, const Vector< T, D > &position, T circumRadius, const Vector< bool, D > &periodic)
void forSpatialLocationsInBlockParticleIntersection(const BlockGeometry< T, DESCRIPTOR::d > &blockGeometry, BlockLattice< T, DESCRIPTOR > &blockLattice, Vector< T, DESCRIPTOR::d > position, T circumRadius, F f)
bool evalSolidVolumeFraction(T output[], T signedDist, T eps) any_platform
Definition sdf.h:453
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
ADf< T, DIM > ceil(const ADf< T, DIM > &a)
Definition aDiff.h:900
ADf< T, DIM > floor(const ADf< T, DIM > &a)
Definition aDiff.h:869
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
Top level namespace for all of OpenLB.
std::unordered_set< int > getLatticeMaterials(const std::vector< SolidBoundary< T, D > > &solidBoundaries)
Get material numbers of multiple solid boundaries in std::vector.
Definition wall.h:64