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-2019 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 #include "syncTools.H"
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 template<class RdeltaTType, class RhoType, class SpType, class SuType>
37 (
38  const RdeltaTType& rDeltaT,
39  const RhoType& rho,
40  volScalarField& psi,
41  const surfaceScalarField& phiCorr,
42  const SpType& Sp,
43  const SuType& Su
44 )
45 {
46  Info<< "MULES: Correcting " << psi.name() << endl;
47 
48  const fvMesh& mesh = psi.mesh();
49 
50  scalarField psiIf(psi.size(), 0);
51  fvc::surfaceIntegrate(psiIf, phiCorr);
52 
53  if (mesh.moving())
54  {
55  psi.primitiveFieldRef() =
56  (
57  rho.field()*psi.primitiveField()*rDeltaT
58  + Su.field()
59  - psiIf
60  )/(rho.field()*rDeltaT - Sp.field());
61  }
62  else
63  {
64  psi.primitiveFieldRef() =
65  (
66  rho.field()*psi.primitiveField()*rDeltaT
67  + Su.field()
68  - psiIf
69  )/(rho.field()*rDeltaT - Sp.field());
70  }
71 
73 }
74 
75 
76 template<class RhoType>
78 (
79  const RhoType& rho,
80  volScalarField& psi,
81  const surfaceScalarField& phiCorr
82 )
83 {
84  correct(rho, psi, phiCorr, zeroField(), zeroField());
85 }
86 
87 
88 template<class RhoType, class SpType, class SuType>
90 (
91  const RhoType& rho,
92  volScalarField& psi,
93  const surfaceScalarField& phiCorr,
94  const SpType& Sp,
95  const SuType& Su
96 )
97 {
98  const fvMesh& mesh = psi.mesh();
99 
100  if (fv::localEulerDdt::enabled(mesh))
101  {
102  const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
103  correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
104  }
105  else
106  {
107  const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
108  correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
109  }
110 }
111 
112 
113 template<class RhoType, class PsiMaxType, class PsiMinType>
115 (
116  const RhoType& rho,
117  volScalarField& psi,
118  const surfaceScalarField& phi,
119  surfaceScalarField& phiCorr,
120  const PsiMaxType& psiMax,
121  const PsiMinType& psiMin
122 )
123 {
124  correct(rho, psi, phi, phiCorr, zeroField(), zeroField(), psiMax, psiMin);
125 }
126 
127 
128 template
129 <
130  class RhoType,
131  class SpType,
132  class SuType,
133  class PsiMaxType,
134  class PsiMinType
135 >
137 (
138  const RhoType& rho,
139  volScalarField& psi,
140  const surfaceScalarField& phi,
141  surfaceScalarField& phiCorr,
142  const SpType& Sp,
143  const SuType& Su,
144  const PsiMaxType& psiMax,
145  const PsiMinType& psiMin
146 )
147 {
148  const fvMesh& mesh = psi.mesh();
149 
150  if (fv::localEulerDdt::enabled(mesh))
151  {
152  const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
153 
154  limitCorr
155  (
156  rDeltaT,
157  rho,
158  psi,
159  phi,
160  phiCorr,
161  Sp,
162  Su,
163  psiMax,
164  psiMin
165  );
166 
167  correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
168  }
169  else
170  {
171  const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
172 
173  limitCorr
174  (
175  rDeltaT,
176  rho,
177  psi,
178  phi,
179  phiCorr,
180  Sp,
181  Su,
182  psiMax,
183  psiMin
184  );
185 
186  correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
187  }
188 }
189 
190 
191 template
192 <
193  class RdeltaTType,
194  class RhoType,
195  class SpType,
196  class SuType,
197  class PsiMaxType,
198  class PsiMinType
199 >
201 (
202  scalarField& allLambda,
203  const RdeltaTType& rDeltaT,
204  const RhoType& rho,
205  const volScalarField& psi,
206  const surfaceScalarField& phi,
207  const surfaceScalarField& phiCorr,
208  const SpType& Sp,
209  const SuType& Su,
210  const PsiMaxType& psiMax,
211  const PsiMinType& psiMin
212 )
213 {
214  const scalarField& psiIf = psi;
215  const volScalarField::Boundary& psiBf = psi.boundaryField();
216 
217  const fvMesh& mesh = psi.mesh();
218 
219  const dictionary& MULEScontrols = mesh.solverDict(psi.name());
220 
221  const label nLimiterIter
222  (
223  MULEScontrols.lookup<label>("nLimiterIter")
224  );
225 
226  const scalar smoothLimiter
227  (
228  MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
229  );
230 
231  const scalar extremaCoeff
232  (
233  MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0)
234  );
235 
236  const scalar boundaryExtremaCoeff
237  (
238  MULEScontrols.lookupOrDefault<scalar>
239  (
240  "boundaryExtremaCoeff",
241  extremaCoeff
242  )
243  );
244 
245  const scalar boundaryDeltaExtremaCoeff
246  (
247  max(boundaryExtremaCoeff - extremaCoeff, 0)
248  );
249 
250  const labelUList& owner = mesh.owner();
251  const labelUList& neighb = mesh.neighbour();
252  tmp<volScalarField::Internal> tVsc = mesh.Vsc();
253  const scalarField& V = tVsc();
254 
255  const surfaceScalarField::Boundary& phiBf =
256  phi.boundaryField();
257 
258  const scalarField& phiCorrIf = phiCorr;
259  const surfaceScalarField::Boundary& phiCorrBf =
260  phiCorr.boundaryField();
261 
263  (
264  IOobject
265  (
266  "lambda",
267  mesh.time().timeName(),
268  mesh,
269  IOobject::NO_READ,
270  IOobject::NO_WRITE,
271  false
272  ),
273  mesh,
274  dimless,
275  allLambda,
276  false // Use slices for the couples
277  );
278 
279  scalarField& lambdaIf = lambda;
280  surfaceScalarField::Boundary& lambdaBf =
281  lambda.boundaryFieldRef();
282 
283  scalarField psiMaxn(psiIf.size());
284  scalarField psiMinn(psiIf.size());
285 
286  psiMaxn = psiMin;
287  psiMinn = psiMax;
288 
289  scalarField sumPhip(psiIf.size(), 0.0);
290  scalarField mSumPhim(psiIf.size(), 0.0);
291 
292  forAll(phiCorrIf, facei)
293  {
294  const label own = owner[facei];
295  const label nei = neighb[facei];
296 
297  psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
298  psiMinn[own] = min(psiMinn[own], psiIf[nei]);
299 
300  psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
301  psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
302 
303  const scalar phiCorrf = phiCorrIf[facei];
304 
305  if (phiCorrf > 0)
306  {
307  sumPhip[own] += phiCorrf;
308  mSumPhim[nei] += phiCorrf;
309  }
310  else
311  {
312  mSumPhim[own] -= phiCorrf;
313  sumPhip[nei] -= phiCorrf;
314  }
315  }
316 
317  forAll(phiCorrBf, patchi)
318  {
319  const fvPatchScalarField& psiPf = psiBf[patchi];
320  const scalarField& phiCorrPf = phiCorrBf[patchi];
321 
322  const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
323 
324  if (psiPf.coupled())
325  {
326  const scalarField psiPNf(psiPf.patchNeighbourField());
327 
328  forAll(phiCorrPf, pFacei)
329  {
330  label pfCelli = pFaceCells[pFacei];
331 
332  psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
333  psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
334  }
335  }
336  else if (psiPf.fixesValue())
337  {
338  forAll(phiCorrPf, pFacei)
339  {
340  const label pfCelli = pFaceCells[pFacei];
341 
342  psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
343  psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
344  }
345  }
346  else
347  {
348  // Add the optional additional allowed boundary extrema
349  if (boundaryDeltaExtremaCoeff > 0)
350  {
351  forAll(phiCorrPf, pFacei)
352  {
353  const label pfCelli = pFaceCells[pFacei];
354 
355  const scalar extrema =
356  boundaryDeltaExtremaCoeff
357  *(psiMax[pfCelli] - psiMin[pfCelli]);
358 
359  psiMaxn[pfCelli] += extrema;
360  psiMinn[pfCelli] -= extrema;
361  }
362  }
363  }
364 
365  forAll(phiCorrPf, pFacei)
366  {
367  const label pfCelli = pFaceCells[pFacei];
368 
369  const scalar phiCorrf = phiCorrPf[pFacei];
370 
371  if (phiCorrf > 0)
372  {
373  sumPhip[pfCelli] += phiCorrf;
374  }
375  else
376  {
377  mSumPhim[pfCelli] -= phiCorrf;
378  }
379  }
380  }
381 
382  psiMaxn = min(psiMaxn + extremaCoeff*(psiMax - psiMin), psiMax);
383  psiMinn = max(psiMinn - extremaCoeff*(psiMax - psiMin), psiMin);
384 
385  if (smoothLimiter > small)
386  {
387  psiMaxn =
388  min(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMaxn, psiMax);
389  psiMinn =
390  max(smoothLimiter*psiIf + (1.0 - smoothLimiter)*psiMinn, psiMin);
391  }
392 
393  psiMaxn =
394  V
395  *(
396  (rho.field()*rDeltaT - Sp.field())*psiMaxn
397  - Su.field()
398  - rho.field()*psi.primitiveField()*rDeltaT
399  );
400 
401  psiMinn =
402  V
403  *(
404  Su.field()
405  - (rho.field()*rDeltaT - Sp.field())*psiMinn
406  + rho.field()*psi.primitiveField()*rDeltaT
407  );
408 
409  scalarField sumlPhip(psiIf.size());
410  scalarField mSumlPhim(psiIf.size());
411 
412  for (int j=0; j<nLimiterIter; j++)
413  {
414  sumlPhip = 0;
415  mSumlPhim = 0;
416 
417  forAll(lambdaIf, facei)
418  {
419  const label own = owner[facei];
420  const label nei = neighb[facei];
421 
422  const scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
423 
424  if (lambdaPhiCorrf > 0)
425  {
426  sumlPhip[own] += lambdaPhiCorrf;
427  mSumlPhim[nei] += lambdaPhiCorrf;
428  }
429  else
430  {
431  mSumlPhim[own] -= lambdaPhiCorrf;
432  sumlPhip[nei] -= lambdaPhiCorrf;
433  }
434  }
435 
436  forAll(lambdaBf, patchi)
437  {
438  scalarField& lambdaPf = lambdaBf[patchi];
439  const scalarField& phiCorrfPf = phiCorrBf[patchi];
440 
441  const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
442 
443  forAll(lambdaPf, pFacei)
444  {
445  label pfCelli = pFaceCells[pFacei];
446 
447  scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
448 
449  if (lambdaPhiCorrf > 0)
450  {
451  sumlPhip[pfCelli] += lambdaPhiCorrf;
452  }
453  else
454  {
455  mSumlPhim[pfCelli] -= lambdaPhiCorrf;
456  }
457  }
458  }
459 
460  forAll(sumlPhip, celli)
461  {
462  sumlPhip[celli] =
463  max(min
464  (
465  (sumlPhip[celli] + psiMaxn[celli])
466  /(mSumPhim[celli] + rootVSmall),
467  1.0), 0.0
468  );
469 
470  mSumlPhim[celli] =
471  max(min
472  (
473  (mSumlPhim[celli] + psiMinn[celli])
474  /(sumPhip[celli] + rootVSmall),
475  1.0), 0.0
476  );
477  }
478 
479  const scalarField& lambdam = sumlPhip;
480  const scalarField& lambdap = mSumlPhim;
481 
482  forAll(lambdaIf, facei)
483  {
484  if (phiCorrIf[facei] > 0)
485  {
486  lambdaIf[facei] = min
487  (
488  lambdaIf[facei],
489  min(lambdap[owner[facei]], lambdam[neighb[facei]])
490  );
491  }
492  else
493  {
494  lambdaIf[facei] = min
495  (
496  lambdaIf[facei],
497  min(lambdam[owner[facei]], lambdap[neighb[facei]])
498  );
499  }
500  }
501 
502 
503  forAll(lambdaBf, patchi)
504  {
505  fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
506  const scalarField& phiCorrfPf = phiCorrBf[patchi];
507  const fvPatchScalarField& psiPf = psiBf[patchi];
508 
509  if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
510  {
511  lambdaPf = 0;
512  }
513  else if (psiPf.coupled())
514  {
515  const labelList& pFaceCells =
516  mesh.boundary()[patchi].faceCells();
517 
518  forAll(lambdaPf, pFacei)
519  {
520  const label pfCelli = pFaceCells[pFacei];
521 
522  if (phiCorrfPf[pFacei] > 0)
523  {
524  lambdaPf[pFacei] =
525  min(lambdaPf[pFacei], lambdap[pfCelli]);
526  }
527  else
528  {
529  lambdaPf[pFacei] =
530  min(lambdaPf[pFacei], lambdam[pfCelli]);
531  }
532  }
533  }
534  else
535  {
536  const labelList& pFaceCells =
537  mesh.boundary()[patchi].faceCells();
538  const scalarField& phiPf = phiBf[patchi];
539 
540  forAll(lambdaPf, pFacei)
541  {
542  // Limit outlet faces only
543  if ((phiPf[pFacei] + phiCorrfPf[pFacei]) > small*small)
544  {
545  const label pfCelli = pFaceCells[pFacei];
546 
547  if (phiCorrfPf[pFacei] > 0)
548  {
549  lambdaPf[pFacei] =
550  min(lambdaPf[pFacei], lambdap[pfCelli]);
551  }
552  else
553  {
554  lambdaPf[pFacei] =
555  min(lambdaPf[pFacei], lambdam[pfCelli]);
556  }
557  }
558  }
559  }
560  }
561 
562  syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
563  }
564 }
565 
566 
567 template
568 <
569  class RdeltaTType,
570  class RhoType,
571  class SpType,
572  class SuType,
573  class PsiMaxType,
574  class PsiMinType
575 >
577 (
578  const RdeltaTType& rDeltaT,
579  const RhoType& rho,
580  const volScalarField& psi,
581  const surfaceScalarField& phi,
582  surfaceScalarField& phiCorr,
583  const SpType& Sp,
584  const SuType& Su,
585  const PsiMaxType& psiMax,
586  const PsiMinType& psiMin
587 )
588 {
589  const fvMesh& mesh = psi.mesh();
590 
591  scalarField allLambda(mesh.nFaces(), 1.0);
592 
594  (
595  IOobject
596  (
597  "lambda",
598  mesh.time().timeName(),
599  mesh,
600  IOobject::NO_READ,
601  IOobject::NO_WRITE,
602  false
603  ),
604  mesh,
605  dimless,
606  allLambda,
607  false // Use slices for the couples
608  );
609 
611  (
612  allLambda,
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 // ************************************************************************* //
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField...
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
const word & name() const
Return name.
Definition: IOobject.H:303
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
void surfaceIntegrate(Field< Type > &ivf, const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
bool moving() const
Is mesh moving.
Definition: polyMesh.H:512
virtual bool fixesValue() const
Return true if this patch field fixes a value.
Definition: fvPatchField.H:296
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const Boundary & boundaryField() const
Return const-reference to the boundary field.
label nFaces() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:636
const dimensionSet dimless
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
const dictionary & solverDict(const word &name) const
Return the solver controls dictionary for the given field.
Definition: solution.C:342
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:284
void correct(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su)
tmp< DimensionedField< scalar, volMesh > > Vsc() const
Return sub-cycle cell volumes.
dynamicFvMesh & mesh
A class representing the concept of a field of 0 used to avoid unnecessary manipulations for objects ...
Definition: zeroField.H:50
virtual tmp< Field< Type > > patchNeighbourField() const
Return patchField on the opposite patch of a coupled patch.
Definition: fvPatchField.H:411
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-9/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:68
scalar deltaTValue() const
Return time step value.
Definition: TimeStateI.H:41
virtual bool coupled() const
Return true if this patch field is coupled.
Definition: fvPatchField.H:309
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
const Mesh & mesh() const
Return mesh.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:278
void limiterCorr(scalarField &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)
label patchi
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
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 correctBoundaryConditions()
Correct boundary field.
messageStream Info
const volScalarField & psi
A class for managing temporary objects.
Definition: PtrList.H:53
Specialisation of GeometricField which holds slices of given complete fields in a form that they act ...
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:65
CMULES: Multidimensional universal limiter for explicit corrected implicit solution.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:540