multiphaseEuler: New implicit drag algorithm replaces partial-elimination corrector

The momentum equation central coefficient and drag matrix is formulated,
inverted and used to eliminate the drag terms from each of the phase momentum
equations which are combined for formulate a drag-implicit pressure equation.
This eliminates the lagged drag terms from the previous formulation which
significantly improves convergence for small particle and Euler-VoF high-drag
cases.

It would also be possible to refactor the virtual-mass terms and include the
central coefficients of the phase acceleration terms in the drag matrix before
inversion to further improve the implicitness of the phase momentum-pressure
coupling for bubbly flows.  This work is pending funding.

applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C

@@ -45,6 +45,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
 {
     volScalarField& p(p_);
 
+    const phaseSystem::phaseModelPartialList& movingPhases
+    (
+        fluid.movingPhases()
+    );
+
     // Face volume fractions
     PtrList<surfaceScalarField> alphafs(phases.size());
     forAll(phases, phasei)
@@ -57,21 +62,25 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
     }
 
     // Diagonal coefficients
-    rAUs.clear();
+    rAs.clear();
-    rAUs.setSize(phases.size());
+    if (fluid.implicitPhasePressure())
-
-    PtrList<surfaceScalarField> rAUfs(phases.size());
-
     {
-        PtrList<volScalarField> Kds(fluid.Kds());
+        rAs.setSize(phases.size());
+    }
 
-        forAll(fluid.movingPhases(), movingPhasei)
+    PtrList<PtrList<volScalarField>> invADs;
+    PtrList<PtrList<surfaceScalarField>> invADfs;
     {
-            const phaseModel& phase = fluid.movingPhases()[movingPhasei];
+        PtrList<volScalarField> As(movingPhases.size());
+        PtrList<surfaceScalarField> Afs(movingPhases.size());
+        forAll(movingPhases, movingPhasei)
+        {
+            const phaseModel& phase = movingPhases[movingPhasei];
             const volScalarField& alpha = phase;
 
-            volScalarField AU
+            As.set
             (
+                movingPhasei,
                 UEqns[phase.index()].A()
               + byDt
                 (
@@ -80,80 +89,64 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                 )
             );
 
-            if (Kds.set(phase.index()))
+            Afs.set
-            {
+            (
-                AU += Kds[phase.index()];
+                movingPhasei,
-            }
+                fvc::interpolate(As[movingPhasei])
+            );
 
-            rAUs.set
+            if (fluid.implicitPhasePressure())
+            {
+                rAs.set
                 (
                     phase.index(),
                     new volScalarField
                     (
-                    IOobject::groupName("rAU", phase.name()),
+                        IOobject::groupName("rA", phase.name()),
-                    1/AU
+                        1/As[movingPhasei]
                     )
                 );
-
-            rAUfs.set(phase.index(), 1/fvc::interpolate(AU));
             }
-        fluid.fillFields("rAU", dimTime/dimDensity, rAUs);
-        fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs);
         }
 
-    // Phase diagonal coefficients
+        fluid.invADs(As, invADs, invADfs);
-    PtrList<surfaceScalarField> alpharAUfs(phases.size());
-    forAll(phases, phasei)
-    {
-        const phaseModel& phase = phases[phasei];
-        const volScalarField& alpha = phase;
-
-        alpharAUfs.set
-        (
-            phasei,
-            (
-                fvc::interpolate(max(alpha, phase.residualAlpha()))
-               *rAUfs[phasei]
-            ).ptr()
-        );
     }
 
-    // Explicit force fluxes
-    PtrList<surfaceScalarField> Fs(fluid.Fs());
-
-    // Mass transfer rates
-    PtrList<volScalarField> dmdts(fluid.dmdts());
-    PtrList<volScalarField> d2mdtdps(fluid.d2mdtdps());
-
-    // --- Pressure corrector loop
-    while (pimple.correct())
-    {
     volScalarField rho("rho", fluid.rho());
 
-        // Correct p_rgh for consistency with p and the updated densities
+    // Explicit force fluxes
-        p_rgh = p - rho*buoyancy.gh;
+    PtrList<surfaceScalarField> alphaByADfs;
-
+    PtrList<surfaceScalarField> FgByADfs;
-        // Correct fixed-flux BCs to be consistent with the velocity BCs
-        fluid_.correctBoundaryFlux();
-
-        // Combined buoyancy and force fluxes
-        PtrList<surfaceScalarField> phigFs(phases.size());
     {
+        PtrList<surfaceScalarField> Ffs(fluid.Fs());
+
         const surfaceScalarField ghSnGradRho
         (
             "ghSnGradRho",
             buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf()
         );
 
-            forAll(phases, phasei)
+        PtrList<surfaceScalarField> lalphafs(movingPhases.size());
-            {
+        PtrList<surfaceScalarField> Fgfs(movingPhases.size());
-                const phaseModel& phase = phases[phasei];
 
-                phigFs.set
+        forAll(movingPhases, movingPhasei)
+        {
+            const phaseModel& phase = movingPhases[movingPhasei];
+            const volScalarField& alpha = phase;
+
+            lalphafs.set
             (
-                    phasei,
+                movingPhasei,
+                fvc::interpolate(max(alpha, phase.residualAlpha()))
+            );
+
+            Fgfs.set
             (
+                movingPhasei,
                 (
+                    Ffs[phase.index()]
+                  + lalphafs[movingPhasei]
+                   *(
                        ghSnGradRho
                      - (fvc::interpolate(phase.rho() - rho))
                       *(buoyancy.g & mesh.Sf())
@@ -162,71 +155,123 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                 ).ptr()
             );
         }
+
+        alphaByADfs = invADfs & lalphafs;
+        FgByADfs = invADfs & Fgfs;
     }
 
-        // Predicted velocities and fluxes for each phase
-        PtrList<volVectorField> HbyAs(phases.size());
-        PtrList<surfaceScalarField> phiHbyAs(phases.size());
-        {
-            // Correction force fluxes
-            PtrList<surfaceScalarField> ddtCorrs(fluid.ddtCorrs());
 
-            forAll(fluid.movingPhases(), movingPhasei)
+    // Mass transfer rates
+    PtrList<volScalarField> dmdts(fluid.dmdts());
+    PtrList<volScalarField> d2mdtdps(fluid.d2mdtdps());
+
+    // --- Optional momentum predictor
+    if (predictMomentum)
     {
-                const phaseModel& phase = fluid.movingPhases()[movingPhasei];
+        InfoInFunction << endl;
+        PtrList<volVectorField> HbyADs;
+        {
+            PtrList<volVectorField> Hs(movingPhases.size());
+
+            forAll(movingPhases, movingPhasei)
+            {
+                const phaseModel& phase = movingPhases[movingPhasei];
                 const volScalarField& alpha = phase;
-                const label phasei = phase.index();
 
-                const volVectorField H
+                Hs.set
                 (
-                    constrainH
+                    movingPhasei,
-                    (
+                    UEqns[phase.index()].H()
-                        UEqns[phasei].H()
                   + byDt
                     (
                         max(phase.residualAlpha() - alpha, scalar(0))
                        *phase.rho()
                     )
-                       *phase.U()().oldTime(),
+                   *phase.U()().oldTime()
-                        rAUs[phasei],
-                        phase.U(),
-                        p_rgh
-                    )
-                );
-
-                HbyAs.set(phasei, rAUs[phasei]*H);
-
-                phiHbyAs.set
-                (
-                    phasei,
-                    new surfaceScalarField
-                    (
-                        IOobject::groupName("phiHbyA", phase.name()),
-                        rAUfs[phasei]*(fvc::flux(H) + ddtCorrs[phasei])
-                      - alpharAUfs[phasei]*phigFs[phasei]
-                      - rAUfs[phasei]*Fs[phasei]
-                    )
                 );
             }
-        }
-        fluid.fillFields("HbyA", dimVelocity, HbyAs);
-        fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs);
 
-        // Add explicit drag forces and fluxes
+            HbyADs = invADs & Hs;
-        PtrList<volVectorField> KdUs(fluid.KdUs());
-        PtrList<surfaceScalarField> KdPhis(fluid.KdPhis());
-
-        forAll(phases, phasei)
-        {
-            if (KdUs.set(phasei))
-            {
-                HbyAs[phasei] -= rAUs[phasei]*KdUs[phasei];
         }
 
-            if (KdPhis.set(phasei))
+        forAll(movingPhases, movingPhasei)
         {
-                phiHbyAs[phasei] -= rAUfs[phasei]*KdPhis[phasei];
+            const phaseModel& phase = movingPhases[movingPhasei];
+            constrainHbyA(HbyADs[movingPhasei], phase.U(), p_rgh);
         }
+
+        const surfaceScalarField mSfGradp(-mesh.magSf()*fvc::snGrad(p_rgh));
+
+        forAll(movingPhases, movingPhasei)
+        {
+            phaseModel& phase = fluid_.movingPhases()[movingPhasei];
+
+            phase.URef() =
+                HbyADs[movingPhasei]
+              + fvc::reconstruct
+                (
+                    alphaByADfs[movingPhasei]*mSfGradp
+                  - FgByADfs[movingPhasei]
+                );
+
+            phase.URef().correctBoundaryConditions();
+            fvConstraints().constrain(phase.URef());
+        }
+    }
+
+
+    // --- Pressure corrector loop
+    while (pimple.correct())
+    {
+        // Correct fixed-flux BCs to be consistent with the velocity BCs
+        fluid_.correctBoundaryFlux();
+
+        PtrList<volVectorField> HbyADs;
+        PtrList<surfaceScalarField> phiHbyADs;
+        {
+            // Predicted velocities and fluxes for each phase
+            PtrList<volVectorField> Hs(movingPhases.size());
+            PtrList<surfaceScalarField> phiHs(movingPhases.size());
+
+            // Correction force fluxes
+            PtrList<surfaceScalarField> ddtCorrs(fluid.ddtCorrs());
+
+            forAll(movingPhases, movingPhasei)
+            {
+                const phaseModel& phase = movingPhases[movingPhasei];
+                const volScalarField& alpha = phase;
+
+                Hs.set
+                (
+                    movingPhasei,
+                    UEqns[phase.index()].H()
+                  + byDt
+                    (
+                        max(phase.residualAlpha() - alpha, scalar(0))
+                       *phase.rho()
+                    )
+                   *phase.U()().oldTime()
+                );
+
+                phiHs.set
+                (
+                    movingPhasei,
+                    fvc::flux(Hs[movingPhasei]) + ddtCorrs[phase.index()]
+                );
+            }
+
+            HbyADs = invADs & Hs;
+            phiHbyADs = invADfs & phiHs;
+        }
+
+        forAll(movingPhases, movingPhasei)
+        {
+            const phaseModel& phase = movingPhases[movingPhasei];
+
+            constrainHbyA(HbyADs[movingPhasei], phase.U(), p_rgh);
+            constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh);
+
+            phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei];
         }
 
         // Total predicted flux
@@ -236,28 +281,28 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
             (
                 "phiHbyA",
                 runTime.name(),
-                mesh,
+                mesh
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
             ),
             mesh,
             dimensionedScalar(dimFlux, 0)
         );
 
-        forAll(phases, phasei)
+        forAll(movingPhases, movingPhasei)
         {
-            phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
+            const phaseModel& phase = movingPhases[movingPhasei];
+
+            phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei];
         }
 
         MRF.makeRelative(phiHbyA);
-        fvc::makeRelative(phiHbyA, fluid.movingPhases()[0].U());
+        fvc::makeRelative(phiHbyA, movingPhases[0].U());
 
-        // Construct pressure "diffusivity"
+        // Pressure "diffusivity"
-        surfaceScalarField rAUf
+        surfaceScalarField rAf
         (
             IOobject
             (
-                "rAUf",
+                "rAf",
                 runTime.name(),
                 mesh
             ),
@@ -265,9 +310,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
             dimensionedScalar(dimTime/dimDensity, 0)
         );
 
-        forAll(phases, phasei)
+        forAll(movingPhases, movingPhasei)
         {
-            rAUf += alphafs[phasei]*alpharAUfs[phasei];
+            const phaseModel& phase = movingPhases[movingPhasei];
+
+            rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei];
         }
 
         // Update the fixedFluxPressure BCs to ensure flux consistency
