reactingMultiphaseEulerFoam: Optimize the handling of optional forces

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C

@@ -511,17 +511,12 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
 
 
 template<class BasePhaseSystem>
-Foam::autoPtr<Foam::PtrList<Foam::volVectorField> >
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
+Foam::volVectorField& Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::setF
+(
+    PtrList<volVectorField>& Fs, const label phasei
+) const
 {
-    autoPtr<PtrList<volVectorField> > tFs
-    (
-        new PtrList<volVectorField>(this->phases().size())
-    );
-
-    PtrList<volVectorField>& Fs = tFs();
-
-    forAll(Fs, phasei)
+    if (!Fs.set(phasei))
     {
         Fs.set
         (
@@ -543,6 +538,20 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
         );
     }
 
+    return Fs[phasei];
+}
+
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::PtrList<Foam::volVectorField> >
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
+{
+    autoPtr<PtrList<volVectorField> > tFs
+    (
+        new PtrList<volVectorField>(this->phases().size())
+    );
+    PtrList<volVectorField>& Fs = tFs();
+
     // Add the lift force
     forAllConstIter
     (
@@ -555,8 +564,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
 
         const phasePair& pair(this->phasePairs_[liftModelIter.key()]);
 
-        Fs[pair.phase1().index()] += F;
-        Fs[pair.phase2().index()] -= F;
+        setF(Fs, pair.phase1().index()) += F;
+        setF(Fs, pair.phase2().index()) -= F;
     }
 
     // Add the wall lubrication force
@@ -572,8 +581,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
         const phasePair&
             pair(this->phasePairs_[wallLubricationModelIter.key()]);
 
-        Fs[pair.phase1().index()] += F;
-        Fs[pair.phase2().index()] -= F;
+        setF(Fs, pair.phase1().index()) += F;
+        setF(Fs, pair.phase2().index()) -= F;
     }
 
     return tFs;
@@ -581,20 +590,13 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
 
 
 template<class BasePhaseSystem>
