chtMultiRegionFoam: SIMPLE operation and transonic switch

Multi-region PIMPLE controls have been applied to the chtMultiRegionFoam
solver, and a transonic option has been implemented.

The new PIMPLE controls let the solver operate SIMPLE mode. The
utilisation of library solution and convergence control functionality
has significantly reduced the amount of code in the solver. The
chtMultiRegionSimpleFoam solver has also been made obsolete, and has
therefore been removed.

A few changes will be necessary to convert an existing
chtMultiRegionSimpleFoam case to chtMultiRegionFoam. All the SIMPLE
sub-dictionaries in the system/<regions>/fvSolution will need to be
renamed PIMPLE. The system/fvSolution file will also need an empty
PIMPLE sub-dictionary. In addition, additional "<variable>Final" solver
and relaxation entries will be needed. For a steady case, adding a
wildcard ending, ".*", to the variable names should be sufficient.

Solution parameters appropriate for a steady case are shown below:

    solvers
    {
        "p_rgh.*"
        {
            solver           GAMG;
            tolerance        1e-7;
            relTol           0.01;
            smoother         DIC;
            maxIter          10;
        }

        "(U|h|e|k|epsilon).*"
        {
            solver           PBiCGStab;
            preconditioner   DILU;
            tolerance        1e-7;
            relTol           0.1;
        }
    }

    PIMPLE
    {
        // ...
    }

    relaxationFactors
    {
        fields
        {
            "p_rgh.*"       0.7;
        }
        equations
        {
            "U.*"           0.5;
            "(h|e).*"       0.3;
            "(k|epsilon).*" 0.2;
        }
    }

This work was supported by Fabian Buelow, at Evonik
Tobias Holzmann provided cases for testing the convergence controls

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -8,6 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \

applications/solvers/heatTransfer/chtMultiRegionFoam/checkResidualControls.H