@@ -279,11 +326,12 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
             );
             phib = 0;
 
-            forAll(phases, phasei)
+            forAll(movingPhases, movingPhasei)
             {
-                const phaseModel& phase = phases[phasei];
+                phaseModel& phase = fluid_.movingPhases()[movingPhasei];
+
                 phib +=
-                    alphafs[phasei].boundaryField()
+                    alphafs[phase.index()].boundaryField()
                    *phase.phi()().boundaryField();
             }
 
@@ -292,7 +340,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                 p_rgh.boundaryFieldRef(),
                 (
                     phiHbyA.boundaryField() - phib
-                )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+                )/(mesh.magSf().boundaryField()*rAf.boundaryField())
             );
         }
 
@@ -309,7 +357,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
             fvScalarMatrix pEqnIncomp
             (
                 fvc::div(phiHbyA)
-              - fvm::laplacian(rAUf, p_rgh)
+              - fvm::laplacian(rAf, p_rgh)
             );
 
             // Solve
@@ -340,86 +388,76 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
             {
                 phi_ = phiHbyA + pEqnIncomp.flux();
 
-                surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
+                surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf);
 
-                forAll(fluid.movingPhases(), movingPhasei)
+                forAll(movingPhases, movingPhasei)
                 {
                     phaseModel& phase = fluid_.movingPhases()[movingPhasei];
 
                     phase.phiRef() =
-                        phiHbyAs[phase.index()]
+                        phiHbyADs[movingPhasei]
-                      + alpharAUfs[phase.index()]*mSfGradp;
+                      + alphaByADfs[movingPhasei]*mSfGradp;
 
                     // Set the phase dilatation rate
                     phase.divU(-pEqnComps[phase.index()] & p_rgh);
+
+                    MRF.makeRelative(phase.phiRef());
+                    fvc::makeRelative(phase.phiRef(), phase.U());
                 }
 
                 // Optionally relax pressure for velocity correction
                 p_rgh.relax();
 
-                mSfGradp = pEqnIncomp.flux()/rAUf;
+                mSfGradp = pEqnIncomp.flux()/rAf;
 
-                if (!dragCorrection)
+                if (dragCorrection)
                 {
-                    forAll(fluid.movingPhases(), movingPhasei)
+                    PtrList<volVectorField> dragCorrs(movingPhases.size());
-                    {
+                    PtrList<surfaceScalarField> dragCorrfs(movingPhases.size());
-                        phaseModel& phase = fluid_.movingPhases()[movingPhasei];
-                        const label phasei = phase.index();
-
-                        phase.URef() =
-                            HbyAs[phasei]
-                          + fvc::reconstruct
-                            (
-                                alpharAUfs[phasei]*(mSfGradp - phigFs[phasei])
-                              - rAUfs[phasei]*Fs[phasei]
-                            );
-                    }
-                }
-                else
-                {
-                    PtrList<volVectorField> dragCorrs(phases.size());
-                    PtrList<surfaceScalarField> dragCorrfs(phases.size());
                     fluid.dragCorrs(dragCorrs, dragCorrfs);
 
-                    forAll(fluid.movingPhases(), movingPhasei)
+                    PtrList<volVectorField> dragCorrByADs
+                    (
+                        invADs & dragCorrs
+                    );
+
+                    PtrList<surfaceScalarField> dragCorrByADfs
+                    (
+                        invADfs & dragCorrfs
+                    );
+
+                    forAll(movingPhases, movingPhasei)
                     {
                         phaseModel& phase = fluid_.movingPhases()[movingPhasei];
-                        const label phasei = phase.index();
 
                         phase.URef() =
-                            HbyAs[phasei]
+                            HbyADs[movingPhasei]
                           + fvc::reconstruct
                             (
-                                alpharAUfs[phasei]*(mSfGradp - phigFs[phasei])
+                                alphaByADfs[movingPhasei]*mSfGradp
-                              + rAUfs[phasei]*(dragCorrfs[phasei] - Fs[phasei])
+                              - FgByADfs[movingPhasei]
+                              + dragCorrByADfs[movingPhasei]
                             )
-                        - rAUs[phasei]*dragCorrs[phasei];
+                          - dragCorrByADs[movingPhasei];
                     }
                 }
-
-                if (partialElimination)
-                {
-                    fluid_.partialElimination
-                    (
-                        rAUs,
-                        KdUs,
-                        alphafs,
-                        rAUfs,
-                        KdPhis
-                    );
-                }
                 else
                 {
-                    forAll(fluid.movingPhases(), movingPhasei)
+                    forAll(movingPhases, movingPhasei)
                     {
                         phaseModel& phase = fluid_.movingPhases()[movingPhasei];
 
-                        MRF.makeRelative(phase.phiRef());
+                        phase.URef() =
-                        fvc::makeRelative(phase.phiRef(), phase.U());
+                            HbyADs[movingPhasei]
+                          + fvc::reconstruct
+                            (
+                                alphaByADfs[movingPhasei]*mSfGradp
+                              - FgByADfs[movingPhasei]
+                            );
                     }
                 }
 
-                forAll(fluid.movingPhases(), movingPhasei)
+                forAll(movingPhases, movingPhasei)
                 {
                     phaseModel& phase = fluid_.movingPhases()[movingPhasei];
 
@@ -472,11 +510,6 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
     }
 
     UEqns.clear();
-
-    if (!fluid.implicitPhasePressure())
-    {
-        rAUs.clear();
-    }
 }
 
 

applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C

@@ -45,9 +45,13 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
 {
     volScalarField& p(p_);
 
+    const phaseSystem::phaseModelPartialList& movingPhases
+    (
+        fluid.movingPhases()
+    );
+
     // Face volume fractions
     PtrList<surfaceScalarField> alphafs(phases.size());
-    PtrList<surfaceScalarField> alphaRho0fs(phases.size());
     forAll(phases, phasei)
     {
         const phaseModel& phase = phases[phasei];
@@ -55,65 +59,113 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
 
         alphafs.set(phasei, fvc::interpolate(alpha).ptr());
         alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
-
-        alphaRho0fs.set
-        (
-            phasei,
-            (
-                fvc::interpolate
-                (
-                    max(alpha.oldTime(), phase.residualAlpha())
-                   *phase.rho().oldTime()
-                )
-            ).ptr()
-        );
     }
 
     // Diagonal coefficients
-    rAUfs.clear();
+    rAs.clear();
-    rAUfs.setSize(phases.size());
+    if (fluid.implicitPhasePressure())
     {
-        PtrList<surfaceScalarField> KdVmfs(fluid.KdVmfs());
+        rAs.setSize(phases.size());
+    }
+
+    PtrList<PtrList<surfaceScalarField>> invADfs;
+    {
+        PtrList<surfaceScalarField> Vmfs(fluid.Vmfs());
+        PtrList<surfaceScalarField> Afs(movingPhases.size());
 
         forAll(fluid.movingPhases(), movingPhasei)
         {
             const phaseModel& phase = fluid.movingPhases()[movingPhasei];
+            const volScalarField& alpha = phase;
 
-            rAUfs.set
+            const volScalarField A
+            (
+                byDt
+                (
+                    max(alpha.oldTime(), phase.residualAlpha())
+                   *phase.rho().oldTime()
+                )
+              + UEqns[phase.index()].A()
+            );
+
+            if (fluid.implicitPhasePressure())
+            {
+                rAs.set
                 (
                     phase.index(),
+                    new volScalarField
+                    (
+                        IOobject::groupName("rA", phase.name()),
+                        1/A
+                    )
+                );
+            }
+
+            Afs.set
+            (
+                movingPhasei,
                 new surfaceScalarField
                 (
-                    IOobject::groupName("rAUf", phase.name()),
+                    IOobject::groupName("rAf", phase.name()),
-                    1.0
+                    fvc::interpolate(A)
-                   /(
-                        byDt(alphaRho0fs[phase.index()])
-                      + fvc::interpolate(UEqns[phase.index()].A())
-                      + KdVmfs[phase.index()]
-                    )
                 )
             );
+
+            if (Vmfs.set(phase.index()))
+            {
+                Afs[movingPhasei] += Vmfs[phase.index()];
             }
         }
-    fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs);
+
+        invADfs = fluid.invADfs(Afs);
+    }
+
+    volScalarField rho("rho", fluid.rho());
 
     // Phase diagonal coefficients