-Foam::autoPtr<Foam::PtrList<Foam::surfaceScalarField> >
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
+Foam::surfaceScalarField&
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::setPhiD
 (
-    const PtrList<volScalarField>& rAUs
+    PtrList<surfaceScalarField>& phiDs, const label phasei
 ) const
 {
-    autoPtr<PtrList<surfaceScalarField> > tphiDs
-    (
-        new PtrList<surfaceScalarField>(this->phases().size())
-    );
-
-    PtrList<surfaceScalarField>& phiDs = tphiDs();
-
-    forAll(phiDs, phasei)
+    if (!phiDs.set(phasei))
     {
         phiDs.set
         (
@@ -603,7 +605,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
             (
                 IOobject
                 (
-                    liftModel::typeName + ":F",
+                    turbulentDispersionModel::typeName + ":phiD",
                     this->mesh_.time().timeName(),
                     this->mesh_,
                     IOobject::NO_READ,
@@ -621,6 +623,23 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
         );
     }
 
+    return phiDs[phasei];
+}
+
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::PtrList<Foam::surfaceScalarField> >
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
+(
+    const PtrList<volScalarField>& rAUs
+) const
+{
+    autoPtr<PtrList<surfaceScalarField> > tphiDs
+    (
+        new PtrList<surfaceScalarField>(this->phases().size())
+    );
+    PtrList<surfaceScalarField>& phiDs = tphiDs();
+
     // Add the turbulent dispersion force
     forAllConstIter
     (
@@ -638,9 +657,9 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
             fvc::snGrad(pair.phase1())*this->mesh_.magSf()
         );
 
-        phiDs[pair.phase1().index()] +=
+        setPhiD(phiDs, pair.phase1().index()) +=
             fvc::interpolate(rAUs[pair.phase1().index()]*D)*snGradAlpha1;
-        phiDs[pair.phase2().index()] -=
+        setPhiD(phiDs, pair.phase2().index()) -=
             fvc::interpolate(rAUs[pair.phase2().index()]*D)*snGradAlpha1;
     }
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H

@@ -133,6 +133,20 @@ private:
             //- Turbulent dispersion models
             turbulentDispersionModelTable turbulentDispersionModels_;
 
+            //- Construct element phasei of Fs if not set and return
+            //  Used by Fs()
+            volVectorField& setF
+            (
+                PtrList<volVectorField>& Fs, const label phasei
+            ) const;
+
+            //- Construct element phasei of phiDs if not set and return
+            //  Used by phiDs()
+            surfaceScalarField& setPhiD
+            (
+                PtrList<surfaceScalarField>& phiDs, const label phasei
+            ) const;
+
 
 public:
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H

@@ -38,22 +38,51 @@ forAll(phases, phasei)
 PtrList<surfaceScalarField> phiFs(phases.size());
 {
     autoPtr<PtrList<volVectorField> > Fs = fluid.Fs();
-    autoPtr<PtrList<surfaceScalarField> > phiDs = fluid.phiDs(rAUs);
 
     forAll(phases, phasei)
     {
         phaseModel& phase = phases[phasei];
 
-        phiFs.set
-        (
-            phasei,
-            new surfaceScalarField
+        if (Fs().set(phasei))
+        {
+            phiFs.set
             (
-                IOobject::groupName("phiF", phase.name()),
-                (fvc::interpolate(rAUs[phasei]*Fs()[phasei]) & mesh.Sf())
-              + phiDs()[phasei]
-            )
-        );
+                phasei,
+                new surfaceScalarField
+                (
+                    IOobject::groupName("phiF", phase.name()),
+                    (fvc::interpolate(rAUs[phasei]*Fs()[phasei]) & mesh.Sf())
+                )
+            );
+        }
+    }
+}
+{
+    autoPtr<PtrList<surfaceScalarField> > phiDs = fluid.phiDs(rAUs);
+
+    forAll(phases, phasei)
+    {
+        phaseModel& phase = phases[phasei];
+
+        if (phiDs().set(phasei))
+        {
+            if (phiFs.set(phasei))
+            {
+                phiFs[phasei] += phiDs()[phasei];
+            }
+            else
+            {
+                phiFs.set
+                (
+                    phasei,
+                    new surfaceScalarField
+                    (
+                        IOobject::groupName("phiF", phase.name()),
+                        phiDs()[phasei]
+                    )
+                );
+            }
+        }
     }
 }
 
@@ -101,12 +130,12 @@ while (pimple.correct())
         ghf*fvc::snGrad(rho)*mesh.magSf()
     );
 
-    PtrList<surfaceScalarField> phigs(phases.size());
+    PtrList<surfaceScalarField> phigFs(phases.size());
     forAll(phases, phasei)
     {
         phaseModel& phase = phases[phasei];
 
-        phigs.set
+        phigFs.set
         (
             phasei,
             (
@@ -118,6 +147,11 @@ while (pimple.correct())
                 )
             ).ptr()
         );
+
+        if (phiFs.set(phasei))
+        {
+            phigFs[phasei] += phiFs[phasei];
+        }
     }
 
     PtrList<surfaceScalarField> phiHbyAs(phases.size());
@@ -177,8 +211,7 @@ while (pimple.correct())
                    MRF.absolute(phase.phi().oldTime())
                  - (fvc::interpolate(phase.U().oldTime()) & mesh.Sf())
                 )/runTime.deltaT()
-              - phiFs[phasei]
-              - phigs[phasei]
+              - phigFs[phasei]
             )
         );
 
@@ -389,8 +422,7 @@ while (pimple.correct())
                   + fvc::reconstruct
                     (
                         alpharAUfs[phasei]*mSfGradp
-                      - phiFs[phasei]
-                      - phigs[phasei]
+                      - phigFs[phasei]
                     );
                 phase.U().correctBoundaryConditions();
                 fvOptions.correct(phase.U());

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties

@@ -151,16 +151,13 @@ heatTransfer
 );
 
 lift
-(
-);
+();
 
 wallLubrication
-(
-);
+();
 
 turbulentDispersion
-(
-);
+();
 
 // Minimum allowable pressure
 pMin            10000;