37 template<
typename Type>
38 Type Foam::isoAdvection::faceValue
40 const GeometricField<Type, fvsPatchField, surfaceMesh>&
f,
46 return f.primitiveField()[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();
92 const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
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);
110 f.boundaryFieldRef()[patchi][patchFacei] = value;
115 template<
class SpType,
class SuType>
116 void Foam::isoAdvection::limitFluxes
125 const scalar aTol = 100*SMALL;
126 scalar maxAlphaMinus1 =
gMax(alpha1_) - 1;
127 scalar minAlpha =
gMin(alpha1_);
128 const label nOvershoots = 20;
130 const labelList& owner = mesh_.faceOwner();
131 const labelList& neighbour = mesh_.faceNeighbour();
134 <<
"isoAdvection: Before conservative bounding: min(alpha) = " 135 << minAlpha <<
", max(alpha) = 1 + " << maxAlphaMinus1 <<
endl;
139 bitSet needBounding(mesh_.nCells(),
false);
140 needBounding.set(surfCells_);
142 extendMarkedCells(needBounding);
145 for (label
n = 0;
n < nAlphaBounds_;
n++)
147 if (maxAlphaMinus1 > aTol || minAlpha < -aTol)
151 DynamicList<label> correctedFaces(3*nOvershoots);
153 boundFlux(needBounding, dVfcorrectionValues, correctedFaces,
Sp,
Su);
155 correctedFaces.append
157 syncProcPatches(dVfcorrectionValues, phi_,
true)
161 for (
const label facei : correctedFaces)
163 if (alreadyUpdated.insert(facei))
165 checkIfOnProcPatch(facei);
166 const label own = owner[facei];
167 scalar Vown = mesh_.V()[own];
168 if (porosityEnabled_)
170 Vown *= porosityPtr_->primitiveField()[own];
172 alpha1_[own] -= faceValue(dVfcorrectionValues, facei)/Vown;
174 if (mesh_.isInternalFace(facei))
176 const label nei = neighbour[facei];
177 scalar Vnei = mesh_.V()[nei];
178 if (porosityEnabled_)
180 Vnei *= porosityPtr_->primitiveField()[nei];
183 faceValue(dVfcorrectionValues, facei)/Vnei;
187 const scalar corrVf =
188 faceValue(dVf_, facei)
189 + faceValue(dVfcorrectionValues, facei);
191 setFaceValue(dVf_, facei, corrVf);
194 syncProcPatches(dVf_, phi_);
201 maxAlphaMinus1 =
gMax(alpha1_) - 1;
202 minAlpha =
gMin(alpha1_);
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));
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;
221 alpha1_.correctBoundaryConditions();
226 template<
class SpType,
class SuType>
227 void Foam::isoAdvection::boundFlux
229 const bitSet& nextToInterface,
231 DynamicList<label>& correctedFaces,
238 scalar rDeltaT = 1/mesh_.time().deltaTValue();
240 correctedFaces.clear();
241 const scalar aTol = 100*SMALL;
244 const scalar dt = mesh_.time().deltaTValue();
246 DynamicList<label> downwindFaces(10);
247 DynamicList<label> facesToPassFluidThrough(downwindFaces.size());
248 DynamicList<scalar> dVfmax(downwindFaces.size());
249 DynamicList<scalar>
phi(downwindFaces.size());
253 for (label celli: nextToInterface)
255 if (alpha1_[celli] < -aTol || alpha1_[celli] > 1 + aTol)
257 scalar Vi = meshV[celli];
258 if (porosityEnabled_)
263 scalar alphaOvershoot =
264 pos0(alpha1_[celli] - 1)*(alpha1_[celli] - 1)
265 +
neg0(alpha1_[celli])*alpha1_[celli];
267 scalar fluidToPassOn = alphaOvershoot*Vi;
268 label nFacesToPassFluidThrough = 1;
270 bool firstLoop =
true;
274 for (label iter=0; iter<10; iter++)
276 if (
mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
282 <<
"\n\nBounding cell " << celli
283 <<
" with alpha overshooting " << alphaOvershoot
286 facesToPassFluidThrough.clear();
291 setDownwindFaces(celli, downwindFaces);
294 nFacesToPassFluidThrough = 0;
296 for (
const label facei : downwindFaces)
298 const scalar phif = faceValue(phi_, facei);
301 faceValue(dVf_, facei)
302 + faceValue(dVfcorrectionValues, facei);
304 const scalar maxExtraFaceFluidTrans =
305 mag(
pos0(fluidToPassOn)*phif*dt - dVff);
314 <<
"downwindFace " << facei
315 <<
" has maxExtraFaceFluidTrans = " 316 << maxExtraFaceFluidTrans <<
endl;
318 if (maxExtraFaceFluidTrans/Vi > aTol)
324 facesToPassFluidThrough.append(facei);
326 dVfmax.append(maxExtraFaceFluidTrans);
327 dVftot +=
mag(phif*dt);
332 <<
"\nfacesToPassFluidThrough: " 333 << facesToPassFluidThrough <<
", dVftot = " 334 << dVftot <<
" m3 corresponding to dalpha = " 335 << dVftot/Vi <<
endl;
337 forAll(facesToPassFluidThrough, fi)
339 const label facei = facesToPassFluidThrough[fi];
340 scalar fluidToPassThroughFace =
341 mag(fluidToPassOn)*
mag(
phi[fi]*dt)/dVftot;
343 nFacesToPassFluidThrough +=
344 pos0(dVfmax[fi] - fluidToPassThroughFace);
346 fluidToPassThroughFace =
347 min(fluidToPassThroughFace, dVfmax[fi]);
349 scalar dVff = faceValue(dVfcorrectionValues, facei);
352 sign(
phi[fi])*fluidToPassThroughFace*
sign(fluidToPassOn);
354 setFaceValue(dVfcorrectionValues, facei, dVff);
358 checkIfOnProcPatch(facei);
359 correctedFaces.append(facei);
367 alphaOld[celli]*rDeltaT +
Su[celli]
368 - netFlux(dVf_, celli)/Vi*rDeltaT
369 - netFlux(dVfcorrectionValues, celli)/Vi*rDeltaT
372 (rDeltaT -
Sp[celli]);
375 pos0(alpha1New - 1)*(alpha1New - 1)
376 +
neg0(alpha1New)*alpha1New;
378 fluidToPassOn = alphaOvershoot*Vi;
381 <<
"\nNew alpha for cell " << celli <<
": " 382 << alpha1New <<
endl;
391 template<
class SpType,
class SuType>
397 if (mesh_.topoChanging())
399 setProcessorPatches();
402 scalar advectionStartTime = mesh_.time().elapsedCpuTime();
404 const scalar rDeltaT = 1/mesh_.time().deltaTValue();
407 surf_->reconstruct();
409 if (timeIndex_ < mesh_.time().timeIndex())
411 timeIndex_= mesh_.time().timeIndex();
412 surf_->normal().oldTime() = surf_->normal();
413 surf_->centre().oldTime() = surf_->centre();
418 dVf_ = upwind<scalar>(mesh_, phi_).
flux(alpha1_)*mesh_.time().deltaT();
421 timeIntegratedFlux();
426 alpha1In_ *= (mesh_.Vsc0()/mesh_.Vsc());
430 if (porosityEnabled_)
433 alpha1_.primitiveFieldRef() =
435 alpha1_.oldTime().primitiveField()*rDeltaT +
Su.
field()
437 /porosityPtr_->primitiveField()
442 alpha1_.primitiveFieldRef() =
444 alpha1_.oldTime().primitiveField()*rDeltaT
450 alpha1_.correctBoundaryConditions();
455 scalar maxAlphaMinus1 =
gMax(alpha1In_) - 1;
456 scalar minAlpha =
gMin(alpha1In_);
458 Info<<
"isoAdvection: After conservative bounding: min(alpha) = " 459 << minAlpha <<
", max(alpha) = 1 + " << maxAlphaMinus1 <<
endl;
463 applyBruteForceBounding();
468 advectionTime_ += (mesh_.time().elapsedCpuTime() - advectionStartTime);
470 <<
"isoAdvection: time consumption = " 471 << label(100*advectionTime_/(mesh_.time().elapsedCpuTime() + SMALL))
474 alphaPhi_ = dVf_/mesh_.time().deltaT();
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
dimensionedScalar sign(const dimensionedScalar &ds)
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField...
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.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
#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))
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
#define DebugInFunction
Report an information message using Foam::Info.
dimensionedScalar neg0(const dimensionedScalar &ds)
label nInternalFaces() const noexcept
Number of internal faces.
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