-    PtrList<surfaceScalarField> alpharAUfs(phases.size());
+    PtrList<surfaceScalarField> alphaByADfs;
-    forAll(phases, phasei)
+    PtrList<surfaceScalarField> FgByADfs;
     {
-        const phaseModel& phase = phases[phasei];
+        // Explicit force fluxes
-        alpharAUfs.set
+        PtrList<surfaceScalarField> Ffs(fluid.Ffs());
+
+        const surfaceScalarField ghSnGradRho
         (
-            phase.index(),
+            "ghSnGradRho",
+            buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf()
+        );
+
+        PtrList<surfaceScalarField> lalphafs(movingPhases.size());
+        PtrList<surfaceScalarField> Fgfs(movingPhases.size());
+
+        forAll(movingPhases, movingPhasei)
+        {
+            const phaseModel& phase = movingPhases[movingPhasei];
+
+            lalphafs.set
             (
+                movingPhasei,
                 max(alphafs[phase.index()], phase.residualAlpha())
-               *rAUfs[phase.index()]
+            );
-            ).ptr()
+
+            Fgfs.set
+            (
+                movingPhasei,
+                Ffs[phase.index()]
+              + lalphafs[movingPhasei]
+               *(
+                    ghSnGradRho
+                  - (fvc::interpolate(phase.rho() - rho))
+                   *(buoyancy.g & mesh.Sf())
+                  - fluid.surfaceTension(phase)*mesh.magSf()
+                )
             );
         }
 
-    // Explicit force fluxes
+        alphaByADfs = invADfs & lalphafs;
-    PtrList<surfaceScalarField> Ffs(fluid.Ffs());
+        FgByADfs = invADfs & Fgfs;
+    }
+
 
     // Mass transfer rates
     PtrList<volScalarField> dmdts(fluid.dmdts());
@@ -122,87 +174,49 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
     // --- Pressure corrector loop
     while (pimple.correct())
     {
-        volScalarField rho("rho", fluid.rho());
-
-        // Correct p_rgh for consistency with p and the updated densities
-        p_rgh = p - rho*buoyancy.gh;
-
         // Correct fixed-flux BCs to be consistent with the velocity BCs
         fluid_.correctBoundaryFlux();
 
-        // Combined buoyancy and force fluxes
-        PtrList<surfaceScalarField> phigFs(phases.size());
-        {
-            const surfaceScalarField ghSnGradRho
-            (
-                "ghSnGradRho",
-                buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf()
-            );
-
-            forAll(phases, phasei)
-            {
-                const phaseModel& phase = phases[phasei];
-
-                phigFs.set
-                (
-                    phasei,
-                    (
-                        alpharAUfs[phasei]
-                       *(
-                           ghSnGradRho
-                         - (fvc::interpolate(phase.rho() - rho))
-                          *(buoyancy.g & mesh.Sf())
-                         - fluid.surfaceTension(phase)*mesh.magSf()
-                        )
-                    ).ptr()
-                );
-
-                if (Ffs.set(phasei))
-                {
-                    phigFs[phasei] += rAUfs[phasei]*Ffs[phasei];
-                }
-            }
-        }
-
         // Predicted fluxes for each phase
-        PtrList<surfaceScalarField> phiHbyAs(phases.size());
+        PtrList<surfaceScalarField> phiHbyADs;
+        {
+            PtrList<surfaceScalarField> phiHs(movingPhases.size());
+
             forAll(fluid.movingPhases(), movingPhasei)
             {
                 const phaseModel& phase = fluid.movingPhases()[movingPhasei];
+                const volScalarField& alpha = phase;
 
-            phiHbyAs.set
+                phiHs.set
                 (
-                phase.index(),
+                    movingPhasei,
-                constrainPhiHbyA
                     (
-                    rAUfs[phase.index()]
+                       fvc::interpolate
-                   *(
+                       (
-                       fvc::flux(UEqns[phase.index()].H())
+                           max(alpha.oldTime(), phase.residualAlpha())
-                     + alphaRho0fs[phase.index()]
+                          *phase.rho().oldTime()
+                       )
                       *byDt
                        (
                            phase.Uf().valid()
                          ? (mesh.Sf() & phase.Uf()().oldTime())
                          : MRF.absolute(phase.phi()().oldTime())
                        )
-                    )
+                     + fvc::flux(UEqns[phase.index()].H())
-                  - phigFs[phase.index()],
-                    phase.U(),
-                    p_rgh
                     )
                 );
             }
-        fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs);
 
-        // Add explicit drag forces and fluxes
+            phiHbyADs = invADfs & phiHs;
-        PtrList<surfaceScalarField> KdPhifs(fluid.KdPhifs());
-
-        forAll(phases, phasei)
-        {
-            if (KdPhifs.set(phasei))
-            {
-                phiHbyAs[phasei] -= rAUfs[phasei]*KdPhifs[phasei];
         }
+
+        forAll(movingPhases, movingPhasei)
+        {
+            const phaseModel& phase = movingPhases[movingPhasei];
+
+            constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh);
+
+            phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei];
         }
 
         // Total predicted flux
@@ -220,20 +234,22 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
             dimensionedScalar(dimFlux, 0)
         );
 
-        forAll(phases, phasei)
+        forAll(movingPhases, movingPhasei)
         {
-            phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
+            const phaseModel& phase = movingPhases[movingPhasei];
+
+            phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei];
         }
 
         MRF.makeRelative(phiHbyA);
         fvc::makeRelative(phiHbyA, phases[0].U());
 
         // Construct pressure "diffusivity"
-        surfaceScalarField rAUf
+        surfaceScalarField rAf
         (
             IOobject
             (
-                "rAUf",
+                "rAf",
                 runTime.name(),
                 mesh
             ),
@@ -241,12 +257,12 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
             dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), 0)
         );
 
-        forAll(phases, phasei)
+        forAll(movingPhases, movingPhasei)
         {
-            rAUf += alphafs[phasei]*alpharAUfs[phasei];
+            const phaseModel& phase = movingPhases[movingPhasei];
-        }
 
-        rAUf = mag(rAUf);
+            rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei];
+        }
 
         // Update the fixedFluxPressure BCs to ensure flux consistency
         {
@@ -257,11 +273,12 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
             );
             phib = 0;
 
-            forAll(phases, phasei)
+            forAll(movingPhases, movingPhasei)
             {
-                const phaseModel& phase = phases[phasei];
+                phaseModel& phase = fluid_.movingPhases()[movingPhasei];
+
                 phib +=
-                    alphafs[phasei].boundaryField()
+                    alphafs[phase.index()].boundaryField()
                    *phase.phi()().boundaryField();
             }
 
@@ -270,7 +287,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
                 p_rgh.boundaryFieldRef(),
                 (
                     phiHbyA.boundaryField() - phib
-                )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+                )/(mesh.magSf().boundaryField()*rAf.boundaryField())
             );
         }
 
@@ -287,7 +304,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
             fvScalarMatrix pEqnIncomp
             (
                 fvc::div(phiHbyA)
-              - fvm::laplacian(rAUf, p_rgh)
+              - fvm::laplacian(rAf, p_rgh)
             );
 
             // Solve
@@ -318,38 +335,27 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
             {
                 phi_ = phiHbyA + pEqnIncomp.flux();
 
-                surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
+                surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf);
 
                 forAll(fluid.movingPhases(), movingPhasei)
                 {
                     phaseModel& phase = fluid_.movingPhases()[movingPhasei];
+                    const label phasei = phase.index();
 
                     phase.phiRef() =
-                        phiHbyAs[phase.index()]
+                        phiHbyADs[movingPhasei]
-                      + alpharAUfs[phase.index()]*mSfGradp;
+                      + alphaByADfs[movingPhasei]*mSfGradp;
 
                     // Set the phase dilatation rate
-                    phase.divU(-pEqnComps[phase.index()] & p_rgh);
+                    phase.divU(-pEqnComps[phasei] & p_rgh);
                 }
 
-                if (partialElimination)
-                {
-                    fluid_.partialEliminationf(rAUfs, alphafs, KdPhifs);
-                }
-                else
-                {
                 forAll(fluid.movingPhases(), movingPhasei)
                 {
                     phaseModel& phase = fluid_.movingPhases()[movingPhasei];
 
                     MRF.makeRelative(phase.phiRef());
                     fvc::makeRelative(phase.phiRef(), phase.U());
-                    }
-                }
-
-                forAll(fluid.movingPhases(), movingPhasei)
-                {
-                    phaseModel& phase = fluid_.movingPhases()[movingPhasei];
 
                     phase.URef() = fvc::reconstruct
                     (
@@ -405,11 +411,6 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
     }
 
     UEqns.clear();
-
-    if (!fluid.implicitPhasePressure())
-    {
-        rAUfs.clear();
-    }
 }
 
 

applications/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C

@@ -40,8 +40,6 @@ void Foam::solvers::multiphaseEuler::cellMomentumPredictor()
     phaseSystem::momentumTransferTable&
         momentumTransfer(momentumTransferPtr());
 
-    const PtrList<volScalarField> Kds(fluid.Kds());
-
     forAll(fluid.movingPhases(), movingPhasei)
     {
         phaseModel& phase = fluid.movingPhases()[movingPhasei];
@@ -59,7 +57,6 @@ void Foam::solvers::multiphaseEuler::cellMomentumPredictor()
              ==
                *momentumTransfer[phase.name()]
               + fvModels().source(alpha, rho, U)
-                // - fvm::Sp(Kds[phase.index()], U)
             )
         );
 

applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C

