From 38177526ffd9719f34fcc0cec9339c301b1f461d Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Fri, 29 May 2015 11:30:40 +0100
Subject: [PATCH] simpleFoam/SRFSimpleFoam: Added support for SIMPLEC SIMPLEC
 (SIMPLE-consistent) is selected by setting "consistent" option true/yes:

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent yes;
}

which relaxes the pressure in a "consistent" manner and additional
relaxation of the pressure is not generally necessary.  In addition
convergence of the p-U system is better and reliable with less
aggressive relaxation of the momentum equation, e.g. for the motorbike
tutorial:

relaxationFactors
{
    equations
    {
        U               0.9;
        k               0.7;
        omega           0.7;
    }
}

The cost per iteration is marginally higher but the convergence rate is
better so the number of iterations can be reduced.

The SIMPLEC algorithm also provides benefit for cases with large
body-forces, e.g. SRF, see tutorials/incompressible/SRFSimpleFoam/mixer
and feature request http://www.openfoam.org/mantisbt/view.php?id=1714
---
 .../simpleFoam/SRFSimpleFoam/pEqn.H           | 17 +++++++++--
 .../solvers/incompressible/simpleFoam/pEqn.H  | 19 ++++++++----
 .../solutionControl/solutionControl.C         | 18 +++++++----
 .../solutionControl/solutionControl.H         | 10 ++++++-
 .../solutionControl/solutionControlI.H        |  8 ++++-
 .../SRFSimpleFoam/mixer/system/fvSolution     | 12 ++------
 .../SRFSimpleFoam/mixer/turbulenceProperties  | 30 -------------------
 .../simpleFoam/motorBike/system/fvSolution    |  7 ++---
 8 files changed, 60 insertions(+), 61 deletions(-)
 delete mode 100644 tutorials/incompressible/SRFSimpleFoam/mixer/turbulenceProperties

diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
index b9fa52c0ff..50691c0154 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
@@ -2,17 +2,28 @@
     volScalarField rAUrel(1.0/UrelEqn().A());
     volVectorField HbyA("HbyA", Urel);
     HbyA = rAUrel*UrelEqn().H();
-    UrelEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
     adjustPhi(phiHbyA, Urel, p);
 
+    tmp<volScalarField> rAtUrel(rAUrel);
+
+    if (simple.consistent())
+    {
+        rAtUrel = 1.0/(1.0/rAUrel - UrelEqn().H1());
+        phiHbyA +=
+            fvc::interpolate(rAtUrel() - rAUrel)*fvc::snGrad(p)*mesh.magSf();
+        HbyA -= (rAUrel - rAtUrel())*fvc::grad(p);
+    }
+
+    UrelEqn.clear();
+
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAUrel, p) == fvc::div(phiHbyA)
+            fvm::laplacian(rAtUrel(), p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -31,7 +42,7 @@
     p.relax();
 
     // Momentum corrector
-    Urel = HbyA - rAUrel*fvc::grad(p);
+    Urel = HbyA - rAtUrel()*fvc::grad(p);
     Urel.correctBoundaryConditions();
     fvOptions.correct(Urel);
 }
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 4d669c8db2..fddb3b05d4 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -2,20 +2,29 @@
     volScalarField rAU(1.0/UEqn().A());
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
-    UEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
-
     fvOptions.makeRelative(phiHbyA);
-
     adjustPhi(phiHbyA, U, p);
 
+    tmp<volScalarField> rAtU(rAU);
+
+    if (simple.consistent())
+    {
+        rAtU = 1.0/(1.0/rAU - UEqn().H1());
+        phiHbyA +=
+            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
+        HbyA -= (rAU - rAtU())*fvc::grad(p);
+    }
+
+    UEqn.clear();
+
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -34,7 +43,7 @@
     p.relax();
 
     // Momentum corrector
