36 template<
class RdeltaTType,
class RhoType,
class SpType,
class SuType>
39 const RdeltaTType& rDeltaT,
47 Info<<
"MULES: Solving for " << psi.name() <<
endl;
49 const fvMesh& mesh = psi.mesh();
61 mesh.Vsc0()().field()*rho.oldTime().field()
62 *psi0*rDeltaT/mesh.Vsc()().field()
65 )/(rho.field()*rDeltaT - Sp.field());
71 rho.oldTime().field()*psi0*rDeltaT
74 )/(rho.field()*rDeltaT - Sp.field());
77 psi.correctBoundaryConditions();
81 template<
class RhoType,
class SpType,
class SuType>
91 const fvMesh& mesh = psi.mesh();
93 if (fv::localEulerDdt::enabled(mesh))
95 const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
100 const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
106 template<
class RhoType,
class SpType,
class SuType>
119 const fvMesh& mesh = psi.mesh();
121 psi.correctBoundaryConditions();
123 if (fv::localEulerDdt::enabled(mesh))
125 const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
145 const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
166 template<
class RdeltaTType,
class RhoType,
class SpType,
class SuType>
170 const RdeltaTType& rDeltaT,
182 const volScalarField::Boundary& psiBf = psi.boundaryField();
184 const fvMesh& mesh = psi.mesh();
186 const dictionary& MULEScontrols = mesh.solverDict(psi.name());
190 MULEScontrols.lookupOrDefault<
label>(
"nLimiterIter", 3)
195 MULEScontrols.lookupOrDefault<scalar>(
"smoothLimiter", 0)
202 tmp<volScalarField::Internal> tVsc = mesh.Vsc();
206 const surfaceScalarField::Boundary& phiBDBf =
207 phiBD.boundaryField();
210 const surfaceScalarField::Boundary& phiCorrBf =
211 phiCorr.boundaryField();
218 mesh.time().timeName(),
231 surfaceScalarField::Boundary& lambdaBf =
232 lambda.boundaryFieldRef();
244 label own = owner[facei];
245 label nei = neighb[facei];
247 psiMaxn[own] =
max(psiMaxn[own], psiIf[nei]);
248 psiMinn[own] =
min(psiMinn[own], psiIf[nei]);
250 psiMaxn[nei] =
max(psiMaxn[nei], psiIf[own]);
251 psiMinn[nei] =
min(psiMinn[nei], psiIf[own]);
253 sumPhiBD[own] += phiBDIf[facei];
254 sumPhiBD[nei] -= phiBDIf[facei];
256 scalar phiCorrf = phiCorrIf[facei];
260 sumPhip[own] += phiCorrf;
261 mSumPhim[nei] += phiCorrf;
265 mSumPhim[own] -= phiCorrf;
266 sumPhip[nei] -= phiCorrf;
280 const scalarField psiPNf(psiPf.patchNeighbourField());
284 label pfCelli = pFaceCells[pFacei];
286 psiMaxn[pfCelli] =
max(psiMaxn[pfCelli], psiPNf[pFacei]);
287 psiMinn[pfCelli] =
min(psiMinn[pfCelli], psiPNf[pFacei]);
294 label pfCelli = pFaceCells[pFacei];
296 psiMaxn[pfCelli] =
max(psiMaxn[pfCelli], psiPf[pFacei]);
297 psiMinn[pfCelli] =
min(psiMinn[pfCelli], psiPf[pFacei]);
303 label pfCelli = pFaceCells[pFacei];
305 sumPhiBD[pfCelli] += phiBDPf[pFacei];
307 scalar phiCorrf = phiCorrPf[pFacei];
311 sumPhip[pfCelli] += phiCorrf;
315 mSumPhim[pfCelli] -= phiCorrf;
320 psiMaxn =
min(psiMaxn, psiMax);
321 psiMinn =
max(psiMinn, psiMin);
323 if (smoothLimiter > SMALL)
326 min(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMaxn, psiMax);
328 max(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMinn, psiMin);
333 tmp<volScalarField::Internal> V0 = mesh.Vsc0();
338 (rho.field()*rDeltaT - Sp.field())*psiMaxn
341 - (V0().field()*rDeltaT)*rho.oldTime().field()*psi0
348 - (rho.field()*rDeltaT - Sp.field())*psiMinn
350 + (V0().field()*rDeltaT)*rho.oldTime().field()*psi0
358 (rho.field()*rDeltaT - Sp.field())*psiMaxn
360 - (rho.oldTime().field()*rDeltaT)*psi0
368 - (rho.field()*rDeltaT - Sp.field())*psiMinn
369 + (rho.oldTime().field()*rDeltaT)*psi0
377 for (
int j=0; j<nLimiterIter; j++)
384 label own = owner[facei];
385 label nei = neighb[facei];
387 scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
389 if (lambdaPhiCorrf > 0.0)
391 sumlPhip[own] += lambdaPhiCorrf;
392 mSumlPhim[nei] += lambdaPhiCorrf;
396 mSumlPhim[own] -= lambdaPhiCorrf;
397 sumlPhip[nei] -= lambdaPhiCorrf;
410 label pfCelli = pFaceCells[pFacei];
412 scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
414 if (lambdaPhiCorrf > 0.0)
416 sumlPhip[pfCelli] += lambdaPhiCorrf;
420 mSumlPhim[pfCelli] -= lambdaPhiCorrf;
430 (sumlPhip[celli] + psiMaxn[celli])
431 /(mSumPhim[celli] - SMALL),
438 (mSumlPhim[celli] + psiMinn[celli])
439 /(sumPhip[celli] + SMALL),
449 if (phiCorrIf[facei] > 0.0)
451 lambdaIf[facei] =
min 454 min(lambdap[owner[facei]], lambdam[neighb[facei]])
459 lambdaIf[facei] =
min 462 min(lambdam[owner[facei]], lambdap[neighb[facei]])
474 if (isA<wedgeFvPatch>(mesh.boundary()[
patchi]))
478 else if (psiPf.coupled())
481 mesh.boundary()[
patchi].faceCells();
485 label pfCelli = pFaceCells[pFacei];
487 if (phiCorrfPf[pFacei] > 0.0)
490 min(lambdaPf[pFacei], lambdap[pfCelli]);
495 min(lambdaPf[pFacei], lambdam[pfCelli]);
502 mesh.boundary()[
patchi].faceCells();
509 if ((phiBDPf[pFacei] + phiCorrPf[pFacei]) > SMALL*SMALL)
511 label pfCelli = pFaceCells[pFacei];
513 if (phiCorrfPf[pFacei] > 0.0)
516 min(lambdaPf[pFacei], lambdap[pfCelli]);
521 min(lambdaPf[pFacei], lambdam[pfCelli]);
528 syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
533 template<
class RdeltaTType,
class RhoType,
class SpType,
class SuType>
536 const RdeltaTType& rDeltaT,
545 const bool returnCorr
548 const fvMesh& mesh = psi.mesh();
562 mesh.time().timeName(),
594 phiPsi = phiBD +
lambda*phiCorr;
599 template<
class SurfaceScalarFieldList>
612 const surfaceScalarField::Boundary& bfld =
613 phiPsiCorrs[0].boundaryField();
617 if (bfld[
patchi].coupled())
fvsPatchField< scalar > fvsPatchScalarField
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField...
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
void limiter(scalarField &allLambda, const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phiBD, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su, const scalar psiMax, const scalar psiMin)
void surfaceIntegrate(Field< Type > &ivf, const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Ostream & endl(Ostream &os)
Add newline and flush stream.
UList< label > labelUList
GeometricField< scalar, fvPatchField, volMesh > volScalarField
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
fvPatchField< scalar > fvPatchScalarField
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
List< label > labelList
A List of labels.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
volScalarField scalarField(fieldObject, mesh)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
MULES: Multidimensional universal limiter for explicit solution.
const volScalarField & psi
void limit(const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phi, surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su, const scalar psiMax, const scalar psiMin, const bool returnCorr)
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
SlicedGeometricField< scalar, fvsPatchField, slicedFvsPatchField, surfaceMesh > slicedSurfaceScalarField