@@ -51,15 +51,15 @@ void Foam::solvers::multiphaseEuler::readControls(const bool construct)
 
     if (construct || mesh.solution().modified())
     {
+        predictMomentum =
+            pimple.dict().lookupOrDefault<bool>("momentumPredictor", false);
+
         faceMomentum =
             pimple.dict().lookupOrDefault<Switch>("faceMomentum", false);
 
         dragCorrection =
             pimple.dict().lookupOrDefault<Switch>("dragCorrection", false);
 
-        partialElimination =
-            pimple.dict().lookupOrDefault<Switch>("partialElimination", false);
-
         nEnergyCorrectors =
             pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1);
     }
@@ -98,6 +98,11 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh)
 :
     fluidSolver(mesh),
 
+    predictMomentum
+    (
+        pimple.dict().lookupOrDefault<Switch>("momentumPredictor", false)
+    ),
+
     faceMomentum
     (
         pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
@@ -108,11 +113,6 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh)
         pimple.dict().lookupOrDefault<Switch>("dragCorrection", false)
     ),
 
-    partialElimination
-    (
-        pimple.dict().lookupOrDefault<Switch>("partialElimination", false)
-    ),
-
     nEnergyCorrectors
     (
         pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1)
@@ -247,7 +247,7 @@ void Foam::solvers::multiphaseEuler::prePredictor()
 {
     if (pimple.thermophysics() || pimple.flow())
     {
-        fluid_.solve(rAUs, rAUfs);
+        fluid_.solve(rAs);
         fluid_.correct();
         fluid_.correctContinuityError();
     }

applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H

@@ -77,6 +77,10 @@ protected:
 
     // Controls
 
+        //- Momentum equation predictor switch
+        //  Defaults to false
+        Switch predictMomentum;
+
         //- Cell/face momentum equation switch
         //  Defaults to false, i.e. uses the cell momentum equation
         Switch faceMomentum;
@@ -85,10 +89,6 @@ protected:
         //  Defaults to false
         Switch dragCorrection;
 
-        //- Partial elimination drag contribution optimisation
-        //  Defaults to false
-        Switch partialElimination;
-
         //- Number of energy correctors
         //  Used to improve stability of phase-change sibulations
         //  Defaults to 1
@@ -144,11 +144,7 @@ protected:
 
         //- Temporary storage for the reciprocal momentum equation diagonal
         //  Used by the phase-fraction predictor and pressure corrector
-        PtrList<volScalarField> rAUs;
+        PtrList<volScalarField> rAs;
-
-        //- Temporary storage for the reciprocal momentum equation diagonal
-        //  Used by the phase-fraction predictor and face pressure corrector
-        PtrList<surfaceScalarField> rAUfs;
 
         //- Stored divU from the previous mesh so that it can be
         //  mapped and used in correctPhi to ensure the corrected phi

applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C

@@ -31,6 +31,8 @@ License
 #include "wallLubricationModel.H"
 #include "turbulentDispersionModel.H"
 
+#include "scalarMatrices.H"
+
 #include "fvmDdt.H"
 #include "fvmDiv.H"
 #include "fvmSup.H"
@@ -474,36 +476,9 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf()
 
 template<class BasePhaseSystem>
 Foam::PtrList<Foam::surfaceScalarField>
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdVmfs() const
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vmfs() const
 {
-    PtrList<surfaceScalarField> KdVmfs(this->phaseModels_.size());
+    PtrList<surfaceScalarField> Vmfs(this->phaseModels_.size());
-
-    // Add the implicit part of the drag force
-    forAllConstIter(KdfTable, Kdfs_, KdfIter)
-    {
-        const surfaceScalarField& Kf(*KdfIter());
-        const phaseInterface interface(*this, KdfIter.key());
-
-        forAllConstIter(phaseInterface, interface, iter)
-        {
-            const phaseModel& phase = iter();
-            const phaseModel& otherPhase = iter.otherPhase();
-
-            const surfaceScalarField alphaf
-            (
-                fvc::interpolate(otherPhase)
-            );
-
-            addField
-            (
-                phase,
-                "KdVmf",
-                (alphaf/max(alphaf, otherPhase.residualAlpha()))
-               *Kf,
-                KdVmfs
-            );
-        }
-    }
 
     // Add the implicit part of the virtual mass force
     forAllConstIter(VmTable, Vms_, VmIter)
@@ -522,13 +497,11 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdVmfs() const
                *Vm
             );
 
-            addField(phase, "KdVmf", byDt(fvc::interpolate(VmPhase)), KdVmfs);
+            addField(phase, "Vmf", byDt(fvc::interpolate(VmPhase)), Vmfs);
         }
     }
 
-    this->fillFields("KdVmf", dimDensity/dimTime, KdVmfs);
+    return Vmfs;
-
-    return KdVmfs;
 }
 
 
@@ -650,11 +623,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
         addField(interface.phase2(), "F", Df*snGradAlpha2By12, Fs);
     }
 
-    if (this->fillFields_)
-    {
-        this->fillFields("F", dimArea*dimDensity*dimAcceleration, Fs);
-    }
-
     return Fs;
 }
 
@@ -813,22 +781,121 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Ffs() const
         addField(interface.phase2(), "F", Df*snGradAlpha2By12, Ffs);
     }
 
-    if (this->fillFields_)
-    {
-        this->fillFields("Ff", dimArea*dimDensity*dimAcceleration, Ffs);
-    }
-
     return Ffs;
 }
 
 
 template<class BasePhaseSystem>
-Foam::PtrList<Foam::surfaceScalarField>
+void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::invADs
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdPhis() const
+(
+    List<UPtrList<scalarField>>& ADs
+) const
 {
-    PtrList<surfaceScalarField> KdPhis(this->phaseModels_.size());
+    const label n = ADs.size();
+
+    scalarSquareMatrix AD(n);
+    scalarField source(n);
+    labelList pivotIndices(n);
+
+    forAll(ADs[0][0], ci)
+    {
+        for (label i=0; i<n; i++)
+        {
+            for (label j=0; j<n; j++)
+            {
+                AD(i, j) = ADs[i][j][ci];
+            }
+        }
+
+        // Calculate the inverse of AD using LD decomposition
+        // and back-substitution
+        LUDecompose(AD, pivotIndices);
+
+        for (label j=0; j<n; j++)
+        {
+            source = Zero;
+            source[j] = 1;
+
+            LUBacksubstitute(AD, pivotIndices, source);
+
+            for (label i=0; i<n; i++)
+            {
+                ADs[i][j][ci] = source[i];
+            }
+        }
+    }
+}
+
+
+template<class BasePhaseSystem>
+template<template<class> class PatchField, class GeoMesh>
+void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::invADs
+(
+    PtrList<PtrList<GeometricField<scalar, PatchField, GeoMesh>>>& ADs
+) const
+{
+    const label n = ADs.size();
+
+    List<UPtrList<scalarField>> ADps(n);
+
+    for (label i=0; i<n; i++)
+    {
+        ADps[i].setSize(n);
+
+        for (label j=0; j<n; j++)
+        {
+            ADps[i].set(j, &ADs[i][j]);
+
+            ADs[i][j].dimensions().reset(dimless/ADs[i][j].dimensions());
+        }
+    }
+
+    invADs(ADps);
+
+    forAll(ADs[0][0].boundaryField(), patchi)
+    {
+        for (label i=0; i<n; i++)
+        {
+            for (label j=0; j<n; j++)
+            {
+                ADps[i].set(j, &ADs[i][j].boundaryFieldRef()[patchi]);
+            }
+        }
+
+        invADs(ADps);
+    }
+}
+
+
+template<class BasePhaseSystem>
+void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::invADs
+(
+    const PtrList<volScalarField>& As,
+    PtrList<PtrList<volScalarField>>& invADs,
+    PtrList<PtrList<surfaceScalarField>>& invADfs
+) const
+{
+    const label n = As.size();
+
+    invADs.setSize(n);
+    invADfs.setSize(n);
+
+    forAll(invADs, i)
+    {
+        invADs.set(i, new PtrList<volScalarField>(n));
+        invADs[i].set(i, As[i].clone());
+
+        invADfs.set(i, new PtrList<surfaceScalarField>(n));
+        invADfs[i].set(i, fvc::interpolate(As[i]));
+    }
+
+    labelList movingPhases(this->phases().size(), -1);
+
+    forAll(this->movingPhases(), movingPhasei)
+    {
+        movingPhases[this->movingPhases()[movingPhasei].index()] = movingPhasei;
+    }
 
-    // Add the explicit part of the drag force
     forAllConstIter(KdTable, Kds_, KdIter)
     {
         const volScalarField& K(*KdIter());
@@ -839,46 +906,91 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdPhis() const
             const phaseModel& phase = iter();
             const phaseModel& otherPhase = iter.otherPhase();
 
-            addField
+            const label i = movingPhases[phase.index()];
-            (
-                phase,
-                "KdPhi",
-                fvc::interpolate
-                (
-                  -(otherPhase/max(otherPhase, otherPhase.residualAlpha()))
-                  *K
-                )
-               *fvc::absolute
-                (
-                    this->MRF().absolute(otherPhase.phi()),
-                    otherPhase.U()
-                ),
-                KdPhis
-            );
-        }
-    }
 
-    if (this->fillFields_)
+            if (i != -1)
             {
-        this->fillFields
+                const volScalarField Kij
                 (
-            "KdPhi",
+                    (otherPhase/max(otherPhase, otherPhase.residualAlpha()))*K
-            dimArea*dimDensity*dimAcceleration,
-            KdPhis
                 );
+
+                const surfaceScalarField Kijf(fvc::interpolate(Kij));
+
+                invADs[i][i] += Kij;
+                invADfs[i][i] += Kijf;
+
+                const label j = movingPhases[otherPhase.index()];
+
+                if (j != -1)
+                {
+                    invADs[i].set(j, -Kij);
+                    invADfs[i].set(j, -Kijf);
+                }
+            }
+        }
     }
 