-    U = HbyA - rAU*fvc::grad(p);
+    U = HbyA - rAtU()*fvc::grad(p);
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 }
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
index 8448153a69..8e3dd89075 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,16 +37,21 @@ namespace Foam
 
 void Foam::solutionControl::read(const bool absTolOnly)
 {
-    const dictionary& solnDict = this->dict();
+    const dictionary& solutionDict = this->dict();
 
     // Read solution controls
     nNonOrthCorr_ =
-        solnDict.lookupOrDefault<label>("nNonOrthogonalCorrectors", 0);
-    momentumPredictor_ = solnDict.lookupOrDefault("momentumPredictor", true);
-    transonic_ = solnDict.lookupOrDefault("transonic", false);
+        solutionDict.lookupOrDefault<label>("nNonOrthogonalCorrectors", 0);
+    momentumPredictor_ =
+        solutionDict.lookupOrDefault("momentumPredictor", true);
+    transonic_ = solutionDict.lookupOrDefault("transonic", false);
+    consistent_ = solutionDict.lookupOrDefault("consistent", false);
 
     // Read residual information
-    const dictionary residualDict(solnDict.subOrEmptyDict("residualControl"));
+    const dictionary residualDict
+    (
+        solutionDict.subOrEmptyDict("residualControl")
+    );
 
     DynamicList<fieldData> data(residualControl_);
 
@@ -177,6 +182,7 @@ Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
     nNonOrthCorr_(0),
     momentumPredictor_(true),
     transonic_(false),
+    consistent_(false),
     corr_(0),
     corrNonOrtho_(0)
 {}
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
index e701813d47..9cd6cbbf16 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
@@ -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
@@ -83,6 +83,10 @@ protected:
             //- Flag to indicate to solve using transonic algorithm
             bool transonic_;
 
+            //- Flag to indicate to relax pressure using the
+            //  "consistent" approach of SIMPLEC
+            bool consistent_;
+
 
         // Evolution
 
@@ -169,6 +173,10 @@ public:
             //- Flag to indicate to solve using transonic algorithm
             inline bool transonic() const;
 
+            //- Flag to indicate to relax pressure using the
+            //  "consistent" approach of SIMPLEC
+            inline bool consistent() const;
+
 
         // Evolution
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
index b9d975cfba..61a647e408 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,6 +67,12 @@ inline bool Foam::solutionControl::transonic() const
 }
 
 
+inline bool Foam::solutionControl::consistent() const
+{
+    return consistent_;
+}
+
+
 inline bool Foam::solutionControl::correctNonOrthogonal()
 {
     corrNonOrtho_++;
diff --git a/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSolution b/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSolution
index 99306fab01..e092956f32 100644
--- a/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSolution
+++ b/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSolution
@@ -37,22 +37,14 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
+    consistent          yes;
 }
 
 relaxationFactors
 {
-    fields
-    {
-        p               0.3;
-    }
     equations
     {
-        Urel            0.7;
-        k               0.7;
-        epsilon         0.7;
-        omega           0.7;
-        R               0.7;
-        nuTilda         0.7;
+        ".*"            0.9;
     }
 }
 
diff --git a/tutorials/incompressible/SRFSimpleFoam/mixer/turbulenceProperties b/tutorials/incompressible/SRFSimpleFoam/mixer/turbulenceProperties
deleted file mode 100644
index 4406361c66..0000000000
--- a/tutorials/incompressible/SRFSimpleFoam/mixer/turbulenceProperties
+++ /dev/null
@@ -1,30 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      turbulenceProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-simulationType RAS;
-
-RAS
-{
-    RASModel        kOmegaSST;
-
-    turbulence      on;
-
-    printCoeffs     on;
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
index 4f0f515aca..6f9e2f67b8 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
@@ -66,6 +66,7 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
+    consistent yes;
 }
 
 potentialFlow
@@ -75,13 +76,9 @@ potentialFlow
 
 relaxationFactors
 {
-    fields
-    {
-        p               0.3;
-    }
     equations
     {
-        U               0.7;
+        U               0.9;
         k               0.7;
         omega           0.7;
     }