gradAlpha.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 -------------------------------------------------------------------------------
10 License
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 "gradAlpha.H"
29 #include "fvc.H"
30 #include "leastSquareGrad.H"
32 #include "profiling.H"
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 namespace reconstruction
38 {
39  defineTypeNameAndDebug(gradAlpha, 0);
40  addToRunTimeSelectionTable(reconstructionSchemes, gradAlpha, components);
41 }
42 }
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 void Foam::reconstruction::gradAlpha::gradSurf(const volScalarField& phi)
48 {
49  addProfilingInFunction(geometricVoF);
50  leastSquareGrad<scalar> lsGrad("polyDegree1",mesh_.geometricD());
51 
52  zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
53 
54  exchangeFields.setUpCommforZone(interfaceCell_,true);
55 
56  Map<vector> mapCC
57  (
58  exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
59  );
60  Map<scalar> mapPhi
61  (
62  exchangeFields.getDatafromOtherProc(interfaceCell_, phi)
63  );
64 
65  DynamicField<vector> cellCentre(100);
66  DynamicField<scalar> phiValues(100);
67 
68  const labelListList& stencil = exchangeFields.getStencil();
69 
71  {
72  const label celli = interfaceLabels_[i];
73 
74  cellCentre.clear();
75  phiValues.clear();
76 
77  for (const label gblIdx : stencil[celli])
78  {
79  cellCentre.append
80  (
81  exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
82  );
83  phiValues.append
84  (
85  exchangeFields.getValue(phi, mapPhi, gblIdx)
86  );
87  }
88 
89  cellCentre -= mesh_.C()[celli];
90  interfaceNormal_[i] = lsGrad.grad(cellCentre, phiValues);
91  }
92 }
93 
94 
95 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
96 
97 Foam::reconstruction::gradAlpha::gradAlpha
98 (
100  const surfaceScalarField& phi,
101  const volVectorField& U,
102  const dictionary& dict
103 )
104 :
106  (
107  typeName,
108  alpha1,
109  phi,
110  U,
111  dict
112  ),
113  mesh_(alpha1.mesh()),
114  interfaceNormal_(fvc::grad(alpha1)),
115  isoFaceTol_(modelDict().getOrDefault<scalar>("isoFaceTol", 1e-8)),
116  surfCellTol_(modelDict().getOrDefault<scalar>("surfCellTol", 1e-8)),
117  sIterPLIC_(mesh_,surfCellTol_)
118 {
119  reconstruct();
120 }
121 
122 
123 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
124 
126 {
127  addProfilingInFunction(geometricVoF);
128  const bool uptodate = alreadyReconstructed(forceUpdate);
129 
130  if (uptodate && !forceUpdate)
131  {
132  return;
133  }
134 
135  if (mesh_.topoChanging())
136  {
137  // Introduced resizing to cope with changing meshes
138  if (interfaceCell_.size() != mesh_.nCells())
139  {
140  interfaceCell_.resize(mesh_.nCells());
141  }
142  }
143  interfaceCell_ = false;
144 
145  interfaceLabels_.clear();
146 
147  forAll(alpha1_, celli)
148  {
149  if (sIterPLIC_.isASurfaceCell(alpha1_[celli]))
150  {
151  interfaceCell_[celli] = true; // is set to false earlier
152  interfaceLabels_.append(celli);
153  }
154  }
155  interfaceNormal_.resize(interfaceLabels_.size());
156  centre_ = dimensionedVector("centre", dimLength, Zero);
157  normal_ = dimensionedVector("normal", dimArea, Zero);
158 
159  gradSurf(alpha1_);
160 
161  forAll(interfaceLabels_, i)
162  {
163  const label celli = interfaceLabels_[i];
164  if (mag(interfaceNormal_[i]) == 0)
165  {
166  continue;
167  }
168 
169  sIterPLIC_.vofCutCell
170  (
171  celli,
172  alpha1_[celli],
173  isoFaceTol_,
174  100,
175  interfaceNormal_[i]
176  );
177 
178  if (sIterPLIC_.cellStatus() == 0)
179  {
180  normal_[celli] = sIterPLIC_.surfaceArea();
181  centre_[celli] = sIterPLIC_.surfaceCentre();
182  if (mag(normal_[celli]) == 0)
183  {
184  normal_[celli] = Zero;
185  centre_[celli] = Zero;
186  }
187  }
188  else
189  {
190  normal_[celli] = Zero;
191  centre_[celli] = Zero;
192  }
193  }
194 }
195 
196 
198 {
199  // Without this line, we seem to get a race condition
200  mesh_.C();
201 
202  cutCellPLIC cutCell(mesh_);
203 
204  forAll(normal_, celli)
205  {
206  if (mag(normal_[celli]) != 0)
207  {
208  vector n = normal_[celli]/mag(normal_[celli]);
209  scalar cutValue = (centre_[celli] - mesh_.C()[celli]) & (n);
210  cutCell.calcSubCell
211  (
212  celli,
213  cutValue,
214  n
215  );
216  alpha1_[celli] = cutCell.VolumeOfFluid();
217  }
218  }
219 
220  alpha1_.correctBoundaryConditions();
221  alpha1_.oldTime () = alpha1_;
222  alpha1_.oldTime().correctBoundaryConditions();
223 }
224 
225 
226 // ************************************************************************* //
profiling.H
Foam::fvc::reconstruct
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcReconstruct.C:56
Foam::reconstruction::defineTypeNameAndDebug
defineTypeNameAndDebug(isoAlpha, 0)
Foam::fvc::grad
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
Foam::dimLength
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
leastSquareGrad.H
Foam::cutCellPLIC
Class for cutting a cell, cellI, of an fvMesh, mesh_, at its intersection with an surface defined by ...
Definition: cutCellPLIC.H:71
alpha1
const volScalarField & alpha1
Definition: setRegionFluidFields.H:8
Foam::polyMesh::geometricD
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:858
Foam::reconstructionSchemes::interfaceCell_
boolList interfaceCell_
Is interface cell?
Definition: reconstructionSchemes.H:88
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::dimensionedVector
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Definition: dimensionedVector.H:50
Foam::reconstruction::gradAlpha::mapAlphaField
virtual void mapAlphaField() const
map VoF Field in case of refinement
Definition: gradAlpha.C:197
Foam::reconstruction::gradAlpha::reconstruct
virtual void reconstruct(bool forceUpdate=true)
Reconstruct interface.
Definition: gradAlpha.C:125
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Foam::fvMesh::C
const volVectorField & C() const
Return cell centres as volVectorField.
Definition: fvMeshGeometry.C:341
Foam::cutCell
Service routines for cutting a cell, celli, of an fvMesh, mesh_, at its intersection with a surface.
Definition: cutCell.H:59
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
gradAlpha.H
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
addProfilingInFunction
#define addProfilingInFunction(name)
Definition: profilingTrigger.H:124
Foam::reconstruction::addToRunTimeSelectionTable
addToRunTimeSelectionTable(reconstructionSchemes, isoAlpha, components)
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::zoneDistribute::New
static zoneDistribute & New(const fvMesh &)
Definition: zoneDistribute.C:105
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
U
U
Definition: pEqn.H:72
Foam::Vector< scalar >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
fvc.H
Foam::reconstructionSchemes
Original code supplied by Henning Scheufler, DLR (2019)
Definition: reconstructionSchemes.H:58
Foam::GeometricField< scalar, fvPatchField, volMesh >
Foam::DynamicField::clear
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicFieldI.H:431
Foam::reconstructionSchemes::interfaceLabels_
DynamicField< label > interfaceLabels_
List of cell labels that have an interface.
Definition: reconstructionSchemes.H:91