applications/solvers/combustion/reactingFoam/EEqn.H

@@ -17,8 +17,7 @@
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        rho*(U&g)
-      + reaction->Sh()
+        reaction->Sh()
       + fvOptions(rho, he)
     );
 

applications/solvers/combustion/reactingFoam/Make/options

 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/combustion/reactingFoam/UEqn.H

-    MRF.correctBoundaryVelocity(U);
+// Solve the Momentum equation
 
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + MRF.DDt(rho, U)
-      + turbulence->divDevRhoReff(U)
-     ==
-        rho*g
-      + fvOptions(rho, U)
-    );
+MRF.correctBoundaryVelocity(U);
 
-    UEqn.relax();
+tmp<fvVectorMatrix> UEqn
+(
+    fvm::ddt(rho, U) + fvm::div(phi, U)
+  + MRF.DDt(rho, U)
+  + turbulence->divDevRhoReff(U)
+ ==
+    fvOptions(rho, U)
+);
 
-    fvOptions.constrain(UEqn);
+UEqn().relax();
 
-    if (pimple.momentumPredictor())
-    {
-        solve(UEqn == -fvc::grad(p));
+fvOptions.constrain(UEqn());
 
-        fvOptions.correct(U);
-        K = 0.5*magSqr(U);
-    }
+if (pimple.momentumPredictor())
+{
+    solve(UEqn() == -fvc::grad(p));
+
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}

applications/solvers/combustion/reactingFoam/createFields.H