-    return KdPhis;
+    for (label i=0; i<n; i++)
+    {
+        for (label j=0; j<n; j++)
+        {
+            if (!invADs[i].set(j))
+            {
+                invADs[i].set
+                (
+                    j,
+                    volScalarField::New
+                    (
+                        "0",
+                        this->mesh(),
+                        dimensionedScalar(As[0].dimensions(), 0)
+                    )
+                );
+
+                invADfs[i].set
+                (
+                    j,
+                    surfaceScalarField::New
+                    (
+                        "0",
+                        this->mesh(),
+                        dimensionedScalar(As[0].dimensions(), 0)
+                    )
+                );
+            }
+        }
+    }
+
+    MomentumTransferPhaseSystem<BasePhaseSystem>::invADs(invADs);
+    MomentumTransferPhaseSystem<BasePhaseSystem>::invADs(invADfs);
 }
 
 
 template<class BasePhaseSystem>
-Foam::PtrList<Foam::surfaceScalarField>
+Foam::PtrList<Foam::PtrList<Foam::surfaceScalarField>>
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdPhifs() const
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::invADfs
+(
+    const PtrList<surfaceScalarField>& As
+) const
 {
-    PtrList<surfaceScalarField> KdPhifs(this->phaseModels_.size());
+    const label n = As.size();
+
+    PtrList<PtrList<surfaceScalarField>> invADfs(n);
+
+    forAll(invADfs, i)
+    {
+        invADfs.set(i, new PtrList<surfaceScalarField>(n));
+        invADfs[i].set(i, As[i].clone());
+    }
+
+    labelList movingPhases(this->phases().size(), -1);
+
+    forAll(this->movingPhases(), movingPhasei)
+    {
+        movingPhases[this->movingPhases()[movingPhasei].index()] = movingPhasei;
+    }
 
-    // Add the explicit part of the drag force
     forAllConstIter(KdfTable, Kdfs_, KdfIter)
     {
         const surfaceScalarField& Kf(*KdfIter());
@@ -889,111 +1001,55 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdPhifs() const
             const phaseModel& phase = iter();
             const phaseModel& otherPhase = iter.otherPhase();
 
+            const label i = movingPhases[phase.index()];
+
+            if (i != -1)
+            {
                 const surfaceScalarField alphaf
                 (
                     fvc::interpolate(otherPhase)
                 );
 
-            addField
+                const surfaceScalarField Kfij
                 (
-                phase,
+                    (alphaf/max(alphaf, otherPhase.residualAlpha()))*Kf
-                "KdPhif",
+                );
-               -(alphaf/max(alphaf, otherPhase.residualAlpha()))
+
-               *Kf
+                invADfs[i][i] += Kfij;
-               *fvc::absolute
+
+                const label j = movingPhases[otherPhase.index()];
+
+                if (j != -1)
+                {
+                    invADfs[i].set(j, -Kfij);
+                }
+            }
+        }
+    }
+
+    for (label i=0; i<n; i++)
+    {
+        for (label j=0; j<n; j++)
+        {
+            if (!invADfs[i].set(j))
+            {
+                invADfs[i].set
                 (
-                    this->MRF().absolute(otherPhase.phi()),
+                    j,
-                    otherPhase.U()
+                    surfaceScalarField::New
-                ),
+                    (
-                KdPhifs
+                        "0",
+                        this->mesh(),
+                        dimensionedScalar(As[0].dimensions(), 0)
+                    )
                 );
             }
         }
-
-    if (this->fillFields_)
-    {
-        this->fillFields
-        (
-            "KdPhif",
-            dimArea*dimDensity*dimAcceleration,
-            KdPhifs
-        );
     }
 
-    return KdPhifs;
+    invADs(invADfs);
-}
 
-
+    return invADfs;
-template<class BasePhaseSystem>
-Foam::PtrList<Foam::volScalarField>
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kds() const
-{
-    PtrList<volScalarField> Kds(this->phaseModels_.size());
-
-    forAllConstIter(KdTable, Kds_, KdIter)
-    {
-        const volScalarField& K(*KdIter());
-        const phaseInterface interface(*this, KdIter.key());
-
-        forAllConstIter(phaseInterface, interface, iter)
-        {
-            const phaseModel& phase = iter();
-            const phaseModel& otherPhase = iter.otherPhase();
-
-            addField
-            (
-                phase,
-                "Kd",
-                (otherPhase/max(otherPhase, otherPhase.residualAlpha()))
-               *K,
-                Kds
-            );
-        }
-    }
-
-    if (this->fillFields_)
-    {
-        this->fillFields("Kd", dimDensity/dimTime, Kds);
-    }
-
-    return Kds;
-}
-
-
-template<class BasePhaseSystem>
-Foam::PtrList<Foam::volVectorField>
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::KdUs() const
-{
-    PtrList<volVectorField> KdUs(this->phaseModels_.size());
-
-    // Add the explicit part of the drag force
-    forAllConstIter(KdTable, Kds_, KdIter)
-    {
-        const volScalarField& K(*KdIter());
-        const phaseInterface interface(*this, KdIter.key());
-
-        forAllConstIter(phaseInterface, interface, iter)
-        {
-            const phaseModel& phase = iter();
-            const phaseModel& otherPhase = iter.otherPhase();
-
-            addField
-            (
-                phase,
-                "KdU",
-               -(otherPhase/max(otherPhase, otherPhase.residualAlpha()))
-               *K*otherPhase.U(),
-                KdUs
-            );
-        }
-    }
-
-    if (this->fillFields_)
-    {
-        this->fillFields("KdU", dimDensity*dimAcceleration, KdUs);
-    }
-
-    return KdUs;
 }
 
 
@@ -1032,8 +1088,7 @@ template<class BasePhaseSystem>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::alphaDByAf
 (
-    const PtrList<volScalarField>& rAUs,
+    const PtrList<volScalarField>& rAs
-    const PtrList<surfaceScalarField>& rAUfs
 ) const
 {
     tmp<surfaceScalarField> alphaDByAf;
@@ -1049,15 +1104,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::alphaDByAf
         (
             alphaDByAf,
             fvc::interpolate(max(phase, scalar(0)))
-           *(
+           *fvc::interpolate(rAs[phasei]*pPrime(), pPrime().name())
-                rAUfs.size()
-
-                // Face-momentum form
-              ? rAUfs[phasei]*fvc::interpolate(pPrime())
-
-                // Cell-momentum form
-              : fvc::interpolate(rAUs[phasei]*pPrime(), pPrime().name())
-            )
         );
     }
 
@@ -1087,28 +1134,15 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::alphaDByAf
             alphaDByAf,
             alpha1f*alpha2f
            /max(alpha1f + alpha2f, interface.phase1().residualAlpha())
-           *(
+           *fvc::interpolate
-                rAUfs.size()
-
-                // Face-momentum form
-              ? max
-                (
-                    rAUfs[interface.phase1().index()],
-                    rAUfs[interface.phase2().index()]
-                )
-               *fvc::interpolate(turbulentDispersionModelIter()->D())
-
-                // Cell-momentum form
-              : fvc::interpolate
             (
                 max
                 (
-                        rAUs[interface.phase1().index()],
+                    rAs[interface.phase1().index()],
-                        rAUs[interface.phase2().index()]
+                    rAs[interface.phase2().index()]
                 )
                *turbulentDispersionModelIter()->D()
             )
-            )
         );
     }
 
