36 template<
typename Type>
37 Type Foam::isoAdvection::faceValue
39 const GeometricField<Type, fvsPatchField, surfaceMesh>&
f,
45 return f.primitiveField()[facei];
54 if (patchi < 0 || patchi >= pbm.size())
57 <<
"Cannot find patch for face " << facei
62 const polyPatch& pp = pbm[patchi];
63 if (isA<emptyPolyPatch>(pp) || pp.empty())
65 return pTraits<Type>::zero;
68 const label patchFacei = pp.whichFace(facei);
69 return f.boundaryField()[patchi][patchFacei];
74 template<
typename Type>
75 void Foam::isoAdvection::setFaceValue
77 GeometricField<Type, fvsPatchField, surfaceMesh>&
f,
82 if (mesh_.isInternalFace(facei))
84 f.primitiveFieldRef()[facei] = value;
88 const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
91 const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
93 if (patchi < 0 || patchi >= pbm.size())
96 <<
"Cannot find patch for face " << facei
101 const polyPatch& pp = pbm[patchi];
102 if (isA<emptyPolyPatch>(pp) || pp.empty())
107 const label patchFacei = pp.whichFace(facei);
109 f.boundaryFieldRef()[patchi][patchFacei] = value;
114 template<
class SpType,
class SuType>
115 void Foam::isoAdvection::limitFluxes
123 const scalar aTol = 1.0e-12;
124 scalar maxAlphaMinus1 =
gMax(alpha1_) - 1;
125 scalar minAlpha =
gMin(alpha1_);
126 const label nOvershoots = 20;
128 const labelList& owner = mesh_.faceOwner();
129 const labelList& neighbour = mesh_.faceNeighbour();
132 <<
"isoAdvection: Before conservative bounding: min(alpha) = "
133 << minAlpha <<
", max(alpha) = 1 + " << maxAlphaMinus1 <<
endl;
139 for (label
n = 0;
n < nAlphaBounds_;
n++)
141 if (maxAlphaMinus1 > aTol || minAlpha < -aTol)
145 DynamicList<label> correctedFaces(3*nOvershoots);
147 boundFlux(alpha1In_, dVfcorrectionValues, correctedFaces,
Sp,
Su);
151 syncProcPatches(dVfcorrectionValues, phi_,
true)
155 forAll(correctedFaces, fi)
157 label facei = correctedFaces[fi];
158 if (alreadyUpdated.insert(facei))
160 checkIfOnProcPatch(facei);
161 const label own = owner[facei];
164 -faceValue(dVfcorrectionValues, facei)/mesh_.V()[own];
165 if (mesh_.isInternalFace(facei))
167 const label nei = neighbour[facei];
169 -faceValue(dVfcorrectionValues, facei)/mesh_.V()[nei];
174 faceValue(dVf_, facei)
175 + faceValue(dVfcorrectionValues, facei);
177 setFaceValue(dVf_, facei, corrVf);
180 syncProcPatches(dVf_, phi_);
187 maxAlphaMinus1 =
gMax(alpha1_) - 1;
188 minAlpha =
gMin(alpha1_);
194 label maxAlphaMinus1 =
max(alpha1_.primitiveField() - 1);
195 scalar minAlpha =
min(alpha1_.primitiveField());
196 label nUndershoots =
sum(
neg0(alpha1_.primitiveField() + aTol));
197 label nOvershoots =
sum(
pos0(alpha1_.primitiveField() - 1 - aTol));
199 Info<<
"After bounding number " <<
n + 1 <<
" of time "
200 << mesh_.time().value() <<
":" <<
nl
201 <<
"nOvershoots = " << nOvershoots <<
" with max(alpha1_-1) = "
202 << maxAlphaMinus1 <<
" and nUndershoots = " << nUndershoots
203 <<
" with min(alpha1_) = " << minAlpha <<
endl;
207 alpha1_.correctBoundaryConditions();
212 template<
class SpType,
class SuType>
213 void Foam::isoAdvection::boundFlux
223 scalar rDeltaT = 1/mesh_.time().deltaTValue();
225 correctedFaces.
clear();
226 scalar aTol = 10*SMALL;
229 const scalar dt = mesh_.time().deltaTValue();
242 const scalar Vi = meshV[celli];
245 scalar fluidToPassOn = alphaOvershoot*Vi;
246 label nFacesToPassFluidThrough = 1;
248 bool firstLoop =
true;
252 while (
mag(alphaOvershoot) > aTol && nFacesToPassFluidThrough > 0)
255 <<
"\n\nBounding cell " << celli
256 <<
" with alpha overshooting " << alphaOvershoot
259 facesToPassFluidThrough.clear();
264 setDownwindFaces(celli, downwindFaces);
267 nFacesToPassFluidThrough = 0;
271 const label facei = downwindFaces[fi];
272 const scalar phif = faceValue(phi_, facei);
275 faceValue(dVf_, facei)
276 + faceValue(dVfcorrectionValues, facei);
278 const scalar maxExtraFaceFluidTrans =
279 mag(
pos0(fluidToPassOn)*phif*dt - dVff);
288 <<
"downwindFace " << facei
289 <<
" has maxExtraFaceFluidTrans = "
290 << maxExtraFaceFluidTrans <<
endl;
292 if (maxExtraFaceFluidTrans/Vi > aTol)
298 facesToPassFluidThrough.append(facei);
300 dVfmax.append(maxExtraFaceFluidTrans);
301 dVftot +=
mag(phif*dt);
306 <<
"\nfacesToPassFluidThrough: "
307 << facesToPassFluidThrough <<
", dVftot = "
308 << dVftot <<
" m3 corresponding to dalpha = "
309 << dVftot/Vi <<
endl;
311 forAll(facesToPassFluidThrough, fi)
313 const label facei = facesToPassFluidThrough[fi];
314 scalar fluidToPassThroughFace =
315 mag(fluidToPassOn)*
mag(
phi[fi]*dt)/dVftot;
317 nFacesToPassFluidThrough +=
318 pos0(dVfmax[fi] - fluidToPassThroughFace);
320 fluidToPassThroughFace =
321 min(fluidToPassThroughFace, dVfmax[fi]);
323 scalar dVff = faceValue(dVfcorrectionValues, facei);
326 sign(
phi[fi])*fluidToPassThroughFace*
sign(fluidToPassOn);
328 setFaceValue(dVfcorrectionValues, facei, dVff);
332 checkIfOnProcPatch(facei);
333 correctedFaces.
append(facei);
341 alphaOld[celli]*rDeltaT +
Su[celli]
342 - netFlux(dVf_, celli)/Vi*rDeltaT
343 - netFlux(dVfcorrectionValues, celli)/Vi*rDeltaT
346 (rDeltaT -
Sp[celli]);
349 pos0(alpha1New-1)*(alpha1New-1)
350 +
neg0(alpha1New)*alpha1New;
352 fluidToPassOn = alphaOvershoot*Vi;
355 <<
"\nNew alpha for cell " << celli <<
": "
356 << alpha1New <<
endl;
365 template<
class SpType,
class SuType>
370 scalar advectionStartTime = mesh_.time().elapsedCpuTime();
372 scalar rDeltaT = 1/mesh_.time().deltaTValue();
375 surf_->reconstruct();
382 timeIntegratedFlux();
387 alpha1In_ *= (mesh_.Vsc0()/mesh_.Vsc());
391 alpha1_.primitiveFieldRef() =
393 alpha1_.oldTime().primitiveField()*rDeltaT
398 alpha1_.correctBoundaryConditions();
407 scalar maxAlphaMinus1 =
gMax(alpha1In_) - 1;
408 scalar minAlpha =
gMin(alpha1In_);
410 Info<<
"isoAdvection: After conservative bounding: min(alpha) = "
411 << minAlpha <<
", max(alpha) = 1 + " << maxAlphaMinus1 <<
endl;
415 applyBruteForceBounding();
420 advectionTime_ += (mesh_.time().elapsedCpuTime() - advectionStartTime);
422 <<
"isoAdvection: time consumption = "
423 << label(100*advectionTime_/(mesh_.time().elapsedCpuTime() + SMALL))
426 alphaPhi_ = dVf_/mesh_.time().deltaT();