37 template<
typename Type>
38 Type Foam::isoAdvection::faceValue
40 const GeometricField<Type, fvsPatchField, surfaceMesh>&
f,
46 return f.primitiveField()[facei];
53 const label patchi =
pbm.patchID(facei);
55 if (patchi < 0 || patchi >=
pbm.size())
58 <<
"Cannot find patch for face " << facei
63 const polyPatch&
pp =
pbm[patchi];
64 if (isA<emptyPolyPatch>(
pp) ||
pp.empty())
66 return pTraits<Type>::zero;
69 const label patchFacei =
pp.whichFace(facei);
70 return f.boundaryField()[patchi][patchFacei];
75 template<
typename Type>
76 void Foam::isoAdvection::setFaceValue
78 GeometricField<Type, fvsPatchField, surfaceMesh>&
f,
83 if (mesh_.isInternalFace(facei))
85 f.primitiveFieldRef()[facei] = value;
89 const polyBoundaryMesh&
pbm = mesh_.boundaryMesh();
94 if (patchi < 0 || patchi >=
pbm.
size())
97 <<
"Cannot find patch for face " << facei
102 const polyPatch&
pp =
pbm[patchi];
103 if (isA<emptyPolyPatch>(
pp) ||
pp.empty())
108 const label patchFacei =
pp.whichFace(facei);
109 f.boundaryFieldRef()[patchi][patchFacei] = value;
114 template<
class SpType,
class SuType>
115 void Foam::isoAdvection::limitFluxes
124 const scalar aTol = 100*SMALL;
125 scalar maxAlphaMinus1 =
gMax(alpha1_) - 1;
126 scalar minAlpha =
gMin(alpha1_);
127 const label nOvershoots = 20;
129 const labelList& owner = mesh_.faceOwner();
130 const labelList& neighbour = mesh_.faceNeighbour();
133 <<
"isoAdvection: Before conservative bounding: min(alpha) = " 134 << minAlpha <<
", max(alpha) = 1 + " << maxAlphaMinus1 <<
endl;
138 bitSet needBounding(mesh_.nCells(),
false);
139 needBounding.set(surfCells_);
141 extendMarkedCells(needBounding);
144 for (label
n = 0;
n < nAlphaBounds_;
n++)
146 if (maxAlphaMinus1 > aTol || minAlpha < -aTol)
150 DynamicList<label> correctedFaces(3*nOvershoots);
152 boundFlux(needBounding, dVfcorrectionValues, correctedFaces,
Sp,
Su);
154 correctedFaces.append
156 syncProcPatches(dVfcorrectionValues, phi_,
true)
160 for (
const label facei : correctedFaces)
162 if (alreadyUpdated.insert(facei))
164 checkIfOnProcPatch(facei);
165 const label own = owner[facei];
166 scalar Vown = mesh_.V()[own];
167 if (porosityEnabled_)
169 Vown *= porosityPtr_->primitiveField()[own];
171 alpha1_[own] -= faceValue(dVfcorrectionValues, facei)/Vown;
173 if (mesh_.isInternalFace(facei))
175 const label nei = neighbour[facei];
176 scalar Vnei = mesh_.V()[nei];
177 if (porosityEnabled_)
179 Vnei *= porosityPtr_->primitiveField()[nei];
182 faceValue(dVfcorrectionValues, facei)/Vnei;
186 const scalar corrVf =
187 faceValue(dVf_, facei)
188 + faceValue(dVfcorrectionValues, facei);
190 setFaceValue(dVf_, facei, corrVf);
193 syncProcPatches(dVf_, phi_);
200 maxAlphaMinus1 =
gMax(alpha1_) - 1;
201 minAlpha =
gMin(alpha1_);
207 label maxAlphaMinus1 =
max(alpha1_.primitiveField() - 1);
208 scalar minAlpha =
min(alpha1_.primitiveField());
209 label nUndershoots =
sum(
neg0(alpha1_.primitiveField() + aTol));
210 label nOvershoots =
sum(
pos0(alpha1_.primitiveField() - 1 - aTol));
212 Info<<
"After bounding number " <<
n + 1 <<
" of time " 213 << mesh_.time().value() <<
":" <<
nl 214 <<
"nOvershoots = " << nOvershoots <<
" with max(alpha1_-1) = " 215 << maxAlphaMinus1 <<
" and nUndershoots = " << nUndershoots
216 <<
" with min(alpha1_) = " << minAlpha <<
endl;
220 alpha1_.correctBoundaryConditions();
225 template<
class SpType,
class SuType>
226 void Foam::isoAdvection::boundFlux
228 const bitSet& nextToInterface,
230 DynamicList<label>& correctedFaces,
237 scalar rDeltaT = 1/mesh_.time().deltaTValue();
239 correctedFaces.clear();
240 const scalar aTol = 100*SMALL;
243 const scalar dt = mesh_.time().deltaTValue();
245 DynamicList<label> downwindFaces(10);
246 DynamicList<label> facesToPassFluidThrough(downwindFaces.size());
247 DynamicList<scalar> dVfmax(downwindFaces.size());
248 DynamicList<scalar>
phi(downwindFaces.size());
252 for (label celli: nextToInterface)
254 if (alpha1_[celli] < -aTol || alpha1_[celli] > 1 + aTol)
256 scalar Vi = meshV[celli];
257 if (porosityEnabled_)
262 scalar alphaOvershoot =
263 pos0(alpha1_[celli] - 1)*(alpha1_[celli] - 1)
264 +
neg0(alpha1_[celli])*alpha1_[celli];
266 scalar fluidToPassOn = alphaOvershoot*Vi;
267 label nFacesToPassFluidThrough = 1;
269 bool firstLoop =
true;
273 for (label iter=0; iter<10; iter++)
275 if (
mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
281 <<
"\n\nBounding cell " << celli
282 <<
" with alpha overshooting " << alphaOvershoot
285 facesToPassFluidThrough.clear();
290 setDownwindFaces(celli, downwindFaces);
293 nFacesToPassFluidThrough = 0;
295 for (
const label facei : downwindFaces)
297 const scalar phif = faceValue(phi_, facei);
300 faceValue(dVf_, facei)
301 + faceValue(dVfcorrectionValues, facei);
303 const scalar maxExtraFaceFluidTrans =
304 mag(
pos0(fluidToPassOn)*phif*dt - dVff);
313 <<
"downwindFace " << facei
314 <<
" has maxExtraFaceFluidTrans = " 315 << maxExtraFaceFluidTrans <<
endl;
317 if (maxExtraFaceFluidTrans/Vi > aTol)
323 facesToPassFluidThrough.append(facei);
325 dVfmax.append(maxExtraFaceFluidTrans);
326 dVftot +=
mag(phif*dt);
331 <<
"\nfacesToPassFluidThrough: " 332 << facesToPassFluidThrough <<
", dVftot = " 333 << dVftot <<
" m3 corresponding to dalpha = " 334 << dVftot/Vi <<
endl;
336 forAll(facesToPassFluidThrough, fi)
338 const label facei = facesToPassFluidThrough[fi];
339 scalar fluidToPassThroughFace =
340 mag(fluidToPassOn)*
mag(
phi[fi]*dt)/dVftot;
342 nFacesToPassFluidThrough +=
343 pos0(dVfmax[fi] - fluidToPassThroughFace);
345 fluidToPassThroughFace =
346 min(fluidToPassThroughFace, dVfmax[fi]);
348 scalar dVff = faceValue(dVfcorrectionValues, facei);
351 sign(
phi[fi])*fluidToPassThroughFace*
sign(fluidToPassOn);
353 setFaceValue(dVfcorrectionValues, facei, dVff);
357 checkIfOnProcPatch(facei);
358 correctedFaces.append(facei);
366 alphaOld[celli]*rDeltaT +
Su[celli]
367 - netFlux(dVf_, celli)/Vi*rDeltaT
368 - netFlux(dVfcorrectionValues, celli)/Vi*rDeltaT
371 (rDeltaT -
Sp[celli]);
374 pos0(alpha1New - 1)*(alpha1New - 1)
375 +
neg0(alpha1New)*alpha1New;
377 fluidToPassOn = alphaOvershoot*Vi;
380 <<
"\nNew alpha for cell " << celli <<
": " 381 << alpha1New <<
endl;
390 template<
class SpType,
class SuType>
396 if (mesh_.topoChanging())
398 setProcessorPatches();
401 scalar advectionStartTime = mesh_.time().elapsedCpuTime();
403 const scalar rDeltaT = 1/mesh_.time().deltaTValue();
406 surf_->reconstruct();
408 if (timeIndex_ < mesh_.time().timeIndex())
410 timeIndex_= mesh_.time().timeIndex();
411 surf_->normal().oldTime() = surf_->normal();
412 surf_->centre().oldTime() = surf_->centre();
417 dVf_ = upwind<scalar>(mesh_, phi_).
flux(alpha1_)*mesh_.time().deltaT();
420 timeIntegratedFlux();
425 alpha1In_ *= (mesh_.Vsc0()/mesh_.Vsc());
429 if (porosityEnabled_)
432 alpha1_.primitiveFieldRef() =
434 alpha1_.oldTime().primitiveField()*rDeltaT +
Su.
field()
436 /porosityPtr_->primitiveField()
441 alpha1_.primitiveFieldRef() =
443 alpha1_.oldTime().primitiveField()*rDeltaT
449 alpha1_.correctBoundaryConditions();
454 scalar maxAlphaMinus1 =
gMax(alpha1In_) - 1;
455 scalar minAlpha =
gMin(alpha1In_);
457 Info<<
"isoAdvection: After conservative bounding: min(alpha) = " 458 << minAlpha <<
", max(alpha) = 1 + " << maxAlphaMinus1 <<
endl;
462 applyBruteForceBounding();
467 advectionTime_ += (mesh_.time().elapsedCpuTime() - advectionStartTime);
469 <<
"isoAdvection: time consumption = " 470 << label(100*advectionTime_/(mesh_.time().elapsedCpuTime() + SMALL))
473 alphaPhi_ = dVf_/mesh_.time().deltaT();
dimensionedScalar sign(const dimensionedScalar &ds)
const polyBoundaryMesh & pbm
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField...
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
const Internal::FieldType & primitiveField() const noexcept
Return a const-reference to the internal field values.
void surfaceIntegrate(Field< Type > &ivf, const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Type gMin(const FieldField< Field, Type > &f)
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
const labelList & patchID() const
Per boundary face label the patch index.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
const dimensionSet dimVolume(pow3(dimLength))
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
#define DebugInFunction
Report an information message using Foam::Info.
label size() const noexcept
The number of entries in the list.
dimensionedScalar neg0(const dimensionedScalar &ds)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
#define DebugInfo
Report an information message using Foam::Info.
dimensionedScalar pos0(const dimensionedScalar &ds)
int debug
Static debugging option.
Type gMax(const FieldField< Field, Type > &f)
#define addProfilingInFunction(name)
Define profiling trigger with specified name and description corresponding to the compiler-defined fu...
zeroField field() const noexcept
Method name compatibility with DimensionedField.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
void advect(const SpType &Sp, const SuType &Su)
Advect the free surface. Updates alpha field, taking into account.
messageStream Info
Information stream (stdout output on master, null elsewhere)
List< label > labelList
A List of labels.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())