@@ -1,66 +0,0 @@
-// Residual control used?
-bool resControlUsed = false;
-int nFluidControlled = fluidRegions.size();
-int nSolidControlled = solidRegions.size();
-
-// Check wheater there is a single regions that uses residual control
-forAll(fluidRegions, i)
-{
-    if (residualControlUsedFluid[i])
-    {
-        resControlUsed = true;
-    }
-}
-
-forAll(solidRegions, i)
-{
-    if(residualControlUsedSolid[i])
-    {
-        resControlUsed = true;
-    }
-}
-
-if (resControlUsed)
-{
-    int nFluidConv = 0;
-    int nSolidConv = 0;
-
-    // Sum of all converged regions (Note: if no residual control is used
-    // the residualReached* flag is already set to true)
-    forAll(fluidRegions, i)
-    {
-        if (residualReachedFluid[i])
-        {
-            nFluidConv++;
-        }
-    }
-
-    forAll(solidRegions, i)
-    {
-        if (residualReachedSolid[i])
-        {
-            nSolidConv++;
-        }
-    }
-
-    if (nFluidConv == nFluidControlled && nSolidConv == nSolidControlled)
-    {
-        // Activate flag to go to the 'Final' loop using the 'Final'
-        // relaxation factors
-        allRegionsConverged = true;
-    }
-}
-
-if (finalIter && resControlUsed && !allRegionsConverged)
-{
-    Info<< "\nRegions not converged after " << nOuterCorr
-        << " outer correctors" << endl;
-}
-else if (finalIter && resControlUsed && allRegionsConverged)
-{
-    Info<< "\nRegions converged after " << oCorr
-        << " outer correctors" << endl;
-
-    // Leave PIMPLE loop
-    break;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,8 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "coordinateSystem.H"
+#include "pimpleMultiRegionControl.H"
+#include "pressureControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,17 +61,18 @@ int main(int argc, char *argv[])
     #include "createMeshes.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
+    pimpleMultiRegionControl pimples(fluidRegions, solidRegions);
+    #include "createFluidPressureControls.H"
     #include "createTimeControls.H"
     #include "readSolidTimeControls.H"
     #include "compressibleMultiRegionCourantNo.H"
     #include "solidRegionDiffusionNo.H"
     #include "setInitialMultiRegionDeltaT.H"
 
-    while (runTime.run())
+    while (pimples.run(runTime))
     {
         #include "readTimeControls.H"
         #include "readSolidTimeControls.H"
-        #include "readPIMPLEControls.H"
 
         #include "compressibleMultiRegionCourantNo.H"
         #include "solidRegionDiffusionNo.H"
@@ -79,36 +82,15 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        if (nOuterCorr != 1)
-        {
-            forAll(fluidRegions, i)
-            {
-                #include "storeOldFluidFields.H"
-            }
-        }
-
-        bool allRegionsConverged = false;
-        bool finalIter = false;
-
         // --- PIMPLE loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        while (pimples.loop())
         {
-            Info<< nl << "Pimple iteration " << oCorr;
-
-            if (oCorr == nOuterCorr-1 || allRegionsConverged)
-            {
-                finalIter = true;
-            }
-
             forAll(fluidRegions, i)
             {
                 Info<< "\nSolving for fluid region "
                     << fluidRegions[i].name() << endl;
                 #include "setRegionFluidFields.H"
-                #include "readFluidMultiRegionPIMPLEControls.H"
-                #include "readFluidMultiRegionResidualControls.H"
                 #include "solveFluid.H"
-                #include "residualControlsFluid.H"
             }
 
             forAll(solidRegions, i)
@@ -116,13 +98,8 @@ int main(int argc, char *argv[])
                 Info<< "\nSolving for solid region "
                     << solidRegions[i].name() << endl;
                 #include "setRegionSolidFields.H"
-                #include "readSolidMultiRegionPIMPLEControls.H"
-                #include "readSolidMultiRegionResidualControls.H"
                 #include "solveSolid.H"
-                #include "residualControlsSolid.H"
             }
-
-            #include "checkResidualControls.H"
         }
 
         runTime.write();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/files

@@ -1,4 +0,0 @@
-chtMultiRegionSimpleFoam.C
-
-EXE = $(FOAM_APPBIN)/chtMultiRegionSimpleFoam
-

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options

@@ -1,30 +0,0 @@
-EXE_INC = \
-    -I. \
-    -I.. \
-    -Ifluid \
-    -Isolid \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/regionModels/regionModel/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lsolidThermo \
-    -lspecie \
-    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lradiationModels \
-    -lfvOptions \
-    -lregionModels \
-    -lsampling

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/checkResidualControls.H

@@ -1,68 +0,0 @@
-// Residual control used?
-bool resControlUsed = false;
-int nFluidControlled = fluidRegions.size();
-int nSolidControlled = solidRegions.size();
-
-// Check wheater there is a single regions that uses residual control
-forAll(fluidRegions, i)
-{
-    if (residualControlUsedFluid[i])
-    {
-        resControlUsed = true;
-        break;
-    }
-}
-
-forAll(solidRegions, i)
-{
-    if(residualControlUsedSolid[i])
-    {
-        resControlUsed = true;
-        break;
-    }
-}
-
-if (resControlUsed)
-{
-    int nFluidConv = 0;
-    int nSolidConv = 0;
-
-    // Sum of all converged regions (Note: if no residual control is used
-    // the residualReached* flag is already set to true)
-    forAll(fluidRegions, i)
-    {
-        if (residualReachedFluid[i])
-        {
-            nFluidConv++;
-        }
-    }
-
-    forAll(solidRegions, i)
-    {
-        if (residualReachedSolid[i])
-        {
-            nSolidConv++;
-        }
-    }
-
-    if (nFluidConv == nFluidControlled && nSolidConv == nSolidControlled)
-    {
-        // Activate flag to go to the 'Final' loop using the 'Final'
-        // relaxation factors
-        allRegionsConverged = true;
-    }
-}
-
-if (!runTime.running())
-{
-    Info<< "\nRegions not converged after " << runTime.timeName()
-        << " iterations" << endl;
-}
-else if (runTime.running() && resControlUsed && allRegionsConverged)
-{
-    Info<< "\nRegions converged after " << runTime.timeName()
-        << " iterations" << endl;
-
-    // Save converged time step and end the run
-    runTime.writeAndEnd();
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -1,103 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 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/>.
-
-Application
-    chtMultiRegionSimpleFoam
-
-Description
-    Steady-state solver for buoyant, turbulent fluid flow and solid heat
-    conduction with conjugate heat transfer between solid and fluid regions.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "rhoThermo.H"
-#include "turbulentFluidThermoModel.H"
-#include "fixedGradientFvPatchFields.H"
-#include "regionProperties.H"
-#include "solidThermo.H"
-#include "radiationModel.H"
-#include "fvOptions.H"
-#include "coordinateSystem.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #define NO_CONTROL
-    #define CREATE_MESH createMeshesPostProcess.H
-    #include "postProcess.H"
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createMeshes.H"
-    #include "createFields.H"
-    #include "initContinuityErrs.H"
-
-    // Residual control
-    bool allRegionsConverged = false;
-
-    while (runTime.run())
-    {
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        forAll(fluidRegions, i)
-        {
-            Info<< "\nSolving for fluid region "
-                << fluidRegions[i].name() << endl;
-            #include "setRegionFluidFields.H"
-            #include "readFluidMultiRegionSIMPLEControls.H"
-            #include "readFluidMultiRegionResidualControls.H"
-            #include "solveFluid.H"
-            #include "residualControlsFluid.H"
-        }
-
-        forAll(solidRegions, i)
-        {
-            Info<< "\nSolving for solid region "
-                << solidRegions[i].name() << endl;
-            #include "setRegionSolidFields.H"
-            #include "readSolidMultiRegionSIMPLEControls.H"
-            #include "readSolidMultiRegionResidualControls.H"
-            #include "solveSolid.H"
-            #include "residualControlsSolid.H"
-        }
-
-        #include "checkResidualControls.H"
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H

@@ -1,32 +0,0 @@
-{
-    volScalarField& he = thermo.he();
-
-    fvScalarMatrix EEqn
-    (
-        fvm::div(phi, he)
-      + (
-            he.name() == "e"
-          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
-          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
-        )
-      - fvm::laplacian(turb.alphaEff(), he)
-     ==
-        rho*(U&g)
-      + rad.Sh(thermo, he)
-      + fvOptions(rho, he)
-    );
-
-    EEqn.relax();
-
-    fvOptions.constrain(EEqn);
-
-    solvPerfh = EEqn.solve();
-
-    fvOptions.correct(he);
-
-    thermo.correct();
-    rad.correct();
-
-    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-        << max(thermo.T()).value() << endl;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H

@@ -1,35 +0,0 @@
-    // Solve the Momentum equation
-
-    MRF.correctBoundaryVelocity(U);
-
-    tmp<fvVectorMatrix> tUEqn
-    (
-        fvm::div(phi, U)
-      + MRF.DDt(rho, U)
-      + turb.divDevRhoReff(U)
-     ==
-        fvOptions(rho, U)
-    );
-    fvVectorMatrix& UEqn = tUEqn.ref();
-
-    UEqn.relax();
-
-    fvOptions.constrain(UEqn);
-
-    if (momentumPredictor)
-    {
-        solvPerfU = solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-                (
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                )*mesh.magSf()
-            )
-        );
-    }
-
-    fvOptions.correct(U);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/compressibleContinuityErrors.H

@@ -1,21 +0,0 @@
-{
-    dimensionedScalar totalMass = fvc::domainIntegrate(rho);
-
-    scalar sumLocalContErr =
-    (
-        fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass
-    ).value();
-
-    scalar globalContErr =
-    (
-        fvc::domainIntegrate(rho - thermo.rho())/totalMass
-    ).value();
-
-    cumulativeContErr[i] += globalContErr;
-
-    Info<< "time step continuity errors (" << mesh.name() << ")"
-        << ": sum local = " << sumLocalContErr
-        << ", global = " << globalContErr
-        << ", cumulative = " << cumulativeContErr[i]
-        << endl;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -1,15 +0,0 @@
-    scalar CoNum = -great;
-    forAll(fluidRegions, regionI)
-    {
-        CoNum = max
-        (
-            compressibleCourantNo
-            (
-                fluidRegions[regionI],
-                runTime,
-                rhoFluid[regionI],
-                phiFluid[regionI]
-            ),
-            CoNum
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H

@@ -1,277 +0,0 @@
-// Initialise fluid field pointer lists
-PtrList<rhoThermo> thermoFluid(fluidRegions.size());
-PtrList<volScalarField> rhoFluid(fluidRegions.size());
-PtrList<volVectorField> UFluid(fluidRegions.size());
-PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
-PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
-PtrList<uniformDimensionedScalarField> hRefFluid(fluidRegions.size());
-PtrList<volScalarField> ghFluid(fluidRegions.size());
-PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
-PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
-PtrList<volScalarField> p_rghFluid(fluidRegions.size());
-PtrList<radiation::radiationModel> radiation(fluidRegions.size());
-
-List<scalar> initialMassFluid(fluidRegions.size());
-List<label> pRefCellFluid(fluidRegions.size(), 0);
-List<scalar> pRefValueFluid(fluidRegions.size(), 0.0);
-List<bool> frozenFlowFluid(fluidRegions.size(), false);
-List<bool> residualReachedFluid(fluidRegions.size(), true);
-List<bool> residualControlUsedFluid(fluidRegions.size(), false);
-List<bool> firstIterationFluid(fluidRegions.size(), true);
-
-PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
-PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
-
-PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
-PtrList<fv::options> fluidFvOptions(fluidRegions.size());
-
-// Populate fluid field pointer lists
-forAll(fluidRegions, i)
-{
-    Info<< "*** Reading fluid mesh thermophysical properties for region "
-        << fluidRegions[i].name() << nl << endl;
-
-    Info<< "    Adding to thermoFluid\n" << endl;
-
-    thermoFluid.set
-    (
-        i,
-        rhoThermo::New(fluidRegions[i]).ptr()
-    );
-
-    Info<< "    Adding to rhoFluid\n" << endl;
-    rhoFluid.set
-    (
-        i,
-        new volScalarField
-        (
-            IOobject
-            (
-                "rho",
-                runTime.timeName(),
-                fluidRegions[i],
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            thermoFluid[i].rho()
-        )
-    );
-
-    Info<< "    Adding to UFluid\n" << endl;
-    UFluid.set
-    (
-        i,
-        new volVectorField
-        (
-            IOobject
-            (
-                "U",
-                runTime.timeName(),
-                fluidRegions[i],
-                IOobject::MUST_READ,
-                IOobject::AUTO_WRITE
-            ),
-            fluidRegions[i]
-        )
-    );
-
-    Info<< "    Adding to phiFluid\n" << endl;
-    phiFluid.set
-    (
-        i,
-        new surfaceScalarField
-        (
-            IOobject
-            (
-                "phi",
-                runTime.timeName(),
-                fluidRegions[i],
-                IOobject::READ_IF_PRESENT,
-                IOobject::AUTO_WRITE
-            ),
-            linearInterpolate(rhoFluid[i]*UFluid[i])
-                & fluidRegions[i].Sf()
-        )
-    );
-
-    Info<< "    Adding to gFluid\n" << endl;
-    gFluid.set
-    (
-        i,
-        new uniformDimensionedVectorField
-        (
-            IOobject
-            (
-                "g",
-                runTime.constant(),
-                fluidRegions[i],
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE
-            )
-        )
-    );
-
-    Info<< "    Adding to hRefFluid\n" << endl;
-    hRefFluid.set
-    (
-        i,
-        new uniformDimensionedScalarField
-        (
-            IOobject
-            (
-                "hRef",
-                runTime.constant(),
-                fluidRegions[i],
-                IOobject::READ_IF_PRESENT,
-                IOobject::NO_WRITE
-            ),
-            dimensionedScalar("hRef", dimLength, 0)
-        )
-    );
-
-    dimensionedScalar ghRef
-    (
-        mag(gFluid[i].value()) > small
-      ? gFluid[i]
-          & (cmptMag(gFluid[i].value())/mag(gFluid[i].value()))*hRefFluid[i]
-      : dimensionedScalar("ghRef", gFluid[i].dimensions()*dimLength, 0)
-    );
-
-    Info<< "    Adding to ghFluid\n" << endl;
-    ghFluid.set
-    (
-        i,
-        new volScalarField
-        (
-            "gh",
-            (gFluid[i] & fluidRegions[i].C()) - ghRef
-        )
-    );
-
-    Info<< "    Adding to ghfFluid\n" << endl;
-    ghfFluid.set
-    (
-        i,
-        new surfaceScalarField
-        (
-            "ghf",
-            (gFluid[i] & fluidRegions[i].Cf()) - ghRef
-        )
-    );
-
-    Info<< "    Adding to turbulence\n" << endl;
-    turbulence.set
-    (
-        i,
-        compressible::turbulenceModel::New
-        (
-            rhoFluid[i],
-            UFluid[i],
-            phiFluid[i],
-            thermoFluid[i]
-        ).ptr()
-    );
-
-    p_rghFluid.set
-    (
-        i,
-        new volScalarField
-        (
-            IOobject
-            (
-                "p_rgh",
-                runTime.timeName(),
-                fluidRegions[i],
-                IOobject::MUST_READ,
-                IOobject::AUTO_WRITE
-            ),
-            fluidRegions[i]
-        )
-    );
-
-    // Force p_rgh to be consistent with p
-    p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
-
-    fluidRegions[i].setFluxRequired(p_rghFluid[i].name());
-
-    radiation.set
-    (
-        i,
-        radiation::radiationModel::New(thermoFluid[i].T())
-    );
-
-    initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
-
-    const dictionary& simpleDict =
-        fluidRegions[i].solutionDict().subDict("SIMPLE");
-
-    setRefCell
-    (
-        thermoFluid[i].p(),
-        p_rghFluid[i],
-        simpleDict,
-        pRefCellFluid[i],
-        pRefValueFluid[i]
-    );
-
-    simpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
-
-    rhoMax.set
-    (
-        i,
-        new dimensionedScalar
-        (
-            dimensionedScalar::lookupOrDefault
-            (
-                "rhoMax",
-                simpleDict,
-                dimDensity,
-                great
-            )
-        )
-    );
-
-    rhoMin.set
-    (
-        i,
-        new dimensionedScalar
-        (
-            dimensionedScalar::lookupOrDefault
-            (
-                "rhoMin",
-                simpleDict,
-                dimDensity,
-                0
-            )
-        )
-    );
-
-    Info<< "    Adding MRF\n" << endl;
-    MRFfluid.set
-    (
-        i,
-        new IOMRFZoneList(fluidRegions[i])
-    );
-
-    Info<< "    Adding fvOptions\n" << endl;
-    fluidFvOptions.set
-    (
-        i,
-        new fv::options(fluidRegions[i])
-    );
-
-    turbulence[i].validate();
-
-    if (simpleDict.isDict("residualControl"))
-    {
-
-        Info<< "    Reading residual controls\n" << endl;
-
-        bool& residualControlUsed = residualControlUsedFluid[i];
-        bool& resReachedFluid = residualReachedFluid[i];
-        bool& firstIteration = firstIterationFluid[i];
-        const fvMesh& mesh = fluidRegions[i];
-
-        #include "readFluidMultiRegionResidualControls.H"
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidMeshes.H

@@ -1,24 +0,0 @@
-    const wordList fluidNames(rp["fluid"]);
-
-    PtrList<fvMesh> fluidRegions(fluidNames.size());
-
-    forAll(fluidNames, i)
-    {
-        Info<< "Create fluid mesh for region " << fluidNames[i]
-            << " for time = " << runTime.timeName() << nl << endl;
-
-        fluidRegions.set
-        (
-            i,
-            new fvMesh
-            (
-                IOobject
-                (
-                    fluidNames[i],
-                    runTime.timeName(),
-                    runTime,
-                    IOobject::MUST_READ
-                )
-            )
-        );
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -1,84 +0,0 @@
-{
-    rho = thermo.rho();
-    rho = max(rho, rhoMin[i]);
-    rho = min(rho, rhoMax[i]);
-    rho.relax();
-
-    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    tUEqn.clear();
-
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::flux(rho*HbyA)
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
-
-    phiHbyA += phig;
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > small);
-
-    // Solve pressure
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        fvScalarMatrix p_rghEqn
-        (
-            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
-        );
-
-        p_rghEqn.setReference
-        (
-            pRefCell,
-            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
-        );
-
-        solvPerfp_rgh = p_rghEqn.solve();
-
-        if (nonOrth == nNonOrthCorr)
-        {
-            // 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);
-        }
-    }
-
-    p = p_rgh + rho*gh;
-
-    #include "continuityErrs.H"
-
-    // For closed-volume cases adjust the pressure level
-    // to obey overall mass continuity
-    if (closedVolume && compressible)
-    {
-        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-            /compressibility;
-        p_rgh = p - rho*gh;
-    }
-
-    rho = thermo.rho();
-    rho = max(rho, rhoMin[i]);
-    rho = min(rho, rhoMax[i]);
-    rho.relax();
-
-    Info<< "Min/max rho:" << min(rho).value() << ' '
-        << max(rho).value() << endl;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionResidualControls.H

@@ -1,60 +0,0 @@
-const dictionary& residualControl =
-    mesh.solutionDict().subDict("SIMPLE").subOrEmptyDict("residualControl");
-
-scalar UTol = -1;
-scalar hTol = -1;
-scalar p_rghTol = -1;
-
-if (!residualControl.empty())
-{
-    if (!residualControl.subOrEmptyDict("U").empty())
-    {
-        UTol = readScalar(residualControl.subDict("U").lookup("tolerance"));
-    }
-
-    if (!residualControl.subOrEmptyDict("p_rgh").empty())
-    {
-        p_rghTol =
-            readScalar
-            (
-                residualControl.subDict("p_rgh").lookup("tolerance")
-            );
-    }
-
-    if (!residualControl.subOrEmptyDict("h").empty())
-    {
-        hTol = readScalar(residualControl.subDict("h").lookup("tolerance"));
-    }
-
-    //- Residual control used?
-    if (UTol != -1 || hTol != -1 || p_rghTol != -1)
-    {
-        residualControlUsed = true;
-        resReachedFluid = false;
-    }
-
-    if (firstIteration)
-    {
-        firstIteration = false;
-
-        if (UTol != -1 || hTol != -1 || p_rghTol != -1)
-        {
-            Info<< "SIMPLE convergence criteria\n";
-
-            if (UTol != -1)
-            {
-                Info<< "   field U      tolerance " << UTol << "\n";
-            }
-            if (hTol != -1)
-            {
-                Info<< "   field h      tolerance " << hTol << "\n";
-            }
-            if (p_rghTol != -1)
-            {
-                Info<< "   field p_rgh  tolerance " << p_rghTol << "\n";
-            }
-
-            Info<< endl;
-        }
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H

@@ -1,7 +0,0 @@
-    const dictionary& simple = fluidRegions[i].solutionDict().subDict("SIMPLE");
-
-    const int nNonOrthCorr =
-        simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    const bool momentumPredictor =
-        simple.lookupOrDefault("momentumPredictor", true);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/residualControlsFluid.H

@@ -1,56 +0,0 @@
-// Residual control used
-if (residualControlUsed)
-{
-    bool UConv = false;
-    bool p_rghConv = false;
-    bool hConv = false;
-
-    // Check which field is not used for control
-    {
-        if (UTol == -1 || !momentumPredictor)
-        {
-            UConv =  true;
-        }
-
-        if (p_rghTol == -1)
-        {
-            p_rghConv = true;
-        }
-
-        if (hTol == -1)
-        {
-            hConv = true;
-        }
-    }
-
-    // Get the last initial residual of the solvers
-    if (momentumPredictor && !UConv)
-    {
-        if (UTol > cmptMax(solvPerfU.initialResidual()))
-        {
-            UConv = true;
-        }
-    }
-
-    if (!p_rghConv)
-    {
-        if (p_rghTol > solvPerfp_rgh.initialResidual())
-        {
-            p_rghConv = true;
-        }
-    }
-
-    if (!hConv)
-    {
-        if (hTol > solvPerfh.initialResidual())
-        {
-            hConv = true;
-        }
-    }
-
-    // Check if each field is converged
-    if (UConv && p_rghConv && hConv)
-    {
-        resReachedFluid = true;
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -1,39 +0,0 @@
-    const fvMesh& mesh = fluidRegions[i];
-
-    rhoThermo& thermo = thermoFluid[i];
-    thermo.validate(args.executable(), "h", "e");
-
-    volScalarField& rho = rhoFluid[i];
-    volVectorField& U = UFluid[i];
-    surfaceScalarField& phi = phiFluid[i];
-
-    compressible::turbulenceModel& turb = turbulence[i];
-
-    volScalarField& p = thermo.p();
-    const volScalarField& psi = thermo.psi();
-
-    IOMRFZoneList& MRF = MRFfluid[i];
-    fv::options& fvOptions = fluidFvOptions[i];
-
-    const dimensionedScalar initialMass
-    (
-        "initialMass",
-        dimMass,
-        initialMassFluid[i]
-    );
-
-    radiation::radiationModel& rad = radiation[i];
-
-    const label pRefCell = pRefCellFluid[i];
-    const scalar pRefValue = pRefValueFluid[i];
-    const bool frozenFlow = frozenFlowFluid[i];
-
-    volScalarField& p_rgh = p_rghFluid[i];
-
-    const dimensionedVector& g = gFluid[i];
-    const volScalarField& gh = ghFluid[i];
-    const surfaceScalarField& ghf = ghfFluid[i];
-
-    bool& resReachedFluid = residualReachedFluid[i];
-    bool& residualControlUsed = residualControlUsedFluid[i];
-    bool& firstIteration = firstIterationFluid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H

@@ -1,19 +0,0 @@
-//  Pressure-velocity SIMPLE corrector
-
-    {
-        if (frozenFlow)
-        {
-            #include "EEqn.H"
-        }
-        else
-        {
-            p_rgh.storePrevIter();
-            rho.storePrevIter();
-
-            #include "UEqn.H"
-            #include "EEqn.H"
-            #include "pEqn.H"
-
-            turb.correct();
-        }
-    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/createSolidFields.H

@@ -1,126 +0,0 @@
-// Initialise solid field pointer lists
-PtrList<coordinateSystem> coordinates(solidRegions.size());
-PtrList<solidThermo> thermos(solidRegions.size());
-PtrList<radiation::radiationModel> radiations(solidRegions.size());
-PtrList<fv::options> solidHeatSources(solidRegions.size());
-PtrList<volScalarField> betavSolid(solidRegions.size());
-PtrList<volSymmTensorField> aniAlphas(solidRegions.size());
-
-List<bool> residualReachedSolid(solidRegions.size(), true);
-List<bool> residualControlUsedSolid(solidRegions.size(), false);
-List<bool> firstIterationSolid(solidRegions.size(), true);
-
-// Populate solid field pointer lists
-forAll(solidRegions, i)
-{
-    Info<< "*** Reading solid mesh thermophysical properties for region "
-        << solidRegions[i].name() << nl << endl;
-
-    Info<< "    Adding to thermos\n" << endl;
-    thermos.set(i, solidThermo::New(solidRegions[i]));
-
-    Info<< "    Adding to radiations\n" << endl;
-    radiations.set(i, radiation::radiationModel::New(thermos[i].T()));
-
-    Info<< "    Adding fvOptions\n" << endl;
-    solidHeatSources.set
-    (
-        i,
-        new fv::options(solidRegions[i])
-    );
-
-    if (!thermos[i].isotropic())
-    {
-        Info<< "    Adding coordinateSystems\n" << endl;
-        coordinates.set
-        (
-            i,
-            coordinateSystem::New(solidRegions[i], thermos[i])
-        );
-
-        tmp<volVectorField> tkappaByCp =
-            thermos[i].Kappa()/thermos[i].Cp();
-
-        aniAlphas.set
-        (
-            i,
-            new volSymmTensorField
-            (
-                IOobject
-                (
-                    "Anialpha",
-                    runTime.timeName(),
-                    solidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                solidRegions[i],
-                dimensionedSymmTensor
-                (
-                    "zero",
-                    tkappaByCp().dimensions(),
-                    Zero
-                ),
-                zeroGradientFvPatchSymmTensorField::typeName
-            )
-        );
-
-        aniAlphas[i].primitiveFieldRef() =
-            coordinates[i].R().transformVector(tkappaByCp());
-        aniAlphas[i].correctBoundaryConditions();
-
-    }
-
-    IOobject betavSolidIO
-    (
-        "betavSolid",
-        runTime.timeName(),
-        solidRegions[i],
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    );
-
-    if (betavSolidIO.typeHeaderOk<volScalarField>(true))
-    {
-        betavSolid.set
-        (
-            i,
-            new volScalarField(betavSolidIO, solidRegions[i])
-        );
-    }
-    else
-    {
-        betavSolid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "betavSolid",
-                    runTime.timeName(),
-                    solidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                solidRegions[i],
-                dimensionedScalar("1", dimless, scalar(1))
-            )
-        );
-    }
-
-    const dictionary& simpleDict =
-        solidRegions[i].solutionDict().subDict("SIMPLE");
-
-    if (simpleDict.isDict("residualControl"))
-    {
-        Info<< "    Reading residual controls\n" << endl;
-
-        bool& residualControlUsed = residualControlUsedSolid[i];
-        bool& resReachedSolid = residualReachedSolid[i];
-        bool& firstIteration = firstIterationSolid[i];
-        const fvMesh& mesh = solidRegions[i];
-
-        #include "readSolidMultiRegionResidualControls.H"
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/createSolidMeshes.H

@@ -1 +0,0 @@
-../../solid/createSolidMeshes.H

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/readSolidMultiRegionResidualControls.H

@@ -1,36 +0,0 @@
-const dictionary& residualControl =
-    mesh.solutionDict().subDict("SIMPLE").subOrEmptyDict("residualControl");
-
-scalar hTol = -1;
-
-if (!residualControl.empty())
-{
-    if (!residualControl.subOrEmptyDict("h").empty())
-    {
-        hTol = readScalar(residualControl.subDict("h").lookup("tolerance"));
-
-        //- Used residual control for actual solid region
-        if (hTol != -1)
-        {
-            residualControlUsed = true;
-            resReachedSolid = false;
-        }
-    }
-
-    if (firstIteration)
-    {
-        firstIteration = false;
-
-        if (hTol != -1)
-        {
-            Info<< "SIMPLE convergence criteria\n";
-
-            if (hTol != -1)
-            {
-                Info<< "   field h      tolerance " << hTol << "\n";
-            }
-
-            Info<< endl;
-        }
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/readSolidMultiRegionSIMPLEControls.H

@@ -1,4 +0,0 @@
-const dictionary& simple = mesh.solutionDict().subDict("SIMPLE");
-
-const int nNonOrthCorr =
-    simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/residualControlsSolid.H

@@ -1 +0,0 @@
-../../solid/residualControlsSolid.H

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H

@@ -1 +0,0 @@
-../../solid/setRegionSolidFields.H

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H

@@ -1,28 +0,0 @@
-{
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        fvScalarMatrix hEqn
-        (
-            (
-                thermo.isotropic()
-              ? -fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
-              : -fvm::laplacian(betav*taniAlpha(), h, "laplacian(alpha,h)")
-            )
-          ==
-            fvOptions(rho, h)
-        );
-
-        hEqn.relax();
-
-        fvOptions.constrain(hEqn);
-
-        solvPerfh = hEqn.solve();
-
-        fvOptions.correct(h);
-    }
-}
-
-thermo.correct();
-
-Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-    << max(thermo.T()).value() << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/createFields.H

@@ -1,7 +1,2 @@
 #include "createFluidFields.H"
 #include "createSolidFields.H"
-
-SolverPerformance<vector> solvPerfU;
-SolverPerformance<scalar> solvPerfh;
-SolverPerformance<scalar> solvPerfhe;
-SolverPerformance<scalar> solvPerfp_rgh;

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H

@@ -27,7 +27,7 @@
 
     fvOptions.constrain(EEqn);
 
-    solvPerfhe = EEqn.solve(mesh.solver(he.select(finalIter)));
+    EEqn.solve();
 
     fvOptions.correct(he);
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H

@@ -16,9 +16,8 @@
 
     fvOptions.constrain(UEqn);
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
-        solvPerfU =
         solve
         (
             UEqn
@@ -29,10 +28,11 @@
                   - ghf*fvc::snGrad(rho)
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
-                ),
+            )
-                mesh.solver(U.select(finalIter))
         );
 
         fvOptions.correct(U);
         K = 0.5*magSqr(U);
     }
+
+    fvOptions.correct(U);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.H

@@ -1,21 +0,0 @@
-{
-    dimensionedScalar totalMass = fvc::domainIntegrate(rho);
-
-    scalar sumLocalContErr =
-    (
-        fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass
-    ).value();
-
-    scalar globalContErr =
-    (
-        fvc::domainIntegrate(rho - thermo.rho())/totalMass
-    ).value();
-
-    cumulativeContErr[i] += globalContErr;
-
-    Info<< "time step continuity errors (" << mesh.name() << ")"
-        << ": sum local = " << sumLocalContErr
-        << ", global = " << globalContErr
-        << ", cumulative = " << cumulativeContErr[i]
-        << endl;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -14,19 +14,3 @@
             CoNum
         );
     }
-/*
-    forAll(porousFluidRegions, porousI)
-    {
-        CoNum = max
-        (
-            compressibleCourantNo
-            (
-                porousFluidRegions[porousI],
-                runTime,
-                rhoPorous[porousI],
-                phiPorous[porousI]
-            ),
-            CoNum
-        );
-    }
-*/

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -18,10 +18,6 @@ PtrList<multivariateSurfaceInterpolationScheme<scalar>::fieldTable>
 PtrList<volScalarField> QdotFluid(fluidRegions.size());
 
 List<scalar> initialMassFluid(fluidRegions.size());
-List<bool> frozenFlowFluid(fluidRegions.size(), false);
-List<bool> residualReachedFluid(fluidRegions.size(), true);
-List<bool> residualControlUsedFluid(fluidRegions.size(), false);
-List<bool> firstIterationFluid(fluidRegions.size(), true);
 
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());
@@ -274,10 +270,6 @@ forAll(fluidRegions, i)
         )
     );
 
-    const dictionary& pimpleDict =
-        fluidRegions[i].solutionDict().subDict("PIMPLE");
-    pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
-
     Info<< "    Adding MRF\n" << endl;
     MRFfluid.set
     (
@@ -293,16 +285,4 @@ forAll(fluidRegions, i)
     );
 
     turbulenceFluid[i].validate();
-
-    if (pimpleDict.isDict("residualControl"))
-    {
-        Info<< "    Reading residual controls\n" << endl;
-
-        bool& residualControlUsed = residualControlUsedFluid[i];
-        bool& resReachedFluid = residualReachedFluid[i];
-        bool& firstIteration = firstIterationFluid[i];
-        const fvMesh& mesh = fluidRegions[i];
-
-        #include "readFluidMultiRegionResidualControls.H"
-    }
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidPressureControls.H

@@ -0,0 +1,16 @@
+PtrList<pressureControl> pressureControlFluid(fluidRegions.size());
+
+forAll(fluidRegions, i)
+{
+    pressureControlFluid.set
+    (
+        i,
+        new pressureControl
+        (
+            p_rghFluid[i],
+            rhoFluid[i],
+            pimples.pimple(i).dict(),
+            false
+        )
+    );
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H

@@ -1 +1 @@
-List<scalar> cumulativeContErr(fluidRegions.size(), 0.0);
+List<scalar> cumulativeContErrs(fluidRegions.size(), 0.0);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -1,69 +1,87 @@
-bool closedVolume = p_rgh.needReference();
+if (!mesh.steady() && !pimple.SIMPLErho())
-dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+{
-bool compressible = (compressibility.value() > small);
+    rho = thermo.rho();
-
+}
-rho = thermo.rho();
 
 volScalarField rAU("rAU", 1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+if (pimple.nCorrPISO() <= 1)
+{
+    tUEqn.clear();
+}
 
 surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    (
+    fvc::flux(rho*HbyA) + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, phi)
-    )
-  + phig
 );
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+const bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
+const bool adjustMass = closedVolume && !thermo.incompressible();
+
+phiHbyA += phig;
+
 // Update the pressure BCs to ensure flux consistency
 constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
 {
-    fvScalarMatrix p_rghDDtEqn
+    fvScalarMatrix p_rghEqnComp
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phiHbyA)
     );
 
+    if (pimple.transonic())
+    {
+        surfaceScalarField phid
+        (
+            "phid",
+            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+        );
+
+        phiHbyA -= fvc::interpolate(psi*p_rgh)*phiHbyA/fvc::interpolate(rho);
+
+        p_rghEqnComp += fvm::div(phid, p_rgh);
+    }
+
     // Thermodynamic density needs to be updated by psi*d(p) after the
     // pressure solution
-    const volScalarField psip0(psi*p);
+    tmp<volScalarField> psip0(mesh.steady() ? tmp<volScalarField>() : psi*p);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix p_rghEqn
+        fvScalarMatrix p_rghEqnIncomp
         (
-            p_rghDDtEqn
+            fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
-        solvPerfp_rgh = p_rghEqn.solve
+        fvScalarMatrix p_rghEqn(p_rghEqnComp + p_rghEqnIncomp);
+
+        p_rghEqn.setReference
         (
-            mesh.solver
+            pressureControl.refCell(),
-            (
+            pressureControl.refValue()
-                p_rgh.select
-                (
-                    (
-                        oCorr == nOuterCorr-1
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                )
-            )
         );
 
-        if (nonOrth == nNonOrthCorr)
+        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);
+                + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rhorAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
             K = 0.5*magSqr(U);
@@ -73,7 +91,10 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
     p = p_rgh + rho*gh;
 
     // Thermodynamic density update
+    if (!mesh.steady())
+    {
         thermo.correctRho(psi*p - psip0);
+    }
 }
 
 // Update pressure time derivative if needed
@@ -83,17 +104,43 @@ if (thermo.dpdt())
 }
 
 // Solve continuity
-#include "rhoEqn.H"
+if (!mesh.steady())
+{
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
+}
+else
+{
+    #include "incompressible/continuityErrs.H"
+}
 
-// Update continuity errors
+// Pressure limiting
-#include "compressibleContinuityErrors.H"
+const bool pLimited = pressureControl.limit(p);
 
-// For closed-volume cases adjust the pressure and density levels
+// For closed-volume compressible cases adjust the pressure level
 // to obey overall mass continuity
-if (closedVolume && compressible)
+if (adjustMass)
 {
     p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-        /compressibility;
+        /fvc::domainIntegrate(psi);
-    rho = thermo.rho();
     p_rgh = p - rho*gh;
 }
+
+if (adjustMass || pLimited)
+{
+    p.correctBoundaryConditions();
+}
+
+// Density updates
+if (adjustMass || pLimited || mesh.steady() || pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+if (mesh.steady() && !pimple.transonic())
+{
+    rho.relax();
+}
+
+Info<< "Min/max rho:" << min(rho).value() << ' '
+    << max(rho).value() << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H

@@ -1,10 +0,0 @@
-    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
-
-    const int nCorr =
-        pimple.lookupOrDefault<int>("nCorrectors", 1);
-
-    const int nNonOrthCorr =
-        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    const bool momentumPredictor =
-        pimple.lookupOrDefault("momentumPredictor", true);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionResidualControls.H

@@ -1,60 +0,0 @@
-const dictionary& residualControl =
-    mesh.solutionDict().subDict("PIMPLE").subOrEmptyDict("residualControl");
-
-scalar UTol = -1;
-scalar heTol = -1;
-scalar p_rghTol = -1;
-
-if (!residualControl.empty())
-{
-    if (!residualControl.subOrEmptyDict("U").empty())
-    {
-        UTol = readScalar(residualControl.subDict("U").lookup("tolerance"));
-    }
-
-    if (!residualControl.subOrEmptyDict("p_rgh").empty())
-    {
-        p_rghTol =
-            readScalar
-            (
-                residualControl.subDict("p_rgh").lookup("tolerance")
-            );
-    }
-
-    if (!residualControl.subOrEmptyDict("h").empty())
-    {
-        heTol = readScalar(residualControl.subDict("h").lookup("tolerance"));
-    }
-
-    // Residual control used?
-    if (UTol != -1 || heTol != -1 || p_rghTol != -1)
-    {
-        residualControlUsed = true;
-        resReachedFluid = false;
-    }
-
-    if (firstIteration)
-    {
-        firstIteration = false;
-
-        if (UTol != -1 || heTol != -1 || p_rghTol != -1)
-        {
-            Info<< "PIMPLE convergence criteria\n";
-
-            if (UTol != -1)
-            {
-                Info<< "   field U      tolerance " << UTol << "\n";
-            }
-            if (heTol != -1)
-            {
-                Info<< "   field h      tolerance " << heTol << "\n";
-            }
-            if (p_rghTol != -1)
-            {
-                Info<< "   field p_rgh  tolerance " << p_rghTol << "\n";
-            }
-
-            Info<< endl;
-        }
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/residualControlsFluid.H

@@ -1,55 +0,0 @@
-if (residualControlUsed)
-{
-    bool UConv = false;
-    bool p_rghConv = false;
-    bool heConv = false;
-
-    // Check which field is not used for control
-    {
-        if (UTol == -1 || !momentumPredictor)
-        {
-            UConv =  true;
-        }
-
-        if (p_rghTol == -1)
-        {
-            p_rghConv = true;
-        }
-
-        if (heTol == -1)
-        {
-            heConv = true;
-        }
-    }
-
-    // Get the last initial residual of the solvers
-    if (momentumPredictor && !UConv)
-    {
-        if (UTol > cmptMax(solvPerfU.initialResidual()))
-        {
-            UConv = true;
-        }
-    }
-
-    if (!p_rghConv)
-    {
-        if (p_rghTol > solvPerfp_rgh.initialResidual())
-        {
-            p_rghConv = true;
-        }
-    }
-
-    if (!heConv)
-    {
-        if (heTol > solvPerfhe.initialResidual())
-        {
-            heConv = true;
-        }
-    }
-
-    // Check if each field is converged
-    if (UConv && p_rghConv && heConv)
-    {
-        resReachedFluid = true;
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -1,8 +1,8 @@
-    fvMesh& mesh = fluidRegions[i];
+    const fvMesh& mesh = fluidRegions[i];
 
     CombustionModel<rhoReactionThermo>& reaction = reactionFluid[i];
 
-    rhoReactionThermo& thermo = reaction.thermo();
+    rhoReactionThermo& thermo = thermoFluid[i];
     thermo.validate(args.executable(), "h", "e");
 
     basicMultiComponentMixture& composition = thermo.composition();
@@ -57,8 +57,8 @@
         initialMassFluid[i]
     );
 
-    const bool frozenFlow = frozenFlowFluid[i];
+    pimpleNoLoopControl& pimple = pimples.pimple(i);
 
-    bool& resReachedFluid = residualReachedFluid[i];
+    pressureControl& pressureControl = pressureControlFluid[i];
-    bool& residualControlUsed = residualControlUsedFluid[i];
+
-    bool& firstIteration = firstIterationFluid[i];
+    scalar cumulativeContErr = cumulativeContErrs[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H

@@ -1,15 +1,10 @@
-if (finalIter)
+if (!pimple.solveFlow())
-{
-    mesh.data::add("finalIteration", true);
-}
-
-if (frozenFlow)
 {
     #include "EEqn.H"
 }
 else
 {
-    if (oCorr == 0)
+    if (!mesh.steady() && pimples.nCorrPIMPLE() <= 1)
     {
         #include "rhoEqn.H"
     }
@@ -19,17 +14,18 @@ else
     #include "EEqn.H"
 
     // --- PISO loop
-    for (int corr=0; corr<nCorr; corr++)
+    while (pimple.correct())
     {
         #include "pEqn.H"
     }
 
+    if (pimples.pimpleTurbCorr(i))
+    {
         turbulence.correct();
+    }
 
+    if (!mesh.steady() && pimples.finalIter())
+    {
         rho = thermo.rho();
-}
+    }
-
-if (finalIter)
-{
-    mesh.data::remove("finalIteration");
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H

@@ -1,2 +0,0 @@
-p_rghFluid[i].storePrevIter();
-rhoFluid[i].storePrevIter();

applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H

@@ -1,8 +0,0 @@
-    // We do not have a top-level mesh. Construct the fvSolution for
-    // the runTime instead.
-    fvSolution solutionDict(runTime);
-
-    const dictionary& pimple = solutionDict.subDict("PIMPLE");
-
-    const int nOuterCorr =
-        pimple.lookupOrDefault<int>("nOuterCorrectors", 1);

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -6,10 +6,6 @@ PtrList<fv::options> solidHeatSources(solidRegions.size());
 PtrList<volScalarField> betavSolid(solidRegions.size());
 PtrList<volSymmTensorField> aniAlphas(solidRegions.size());
 
-List<bool> residualReachedSolid(solidRegions.size(), true);
-List<bool> residualControlUsedSolid(solidRegions.size(), false);
-List<bool> firstIterationSolid(solidRegions.size(), true);
-
 // Populate solid field pointer lists
 forAll(solidRegions, i)
 {
@@ -108,19 +104,4 @@ forAll(solidRegions, i)
             )
         );
     }
-
-    const dictionary& pimpleDict =
-        solidRegions[i].solutionDict().subDict("PIMPLE");
-
-    if (pimpleDict.isDict("residualControl"))
-    {
-        Info<< "    Reading residual controls\n" << endl;
-
-        bool& residualControlUsed = residualControlUsedSolid[i];
-        bool& resReachedSolid = residualReachedSolid[i];
-        bool& firstIteration = firstIterationSolid[i];
-        const fvMesh& mesh = solidRegions[i];
-
-        #include "readSolidMultiRegionResidualControls.H"
-    }
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H

@@ -1,4 +0,0 @@
-const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
-
-int nNonOrthCorr =
-    pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionResidualControls.H

@@ -1,36 +0,0 @@
-const dictionary& residualControl =
-    mesh.solutionDict().subDict("PIMPLE").subOrEmptyDict("residualControl");
-
-scalar hTol = -1;
-
-if (!residualControl.empty())
-{
-    if (!residualControl.subOrEmptyDict("h").empty())
-    {
-        hTol = readScalar(residualControl.subDict("h").lookup("tolerance"));
-
-        // Used residual control for actual solid region
-        if (hTol != -1)
-        {
-            residualControlUsed = true;
-            resReachedSolid = false;
-        }
-    }
-
-    if (firstIteration)
-    {
-        firstIteration = false;
-
-        if (hTol != -1)
-        {
-            Info<< "PIMPLE convergence criteria\n";
-
-            if (hTol != -1)
-            {
-                Info<< "   field h      tolerance " << hTol << "\n";
-            }
-
-            Info<< endl;
-        }
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/residualControlsSolid.H

@@ -1,8 +0,0 @@
-// Residual control used
-if (residualControlUsed)
-{
-    if (hTol > solvPerfh.initialResidual())
-    {
-        resReachedSolid = true;
-    }
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -31,6 +31,4 @@ const volScalarField& betav = betavSolid[i];
 
 fv::options& fvOptions = solidHeatSources[i];
 
-bool& resReachedSolid = residualReachedSolid[i];
+solidNoLoopControl& pimple = pimples.solid(i);
-bool& residualControlUsed = residualControlUsedSolid[i];
-bool& firstIteration = firstIterationSolid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -1,10 +1,5 @@
-if (finalIter)
 {
-    mesh.data::add("finalIteration", true);
+    while (pimple.correctNonOrthogonal())
-}
-
-{
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix hEqn
         (
@@ -22,7 +17,7 @@ if (finalIter)
 
         fvOptions.constrain(hEqn);
 
-        solvPerfh = hEqn.solve(mesh.solver(h.select(finalIter)));
+        hEqn.solve(mesh.solver(h.select(pimples.finalIter())));
 
         fvOptions.correct(h);
     }
@@ -32,8 +27,3 @@ thermo.correct();
 
 Info<< "Min/max T:" << min(thermo.T()).value() << ' '
     << max(thermo.T()).value() << endl;
-
-if (finalIter)
-{
-    mesh.data::remove("finalIteration");
-}

tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/system/air/fvSolution

@@ -10,26 +10,23 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
+    location    "system/air";
     object      fvSolution;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 solvers
 {
-    p_rgh
+    "p_rgh.*"
     {
         solver           GAMG;
         tolerance        1e-7;
         relTol           0.01;
-
         smoother         DIC;
-
-
         maxIter          100;
     }
 
-    "(U|h|e|k|epsilon)"
+    "(U|h|e|k|epsilon).*"
     {
         solver           PBiCGStab;
         preconditioner   DILU;
@@ -38,27 +35,24 @@ solvers
     }
 }
 
-SIMPLE
+PIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    rhoMax          2;
-    rhoMin          1;
 }
 
-
 relaxationFactors
 {
     fields
     {
-        rho             1;
+        "rho.*"         1;
-        p_rgh           0.7;
+        "p_rgh.*"       0.7;
     }
     equations
     {
-        U               0.3;
+        "U.*"           0.3;
-        "(h|e)"         0.3;
+        "(h|e).*"       0.3;
-        k               0.3;
+        "k.*"           0.3;
-        epsilon         0.3;
+        "epsilon.*"     0.3;
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/system/controlDict

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     chtMultiRegionSimpleFoam;
+application     chtMultiRegionFoam;
 
 startFrom       latestTime;
 

tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/system/porous/fvSolution

@@ -10,25 +10,23 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
+    location    "system/porous";
     object      fvSolution;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 solvers
 {
-    p_rgh
+    "p_rgh.*"
     {
         solver           GAMG;
         tolerance        1e-7;
         relTol           0.01;
-
         smoother         DIC;
-
-
         maxIter          10;
     }
 
-    "(U|h|e|k|epsilon)"
+    "(U|h|e|k|epsilon).*"
     {
         solver           PBiCGStab;
         preconditioner   DILU;
@@ -37,26 +35,24 @@ solvers
     }
 }
 
-SIMPLE
+PIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    rhoMax          1100;
-    rhoMin          900;
 }
 
 relaxationFactors
 {
     fields
     {
-        rho             1;
+        "rho.*"         1;
-        p_rgh           0.7;
+        "p_rgh.*"       0.7;
     }
     equations
     {
-        U               0.3;
+        "U.*"           0.3;
-        "(h|e)"         0.3;
+        "(h|e).*"       0.3;
-        k               0.3;
+        "k.*"           0.3;
-        epsilon         0.3;
+        "epsilon.*"     0.3;
     }
 }