plicRDF.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2019-2020 DLR
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "plicRDF.H"
30#include "fvc.H"
31#include "leastSquareGrad.H"
32#include "zoneDistribute.H"
34#include "profiling.H"
35
36// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37
38namespace Foam
39{
40namespace reconstruction
41{
44}
45}
46
47
48// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49
50void Foam::reconstruction::plicRDF::interpolateNormal()
51{
52 addProfilingInFunction(geometricVoF);
53 scalar dt = mesh_.time().deltaTValue();
54
55 leastSquareGrad<scalar> lsGrad("polyDegree1",mesh_.geometricD());
56
57 zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
58 exchangeFields.setUpCommforZone(interfaceCell_, false);
59
60 Map<vector> mapCentre
61 (
62 exchangeFields.getDatafromOtherProc(interfaceCell_, centre_)
63 );
64 Map<vector> mapNormal
65 (
66 exchangeFields.getDatafromOtherProc(interfaceCell_, normal_)
67 );
68
69 Map<vector> mapCC
70 (
71 exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
72 );
73 Map<scalar> mapAlpha
74 (
75 exchangeFields.getDatafromOtherProc(interfaceCell_, alpha1_)
76 );
77
78 DynamicField<vector> cellCentre(100);
79 DynamicField<scalar> alphaValues(100);
80
81 DynamicList<vector> foundNormals(30);
82
83 const labelListList& stencil = exchangeFields.getStencil();
84
86 {
87 const label celli = interfaceLabels_[i];
88 vector estimatedNormal{Zero};
89 scalar weight{0};
90 foundNormals.clear();
91 for (const label gblIdx : stencil[celli])
92 {
93 vector n =
94 exchangeFields.getValue(normal_, mapNormal, gblIdx);
95 point p = mesh_.C()[celli]-U_[celli]*dt;
96 if (mag(n) != 0)
97 {
98 n /= mag(n);
100 exchangeFields.getValue(centre_, mapCentre, gblIdx);
101 vector distanceToIntSeg = (tensor::I- n*n) & (p - centre);
102 estimatedNormal += n /max(mag(distanceToIntSeg), SMALL);
103 weight += 1/max(mag(distanceToIntSeg), SMALL);
104 foundNormals.append(n);
105 }
106 }
107
108 if (weight != 0 && mag(estimatedNormal) != 0)
109 {
110 estimatedNormal /= weight;
111 estimatedNormal /= mag(estimatedNormal);
112 }
113
114 bool tooCoarse = false;
115
116 if (foundNormals.size() > 1 && mag(estimatedNormal) != 0)
117 {
118 forAll(foundNormals, i)
119 {
120 // all have the length of 1
121 // to coarse if normal angle is bigger than 10 deg
122 if ((estimatedNormal & foundNormals[i]) <= 0.98)
123 {
124 tooCoarse = true;
125 }
126 }
127 }
128 else
129 {
130 tooCoarse = true;
131 }
132
133 // if a normal was found and the interface is fine enough
134 // smallDist is always smallDist
135 if (mag(estimatedNormal) != 0 && !tooCoarse)
136 {
137 interfaceNormal_[i] = estimatedNormal;
138 }
139 else
140 {
141 cellCentre.clear();
142 alphaValues.clear();
143
144 forAll(stencil[celli],i)
145 {
146 const label gblIdx = stencil[celli][i];
147 cellCentre.append
148 (
149 exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
150 );
151 alphaValues.append
152 (
153 exchangeFields.getValue(alpha1_, mapAlpha, gblIdx)
154 );
155 }
156 cellCentre -= mesh_.C()[celli];
157 interfaceNormal_[i] = lsGrad.grad(cellCentre, alphaValues);
158 }
159
160 }
161}
162
163void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
164{
165 addProfilingInFunction(geometricVoF);
166 leastSquareGrad<scalar> lsGrad("polyDegree1", mesh_.geometricD());
167 zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
168
169 exchangeFields.setUpCommforZone(interfaceCell_, false);
170
171 Map<vector> mapCC
172 (
173 exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
174 );
175 Map<scalar> mapPhi
176 (
177 exchangeFields.getDatafromOtherProc(interfaceCell_, phi)
178 );
179
180 DynamicField<vector> cellCentre(100);
181 DynamicField<scalar> phiValues(100);
182
183 const labelListList& stencil = exchangeFields.getStencil();
184
185 forAll(interfaceLabels_, i)
186 {
187 const label celli = interfaceLabels_[i];
188
189 cellCentre.clear();
190 phiValues.clear();
191
192 for (const label gblIdx : stencil[celli])
193 {
194 cellCentre.append
195 (
196 exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
197 );
198 phiValues.append
199 (
200 exchangeFields.getValue(phi, mapPhi, gblIdx)
201 );
202 }
203
204 cellCentre -= mesh_.C()[celli];
205 interfaceNormal_[i] = lsGrad.grad(cellCentre, phiValues);
206 }
207}
208
209
210void Foam::reconstruction::plicRDF::setInitNormals(bool interpolate)
211{
212 addProfilingInFunction(geometricVoF);
213 zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
214
215 interfaceLabels_.clear();
216
217 forAll(alpha1_, celli)
218 {
219 if (sIterPLIC_.isASurfaceCell(alpha1_[celli]))
220 {
221 interfaceCell_[celli] = true; // is set to false earlier
222 interfaceLabels_.append(celli);
223 }
224 }
225 interfaceNormal_.setSize(interfaceLabels_.size());
226
227 RDF_.markCellsNearSurf(interfaceCell_, 1);
228 const boolList& nextToInterface_ = RDF_.nextToInterface();
229 exchangeFields.updateStencil(nextToInterface_);
230
231 if (interpolate)
232 {
233 interpolateNormal();
234 }
235 else
236 {
237 gradSurf(alpha1_);
238 }
239}
240
241
242void Foam::reconstruction::plicRDF::calcResidual
243(
244 List<normalRes>& normalResidual
245)
246{
247 addProfilingInFunction(geometricVoF);
248 zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
249 exchangeFields.setUpCommforZone(interfaceCell_,false);
250
251 Map<vector> mapNormal
252 (
253 exchangeFields.getDatafromOtherProc(interfaceCell_, normal_)
254 );
255
256 const labelListList& stencil = exchangeFields.getStencil();
257
258 forAll(interfaceLabels_, i)
259 {
260 const label celli = interfaceLabels_[i];
261 if (mag(normal_[celli]) == 0 || mag(interfaceNormal_[i]) == 0)
262 {
263 normalResidual[i].celli = celli;
264 normalResidual[i].normalResidual = 0;
265 normalResidual[i].avgAngle = 0;
266 continue;
267 }
268
269 scalar avgDiffNormal = 0;
270 scalar maxDiffNormal = GREAT;
271 scalar weight= 0;
272 const vector cellNormal = normal_[celli]/mag(normal_[celli]);
273 forAll(stencil[celli],j)
274 {
275 const label gblIdx = stencil[celli][j];
276 vector normal =
277 exchangeFields.getValue(normal_, mapNormal, gblIdx);
278
279 if (mag(normal) != 0 && j != 0)
280 {
281 vector n = normal/mag(normal);
282 scalar cosAngle = max(min((cellNormal & n), 1), -1);
283 avgDiffNormal += acos(cosAngle) * mag(normal);
284 weight += mag(normal);
285 if (cosAngle < maxDiffNormal)
286 {
287 maxDiffNormal = cosAngle;
288 }
289 }
290 }
291
292 if (weight != 0)
293 {
294 avgDiffNormal /= weight;
295 }
296 else
297 {
298 avgDiffNormal = 0;
299 }
300
301 vector newCellNormal = normalised(interfaceNormal_[i]);
302
303 scalar normalRes = (1 - (cellNormal & newCellNormal));
304 normalResidual[i].celli = celli;
305 normalResidual[i].normalResidual = normalRes;
306 normalResidual[i].avgAngle = avgDiffNormal;
307 }
308}
309
310
311// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
312
314(
316 const surfaceScalarField& phi,
317 const volVectorField& U,
318 const dictionary& dict
319)
320:
322 (
323 typeName,
324 alpha1,
325 phi,
326 U,
327 dict
328 ),
329 mesh_(alpha1.mesh()),
330
331 interfaceNormal_(0.2*mesh_.nCells()),
332
333 isoFaceTol_(modelDict().getOrDefault<scalar>("isoFaceTol", 1e-8)),
334 surfCellTol_(modelDict().getOrDefault<scalar>("surfCellTol", 1e-8)),
335 tol_(modelDict().getOrDefault<scalar>("tol", 1e-6)),
336 relTol_(modelDict().getOrDefault<scalar>("relTol", 0.1)),
337 iteration_(modelDict().getOrDefault<label>("iterations", 5)),
338 interpolateNormal_(modelDict().getOrDefault("interpolateNormal", true)),
339 RDF_(mesh_),
340 sIterPLIC_(mesh_,surfCellTol_)
341{
342 setInitNormals(false);
343
346
348 {
349 const label celli = interfaceLabels_[i];
350 if (mag(interfaceNormal_[i]) == 0)
351 {
352 continue;
353 }
354 sIterPLIC_.vofCutCell
355 (
356 celli,
357 alpha1_[celli],
358 isoFaceTol_,
359 100,
360 interfaceNormal_[i]
361 );
362
363 if (sIterPLIC_.cellStatus() == 0)
364 {
365 normal_[celli] = sIterPLIC_.surfaceArea();
366 centre_[celli] = sIterPLIC_.surfaceCentre();
367 if (mag(normal_[celli]) == 0)
368 {
369 normal_[celli] = Zero;
370 centre_[celli] = Zero;
371 }
372 }
373 else
374 {
375 normal_[celli] = Zero;
376 centre_[celli] = Zero;
377 }
378 }
379}
380
381
382// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
383
385{
386 addProfilingInFunction(geometricVoF);
387 zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
388 const bool uptodate = alreadyReconstructed(forceUpdate);
389
390 if (uptodate && !forceUpdate)
391 {
392 return;
393 }
394
395 if (mesh_.topoChanging())
396 {
397 // Introduced resizing to cope with changing meshes
398 if (interfaceCell_.size() != mesh_.nCells())
399 {
400 interfaceCell_.resize(mesh_.nCells());
401 }
402 }
403 interfaceCell_ = false;
404
405 // Sets interfaceCell_ and interfaceNormal
406 setInitNormals(interpolateNormal_);
407
408 centre_ = dimensionedVector("centre", dimLength, Zero);
409 normal_ = dimensionedVector("normal", dimArea, Zero);
410
411 // nextToInterface is update on setInitNormals
412 const boolList& nextToInterface_ = RDF_.nextToInterface();
413
414 bitSet tooCoarse(mesh_.nCells(),false);
415
416 for (label iter=0; iter<iteration_; ++iter)
417 {
418 forAll(interfaceLabels_, i)
419 {
420 const label celli = interfaceLabels_[i];
421 if (mag(interfaceNormal_[i]) == 0 || tooCoarse.test(celli))
422 {
423 continue;
424 }
425 sIterPLIC_.vofCutCell
426 (
427 celli,
428 alpha1_[celli],
429 isoFaceTol_,
430 100,
431 interfaceNormal_[i]
432 );
433
434 if (sIterPLIC_.cellStatus() == 0)
435 {
436
437 normal_[celli] = sIterPLIC_.surfaceArea();
438 centre_[celli] = sIterPLIC_.surfaceCentre();
439 if (mag(normal_[celli]) == 0)
440 {
441 normal_[celli] = Zero;
442 centre_[celli] = Zero;
443 }
444 }
445 else
446 {
447 normal_[celli] = Zero;
448 centre_[celli] = Zero;
449 }
450 }
451
452 normal_.correctBoundaryConditions();
453 centre_.correctBoundaryConditions();
454 List<normalRes> normalResidual(interfaceLabels_.size());
455
456 surfaceVectorField::Boundary nHatb(mesh_.Sf().boundaryField());
457 nHatb *= 1/(mesh_.magSf().boundaryField());
458
459 {
460 RDF_.constructRDF
461 (
462 nextToInterface_,
463 centre_,
464 normal_,
465 exchangeFields,
466 false
467 );
468 RDF_.updateContactAngle(alpha1_, U_, nHatb);
469 gradSurf(RDF_);
470 calcResidual(normalResidual);
471 }
472
473 label resCounter = 0;
474 scalar avgRes = 0;
475 scalar avgNormRes = 0;
476
477 forAll(normalResidual,i)
478 {
479
480 const label celli = normalResidual[i].celli;
481 const scalar normalRes= normalResidual[i].normalResidual;
482 const scalar avgA = normalResidual[i].avgAngle;
483
484 if (avgA > 0.26 && iter > 0) // 15 deg
485 {
486 tooCoarse.set(celli);
487 }
488 else
489 {
490 avgRes += normalRes;
491 scalar normRes = 0;
492 scalar discreteError = 0.01*sqr(avgA);
493 if (discreteError != 0)
494 {
495 normRes= normalRes/max(discreteError, tol_);
496 }
497 else
498 {
499 normRes= normalRes/tol_;
500 }
501 avgNormRes += normRes;
502 resCounter++;
503
504 }
505 }
506
507 reduce(avgRes,sumOp<scalar>());
508 reduce(avgNormRes,sumOp<scalar>());
509 reduce(resCounter,sumOp<label>());
510
511 if (resCounter == 0) // avoid division by zero and leave loop
512 {
513 resCounter = 1;
514 avgRes = 0;
515 avgNormRes = 0;
516 }
517
518
519 if (iter == 0)
520 {
522 << "initial residual absolute = "
523 << avgRes/resCounter << nl
524 << "initial residual normalized = "
525 << avgNormRes/resCounter << nl;
526 }
527
528 if
529 (
530 (
531 (avgNormRes/resCounter < relTol_ || avgRes/resCounter < tol_)
532 && (iter >= 1 )
533 )
534 || iter + 1 == iteration_
535 )
536 {
538 << "iterations = " << iter << nl
539 << "final residual absolute = "
540 << avgRes/resCounter << nl
541 << "final residual normalized = " << avgNormRes/resCounter
542 << endl;
543
544 break;
545 }
546 }
547}
548
549
551{
552 // without it we seem to get a race condition
553 mesh_.C();
554
555 cutCellPLIC cutCell(mesh_);
556
557 forAll(normal_, celli)
558 {
559 if (mag(normal_[celli]) != 0)
560 {
561 vector n = normal_[celli]/mag(normal_[celli]);
562 scalar cutValue = (centre_[celli] - mesh_.C()[celli]) & (n);
563 cutCell.calcSubCell
564 (
565 celli,
566 cutValue,
567 n
568 );
569 alpha1_[celli] = cutCell.VolumeOfFluid();
570
571 }
572 }
573 alpha1_.correctBoundaryConditions();
574 alpha1_.oldTime () = alpha1_;
575 alpha1_.oldTime().correctBoundaryConditions();
576}
577
578
579// ************************************************************************* //
label n
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
surfaceScalarField & phi
const volScalarField & alpha1
static const Tensor I
Definition: Tensor.H:81
scalar deltaTValue() const noexcept
Return time step value.
Definition: TimeStateI.H:43
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:521
Class for cutting a cell, cellI, of an fvMesh, mesh_, at its intersection with an surface defined by ...
Definition: cutCellPLIC.H:74
Service routines for cutting a cell, celli, of an fvMesh, mesh_, at its intersection with a surface.
Definition: cutCell.H:60
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
const volVectorField & C() const
Return cell centres as volVectorField.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:872
Original code supplied by Henning Scheufler, DLR (2019)
boolList interfaceCell_
Is interface cell?
volVectorField normal_
Interface area normals.
const volVectorField & U_
Reference to the velocity field.
const volVectorField & centre() const noexcept
const-Reference to interface centres
DynamicField< label > interfaceLabels_
List of cell labels that have an interface.
volScalarField & alpha1_
Reference to the VoF Field.
volVectorField centre_
Interface centres.
Reconstructs an interface (centre and normal vector) consisting of planes to match the internal fluid...
Definition: plicRDF.H:76
virtual void reconstruct(bool forceUpdate=true)
Reconstruct interface.
Definition: plicRDF.C:384
virtual void mapAlphaField() const
Map VoF Field in case of refinement.
Definition: plicRDF.C:550
const point & surfaceCentre() const
The centre of cutting isosurface.
label cellStatus()
The cellStatus.
label vofCutCell(const label celli, const scalar alpha1, const scalar tol, const label maxIter, vector normal)
Finds matching cutValue for the given value fraction.
const vector & surfaceArea() const
The area vector of cutting isosurface.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
Class for parallel communication in a narrow band. It either provides a Map with the neighbouring val...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
U
Definition: pEqn.H:72
volScalarField & p
dynamicFvMesh & mesh
#define DebugInfo
Report an information message using Foam::Info.
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
vector point
Point is a vector.
Definition: point.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
List< bool > boolList
A List of bools.
Definition: List.H:64
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
Definition: curveTools.C:75
dimensionedScalar acos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
#define addProfilingInFunction(name)
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333