CMULESTemplates.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2024 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "CMULES.H"
27 #include "fvcSurfaceIntegrate.H"
28 #include "localEulerDdtScheme.H"
29 #include "slicedSurfaceFields.H"
30 #include "wedgeFvPatch.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 template<class RdeltaTType, class RhoType, class SpType, class SuType>
36 (
37  const RdeltaTType& rDeltaT,
38  const RhoType& rho,
40  const surfaceScalarField& phiCorr,
41  const SpType& Sp,
42  const SuType& Su
43 )
44 {
45  Info<< "MULES: Correcting " << psi.name() << endl;
46 
47  const fvMesh& mesh = psi.mesh();
48 
49  scalarField psiIf(psi.size(), 0);
50  fvc::surfaceIntegrate(psiIf, phiCorr);
51 
52  if (mesh.moving())
53  {
54  psi.primitiveFieldRef() =
55  (
56  rho.primitiveField()*psi.primitiveField()*rDeltaT
57  + Su.primitiveField()
58  - psiIf
59  )/(rho.primitiveField()*rDeltaT - Sp.primitiveField());
60  }
61  else
62  {
63  psi.primitiveFieldRef() =
64  (
65  rho.primitiveField()*psi.primitiveField()*rDeltaT
66  + Su.primitiveField()
67  - psiIf
68  )/(rho.primitiveField()*rDeltaT - Sp.primitiveField());
69  }
70 
71  psi.correctBoundaryConditions();
72 }
73 
74 
75 template<class RhoType>
77 (
78  const RhoType& rho,
80  const surfaceScalarField& phiCorr
81 )
82 {
83  correct(rho, psi, phiCorr, zeroField(), zeroField());
84 }
85 
86 
87 template<class RhoType, class SpType, class SuType>
89 (
90  const RhoType& rho,
92  const surfaceScalarField& phiCorr,
93  const SpType& Sp,
94  const SuType& Su
95 )
96 {
97  const fvMesh& mesh = psi.mesh();
98 
99  if (fv::localEulerDdt::enabled(mesh))
100  {
101  const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
102  correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
103  }
104  else
105  {
106  const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
107  correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
108  }
109 }
110 
111 
112 template<class RhoType, class PsiMaxType, class PsiMinType>
114 (
115  const RhoType& rho,
117  const surfaceScalarField& phi,
118  surfaceScalarField& phiCorr,
119  const PsiMaxType& psiMax,
120  const PsiMinType& psiMin
121 )
122 {
123  correct(rho, psi, phi, phiCorr, zeroField(), zeroField(), psiMax, psiMin);
124 }
125 
126 
127 template
128 <
129  class RhoType,
130  class SpType,
131  class SuType,
132  class PsiMaxType,
133  class PsiMinType
134 >
136 (
137  const RhoType& rho,
139  const surfaceScalarField& phi,
140  surfaceScalarField& phiCorr,
141  const SpType& Sp,
142  const SuType& Su,
143  const PsiMaxType& psiMax,
144  const PsiMinType& psiMin
145 )
146 {
147  const fvMesh& mesh = psi.mesh();
148 
149  if (fv::localEulerDdt::enabled(mesh))
150  {
151  const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
152 
153  limitCorr
154  (
155  rDeltaT,
156  rho,
157  psi,
158  phi,
159  phiCorr,
160  Sp,
161  Su,
162  psiMax,
163  psiMin
164  );
165 
166  correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
167  }
168  else
169  {
170  const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
171 
172  limitCorr
173  (
174  rDeltaT,
175  rho,
176  psi,
177  phi,
178  phiCorr,
179  Sp,
180  Su,
181  psiMax,
182  psiMin
183  );
184 
185  correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
186  }
187 }
188 
189 
190 template
191 <
192  class RdeltaTType,
193  class RhoType,
194  class SpType,
195  class SuType,
196  class PsiMaxType,
197  class PsiMinType
198 >
200 (
202  const RdeltaTType& rDeltaT,
203  const RhoType& rho,
204  const volScalarField& psi,
205  const surfaceScalarField& phi,
206  const surfaceScalarField& phiCorr,
207  const SpType& Sp,
208  const SuType& Su,
209  const PsiMaxType& psiMax,
210  const PsiMinType& psiMin
211 )
212 {
213  const scalarField& psiIf = psi;
214  const volScalarField::Boundary& psiBf = psi.boundaryField();
215 
216  const fvMesh& mesh = psi.mesh();
217 
218  const dictionary& MULEScontrols = mesh.solution().solverDict(psi.name());
219 
220  const label nLimiterIter
221  (
222  MULEScontrols.lookup<label>("nLimiterIter")
223  );
224 
225  const scalar smoothLimiter
226  (
227  MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
228  );
229 
230  const scalar extremaCoeff
231  (
232  MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0)
233  );
234 
235  const scalar boundaryExtremaCoeff
236  (
237  MULEScontrols.lookupOrDefault<scalar>
238  (
239  "boundaryExtremaCoeff",
240  extremaCoeff
241  )
242  );
243 
244  const scalar boundaryDeltaExtremaCoeff
245  (
246  max(boundaryExtremaCoeff - extremaCoeff, 0)
247  );
248 
249  const labelUList& owner = mesh.owner();
250  const labelUList& neighb = mesh.neighbour();
251  tmp<volScalarField::Internal> tVsc = mesh.Vsc();
252  const scalarField& V = tVsc();
253 
254  const surfaceScalarField::Boundary& phiBf =
255  phi.boundaryField();
256 
257  const scalarField& phiCorrIf = phiCorr;
258  const surfaceScalarField::Boundary& phiCorrBf =
259  phiCorr.boundaryField();
260 
261  scalarField& lambdaIf = lambda;
262  surfaceScalarField::Boundary& lambdaBf =
263  lambda.boundaryFieldRef();
264 
265  scalarField psiMaxn(psiIf.size());
266  scalarField psiMinn(psiIf.size());
267 
268  psiMaxn = psiMin;
269  psiMinn = psiMax;
270 
271  scalarField sumPhip(psiIf.size(), 0.0);
272  scalarField mSumPhim(psiIf.size(), 0.0);
273 
274  forAll(phiCorrIf, facei)
275  {
276  const label own = owner[facei];
277  const label nei = neighb[facei];
278 
279  psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
280  psiMinn[own] = min(psiMinn[own], psiIf[nei]);
281 
282  psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
283  psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
284 
285  const scalar phiCorrf = phiCorrIf[facei];
286 
287  if (phiCorrf > 0)
288  {
289  sumPhip[own] += phiCorrf;
290  mSumPhim[nei] += phiCorrf;
291  }
292  else
293  {
294  mSumPhim[own] -= phiCorrf;
295  sumPhip[nei] -= phiCorrf;
296  }
297  }
298 
299  forAll(phiCorrBf, patchi)
300  {
301  const fvPatchScalarField& psiPf = psiBf[patchi];
302  const scalarField& phiCorrPf = phiCorrBf[patchi];
303 
304  const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
305 
306  if (psiPf.coupled())
307  {
308  const scalarField psiPNf(psiPf.patchNeighbourField());
309 
310  forAll(phiCorrPf, pFacei)
311  {
312  label pfCelli = pFaceCells[pFacei];
313 
314  psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
315  psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
316  }
317  }
318  else if (psiPf.fixesValue())
319  {
320  forAll(phiCorrPf, pFacei)
321  {
322  const label pfCelli = pFaceCells[pFacei];
323 
324  psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
325  psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
326  }
327  }
328  else
329  {
330  // Add the optional additional allowed boundary extrema
331  if (boundaryDeltaExtremaCoeff > 0)
332  {
333  forAll(phiCorrPf, pFacei)
334  {
335  const label pfCelli = pFaceCells[pFacei];
336 
337  const scalar extrema =
338  boundaryDeltaExtremaCoeff
339  *(psiMax[pfCelli] - psiMin[pfCelli]);
340 
341  psiMaxn[pfCelli] += extrema;
342  psiMinn[pfCelli] -= extrema;
343  }
344  }
345  }
346 
347  forAll(phiCorrPf, pFacei)
348  {
349  const label pfCelli = pFaceCells[pFacei];
350 
351  const scalar phiCorrf = phiCorrPf[pFacei];
352 
353  if (phiCorrf > 0)
354  {
355  sumPhip[pfCelli] += phiCorrf;
356  }
357  else
358  {
359  mSumPhim[pfCelli] -= phiCorrf;
360  }
361  }
362  }
363 
364  if (extremaCoeff > 0)
365  {
366  psiMaxn = min(psiMaxn + extremaCoeff*(psiMax - psiMin), psiMax);
367  psiMinn = max(psiMinn - extremaCoeff*(psiMax - psiMin), psiMin);
368  }
369  else
370  {
371  psiMaxn = min(psiMaxn, psiMax);
372  psiMinn = max(psiMinn, psiMin);
373  }
374 
375  if (smoothLimiter > small)
376  {
377  psiMaxn =
378  min(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMaxn, psiMax);
379  psiMinn =
380  max(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMinn, psiMin);
381  }
382 
383  psiMaxn =
384  V
385  *(
386  (rho.primitiveField()*rDeltaT - Sp.primitiveField())*psiMaxn
387  - Su.primitiveField()
388  - rho.primitiveField()*psi.primitiveField()*rDeltaT
389  );
390 
391  psiMinn =
392  V
393  *(
394  Su.primitiveField()
395  - (rho.primitiveField()*rDeltaT - Sp.primitiveField())*psiMinn
396  + rho.primitiveField()*psi.primitiveField()*rDeltaT
397  );
398 
399  scalarField sumlPhip(psiIf.size());
400  scalarField mSumlPhim(psiIf.size());
401 
402  for (int j=0; j<nLimiterIter; j++)
403  {
404  sumlPhip = 0;
405  mSumlPhim = 0;
406 
407  forAll(lambdaIf, facei)
408  {
409  const label own = owner[facei];
410  const label nei = neighb[facei];
411 
412  const scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
413 
414  if (lambdaPhiCorrf > 0)
415  {
416  sumlPhip[own] += lambdaPhiCorrf;
417  mSumlPhim[nei] += lambdaPhiCorrf;
418  }
419  else
420  {
421  mSumlPhim[own] -= lambdaPhiCorrf;
422  sumlPhip[nei] -= lambdaPhiCorrf;
423  }
424  }
425 
426  forAll(lambdaBf, patchi)
427  {
428  scalarField& lambdaPf = lambdaBf[patchi];
429  const scalarField& phiCorrfPf = phiCorrBf[patchi];
430 
431  const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
432 
433  forAll(lambdaPf, pFacei)
434  {
435  label pfCelli = pFaceCells[pFacei];
436 
437  scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
438 
439  if (lambdaPhiCorrf > 0)
440  {
441  sumlPhip[pfCelli] += lambdaPhiCorrf;
442  }
443  else
444  {
445  mSumlPhim[pfCelli] -= lambdaPhiCorrf;
446  }
447  }
448  }
449 
450  forAll(sumlPhip, celli)
451  {
452  sumlPhip[celli] =
453  max(min
454  (
455  (sumlPhip[celli] + psiMaxn[celli])
456  /(mSumPhim[celli] + rootVSmall),
457  1.0), 0.0
458  );
459 
460  mSumlPhim[celli] =
461  max(min
462  (
463  (mSumlPhim[celli] + psiMinn[celli])
464  /(sumPhip[celli] + rootVSmall),
465  1.0), 0.0
466  );
467  }
468 
469  const scalarField& lambdam = sumlPhip;
470  const scalarField& lambdap = mSumlPhim;
471 
472  forAll(lambdaIf, facei)
473  {
474  if (phiCorrIf[facei] > 0)
475  {
476  lambdaIf[facei] = min
477  (
478  lambdaIf[facei],
479  min(lambdap[owner[facei]], lambdam[neighb[facei]])
480  );
481  }
482  else
483  {
484  lambdaIf[facei] = min
485  (
486  lambdaIf[facei],
487  min(lambdam[owner[facei]], lambdap[neighb[facei]])
488  );
489  }
490  }
491 
492 
493  forAll(lambdaBf, patchi)
494  {
495  fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
496  const scalarField& phiCorrfPf = phiCorrBf[patchi];
497  const fvPatchScalarField& psiPf = psiBf[patchi];
498 
499  if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
500  {
501  lambdaPf = 0;
502  }
503  else if (psiPf.coupled())
504  {
505  const labelList& pFaceCells =
506  mesh.boundary()[patchi].faceCells();
507 
508  forAll(lambdaPf, pFacei)
509  {
510  const label pfCelli = pFaceCells[pFacei];
511 
512  if (phiCorrfPf[pFacei] > 0)
513  {
514  lambdaPf[pFacei] =
515  min(lambdaPf[pFacei], lambdap[pfCelli]);
516  }
517  else
518  {
519  lambdaPf[pFacei] =
520  min(lambdaPf[pFacei], lambdam[pfCelli]);
521  }
522  }
523  }
524  else
525  {
526  const labelList& pFaceCells =
527  mesh.boundary()[patchi].faceCells();
528  const scalarField& phiPf = phiBf[patchi];
529 
530  forAll(lambdaPf, pFacei)
531  {
532  // Limit outlet faces only
533  if ((phiPf[pFacei] + phiCorrfPf[pFacei]) > small*small)
534  {
535  const label pfCelli = pFaceCells[pFacei];
536 
537  if (phiCorrfPf[pFacei] > 0)
538  {
539  lambdaPf[pFacei] =
540  min(lambdaPf[pFacei], lambdap[pfCelli]);
541  }
542  else
543  {
544  lambdaPf[pFacei] =
545  min(lambdaPf[pFacei], lambdam[pfCelli]);
546  }
547  }
548  }
549  }
550  }
551 
552  // Take minimum of value across coupled patches
553  surfaceScalarField::Boundary lambdaNbrBf
554  (
555  surfaceScalarField::Internal::null(),
556  lambdaBf.boundaryNeighbourField()
557  );
558  forAll(lambdaBf, patchi)
559  {
560  fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
561  const fvsPatchScalarField& lambdaNbrPf = lambdaNbrBf[patchi];
562  if (lambdaPf.coupled())
563  {
564  lambdaPf = min(lambdaPf, lambdaNbrPf);
565  }
566  }
567  }
568 }
569 
570 
571 template
572 <
573  class RdeltaTType,
574  class RhoType,
575  class SpType,
576  class SuType,
577  class PsiMaxType,
578  class PsiMinType
579 >
581 (
582  const RdeltaTType& rDeltaT,
583  const RhoType& rho,
584  const volScalarField& psi,
585  const surfaceScalarField& phi,
586  surfaceScalarField& phiCorr,
587  const SpType& Sp,
588  const SuType& Su,
589  const PsiMaxType& psiMax,
590  const PsiMinType& psiMin
591 )
592 {
593  const fvMesh& mesh = psi.mesh();
594 
596  (
597  IOobject
598  (
599  "lambda",
600  mesh.time().name(),
601  mesh,
602  IOobject::NO_READ,
603  IOobject::NO_WRITE,
604  false
605  ),
606  mesh,
608  );
609 
611  (
612  lambda,
613  rDeltaT,
614  rho,
615  psi,
616  phi,
617  phiCorr,
618  Sp,
619  Su,
620  psiMax,
621  psiMin
622  );
623 
624  phiCorr *= lambda;
625 }
626 
627 
628 // ************************************************************************* //
CMULES: Multidimensional universal limiter for explicit corrected implicit solution.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Generic GeometricBoundaryField class.
tmp< GeometricBoundaryField > boundaryNeighbourField() const
Return BoundaryField of the values on the other side of couples.
Generic GeometricField class.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
scalar deltaTValue() const
Return time step value.
Definition: TimeStateI.H:34
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T, if not found return the given default.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:710
const word & name() const
Return const reference to name.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:99
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:418
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:469
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:857
const fvSolution & solution() const
Return the fvSolution.
Definition: fvMesh.C:1762
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:475
tmp< DimensionedField< scalar, volMesh > > Vsc() const
Return sub-cycle cell volumes.
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:88
virtual bool fixesValue() const
Return true if this patch field fixes a value.
Definition: fvPatchField.H:315
virtual bool coupled() const
Return true if this patch field is coupled.
Definition: fvPatchField.H:328
virtual tmp< Field< Type > > patchNeighbourField(const Pstream::commsTypes commsType=Pstream::commsTypes::blocking) const
Return patchField on the opposite patch of a coupled patch.
Definition: fvPatchField.H:433
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:83
virtual bool coupled() const
Return true if this patch field is coupled.
bool moving() const
Is mesh moving.
Definition: polyMesh.H:476
const dictionary & solverDict(const word &name) const
Return the solver controls dictionary for the given field.
Definition: solution.C:348
A class for managing temporary objects.
Definition: tmp.H:55
A class representing the concept of a field of 0 used to avoid unnecessary manipulations for objects ...
Definition: zeroField.H:53
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField.
label patchi
const volScalarField & psi
dimensionedScalar lambda(viscosity->lookup("lambda"))
void limiterCorr(surfaceScalarField &allLambda, const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phi, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su, const PsiMaxType &psiMax, const PsiMinType &psiMin)
void limitCorr(const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phi, surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su, const PsiMaxType &psiMax, const PsiMinType &psiMin)
void correct(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su)
void surfaceIntegrate(Field< Type > &ivf, const SurfaceField< Type > &ssf)
tmp< VolField< Type > > Su(const VolField< Type > &su, const VolField< Type > &vf)
Definition: fvcSup.C:44
tmp< VolField< Type > > Sp(const volScalarField &sp, const VolField< Type > &vf)
Definition: fvcSup.C:67
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
const dimensionSet dimless
messageStream Info
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)