All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
phaseSystemSolve.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-2020 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 "phaseSystem.H"
27 
28 #include "MULES.H"
29 #include "subCycle.H"
30 
31 #include "fvcDdt.H"
32 #include "fvcDiv.H"
33 #include "fvcSnGrad.H"
34 #include "fvcFlux.H"
35 #include "fvcMeshPhi.H"
36 #include "fvcSup.H"
37 
38 #include "fvmDdt.H"
39 #include "fvmLaplacian.H"
40 #include "fvmSup.H"
41 
42 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
43 
45 (
46  const PtrList<volScalarField>& rAUs,
47  const PtrList<surfaceScalarField>& rAUfs
48 )
49 {
50  const dictionary& alphaControls = mesh_.solverDict("alpha");
51 
52  const label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
53  const label nAlphaCorr(alphaControls.lookup<label>("nAlphaCorr"));
54 
55  const bool LTS = fv::localEulerDdt::enabled(mesh_);
56 
57  // Temporary switch for testing and comparing the standard split
58  // and the new un-split phase flux discretisation
59  const bool splitPhaseFlux
60  (
61  alphaControls.lookupOrDefault<Switch>("splitPhaseFlux", false)
62  );
63 
64  // Temporary switch for testing and comparing the standard mean flux
65  // and the new phase flux reference for the phase flux correction
66  const bool meanFluxReference
67  (
68  alphaControls.lookupOrDefault<Switch>("meanFluxReference", false)
69  );
70 
71  // Optional reference phase which is not solved for
72  // but obtained from the sum of the other phases
73  phaseModel* referencePhasePtr = nullptr;
74 
75  // The phases which are solved
76  // i.e. the moving phases less the optional reference phase
77  phaseModelPartialList solvePhases;
78 
80  {
81  referencePhasePtr = &phases()[referencePhaseName_];
82 
83  solvePhases.setSize(movingPhases().size() - 1);
84  label solvePhasesi = 0;
85  forAll(movingPhases(), movingPhasei)
86  {
87  if (&movingPhases()[movingPhasei] != referencePhasePtr)
88  {
89  solvePhases.set(solvePhasesi++, &movingPhases()[movingPhasei]);
90  }
91  }
92  }
93  else
94  {
95  solvePhases = movingPhases();
96  }
97 
98  forAll(phases(), phasei)
99  {
100  phases()[phasei].correctBoundaryConditions();
101  }
102 
103  PtrList<surfaceScalarField> alphaPhiDbyA0s(phases().size());
104  if (implicitPhasePressure() && (rAUs.size() || rAUfs.size()))
105  {
106  const PtrList<surfaceScalarField> DByAfs(this->DByAfs(rAUs, rAUfs));
107 
108  forAll(solvePhases, solvePhasei)
109  {
110  const phaseModel& phase = solvePhases[solvePhasei];
111  const volScalarField& alpha = phase;
112 
113  alphaPhiDbyA0s.set
114  (
115  phase.index(),
116  DByAfs[phase.index()]
117  *fvc::snGrad(alpha, "bounded")*mesh_.magSf()
118  );
119  }
120  }
121 
122  // Calculate the void fraction
123  volScalarField alphaVoid
124  (
125  IOobject
126  (
127  "alphaVoid",
128  mesh_.time().timeName(),
129  mesh_
130  ),
131  mesh_,
133  );
134  forAll(stationaryPhases(), stationaryPhasei)
135  {
136  alphaVoid -= stationaryPhases()[stationaryPhasei];
137  }
138 
139  // Calculate the effective flux of the moving phases
140  tmp<surfaceScalarField> tphiMoving(phi_);
141  if (stationaryPhases().size())
142  {
143  tphiMoving = phi_/upwind<scalar>(mesh_, phi_).interpolate(alphaVoid);
144  }
145  const surfaceScalarField& phiMoving = tphiMoving();
146 
147  bool dilatation = false;
148  forAll(movingPhases(), movingPhasei)
149  {
150  if (movingPhases()[movingPhasei].divU().valid())
151  {
152  dilatation = true;
153  break;
154  }
155  }
156 
157  for (int acorr=0; acorr<nAlphaCorr; acorr++)
158  {
159  PtrList<volScalarField::Internal> Sps(phases().size());
160  PtrList<volScalarField::Internal> Sus(phases().size());
161 
162  forAll(movingPhases(), movingPhasei)
163  {
164  const phaseModel& phase = movingPhases()[movingPhasei];
165  const volScalarField& alpha = phase;
166  const label phasei = phase.index();
167 
168  Sps.set
169  (
170  phasei,
172  (
173  IOobject
174  (
175  "Sp",
176  mesh_.time().timeName(),
177  mesh_
178  ),
179  mesh_,
181  )
182  );
183 
184  Sus.set
185  (
186  phasei,
188  (
189  "Su",
190  min(alpha.v(), scalar(1))
191  *fvc::div(fvc::absolute(phi_, phase.U()))->v()
192  )
193  );
194 
195  if (dilatation)
196  {
197  // Construct the dilatation rate source term
199  (
201  (
202  "dgdt",
203  mesh_,
205  )
206  );
207 
208  forAll(phases(), phasej)
209  {
210  const phaseModel& phase2 = phases()[phasej];
211  const volScalarField& alpha2 = phase2;
212 
213  if (&phase2 != &phase)
214  {
215  if (phase.divU().valid())
216  {
217  dgdt += alpha2()*phase.divU()()();
218  }
219 
220  if (phase2.divU().valid())
221  {
222  dgdt -= alpha()*phase2.divU()()();
223  }
224  }
225  }
226 
227  volScalarField::Internal& Sp = Sps[phasei];
228  volScalarField::Internal& Su = Sus[phasei];
229 
230  forAll(dgdt, celli)
231  {
232  if (dgdt[celli] > 0)
233  {
234  Sp[celli] -= dgdt[celli]/max(1 - alpha[celli], 1e-4);
235  Su[celli] += dgdt[celli]/max(1 - alpha[celli], 1e-4);
236  }
237  else if (dgdt[celli] < 0)
238  {
239  Sp[celli] += dgdt[celli]/max(alpha[celli], 1e-4);
240  }
241  }
242  }
243  }
244 
245  tmp<volScalarField> trSubDeltaT;
246 
247  if (LTS && nAlphaSubCycles > 1)
248  {
249  trSubDeltaT =
251  }
252 
253  List<volScalarField*> alphaPtrs(phases().size());
254  forAll(phases(), phasei)
255  {
256  alphaPtrs[phasei] = &phases()[phasei];
257  }
258 
259  for
260  (
261  subCycle<volScalarField, subCycleFields> alphaSubCycle
262  (
263  alphaPtrs,
265  );
266  !(++alphaSubCycle).end();
267  )
268  {
269  // Create correction fluxes
270  PtrList<surfaceScalarField> alphaPhis(phases().size());
271 
272  forAll(movingPhases(), movingPhasei)
273  {
274  const phaseModel& phase = movingPhases()[movingPhasei];
275  const volScalarField& alpha = phase;
276 
277  alphaPhis.set
278  (
279  phase.index(),
281  (
282  IOobject::groupName("alphaPhiCorr", phase.name()),
283  fvc::flux
284  (
285  splitPhaseFlux ? phi_ : phase.phi()(),
286  alpha,
287  "div(phi," + alpha.name() + ')'
288  )
289  )
290  );
291 
292  surfaceScalarField& alphaPhi = alphaPhis[phase.index()];
293 
294  if (splitPhaseFlux)
295  {
296  forAll(phases(), phasei)
297  {
298  const phaseModel& phase2 = phases()[phasei];
299  const volScalarField& alpha2 = phase2;
300 
301  if (&phase2 == &phase) continue;
302 
303  surfaceScalarField phir(phase.phi() - phase2.phi());
304 
306  (
307  cAlphas_.find
308  (
309  phasePairKey(phase.name(), phase2.name())
310  )
311  );
312 
313  if (cAlpha != cAlphas_.end())
314  {
316  (
317  (mag(phi_) + mag(phir))/mesh_.magSf()
318  );
319 
320  phir +=
321  min(cAlpha()*phic, max(phic))
322  *nHatf(alpha, alpha2);
323  }
324 
325  const word phirScheme
326  (
327  "div(phir,"
328  + alpha2.name() + ',' + alpha.name()
329  + ')'
330  );
331 
332  alphaPhi += fvc::flux
333  (
334  -fvc::flux(-phir, alpha2, phirScheme),
335  alpha,
336  phirScheme
337  );
338  }
339  }
340  else if (!cAlphas_.empty())
341  {
342  forAll(phases(), phasei)
343  {
344  const phaseModel& phase2 = phases()[phasei];
345  const volScalarField& alpha2 = phase2;
346 
347  if (&phase2 == &phase) continue;
348 
350  (
351  cAlphas_.find
352  (
353  phasePairKey(phase.name(), phase2.name())
354  )
355  );
356 
357  if (cAlpha != cAlphas_.end())
358  {
360  (
361  phase.phi() - phase2.phi()
362  );
363 
365  (
366  (mag(phi_) + mag(phir))/mesh_.magSf()
367  );
368 
369  const surfaceScalarField phirc
370  (
371  min(cAlpha()*phic, max(phic))
372  *nHatf(alpha, alpha2)
373  );
374 
375  const word phirScheme
376  (
377  "div(phir,"
378  + alpha2.name() + ',' + alpha.name()
379  + ')'
380  );
381 
382  alphaPhi += fvc::flux
383  (
384  -fvc::flux(-phirc, alpha2, phirScheme),
385  alpha,
386  phirScheme
387  );
388  }
389  }
390  }
391 
392  if (alphaPhiDbyA0s.set(phase.index()))
393  {
394  alphaPhi +=
395  fvc::interpolate(max(alpha, scalar(0)))
396  *fvc::interpolate(max(1 - alpha, scalar(0)))
397  *alphaPhiDbyA0s[phase.index()];
398  }
399 
400  phase.correctInflowOutflow(alphaPhi);
401 
403  (
404  geometricOneField(),
405  alpha,
406  meanFluxReference
407  ? phiMoving // Guarantees boundedness but less accurate
408  : phase.phi()(), // Less robust but more accurate
409  alphaPhi,
410  Sps[phase.index()],
411  Sus[phase.index()],
412  min(alphaVoid.primitiveField(), phase.alphaMax())(),
413  zeroField(),
414  false
415  );
416  }
417 
418  // Limit the flux corrections to ensure the phase fractions sum to 1
419  {
420  // Generate alphas for the moving phases
421  UPtrList<const volScalarField> alphasMoving
422  (
423  movingPhases().size()
424  );
425 
426  UPtrList<surfaceScalarField> alphaPhisMoving
427  (
428  movingPhases().size()
429  );
430 
431  forAll(movingPhases(), movingPhasei)
432  {
433  const phaseModel& phase = movingPhases()[movingPhasei];
434 
435  alphasMoving.set(movingPhasei, &phase);
436 
437  alphaPhisMoving.set
438  (
439  movingPhasei,
440  &alphaPhis[phase.index()]
441  );
442  }
443 
444  MULES::limitSum(alphasMoving, alphaPhisMoving, phiMoving);
445  }
446 
447  forAll(solvePhases, solvePhasei)
448  {
449  phaseModel& phase = solvePhases[solvePhasei];
450  volScalarField& alpha = phase;
451 
452  surfaceScalarField& alphaPhi = alphaPhis[phase.index()];
453  phase.correctInflowOutflow(alphaPhi);
454 
456  (
457  geometricOneField(),
458  alpha,
459  alphaPhi,
460  Sps[phase.index()],
461  Sus[phase.index()]
462  );
463 
464  if (alphaSubCycle.index() == 1)
465  {
466  phase.alphaPhiRef() = alphaPhi;
467  }
468  else
469  {
470  phase.alphaPhiRef() += alphaPhi;
471  }
472  }
473 
474  if (implicitPhasePressure() && (rAUs.size() || rAUfs.size()))
475  {
476  const PtrList<surfaceScalarField> DByAfs
477  (
478  this->DByAfs(rAUs, rAUfs)
479  );
480 
481  forAll(solvePhases, solvePhasei)
482  {
483  phaseModel& phase = solvePhases[solvePhasei];
484  volScalarField& alpha = phase;
485 
486  const surfaceScalarField alphaDbyA
487  (
488  fvc::interpolate(max(alpha, scalar(0)))
489  *fvc::interpolate(max(1 - alpha, scalar(0)))
490  *DByAfs[phase.index()]
491  );
492 
493  fvScalarMatrix alphaEqn
494  (
495  fvm::ddt(alpha) - fvc::ddt(alpha)
496  - fvm::laplacian(alphaDbyA, alpha, "bounded")
497  );
498 
499  alphaEqn.solve();
500 
501  phase.alphaPhiRef() += alphaEqn.flux();
502  }
503  }
504 
505  // Report the phase fractions and the phase fraction sum
506  forAll(solvePhases, solvePhasei)
507  {
508  phaseModel& phase = solvePhases[solvePhasei];
509 
510  Info<< phase.name() << " fraction, min, max = "
511  << phase.weightedAverage(mesh_.V()).value()
512  << ' ' << min(phase).value()
513  << ' ' << max(phase).value()
514  << endl;
515  }
516 
517  if (referencePhasePtr)
518  {
519  volScalarField& referenceAlpha = *referencePhasePtr;
520  referenceAlpha = alphaVoid;
521 
522  forAll(solvePhases, solvePhasei)
523  {
524  referenceAlpha -= solvePhases[solvePhasei];
525  }
526  }
527  else
528  {
530  (
531  IOobject
532  (
533  "sumAlphaMoving",
534  mesh_.time().timeName(),
535  mesh_
536  ),
537  mesh_,
539  );
540  forAll(movingPhases(), movingPhasei)
541  {
542  sumAlphaMoving += movingPhases()[movingPhasei];
543  }
544 
545  Info<< "Phase-sum volume fraction, min, max = "
546  << (sumAlphaMoving + 1 - alphaVoid)()
547  .weightedAverage(mesh_.V()).value()
548  << ' ' << min(sumAlphaMoving + 1 - alphaVoid).value()
549  << ' ' << max(sumAlphaMoving + 1 - alphaVoid).value()
550  << endl;
551 
552  // Correct the sum of the phase fractions to avoid drift
553  forAll(movingPhases(), movingPhasei)
554  {
555  movingPhases()[movingPhasei] *= alphaVoid/sumAlphaMoving;
556  }
557  }
558  }
559 
560  if (nAlphaSubCycles > 1)
561  {
562  forAll(solvePhases, solvePhasei)
563  {
564  phaseModel& phase = solvePhases[solvePhasei];
565 
566  phase.alphaPhiRef() /= nAlphaSubCycles;
567  }
568  }
569 
570  forAll(solvePhases, solvePhasei)
571  {
572  phaseModel& phase = solvePhases[solvePhasei];
573 
574  phase.alphaRhoPhiRef() =
575  fvc::interpolate(phase.rho())*phase.alphaPhi();
576 
577  phase.maxMin(0, 1);
578  }
579 
580  if (referencePhasePtr)
581  {
582  phaseModel& referencePhase = *referencePhasePtr;
583 
584  referencePhase.alphaPhiRef() = phi_;
585 
586  forAll(solvePhases, solvePhasei)
587  {
588  phaseModel& phase = solvePhases[solvePhasei];
589  referencePhase.alphaPhiRef() -= phase.alphaPhi();
590  }
591 
592  referencePhase.alphaRhoPhiRef() =
593  fvc::interpolate(referencePhase.rho())
594  *referencePhase.alphaPhi();
595 
596  volScalarField& referenceAlpha = referencePhase;
597  referenceAlpha = alphaVoid;
598 
599  forAll(solvePhases, solvePhasei)
600  {
601  referenceAlpha -= solvePhases[solvePhasei];
602  }
603  }
604  }
605 }
606 
607 
608 // ************************************************************************* //
void limit(const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phi, surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su, const PsiMaxType &psiMax, const PsiMinType &psiMin, const bool returnCorr)
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const phaseModelPartialList & movingPhases() const
Return the models for phases that are moving.
Definition: phaseSystemI.H:49
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
static tmp< DimensionedField< Type, GeoMesh > > New(const word &name, const Mesh &mesh, const dimensionSet &)
Return a temporary field constructed from name, mesh.
volScalarField dgdt(alpha1 *fvc::div(phi))
label phasei
Definition: pEqn.H:27
tmp< surfaceScalarField > nHatf(const volScalarField &alpha1, const volScalarField &alpha2) const
Normal to interface between two phases dotted with face areas.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
Definition: curveTools.C:75
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
static tmp< volScalarField > localRSubDeltaT(const fvMesh &mesh, const label nAlphaSubCycles)
Calculate and return the reciprocal of the local sub-cycling.
Definition: localEulerDdt.C:70
IOobject(const word &name, const fileName &instance, const objectRegistry &registry, readOption r=NO_READ, writeOption w=NO_WRITE, bool registerObject=true)
Construct from name, instance, registry, io options.
Definition: IOobject.C:167
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
bool empty() const
Return true if the hash table is empty.
Definition: HashTableI.H:72
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Calculate the matrix for the laplacian of the field.
Calculate the snGrad of the given volField.
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
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 DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
alpha2
Definition: alphaEqn.H:115
const phaseModelPartialList & stationaryPhases() const
Return the models for phases that are stationary.
Definition: phaseSystemI.H:63
virtual void solve(const PtrList< volScalarField > &rAUs, const PtrList< surfaceScalarField > &rAUfs)
Solve for the phase fractions.
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
dictionary()
Construct top-level dictionary null.
Definition: dictionary.C:440
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Calculate the first temporal derivative.
const phaseModelList & phases() const
Return the phase models.
Definition: phaseSystemI.H:35
const dimensionSet dimTime
static bool enabled(const fvMesh &mesh)
Return true if LTS is enabled.
Definition: localEulerDdt.C:37
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
Definition: MULES.C:30
phic
Definition: correctPhic.H:2
Calculate the face-flux of the given field.
static word groupName(Name name, const word &group)
friend class const_iterator
Declare friendship with the const_iterator.
Definition: HashTable.H:197
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
Calculate the matrix for the first temporal derivative.
word referencePhaseName_
Name of optional reference phase which is not solved for.
Definition: phaseSystem.H:148
static const word null
An empty word.
Definition: word.H:77
UPtrList< phaseModel > phaseModelPartialList
Definition: phaseSystem.H:86
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
const fvMesh & mesh_
Reference to the mesh.
Definition: phaseSystem.H:144
const iterator & end()
Definition: UILList.H:223
Calculate the field for explicit evaluation of implicit and explicit sources.
cAlphaTable cAlphas_
Interface compression coefficients.
Definition: phaseSystem.H:181
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
Calculate the divergence of the given field.
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
bool set(const label) const
Is element set.
Definition: UPtrListI.H:78
tmp< volScalarField > sumAlphaMoving() const
Return the sum of the phase fractions of the moving phases.
const label nAlphaSubCycles(alphaControls.lookup< label >("nAlphaSubCycles"))
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< volScalarField > divU
Definition: alphaSuSp.H:9
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:188
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
virtual PtrList< surfaceScalarField > DByAfs(const PtrList< volScalarField > &rAUs, const PtrList< surfaceScalarField > &rAUfs) const =0
Return the phase diffusivity.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const label nAlphaCorr(alphaControls.lookup< label >("nAlphaCorr"))
MULES: Multidimensional universal limiter for explicit solution.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
surfaceScalarField phi_
Total volumetric flux.
Definition: phaseSystem.H:169
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
virtual bool implicitPhasePressure() const
Returns true if the phase pressure is treated implicitly.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
Calculate the matrix for implicit and explicit sources.