@@ -45,6 +45,28 @@ const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        pimple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        pimple.dict(),
+        dimDensity,
+        0
+    )
+);
+
 mesh.setFluxRequired(p.name());
 
 Info << "Creating turbulence model.\n" << nl;

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());
+volScalarField rAU(1.0/UEqn().A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
+HbyA = rAU*UEqn().H();
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
 
 if (pimple.transonic())
 {
@@ -26,7 +34,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -58,7 +66,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -75,6 +83,17 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
+// 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);

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()));
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
+
+if (pimple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+           /fvc::interpolate(rho)
+        )
+    );
+
+    MRF.makeRelative(fvc::interpolate(psi), phid);
+
+    surfaceScalarField phic
+    (
+        "phic",
+        fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
+    );
+
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          + fvc::div(phic)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi == phic + pEqn.flux();
+        }
+    }
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+        )
+    );
+
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+U = HbyA - rAtU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+
+if (!pimple.transonic())
+{
+    rho.relax();
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -24,6 +24,9 @@ License
 Application
     reactingFoam
 
+Group
+    grpCombustionSolvers
+
 Description
     Solver for combustion with chemical reactions.
 
@@ -34,7 +37,7 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
-#include "fvIOoptionList.H"
+#include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
@@ -51,11 +54,12 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRF.H"
     #include "createFvOptions.H"
 
+    turbulence->validate();
+
     if (!LTS)
     {
         #include "compressibleCourantNo.H"
@@ -95,7 +99,14 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-                #include "pEqn.H"
+                if (pimple.consistent())
+                {
+                    #include "pcEqn.H"
+                }
+                else
+                {
+                    #include "pEqn.H"
+                }
             }
 
             if (pimple.turbCorr())

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options

 EXE_INC = \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C

@@ -24,6 +24,9 @@ License
 Application
     rhoReactingBuoyantFoam
 
+Group
+    grpCombustionSolvers
+
 Description
     Solver for combustion with chemical reactions using density based
     thermodynamics package, using enahanced buoyancy treatment.
@@ -35,7 +38,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
-#include "fvIOoptionList.H"
+#include "fvOptions.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -57,6 +60,8 @@ int main(int argc, char *argv[])
     #include "createMRF.H"
     #include "createFvOptions.H"
 
+    turbulence->validate();
+
     if (!LTS)
     {
         #include "compressibleCourantNo.H"

applications/solvers/combustion/reactingFoam/rhoReactingFoam/Make/options

 EXE_INC = \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

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

@@ -5,11 +5,11 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    volScalarField rAU(1.0/UEqn.A());
+    volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn.H();
+    HbyA = rAU*UEqn().H();
 
     if (pimple.transonic())
     {

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

@@ -24,6 +24,9 @@ License
 Application
     rhoReactingFoam
 
+Group
+    grpCombustionSolvers
+
 Description
     Solver for combustion with chemical reactions using density based
     thermodynamics package.
@@ -35,7 +38,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
-#include "fvIOoptionList.H"
+#include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
@@ -52,11 +55,12 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRF.H"
     #include "createFvOptions.H"
 
+    turbulence->validate();
+
     if (!LTS)
     {
         #include "compressibleCourantNo.H"

applications/solvers/combustion/reactingFoam/setRDeltaT.H

@@ -116,6 +116,9 @@ License
         );
     }
 
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
     Info<< "    Overall     = "
         << gMin(1/rDeltaT.internalField())
         << ", " << gMax(1/rDeltaT.internalField()) << endl;

applications/solvers/compressible/compressibleSolvers.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it under
+    the terms of the GNU General Public License as published by the Free
+    Software Foundation, either version 3 of the License, or (at your option)
+    any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY
+    WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+    FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+    details.
+
+    You should have received a copy of the GNU General Public License along with
+    OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\defgroup grpCompressibleSolvers Compressible flow solvers
+@{
+    \ingroup grpSolvers
+    This group contains compressible flow solvers.
+@}
+
+\*---------------------------------------------------------------------------*/

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C

@@ -94,15 +94,8 @@ Foam::smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
      || mag(accommodationCoeff_) > 2.0
     )
     {
-        FatalIOErrorIn
+        FatalIOErrorInFunction
         (
-            "smoluchowskiJumpTFvPatchScalarField::"
-            "smoluchowskiJumpTFvPatchScalarField"
-            "("
-            "    const fvPatch&,"
-            "    const DimensionedField<scalar, volMesh>&,"
-            "    const dictionary&"
-            ")",
             dict
         )   << "unphysical accommodationCoeff specified"
             << "(0 < accommodationCoeff <= 1)" << endl

applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,14 +96,8 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
      || mag(accommodationCoeff_) > 2.0
     )
     {
-        FatalIOErrorIn
+        FatalIOErrorInFunction
         (
-            "maxwellSlipUFvPatchScalarField::maxwellSlipUFvPatchScalarField"
-            "("
-                "const fvPatch&, "
-                "const DimensionedField<vector, volMesh>&, "
-                "const dictionary&"
-            ")",
             dict
         )   << "unphysical accommodationCoeff_ specified"
             << "(0 < accommodationCoeff_ <= 1)" << endl

applications/solvers/compressible/rhoCentralFoam/Make/options

@@ -11,6 +11,7 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfvOptions \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \

applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H

@@ -7,10 +7,8 @@ if (mesh.schemesDict().readIfPresent("fluxScheme", fluxScheme))
     }
     else
     {
-        FatalErrorIn
-        (
-            "rhoCentralFoam::readFluxScheme"
-        )   << "fluxScheme: " << fluxScheme
+        FatalErrorInFunction
+            << "fluxScheme: " << fluxScheme
             << " is not a valid choice. "
             << "Options are: Tadmor, Kurganov"
             << abort(FatalError);

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options

@@ -13,6 +13,7 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfvOptions \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C

@@ -24,6 +24,9 @@ License
 Application
     rhoCentralDyMFoam
 
+Group
+    grpCompressibleSolvers grpMovingMeshSolvers
+
 Description
     Density-based compressible flow solver based on central-upwind schemes of
     Kurganov and Tadmor with support for mesh-motion and topology changes
@@ -49,6 +52,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createTimeControls.H"
 
+    turbulence->validate();
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     #include "readFluxScheme.H"