@@ -1232,405 +1266,59 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::dragCorrs
     PtrList<surfaceScalarField>& dragCorrfs
 ) const
 {
-    const phaseSystem::phaseModelList& phases = this->phaseModels_;
+    labelList movingPhases(this->phases().size(), -1);
+    PtrList<volVectorField> Uphis(this->movingPhases().size());
 
-    PtrList<volVectorField> Uphis(phases.size());
+    forAll(this->movingPhases(), movingPhasei)
+    {
+        movingPhases[this->movingPhases()[movingPhasei].index()] = movingPhasei;
 
-    forAll(phases, i)
-    {
-        if (!phases[i].stationary())
-        {
         Uphis.set
         (
-                i,
+            movingPhasei,
-                fvc::reconstruct(phases[i].phi())
+            fvc::reconstruct(this->movingPhases()[movingPhasei].phi())
         );
     }
-    }
 
     forAllConstIter(KdTable, Kds_, KdIter)
     {
         const volScalarField& K(*KdIter());
         const phaseInterface interface(*this, KdIter.key());
 
-        const phaseModel& phase1 = interface.phase1();
+        forAllConstIter(phaseInterface, interface, iter)
-        const phaseModel& phase2 = interface.phase2();
+        {
+            const phaseModel& phase = iter();
+            const phaseModel& otherPhase = iter.otherPhase();
 
-        const label phase1i = phase1.index();
+            const label i = movingPhases[phase.index()];
-        const label phase2i = phase2.index();
+            const label j = movingPhases[otherPhase.index()];
 
-        if (!phase1.stationary())
+            if (i != -1)
             {
                 const volScalarField K1
                 (
-                (phase2/max(phase2, phase2.residualAlpha()))*K
+                    (otherPhase/max(otherPhase, otherPhase.residualAlpha()))*K
                 );
 
                 addField
                 (
-                phase1,
+                    i,
                     "dragCorr",
-                K1
+                    K1*(j == -1 ? -Uphis[i] : (Uphis[j] - Uphis[i])),
-               *(
-                    phase2.stationary()
-                 ? -Uphis[phase1i]
-                 :  (Uphis[phase2i] - Uphis[phase1i])
-                ),
                     dragCorrs
                 );
 
                 addField
                 (
-                phase1,
+                    i,
                     "dragCorrf",
                     fvc::interpolate(K1)
-               *(
+                   *(j == -1 ? -phase.phi() : (otherPhase.phi() - phase.phi())),
-                    phase2.stationary()
-                 ? -phase1.phi()
-                 :  (phase2.phi() - phase1.phi())
-                ),
-                dragCorrfs
-            );
-        }
-
-        if (!phase2.stationary())
-        {
-            const volScalarField K2
-            (
-                (phase1/max(phase1, phase1.residualAlpha()))*K
-            );
-
-            addField
-            (
-                phase2,
-                "dragCorr",
-                K2
-               *(
-                    phase1.stationary()
-                 ? -Uphis[phase2i]
-                 : (Uphis[phase1i] - Uphis[phase2i])
-                ),
-                dragCorrs
-            );
-
-            addField
-            (
-                phase2,
-                "dragCorrf",
-                fvc::interpolate(K2)
-               *(
-                    phase1.stationary()
-                 ? -phase2.phi()
-                 :  (phase1.phi() - phase2.phi())
-                ),
                     dragCorrfs
                 );
             }
         }
     }
-
-
-template<class BasePhaseSystem>
-void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::partialElimination
-(
-    const PtrList<volScalarField>& rAUs,
-    const PtrList<volVectorField>& KdUs,
-    const PtrList<surfaceScalarField>& alphafs,
-    const PtrList<surfaceScalarField>& rAUfs,
-    const PtrList<surfaceScalarField>& KdPhis
-)
-{
-    Info<< "Inverting drag systems: ";
-
-    phaseSystem::phaseModelList& phases = this->phaseModels_;
-
-    // Remove the drag contributions from the velocity and flux of the phases
-    // in preparation for the partial elimination of these terms
-    forAll(phases, i)
-    {
-        if (!phases[i].stationary())
-        {
-            phases[i].URef() += rAUs[i]*KdUs[i];
-            phases[i].phiRef() += rAUfs[i]*KdPhis[i];
-        }
-    }
-
-    {
-        // Create drag coefficient matrices
-        PtrList<PtrList<volScalarField>> KdByAs(phases.size());
-        PtrList<PtrList<surfaceScalarField>> KdByAfs(phases.size());
-
-        forAll(phases, i)
-        {
-            KdByAs.set
-            (
-                i,
-                new PtrList<volScalarField>(phases.size())
-            );
-
-            KdByAfs.set
-            (
-                i,
-                new PtrList<surfaceScalarField>(phases.size())
-            );
-        }
-
-        forAllConstIter(KdTable, Kds_, KdIter)
-        {
-            const volScalarField& K(*KdIter());
-            const phaseInterface interface(*this, KdIter.key());
-
-            const label phase1i = interface.phase1().index();
-            const label phase2i = interface.phase2().index();
-
-            const volScalarField K1
-            (
-                interface.phase2()
-               /max(interface.phase2(), interface.phase2().residualAlpha())
-               *K
-            );
-
-            const volScalarField K2
-            (
-                interface.phase1()
-               /max(interface.phase1(), interface.phase1().residualAlpha())
-               *K
-            );
-
-            addField
-            (
-                interface.phase2(),
-                "KdByA",
-               -rAUs[phase1i]*K1,
-                KdByAs[phase1i]
-            );
-            addField
-            (
-                interface.phase1(),
-                "KdByA",
-               -rAUs[phase2i]*K2,
-                KdByAs[phase2i]
-            );
-
-            addField
-            (
-                interface.phase2(),
-                "KdByAf",
-               -rAUfs[phase1i]*fvc::interpolate(K1),
-                KdByAfs[phase1i]
-            );
-            addField
-            (
-                interface.phase1(),
-                "KdByAf",
-               -rAUfs[phase2i]*fvc::interpolate(K2),
-                KdByAfs[phase2i]
-            );
-        }
-
-        forAll(phases, i)
-        {
-            this->fillFields("KdByAs", dimless, KdByAs[i]);
-            this->fillFields("KdByAfs", dimless, KdByAfs[i]);
-
-            KdByAs[i][i] = 1;
-            KdByAfs[i][i] = 1;
-        }
-
-        // Decompose
-        for (label i = 0; i < phases.size(); i++)
-        {
-            for (label j = i + 1; j < phases.size(); j++)
-            {
-                KdByAs[j][i] /= KdByAs[i][i];
-                KdByAfs[j][i] /= KdByAfs[i][i];
-                for (label k = i + 1; k < phases.size(); ++ k)
-                {
-                    KdByAs[j][k] -= KdByAs[j][i]*KdByAs[i][k];
-                    KdByAfs[j][k] -= KdByAfs[j][i]*KdByAfs[i][k];
-                }
-            }
-        }
-
-        {
-            volScalarField detKdByAs(KdByAs[0][0]);
-            surfaceScalarField detPhiKdfs(KdByAfs[0][0]);
-
-            for (label i = 1; i < phases.size(); i++)
-            {
-                detKdByAs *= KdByAs[i][i];
-                detPhiKdfs *= KdByAfs[i][i];
-            }
-
-            Info<< "Min cell/face det = " << gMin(detKdByAs.primitiveField())
-                << "/" << gMin(detPhiKdfs.primitiveField()) << endl;
-        }
-
-        // Solve for the velocities and fluxes
-        for (label i = 1; i < phases.size(); i++)
-        {
-            if (!phases[i].stationary())
-            {
-                for (label j = 0; j < i; j ++)
-                {
-                    phases[i].URef() -= KdByAs[i][j]*phases[j].U();
-                    phases[i].phiRef() -= KdByAfs[i][j]*phases[j].phi();
-                }
-            }
-        }
-        for (label i = phases.size() - 1; i >= 0; i--)
-        {
-            if (!phases[i].stationary())
-            {
-                for (label j = phases.size() - 1; j > i; j--)
-                {
-                    phases[i].URef() -= KdByAs[i][j]*phases[j].U();
-                    phases[i].phiRef() -= KdByAfs[i][j]*phases[j].phi();
-                }
-                phases[i].URef() /= KdByAs[i][i];
-                phases[i].phiRef() /= KdByAfs[i][i];
-            }
-        }
-    }
-
-    // Adjust the phase-fluxes such that the mean flux
-    // obtained from the pressure solution is retained
-    this->setMixturePhi(alphafs, this->phi());
-}
-
-
-template<class BasePhaseSystem>
-void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::partialEliminationf
-(
-    const PtrList<surfaceScalarField>& rAUfs,
-    const PtrList<surfaceScalarField>& alphafs,
-    const PtrList<surfaceScalarField>& KdPhifs
-)
-{
-    Info<< "Inverting drag system: ";
-
-    phaseSystem::phaseModelList& phases = this->phaseModels_;
-
-    // Remove the drag contributions from the flux of the phases
-    // in preparation for the partial elimination of these terms
-    forAll(phases, i)
-    {
-        if (!phases[i].stationary())
-        {
-            phases[i].phiRef() += rAUfs[i]*KdPhifs[i];
-        }
-    }
-
-    {
-        // Create drag coefficient matrix
-        PtrList<PtrList<surfaceScalarField>> phiKdfs(phases.size());
-
-        forAll(phases, phasei)
-        {
-            phiKdfs.set
-            (
-                phasei,
-                new PtrList<surfaceScalarField>(phases.size())
-            );
-        }
-
-        forAllConstIter(KdfTable, Kdfs_, KdfIter)
-        {
-            const surfaceScalarField& Kf(*KdfIter());
-            const phaseInterface interface(*this, KdfIter.key());
-
-            const label phase1i = interface.phase1().index();
-            const label phase2i = interface.phase2().index();
-
-            const surfaceScalarField alpha1f
-            (
-                fvc::interpolate(interface.phase1())
-            );
-
-            const surfaceScalarField alpha2f
-            (
-                fvc::interpolate(interface.phase2())
-            );
-
-            addField
-            (
-                interface.phase2(),
-                "phiKdf",
-               -rAUfs[phase1i]
-               *alpha2f/max(alpha2f, interface.phase2().residualAlpha())
-               *Kf,
-                phiKdfs[phase1i]
-            );
-            addField
-            (
-                interface.phase1(),
-                "phiKdf",
-               -rAUfs[phase2i]
-               *alpha1f/max(alpha1f, interface.phase1().residualAlpha())
-               *Kf,
-                phiKdfs[phase2i]
-            );
-        }
-
-        forAll(phases, phasei)
-        {
-            this->fillFields("phiKdf", dimless, phiKdfs[phasei]);
-
-            phiKdfs[phasei][phasei] = 1;
-        }
-
-        // Decompose
-        for (label i = 0; i < phases.size(); i++)
-        {
-            for (label j = i + 1; j < phases.size(); j++)
-            {
-                phiKdfs[j][i] /= phiKdfs[i][i];
-                for (label k = i + 1; k < phases.size(); ++ k)
-                {
-                    phiKdfs[j][k] -= phiKdfs[j][i]*phiKdfs[i][k];
-                }
-            }
-        }
-
-        {
-            surfaceScalarField detPhiKdfs(phiKdfs[0][0]);
-
-            for (label i = 1; i < phases.size(); i++)
-            {
-                detPhiKdfs *= phiKdfs[i][i];
-            }
-
-            Info<< "Min face det = "
-                << gMin(detPhiKdfs.primitiveField()) << endl;
-        }
-
-        // Solve for the fluxes
-        for (label i = 1; i < phases.size(); i++)
-        {
-            if (!phases[i].stationary())
-            {
-                for (label j = 0; j < i; j ++)
-                {
-                    phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi();
-                }
-            }
-        }
-        for (label i = phases.size() - 1; i >= 0; i--)
-        {
-            if (!phases[i].stationary())
-            {
-                for (label j = phases.size() - 1; j > i; j--)
-                {
-                    phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi();
-                }
-                phases[i].phiRef() /= phiKdfs[i][i];
-            }
-        }
-    }
-
-    // Adjust the phase-fluxes such that the mean flux
-    // obtained from the pressure solution is retained
-    this->setMixturePhi(alphafs, this->phi());
 }
 
 

applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H

@@ -169,6 +169,16 @@ protected:
             const tmp<surfaceScalarField>& field
         ) const;
 
+        //- Invert the ADs matrix inplace
+        void invADs(List<UPtrList<scalarField>>& ADs) const;
+
+        //- Invert the ADs matrix inplace
+        template<template<class> class PatchField, class GeoMesh>
+        void invADs
+        (
+            PtrList<PtrList<GeometricField<scalar, PatchField, GeoMesh>>>& ADs
+        ) const;
+
 
 public:
 
@@ -194,7 +204,7 @@ public:
 
         //- Return implicit force coefficients on the faces, for the face-based
         //  algorithm.
-        virtual PtrList<surfaceScalarField> KdVmfs() const;
+        virtual PtrList<surfaceScalarField> Vmfs() const;
 
         //- Return the explicit force fluxes for the cell-based algorithm, that
         //  do not depend on phase mass/volume fluxes, and can therefore be
@@ -205,22 +215,19 @@ public:
         //- As Fs, but for the face-based algorithm
         virtual PtrList<surfaceScalarField> Ffs() const;
 
-        //- Return the explicit drag force fluxes for the cell-based algorithm.
+        //- Return the inverse of (As + Ds)
-        //  These depend on phase mass/volume fluxes, and must therefore be
+        virtual void invADs
-        //  evaluated inside the corrector loop.
+        (
-        virtual PtrList<surfaceScalarField> KdPhis() const;
+            const PtrList<volScalarField>& As,
+            PtrList<PtrList<volScalarField>>& invADs,
+            PtrList<PtrList<surfaceScalarField>>& invADfs
+        ) const;
 
-        //- As KdPhis, but for the face-based algorithm
+        //- Return the inverse of (Afs + Dfs)
-        virtual PtrList<surfaceScalarField> KdPhifs() const;
+        virtual PtrList<PtrList<surfaceScalarField>> invADfs
-
+        (
-        //- Return the implicit part of the drag force
+            const PtrList<surfaceScalarField>& Afs
-        virtual PtrList<volScalarField> Kds() const;
+        ) const;
-
-        //- Return the explicit part of the drag force for the cell-based
-        //  algorithm. This is the cell-equivalent of KdPhis. These depend on
-        //  phase velocities, and must therefore be evaluated inside the
-        //  corrector loop.
-        virtual PtrList<volVectorField> KdUs() const;
 
         //- Returns true if the phase pressure is treated implicitly
         //  in the phase fraction equation
@@ -234,8 +241,7 @@ public:
         //  divided by the momentum central coefficient
         virtual tmp<surfaceScalarField> alphaDByAf
         (
-            const PtrList<volScalarField>& rAUs,
+            const PtrList<volScalarField>& rAs
-            const PtrList<surfaceScalarField>& rAUfs
         ) const;
 
         //- Return the flux corrections for the cell-based algorithm. These
@@ -250,25 +256,6 @@ public:
             PtrList<surfaceScalarField>& dragCorrf
         ) const;
 
