isoAdvectionTemplates.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) 2016-2017 DHI
9  Modified code Copyright (C) 2016-2017 OpenCFD Ltd.
10  Modified code Copyright (C) 2019-2020 DLR
11  Modified code Copyright (C) 2021 Johan Roenby
12 -------------------------------------------------------------------------------
13 License
14  This file is part of OpenFOAM.
15 
16  OpenFOAM is free software: you can redistribute it and/or modify it
17  under the terms of the GNU General Public License as published by
18  the Free Software Foundation, either version 3 of the License, or
19  (at your option) any later version.
20 
21  OpenFOAM 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  You should have received a copy of the GNU General Public License
27  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
28 
29 \*---------------------------------------------------------------------------*/
30 
31 #include "isoAdvection.H"
32 #include "fvcSurfaceIntegrate.H"
33 #include "upwind.H"
34 
35 // ************************************************************************* //
36 
37 template<typename Type>
38 Type Foam::isoAdvection::faceValue
39 (
40  const GeometricField<Type, fvsPatchField, surfaceMesh>& f,
41  const label facei
42 ) const
43 {
44  if (mesh_.isInternalFace(facei))
45  {
46  return f.primitiveField()[facei];
47  }
48  else
49  {
50  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
51 
52  // Boundary face. Find out which face of which patch
53  const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
54 
55  if (patchi < 0 || patchi >= pbm.size())
56  {
58  << "Cannot find patch for face " << facei
59  << abort(FatalError);
60  }
61 
62  // Handle empty patches
63  const polyPatch& pp = pbm[patchi];
64  if (isA<emptyPolyPatch>(pp) || pp.empty())
65  {
66  return pTraits<Type>::zero;
67  }
68 
69  const label patchFacei = pp.whichFace(facei);
70  return f.boundaryField()[patchi][patchFacei];
71  }
72 }
73 
74 
75 template<typename Type>
76 void Foam::isoAdvection::setFaceValue
77 (
78  GeometricField<Type, fvsPatchField, surfaceMesh>& f,
79  const label facei,
80  const Type& value
81 ) const
82 {
83  if (mesh_.isInternalFace(facei))
84  {
85  f.primitiveFieldRef()[facei] = value;
86  }
87  else
88  {
89  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
90 
91  // Boundary face. Find out which face of which patch
92  const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
93 
94  if (patchi < 0 || patchi >= pbm.size())
95  {
97  << "Cannot find patch for face " << facei
98  << abort(FatalError);
99  }
100 
101  // Handle empty patches
102  const polyPatch& pp = pbm[patchi];
103  if (isA<emptyPolyPatch>(pp) || pp.empty())
104  {
105  return;
106  }
107 
108  const label patchFacei = pp.whichFace(facei);
109 
110  f.boundaryFieldRef()[patchi][patchFacei] = value;
111  }
112 }
113 
114 
115 template<class SpType, class SuType>
116 void Foam::isoAdvection::limitFluxes
117 (
118  const SpType& Sp,
119  const SuType& Su
120 )
121 {
122  addProfilingInFunction(geometricVoF);
124 
125  const scalar aTol = 100*SMALL; // Note: tolerances
126  scalar maxAlphaMinus1 = gMax(alpha1_) - 1; // max(alphaNew - 1);
127  scalar minAlpha = gMin(alpha1_); // min(alphaNew);
128  const label nOvershoots = 20; // sum(pos0(alphaNew - 1 - aTol));
129 
130  const labelList& owner = mesh_.faceOwner();
131  const labelList& neighbour = mesh_.faceNeighbour();
132 
133  DebugInfo
134  << "isoAdvection: Before conservative bounding: min(alpha) = "
135  << minAlpha << ", max(alpha) = 1 + " << maxAlphaMinus1 << endl;
136 
137  surfaceScalarField dVfcorrectionValues("dVfcorrectionValues", dVf_*0.0);
138 
139  bitSet needBounding(mesh_.nCells(),false);
140  needBounding.set(surfCells_);
141 
142  extendMarkedCells(needBounding);
143 
144  // Loop number of bounding steps
145  for (label n = 0; n < nAlphaBounds_; n++)
146  {
147  if (maxAlphaMinus1 > aTol || minAlpha < -aTol) // Note: tolerances
148  {
149  DebugInfo << "boundAlpha... " << endl;
150 
151  DynamicList<label> correctedFaces(3*nOvershoots);
152  dVfcorrectionValues = dimensionedScalar("0",dimVolume,0.0);
153  boundFlux(needBounding, dVfcorrectionValues, correctedFaces,Sp,Su);
154 
155  correctedFaces.append
156  (
157  syncProcPatches(dVfcorrectionValues, phi_,true)
158  );
159 
160  labelHashSet alreadyUpdated;
161  for (const label facei : correctedFaces)
162  {
163  if (alreadyUpdated.insert(facei))
164  {
165  checkIfOnProcPatch(facei);
166  const label own = owner[facei];
167  scalar Vown = mesh_.V()[own];
168  if (porosityEnabled_)
169  {
170  Vown *= porosityPtr_->primitiveField()[own];
171  }
172  alpha1_[own] -= faceValue(dVfcorrectionValues, facei)/Vown;
173 
174  if (mesh_.isInternalFace(facei))
175  {
176  const label nei = neighbour[facei];
177  scalar Vnei = mesh_.V()[nei];
178  if (porosityEnabled_)
179  {
180  Vnei *= porosityPtr_->primitiveField()[nei];
181  }
182  alpha1_[nei] +=
183  faceValue(dVfcorrectionValues, facei)/Vnei;
184  }
185 
186  // Change to treat boundaries consistently
187  const scalar corrVf =
188  faceValue(dVf_, facei)
189  + faceValue(dVfcorrectionValues, facei);
190 
191  setFaceValue(dVf_, facei, corrVf);
192  }
193  }
194  syncProcPatches(dVf_, phi_);
195  }
196  else
197  {
198  break;
199  }
200 
201  maxAlphaMinus1 = gMax(alpha1_) - 1; // max(alphaNew - 1);
202  minAlpha = gMin(alpha1_); // min(alphaNew);
203 
204  if (debug)
205  {
206  // Check if still unbounded
207  //scalarField alphaNew(alpha1In_ - fvc::surfaceIntegrate(dVf_)());
208  label maxAlphaMinus1 = max(alpha1_.primitiveField() - 1);
209  scalar minAlpha = min(alpha1_.primitiveField());
210  label nUndershoots = sum(neg0(alpha1_.primitiveField() + aTol));
211  label nOvershoots = sum(pos0(alpha1_.primitiveField() - 1 - aTol));
212 
213  Info<< "After bounding number " << n + 1 << " of time "
214  << mesh_.time().value() << ":" << nl
215  << "nOvershoots = " << nOvershoots << " with max(alpha1_-1) = "
216  << maxAlphaMinus1 << " and nUndershoots = " << nUndershoots
217  << " with min(alpha1_) = " << minAlpha << endl;
218  }
219  }
220 
221  alpha1_.correctBoundaryConditions();
222 
223 }
224 
225 
226 template<class SpType, class SuType>
227 void Foam::isoAdvection::boundFlux
228 (
229  const bitSet& nextToInterface,
230  surfaceScalarField& dVfcorrectionValues,
231  DynamicList<label>& correctedFaces,
232  const SpType& Sp,
233  const SuType& Su
234 )
235 {
236  addProfilingInFunction(geometricVoF);
238  scalar rDeltaT = 1/mesh_.time().deltaTValue();
239 
240  correctedFaces.clear();
241  const scalar aTol = 100*SMALL; // Note: tolerances
242 
243  const scalarField& meshV = mesh_.cellVolumes();
244  const scalar dt = mesh_.time().deltaTValue();
245 
246  DynamicList<label> downwindFaces(10);
247  DynamicList<label> facesToPassFluidThrough(downwindFaces.size());
248  DynamicList<scalar> dVfmax(downwindFaces.size());
249  DynamicList<scalar> phi(downwindFaces.size());
250  const volScalarField& alphaOld = alpha1_.oldTime();
251 
252  // Loop through alpha cell centred field
253  for (label celli: nextToInterface)
254  {
255  if (alpha1_[celli] < -aTol || alpha1_[celli] > 1 + aTol)
256  {
257  scalar Vi = meshV[celli];
258  if (porosityEnabled_)
259  {
260  Vi *= porosityPtr_->primitiveField()[celli];
261  }
262 
263  scalar alphaOvershoot =
264  pos0(alpha1_[celli] - 1)*(alpha1_[celli] - 1)
265  + neg0(alpha1_[celli])*alpha1_[celli];
266 
267  scalar fluidToPassOn = alphaOvershoot*Vi;
268  label nFacesToPassFluidThrough = 1;
269 
270  bool firstLoop = true;
271 
272  // First try to pass surplus fluid on to neighbour cells that are
273  // not filled and to which dVf < phi*dt
274  for (label iter=0; iter<10; iter++)
275  {
276  if (mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
277  {
278  break;
279  }
280 
281  DebugInfo
282  << "\n\nBounding cell " << celli
283  << " with alpha overshooting " << alphaOvershoot
284  << endl;
285 
286  facesToPassFluidThrough.clear();
287  dVfmax.clear();
288  phi.clear();
289 
290  // Find potential neighbour cells to pass surplus phase to
291  setDownwindFaces(celli, downwindFaces);
292 
293  scalar dVftot = 0;
294  nFacesToPassFluidThrough = 0;
295 
296  for (const label facei : downwindFaces)
297  {
298  const scalar phif = faceValue(phi_, facei);
299 
300  const scalar dVff =
301  faceValue(dVf_, facei)
302  + faceValue(dVfcorrectionValues, facei);
303 
304  const scalar maxExtraFaceFluidTrans =
305  mag(pos0(fluidToPassOn)*phif*dt - dVff);
306 
307  // dVf has same sign as phi and so if phi>0 we have
308  // mag(phi_[facei]*dt) - mag(dVf[facei]) = phi_[facei]*dt
309  // - dVf[facei]
310  // If phi < 0 we have mag(phi_[facei]*dt) -
311  // mag(dVf[facei]) = -phi_[facei]*dt - (-dVf[facei]) > 0
312  // since mag(dVf) < phi*dt
313  DebugInfo
314  << "downwindFace " << facei
315  << " has maxExtraFaceFluidTrans = "
316  << maxExtraFaceFluidTrans << endl;
317 
318  if (maxExtraFaceFluidTrans/Vi > aTol)
319  {
320 // if (maxExtraFaceFluidTrans/Vi > aTol &&
321 // mag(dVfIn[facei])/Vi > aTol) //Last condition may be
322 // important because without this we will flux through uncut
323 // downwind faces
324  facesToPassFluidThrough.append(facei);
325  phi.append(phif);
326  dVfmax.append(maxExtraFaceFluidTrans);
327  dVftot += mag(phif*dt);
328  }
329  }
330 
331  DebugInfo
332  << "\nfacesToPassFluidThrough: "
333  << facesToPassFluidThrough << ", dVftot = "
334  << dVftot << " m3 corresponding to dalpha = "
335  << dVftot/Vi << endl;
336 
337  forAll(facesToPassFluidThrough, fi)
338  {
339  const label facei = facesToPassFluidThrough[fi];
340  scalar fluidToPassThroughFace =
341  mag(fluidToPassOn)*mag(phi[fi]*dt)/dVftot;
342 
343  nFacesToPassFluidThrough +=
344  pos0(dVfmax[fi] - fluidToPassThroughFace);
345 
346  fluidToPassThroughFace =
347  min(fluidToPassThroughFace, dVfmax[fi]);
348 
349  scalar dVff = faceValue(dVfcorrectionValues, facei);
350 
351  dVff +=
352  sign(phi[fi])*fluidToPassThroughFace*sign(fluidToPassOn);
353 
354  setFaceValue(dVfcorrectionValues, facei, dVff);
355 
356  if (firstLoop)
357  {
358  checkIfOnProcPatch(facei);
359  correctedFaces.append(facei);
360  }
361  }
362 
363  firstLoop = false;
364 
365  scalar alpha1New =
366  (
367  alphaOld[celli]*rDeltaT + Su[celli]
368  - netFlux(dVf_, celli)/Vi*rDeltaT
369  - netFlux(dVfcorrectionValues, celli)/Vi*rDeltaT
370  )
371  /
372  (rDeltaT - Sp[celli]);
373 
374  alphaOvershoot =
375  pos0(alpha1New - 1)*(alpha1New - 1)
376  + neg0(alpha1New)*alpha1New;
377 
378  fluidToPassOn = alphaOvershoot*Vi;
379 
380  DebugInfo
381  << "\nNew alpha for cell " << celli << ": "
382  << alpha1New << endl;
383  }
384  }
385  }
386 
387  DebugInfo << "correctedFaces = " << correctedFaces << endl;
388 }
389 
390 
391 template<class SpType, class SuType>
392 void Foam::isoAdvection::advect(const SpType& Sp, const SuType& Su)
393 {
394  addProfilingInFunction(geometricVoF);
396 
397  if (mesh_.topoChanging())
398  {
399  setProcessorPatches();
400  }
401 
402  scalar advectionStartTime = mesh_.time().elapsedCpuTime();
403 
404  const scalar rDeltaT = 1/mesh_.time().deltaTValue();
405 
406  // reconstruct the interface
407  surf_->reconstruct();
408 
409  if (timeIndex_ < mesh_.time().timeIndex())
410  {
411  timeIndex_= mesh_.time().timeIndex();
412  surf_->normal().oldTime() = surf_->normal();
413  surf_->centre().oldTime() = surf_->centre();
414  }
415 
416  // Initialising dVf with upwind values
417  // i.e. phi[facei]*alpha1[upwindCell[facei]]*dt
418  dVf_ = upwind<scalar>(mesh_, phi_).flux(alpha1_)*mesh_.time().deltaT();
419 
420  // Do the isoAdvection on surface cells
421  timeIntegratedFlux();
422 
423  // Adjust alpha for mesh motion
424  if (mesh_.moving())
425  {
426  alpha1In_ *= (mesh_.Vsc0()/mesh_.Vsc());
427  }
428 
429  // Advect the free surface
430  if (porosityEnabled_)
431  {
432  // Should Su and Sp also be divided by porosity?
433  alpha1_.primitiveFieldRef() =
434  (
435  alpha1_.oldTime().primitiveField()*rDeltaT + Su.field()
436  - fvc::surfaceIntegrate(dVf_)().primitiveField()*rDeltaT
437  /porosityPtr_->primitiveField()
438  )/(rDeltaT - Sp.field());
439  }
440  else
441  {
442  alpha1_.primitiveFieldRef() =
443  (
444  alpha1_.oldTime().primitiveField()*rDeltaT
445  + Su.field()
446  - fvc::surfaceIntegrate(dVf_)().primitiveField()*rDeltaT
447  )/(rDeltaT - Sp.field());
448  }
449 
450  alpha1_.correctBoundaryConditions();
451 
452  // Adjust dVf for unbounded cells
453  limitFluxes(Sp, Su);
454 
455  scalar maxAlphaMinus1 = gMax(alpha1In_) - 1;
456  scalar minAlpha = gMin(alpha1In_);
457 
458  Info<< "isoAdvection: After conservative bounding: min(alpha) = "
459  << minAlpha << ", max(alpha) = 1 + " << maxAlphaMinus1 << endl;
460 
461  // Apply non-conservative bounding mechanisms (clipping and snapping)
462  // Note: We should be able to write out alpha before this is done!
463  applyBruteForceBounding();
464 
465  // Write surface cell set and bound cell set if required by user
466  writeSurfaceCells();
467 
468  advectionTime_ += (mesh_.time().elapsedCpuTime() - advectionStartTime);
469  DebugInfo
470  << "isoAdvection: time consumption = "
471  << label(100*advectionTime_/(mesh_.time().elapsedCpuTime() + SMALL))
472  << "%" << endl;
473 
474  alphaPhi_ = dVf_/mesh_.time().deltaT();
475 }
476 
477 
478 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::upwind
Upwind differencing scheme class.
Definition: upwind.H:56
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Definition: GeometricFieldI.H:53
Foam::neg0
dimensionedScalar neg0(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:210
Sp
zeroField Sp
Definition: alphaSuSp.H:2
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:850
Foam::isoAdvection::advect
void advect(const SpType &Sp, const SuType &Su)
Advect the free surface. Updates alpha field, taking into account.
Definition: isoAdvectionTemplates.C:392
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:444
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::pos0
dimensionedScalar pos0(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:188
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:166
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Su
zeroField Su
Definition: alphaSuSp.H:1
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::polyBoundaryMesh::patchID
const labelList & patchID() const
Per boundary face label the patch index.
Definition: polyBoundaryMesh.C:456
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
DebugInFunction
#define DebugInFunction
Report an information message using Foam::Info.
Definition: messageStream.H:388
upwind.H
Foam::zeroField::field
zeroField field() const noexcept
Definition: zeroField.H:67
Foam::limitedSurfaceInterpolationScheme::flux
virtual tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > flux(const GeometricField< Type, fvPatchField, volMesh > &) const
Return the interpolation weighting factors.
Definition: limitedSurfaceInterpolationScheme.C:195
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:42
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
fvcSurfaceIntegrate.H
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField.
Foam::FatalError
error FatalError
addProfilingInFunction
#define addProfilingInFunction(name)
Definition: profilingTrigger.H:124
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
isoAdvection.H
Foam::primitiveMesh::isInternalFace
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
Definition: primitiveMeshI.H:103
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
DebugInfo
#define DebugInfo
Report an information message using Foam::Info.
Definition: messageStream.H:382
Foam::nl
constexpr char nl
Definition: Ostream.H:404
f
labelList f(nPoints)
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition: surfaceFieldsFwd.H:54
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:78
Foam::fvc::surfaceIntegrate
void surfaceIntegrate(Field< Type > &ivf, const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcSurfaceIntegrate.C:46
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
Foam::gMin
Type gMin(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:593
Foam::labelHashSet
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
Foam::dimVolume
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:60
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592