GIT: Resolve conflict associated with cherry-pick of Foundation commit 79ff9135

79ff9135 - rhoPimpleFoam: Improved support for compressible liquids
(2017-05-17 17:05:43 +0100) <Henry Weller>

applications/solvers/combustion/reactingFoam/createFields.H

@@ -52,27 +52,7 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
+pressureControl pressureControl(p, rho, pimple.dict(), false);
 
 mesh.setFluxRequired(p.name());
 

applications/solvers/combustion/reactingFoam/pEqn.H

 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -87,19 +84,17 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);

applications/solvers/combustion/reactingFoam/pcEqn.H

 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
 
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@@ -109,19 +106,13 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-if (thermo.dpdt())
+if (pressureControl.limit(p))
 {
-    dpdt = fvc::ddt(p);
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
 }
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-
-if (!pimple.transonic())
+if (thermo.dpdt())
 {
-    rho.relax();
+    dpdt = fvc::ddt(p);
 }
-
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -114,6 +115,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H

-{
-    rho = thermo.rho();
+rho = thermo.rho();
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
-    surfaceScalarField phiHbyA
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-      + phig
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-    fvScalarMatrix p_rghDDtEqn
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+  + phig
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+  + fvc::div(phiHbyA)
+  ==
+    fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix p_rghEqn
     (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phiHbyA)
-      ==
-        fvOptions(psi, p_rgh, rho.name())
+        p_rghDDtEqn
+      - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    while (pimple.correctNonOrthogonal())
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
     {
-        fvScalarMatrix p_rghEqn
-        (
-            p_rghDDtEqn
-          - fvm::laplacian(rhorAUf, p_rgh)
-        );
-
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            // Calculate the conservative fluxes
-            phi = phiHbyA + p_rghEqn.flux();
-
-            // Explicitly relax pressure for momentum corrector
-            p_rgh.relax();
-
-            // Correct the momentum source with the pressure gradient flux
-            // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-            U.correctBoundaryConditions();
-            fvOptions.correct(U);
-            K = 0.5*magSqr(U);
-        }
+        // Calculate the conservative fluxes
+        phi = phiHbyA + p_rghEqn.flux();
+
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
+
+        // Correct the momentum source with the pressure gradient flux
+        // calculated from the relaxed pressure
+        U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
     }
+}
 
-    p = p_rgh + rho*gh;
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
+p = p_rgh + rho*gh;
 
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
 
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
 }
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"

applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H

@@ -52,6 +52,8 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
 mesh.setFluxRequired(p.name());
 
 

applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H

-{
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-    if (pimple.transonic())
-    {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            (
-                fvc::flux(HbyA)
-              + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-            )
-        );
-
-        MRF.makeRelative(phiHbyA);
-
-        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
-
-        phiHbyA *= fvc::interpolate(rho);
-
-        fvScalarMatrix pDDtEqn
-        (
-            fvc::ddt(rho) + fvc::div(phiHbyA)
-          + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
-        );
-
-        while (pimple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                pDDtEqn
-              - fvm::laplacian(rhorAUf, p)
-             ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            if (pimple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-    else
-    {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            (
-                fvc::flux(rho*HbyA)
-              + rhorAUf*fvc::ddtCorr(rho, U, phi)
-            )
-        );
-
-        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-        // Update the pressure BCs to ensure flux consistency
-        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
-        fvScalarMatrix pDDtEqn
-        (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phiHbyA)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
-
-        while (pimple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                pDDtEqn
-              - fvm::laplacian(rhorAUf, p)
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            if (pimple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-    K = 0.5*magSqr(U);
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-}

applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -100,7 +101,14 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-                #include "pEqn.H"
+                if (pimple.consistent())
+                {
+                    #include "../../../compressible/rhoPimpleFoam/pcEqn.H"
+                }
+                else
+                {
+                    #include "../../../compressible/rhoPimpleFoam/pEqn.H"
+                }
             }
 
             if (pimple.turbCorr())

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -42,6 +42,8 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
+mesh.setFluxRequired(p.name());
+
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -54,8 +56,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-mesh.setFluxRequired(p.name());
-
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -31,17 +38,17 @@ if (pimple.transonic())
 
     phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
@@ -56,16 +63,17 @@ if (pimple.transonic())
 }
 else
 {
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -76,6 +84,9 @@ else
     }
 }
 
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -92,10 +103,9 @@ if (pressureControl.limit(p))
     p.correctBoundaryConditions();
     rho = thermo.rho();
 }
-
-if (!pimple.transonic())
+else if (pimple.SIMPLErho())
 {
-    rho.relax();
+    rho = thermo.rho();
 }
 
 if (thermo.dpdt())

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@@ -39,17 +46,17 @@ if (pimple.transonic())
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAtU, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
@@ -67,16 +74,17 @@ else
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAtU, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -87,6 +95,9 @@ else
     }
 }
 
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -103,10 +114,9 @@ if (pressureControl.limit(p))