51 Foam::isoAdvection::isoAdvection
93 nAlphaBounds_(dict_.getOrDefault<label>(
"nAlphaBounds", 3)),
94 isoFaceTol_(dict_.getOrDefault<scalar>(
"isoFaceTol", 1
e-8)),
95 surfCellTol_(dict_.getOrDefault<scalar>(
"surfCellTol", 1
e-8)),
96 writeIsoFacesToFile_(dict_.getOrDefault(
"writeIsoFaces", false)),
99 surfCells_(label(0.2*mesh_.nCells())),
100 surf_(reconstructionSchemes::
New(alpha1_, phi_, U_, dict_)),
102 bsFaces_(label(0.2*mesh_.nBoundaryFaces())),
103 bsx0_(bsFaces_.size()),
104 bsn0_(bsFaces_.size()),
105 bsUn0_(bsFaces_.size()),
108 procPatchLabels_(mesh_.
boundary().size()),
109 surfaceCellFacesOnProcPatches_(0)
131 surfaceCellFacesOnProcPatches_.resize(
patches.size());
138 isA<processorPolyPatch>(
patches[patchi])
142 procPatchLabels_.
append(patchi);
151 void Foam::isoAdvection::timeIntegratedFlux()
154 const scalar dt = mesh_.time().deltaTValue();
157 interpolationCellPoint<vector> UInterp(U_);
160 label nSurfaceCells = 0;
168 const scalarField& magSfIn = mesh_.magSf().primitiveField();
172 const cellList& cellFaces = mesh_.cells();
173 const labelList& own = mesh_.faceOwner();
177 DynamicList<List<point>> isoFacePts;
178 const DynamicField<label>& interfaceLabels = surf_->interfaceLabels();
181 forAll(interfaceLabels, i)
183 const label celli = interfaceLabels[i];
184 if (
mag(surf_->normal()[celli]) != 0)
189 surfCells_.append(celli);
192 <<
"\n------------ Cell " << celli <<
" with alpha1 = "
193 << alpha1In_[celli] <<
" and 1-alpha1 = "
194 << 1.0 - alpha1In_[celli] <<
" ------------"
199 vector n0 = -surf_->normal()[celli];
204 const scalar Un0 = UInterp.interpolate(x0, celli) & n0;
207 <<
"calcIsoFace gives initial surface: \nx0 = " << x0
208 <<
", \nn0 = " << n0 <<
", \nUn0 = "
214 const cell& celliFaces = cellFaces[celli];
217 const label facei = celliFaces[fi];
219 if (mesh_.isInternalFace(facei))
221 bool isDownwindFace =
false;
223 if (celli == own[facei])
225 if (phiIn[facei] >= 0)
227 isDownwindFace =
true;
232 if (phiIn[facei] < 0)
234 isDownwindFace =
true;
240 dVfIn[facei] = advectFace_.timeIntegratedFaceFlux
255 bsFaces_.append(facei);
268 const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
269 const surfaceScalarField::Boundary&
phib = phi_.boundaryField();
270 const surfaceScalarField::Boundary& magSfb = mesh_.magSf().boundaryField();
271 surfaceScalarField::Boundary& dVfb = dVf_.boundaryFieldRef();
272 const label nInternalFaces = mesh_.nInternalFaces();
278 const label facei = bsFaces_[i];
279 const label patchi = boundaryMesh.patchID()[facei - nInternalFaces];
280 const label start = boundaryMesh[patchi].start();
282 if (
phib[patchi].size())
284 const label patchFacei = facei - start;
285 const scalar phiP =
phib[patchi][patchFacei];
289 const scalar magSf = magSfb[patchi][patchFacei];
291 dVfb[patchi][patchFacei] = advectFace_.timeIntegratedFaceFlux
304 checkIfOnProcPatch(facei);
310 syncProcPatches(dVf_, phi_);
312 writeIsoFaces(isoFacePts);
314 DebugInfo <<
"Number of isoAdvector surface cells = "
319 void Foam::isoAdvection::setDownwindFaces
322 DynamicLabelList& downwindFaces
328 const labelList& own = mesh_.faceOwner();
330 const cell&
c =
cells[celli];
332 downwindFaces.
clear();
338 const label facei =
c[fi];
339 const scalar
phi = faceValue(phi_, facei);
341 if (own[facei] == celli)
345 downwindFaces.append(facei);
350 downwindFaces.append(facei);
354 downwindFaces.shrink();
358 Foam::scalar Foam::isoAdvection::netFlux
367 const cell&
c = mesh_.cells()[celli];
370 const labelList& own = mesh_.faceOwner();
374 const label facei =
c[fi];
375 const scalar dVff = faceValue(dVf, facei);
377 if (own[facei] == celli)
395 bool returnSyncedFaces
398 DynamicLabelList syncedFaces(0);
399 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
406 forAll(procPatchLabels_, i)
408 const label patchi = procPatchLabels_[i];
410 const processorPolyPatch& procPatch =
411 refCast<const processorPolyPatch>(
patches[patchi]);
413 UOPstream toNbr(procPatch.neighbProcNo(), pBufs);
414 const scalarField& pFlux = dVf.boundaryField()[patchi];
416 const List<label>& surfCellFacesOnProcPatch =
417 surfaceCellFacesOnProcPatches_[patchi];
419 const UIndirectList<scalar> dVfPatch
422 surfCellFacesOnProcPatch
425 toNbr << surfCellFacesOnProcPatch << dVfPatch;
428 pBufs.finishedSends();
432 forAll(procPatchLabels_, patchLabeli)
434 const label patchi = procPatchLabels_[patchLabeli];
436 const processorPolyPatch& procPatch =
437 refCast<const processorPolyPatch>(
patches[patchi]);
439 UIPstream fromNeighb(procPatch.neighbProcNo(), pBufs);
441 List<scalar> nbrdVfs;
443 fromNeighb >> faceIDs >> nbrdVfs;
444 if (returnSyncedFaces)
446 List<label> syncedFaceI(faceIDs);
447 for (label& faceI : syncedFaceI)
449 faceI += procPatch.start();
451 syncedFaces.append(syncedFaceI);
456 Pout<<
"Received at time = " << mesh_.time().value()
457 <<
": surfCellFacesOnProcPatch = " << faceIDs <<
nl
458 <<
"Received at time = " << mesh_.time().value()
459 <<
": dVfPatch = " << nbrdVfs <<
endl;
463 scalarField& localFlux = dVf.boundaryFieldRef()[patchi];
467 const label facei = faceIDs[i];
468 localFlux[facei] = - nbrdVfs[i];
469 if (
debug &&
mag(localFlux[facei] + nbrdVfs[i]) > 10*SMALL)
471 Pout<<
"localFlux[facei] = " << localFlux[facei]
472 <<
" and nbrdVfs[i] = " << nbrdVfs[i]
473 <<
" for facei = " << facei <<
endl;
481 forAll(procPatchLabels_, patchLabeli)
483 const label patchi = procPatchLabels_[patchLabeli];
484 const scalarField& localFlux = dVf.boundaryField()[patchi];
485 Pout<<
"time = " << mesh_.time().value() <<
": localFlux = "
486 << localFlux <<
endl;
491 forAll(surfaceCellFacesOnProcPatches_, patchi)
493 surfaceCellFacesOnProcPatches_[patchi].clear();
501 void Foam::isoAdvection::checkIfOnProcPatch(
const label facei)
503 if (!mesh_.isInternalFace(facei))
505 const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
506 const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
508 if (isA<processorPolyPatch>(pbm[patchi]) && pbm[patchi].size())
510 const label patchFacei = pbm[patchi].whichFace(facei);
511 surfaceCellFacesOnProcPatches_[patchi].append(patchFacei);
521 bool alpha1Changed =
false;
523 const scalar snapAlphaTol = dict_.getOrDefault<scalar>(
"snapTol", 0);
524 if (snapAlphaTol > 0)
528 *
pos0(alpha1_ - snapAlphaTol)
529 *
neg0(alpha1_ - (1.0 - snapAlphaTol))
530 +
pos0(alpha1_ - (1.0 - snapAlphaTol));
532 alpha1Changed =
true;
535 if (dict_.getOrDefault(
"clip",
true))
537 alpha1_ =
min(scalar(1),
max(scalar(0), alpha1_));
538 alpha1Changed =
true;
543 alpha1_.correctBoundaryConditions();
550 if (!mesh_.time().writeTime())
return;
552 if (dict_.getOrDefault(
"writeSurfCells",
false))
559 mesh_.time().timeName(),
576 if (!writeIsoFacesToFile_ || !mesh_.time().writeTime())
return;
581 mesh_.time().globalPath()
583 /
word::printf(
"isoFaces_%012d.obj", mesh_.time().timeIndex())
597 Info<<
nl <<
"isoAdvection: writing iso faces to file: "
598 << os.name() <<
nl <<
endl;
601 forAll(allProcFaces, proci)
610 if (facePts.size() !=
f.size())
615 os.write(
f, facePts,
false);
624 Info<<
nl <<
"isoAdvection: writing iso faces to file: "
625 << os.name() <<
nl <<
endl;
632 if (facePts.size() !=
f.size())
637 os.write(
f, facePts,
false);