-        //- Solve the drag system for the velocities and fluxes
-        virtual void partialElimination
-        (
-            const PtrList<volScalarField>& rAUs,
-            const PtrList<volVectorField>& KdUs,
-            const PtrList<surfaceScalarField>& alphafs,
-            const PtrList<surfaceScalarField>& rAUfs,
-            const PtrList<surfaceScalarField>& KdPhis
-        );
-
-        //- As partialElimination, but for the face-based algorithm. Only solves
-        //  for the fluxes.
-        virtual void partialEliminationf
-        (
-            const PtrList<surfaceScalarField>& rAUfs,
-            const PtrList<surfaceScalarField>& alphafs,
-            const PtrList<surfaceScalarField>& KdPhifs
-        );
-
         //- Read base phaseProperties dictionary
         virtual bool read();
 };

applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C

@@ -175,11 +175,10 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::specieTransfer() const
 template<class BasePhaseSystem>
 void Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::solve
 (
-    const PtrList<volScalarField>& rAUs,
+    const PtrList<volScalarField>& rAs
-    const PtrList<surfaceScalarField>& rAUfs
 )
 {
-    BasePhaseSystem::solve(rAUs, rAUfs);
+    BasePhaseSystem::solve(rAs);
 
     forAll(populationBalances_, i)
     {

applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H

@@ -97,11 +97,7 @@ public:
             specieTransfer() const;
 
         //- Solve all population balance equations
-        virtual void solve
+        virtual void solve(const PtrList<volScalarField>& rAs);
-        (
-            const PtrList<volScalarField>& rAUs,
-            const PtrList<surfaceScalarField>& rAUfs
-        );
 
         //- Correct derived properties
         virtual void correct();

applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H

@@ -523,7 +523,7 @@ public:
 
             //- Return the implicit force coefficients for the face-based
             //  algorithm
-            virtual PtrList<surfaceScalarField> KdVmfs() const = 0;
+            virtual PtrList<surfaceScalarField> Vmfs() const = 0;
 
             //- Return the force fluxes for the cell-based algorithm
             virtual PtrList<surfaceScalarField> Fs() const = 0;
@@ -531,17 +531,19 @@ public:
             //- Return the force fluxes for the face-based algorithm
             virtual PtrList<surfaceScalarField> Ffs() const = 0;
 
-            //- Return the force fluxes for the cell-based algorithm
+            //- Return the inverse of (As + Ds)
-            virtual PtrList<surfaceScalarField> KdPhis() const = 0;
+            virtual void invADs
+            (
+                const PtrList<volScalarField>& As,
+                PtrList<PtrList<volScalarField>>& invADs,
+                PtrList<PtrList<surfaceScalarField>>& invADfs
+            ) const = 0;
 
-            //- Return the force fluxes for the face-based algorithm
+            //- Return the inverse of (Afs + Dfs)
-            virtual PtrList<surfaceScalarField> KdPhifs() const = 0;
+            virtual PtrList<PtrList<surfaceScalarField>> invADfs
-
+            (
-            //- Return the implicit part of the drag force
+                const PtrList<surfaceScalarField>& Afs
-            virtual PtrList<volScalarField> Kds() const = 0;
+            ) const = 0;
-
-            //- Return the explicit part of the drag force
-            virtual PtrList<volVectorField> KdUs() const = 0;
 
             //- Returns true if the phase pressure is treated implicitly
             //  in the phase fraction equation
@@ -555,8 +557,7 @@ public:
             //  divided by the momentum central coefficient
             virtual tmp<surfaceScalarField> alphaDByAf
             (
-                const PtrList<volScalarField>& rAUs,
+                const PtrList<volScalarField>& rAs
-                const PtrList<surfaceScalarField>& rAUfs
             ) const = 0;
 
             //- Return the flux corrections for the cell-based algorithm
@@ -569,24 +570,6 @@ public:
                 PtrList<surfaceScalarField>& dragCorrf
             ) const = 0;
 
-            //- Solve the drag system for the new velocities and fluxes
-            virtual void partialElimination
-            (
-                const PtrList<volScalarField>& rAUs,
-                const PtrList<volVectorField>& KdUs,
-                const PtrList<surfaceScalarField>& alphafs,
-                const PtrList<surfaceScalarField>& rAUfs,
-                const PtrList<surfaceScalarField>& KdPhis
-            ) = 0;
-
-            //- Solve the drag system for the new fluxes
-            virtual void partialEliminationf
-            (
-                const PtrList<surfaceScalarField>& rAUfs,
-                const PtrList<surfaceScalarField>& alphafs,
-                const PtrList<surfaceScalarField>& KdPhifs
-            ) = 0;
-
             //- Re-normalise the flux of the phases
             //  around the specified mixture mean
             void setMixturePhi
@@ -611,11 +594,7 @@ public:
         // Evolution
 
             //- Solve for the phase fractions
-            virtual void solve
+            virtual void solve(const PtrList<volScalarField>& rAs);
-            (
-                const PtrList<volScalarField>& rAUs,
-                const PtrList<surfaceScalarField>& rAUfs
-            );
 
             //- Correct the fluid properties other than those listed below
             virtual void correct();

applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H

@@ -23,6 +23,67 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+class wordListAndType
+{
+    public:
+
+        wordList wl;
+        Type t;
+
+        wordListAndType()
+        {}
+
+        wordListAndType(Istream& is)
+        :
+            wl(is),
+            t(is)
+        {}
+};
+
+template<class Type>
+inline Istream& operator>>(Istream& is, wordListAndType<Type>& wlat)
+{
+    return is >> wlat.wl >> wlat.t;
+}
+
+template<class Type>
+inline Ostream& operator<<(Ostream& os, const wordListAndType<Type>& wlat)
+{
+    return os << wlat.wl << wlat.t;
+}
+
+template<class Type>
+inline bool operator==
+(
+    const wordListAndType<Type>& a,
+    const wordListAndType<Type>& b
+)
+{
+    return a.wl == b.wl && a.t == b.t;
+}
+
+template<class Type>
+inline bool operator!=
+(
+    const wordListAndType<Type>& a,
+    const wordListAndType<Type>& b
+)
+{
+    return !(a == b);
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline const Foam::fvMesh& Foam::phaseSystem::mesh() const
@@ -190,4 +251,151 @@ inline const Foam::dimensionedScalar& Foam::phaseSystem::deltaN() const
 }
 
 
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class GeoField>
+inline void addField
+(
+    const label phasei,
+    const word& name,
+    tmp<GeoField> field,
+    PtrList<GeoField>& fieldList
+)
+{
+    if (fieldList.set(phasei))
+    {
+        fieldList[phasei] += field;
+    }
+    else
+    {
+        fieldList.set
+        (
+            phasei,
+            new GeoField
+            (
+                name, // IOobject::groupName(name, group.name()),
+                field
+            )
+        );
+    }
+}
+
+
+template<class GeoField, class Group>
+inline void addField
+(
+    const Group& group,
+    const word& name,
+    tmp<GeoField> field,
+    PtrList<GeoField>& fieldList
+)
+{
+    if (fieldList.set(group.index()))
+    {
+        fieldList[group.index()] += field;
+    }
+    else
+    {
+        fieldList.set
+        (
+            group.index(),
+            new GeoField
+            (
+                IOobject::groupName(name, group.name()),
+                field
+            )
+        );
+    }
+}
+
+
+template<class GeoField, class Group>
+inline void addField
+(
+    const Group& group,
+    const word& name,
+    const GeoField& field,
+    PtrList<GeoField>& fieldList
+)
+{
+    addField(group, name, tmp<GeoField>(field), fieldList);
+}
+
+
+template<class GeoField, class Group>
+inline void addField
+(
+    const Group& group,
+    const word& name,
+    tmp<GeoField> field,
+    HashPtrTable<GeoField>& fieldTable
+)
+{
+    if (fieldTable.found(group.name()))
+    {
+        *fieldTable[group.name()] += field;
+    }
+    else
+    {
+        fieldTable.set
+        (
+            group.name(),
+            new GeoField
+            (
+                IOobject::groupName(name, group.name()),
+                field
+            )
+        );
+    }
+}
+
+
+template<class GeoField, class Group>
+inline void addField
+(
+    const Group& group,
+    const word& name,
+    const GeoField& field,
+    HashPtrTable<GeoField>& fieldTable
+)
+{
+    addField(group, name, tmp<GeoField>(field), fieldTable);
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+PtrList<GeometricField<Type, PatchField, GeoMesh>> operator&
+(
+    const PtrList<PtrList<GeometricField<scalar, PatchField, GeoMesh>>>& As,
+    const PtrList<GeometricField<Type, PatchField, GeoMesh>>& fs
+)
+{
+    PtrList<GeometricField<Type, PatchField, GeoMesh >> Afs(fs.size());
+
+    forAll(Afs, i)
+    {
+        Afs.set(i, As[i][i]*fs[i]);
+
+        forAll(Afs, j)
+        {
+            if (j != i)
+            {
+                Afs[i] += As[i][j]*fs[j];
+            }
+        }
+    }
+
+    return Afs;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 // ************************************************************************* //

applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C

@@ -41,11 +41,7 @@ License
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::phaseSystem::solve
+void Foam::phaseSystem::solve(const PtrList<volScalarField>& rAs)
-(
-    const PtrList<volScalarField>& rAUs,
-    const PtrList<surfaceScalarField>& rAUfs
-)
 {
     const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
 
@@ -251,9 +247,9 @@ void Foam::phaseSystem::solve
             PtrList<surfaceScalarField> alphaPhis(phases().size());
 
             tmp<surfaceScalarField> alphaDByAf;
-            if (implicitPhasePressure() && (rAUs.size() || rAUfs.size()))
+            if (implicitPhasePressure() && (rAs.size()))
             {
-                alphaDByAf = this->alphaDByAf(rAUs, rAUfs);
+                alphaDByAf = this->alphaDByAf(rAs);
             }
 
             forAll(movingPhases(), movingPhasei)
@@ -454,7 +450,7 @@ void Foam::phaseSystem::solve
             if (alphaDByAf.valid())
             {
                 // Update alphaDByAf due to changes in alpha
-                alphaDByAf = this->alphaDByAf(rAUs, rAUfs);
+                alphaDByAf = this->alphaDByAf(rAs);
 
                 forAll(solvePhases, solvePhasei)
                 {

applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C

@@ -28,62 +28,6 @@ License
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-namespace Foam
-{
-
-template<class Type>
-class wordListAndType
-{
-    public:
-
-        wordList wl;
-        Type t;
-
-        wordListAndType()
-        {}
-
-        wordListAndType(Istream& is)
-        :
-            wl(is),
-            t(is)
-        {}
-};
-
-template<class Type>
-inline Istream& operator>>(Istream& is, wordListAndType<Type>& wlat)
-{
-    return is >> wlat.wl >> wlat.t;
-}
-
-template<class Type>
-inline Ostream& operator<<(Ostream& os, const wordListAndType<Type>& wlat)
-{
-    return os << wlat.wl << wlat.t;
-}
-
-template<class Type>
-inline bool operator==
-(
-    const wordListAndType<Type>& a,
-    const wordListAndType<Type>& b
-)
-{
-    return a.wl == b.wl && a.t == b.t;
-}
-
-template<class Type>
-inline bool operator!=
-(
-    const wordListAndType<Type>& a,
-    const wordListAndType<Type>& b
-)
-{
-    return !(a == b);
-}
-
-}
-
-
 template<class Type>
 Foam::dictionary Foam::phaseSystem::interfacialDict(const word& name) const
 {
@@ -362,7 +306,6 @@ void Foam::phaseSystem::generateInterfacialModels
 }
 
 
-
 template<class ModelType>
 void Foam::phaseSystem::generateInterfacialModels
 (
@@ -483,96 +426,4 @@ const ModelType& Foam::phaseSystem::lookupInterfacialModel
 }
 
 
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class GeoField, class Group>
-inline void addField
-(
-    const Group& group,
-    const word& name,
-    tmp<GeoField> field,
-    PtrList<GeoField>& fieldList
-)
-{
-    if (fieldList.set(group.index()))
-    {
-        fieldList[group.index()] += field;
-    }
-    else
-    {
-        fieldList.set
-        (
-            group.index(),
-            new GeoField
-            (
-                IOobject::groupName(name, group.name()),
-                field
-            )
-        );
-    }
-}
-
-
-template<class GeoField, class Group>
-inline void addField
-(
-    const Group& group,
-    const word& name,
-    const GeoField& field,
-    PtrList<GeoField>& fieldList
-)
-{
-    addField(group, name, tmp<GeoField>(field), fieldList);
-}
-
-
-template<class GeoField, class Group>
-inline void addField
-(
-    const Group& group,
-    const word& name,
-    tmp<GeoField> field,
-    HashPtrTable<GeoField>& fieldTable
-)
-{
-    if (fieldTable.found(group.name()))
-    {
-        *fieldTable[group.name()] += field;
-    }
-    else
-    {
-        fieldTable.set
-        (
-            group.name(),
-            new GeoField
-            (
-                IOobject::groupName(name, group.name()),
-                field
-            )
-        );
-    }
-}
-
-
-template<class GeoField, class Group>
-inline void addField
-(
-    const Group& group,
-    const word& name,
-    const GeoField& field,
-    HashPtrTable<GeoField>& fieldTable
-)
-{
-    addField(group, name, tmp<GeoField>(field), fieldTable);
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -90,6 +90,11 @@ public:
             return
                 word("SymmetricSquareMatrix<") + pTraits<Type>::typeName + '>';
         }
+
+
+    // Member Operators
+
+        using Matrix<SymmetricSquareMatrix<Type>, Type>::operator=;
 };
 
 

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleSingleRegionControl/pimpleSingleRegionControl.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,7 +45,7 @@ Foam::pimpleSingleRegionControl::pimpleSingleRegionControl
 {
     pimple_.pimpleLoopPtr_ = this;
 
-    read();
+    pimpleLoop::read();
 
     pimple_.printResidualControls();
 

src/finiteVolume/fvMesh/fvMesh.H

@@ -536,7 +536,7 @@ public:
             //- Return the fvSchemes
             const fvSchemes& schemes() const;
 
-            //- Return the fvSchemes
+            //- Return the fvSolution
             const fvSolution& solution() const;
 
 

tutorials/multiRegion/CHT/wallBoiling/system/fluid/fvSolution

@@ -68,7 +68,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/Grossetete/system/fvSolution

@@ -77,7 +77,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/LBend/system/fvSolution

@@ -92,7 +92,6 @@ PIMPLE
     faceMomentum        yes;
     VmDdtCorrection     no;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bed/system/fvSolution

@@ -61,6 +61,7 @@ PIMPLE
     nCorrectors         1;
     nNonOrthogonalCorrectors 0;
 
+    dragCorrection      yes;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bubbleColumn/system/fvSolution

@@ -73,7 +73,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bubbleColumnEvaporating/system/fvSolution

@@ -73,7 +73,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bubbleColumnEvaporatingDissolving/system/fvSolution

@@ -70,7 +70,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bubbleColumnEvaporatingReacting/system/fvSolution

@@ -82,7 +82,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bubbleColumnIATE/system/fvSolution

@@ -63,7 +63,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bubbleColumnLES/system/fvSolution

@@ -73,7 +73,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bubbleColumnLaminar/system/fvSolution

@@ -64,7 +64,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/bubblePipe/system/fvSolution

@@ -73,7 +73,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/damBreak4phase/system/fvSolution

@@ -79,8 +79,7 @@ PIMPLE
 
     faceMomentum        no;
     VmDdtCorrection     no;
-    dragCorrection      no;
+    dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/fluidisedBed/system/fvSolution

@@ -90,7 +90,6 @@ PIMPLE
 
     faceMomentum        no;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/fluidisedBedLaminar/system/fvSolution

@@ -89,8 +89,7 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 
     faceMomentum        no;
-    dragCorrection      yes;
+    dragCorrection      no;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/hydrofoil/system/fvSolution

@@ -83,7 +83,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     no;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 cache

tutorials/multiphaseEuler/injection/system/fvSolution

@@ -64,7 +64,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/mixerVessel2D/0/alpha.map

@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [];
+dimensions      [0 0 0 0 0 0 0];
 
 internalField   uniform 0.6;
 

tutorials/multiphaseEuler/mixerVessel2D/system/fvSolution

@@ -98,7 +98,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/mixerVessel2DMRF/system/fvSolution

@@ -85,7 +85,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/pipeBend/system/fvSolution

@@ -89,7 +89,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/steamInjection/system/fvSolution

@@ -86,7 +86,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  yes;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/titaniaSynthesis/system/fvSolution

@@ -86,7 +86,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     no;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/titaniaSynthesisSurface/system/fvSolution

@@ -86,7 +86,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     no;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/wallBoilingIATE/system/fvSolution

@@ -68,6 +68,7 @@ solvers
 
 PIMPLE
 {
+    momentumPredictor   no;
     nOuterCorrectors    1;
     nCorrectors         1;
     nNonOrthogonalCorrectors 0;
@@ -76,19 +77,19 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors
 {
     fields
     {
-        thermalPhaseChange:dmdtf 1.0;
+        thermalPhaseChange:dmdtf 1;
     }
 
     equations
     {
         ".*"            1;
+        "U\..*"         1;
         "h\..*"         1;
     }
 }

tutorials/multiphaseEuler/wallBoilingPolydisperse/system/fvSolution

@@ -85,7 +85,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors

tutorials/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/fvSolution

@@ -85,7 +85,6 @@ PIMPLE
     faceMomentum        no;
     VmDdtCorrection     yes;
     dragCorrection      yes;
-    partialElimination  no;
 }
 
 relaxationFactors