diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
index 9d9152930a9c359695a4e969886336ce38466f62..9b95197f541bfa2b20797804bf04054db815eb9f 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
@@ -1,10 +1,6 @@
 regionProperties/regionProperties.C
 
-coupleManager/coupleManager.C
-
 derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
-derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
-derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
 derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
 fluid/compressibleCourantNo.C
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
index a8797ca4564489292b11cf720b6ce80816a67192..12316d6bc7b7f06ed9d1f31316b9fd5c616d7e4e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@@ -2,9 +2,7 @@ EXE_INC = \
     -Ifluid \
     -Isolid \
     -IregionProperties \
-    -IcoupleManager \
-    -IderivedFvPatchFields/solidWallTemperatureCoupled \
-    -IderivedFvPatchFields/solidWallHeatFluxTemperatureCoupled \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index faff0edba042b2976c9fca21ec5fb09a0d8d6bca..7ad940bcff1bf5910c2c78253d5068b78a80f768 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -67,6 +67,7 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
+#       include "readPIMPLEControls.H"
 
         if (fluidRegions.size())
         {
@@ -78,22 +79,36 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        forAll(fluidRegions, i)
+        if (nOuterCorr != 1)
         {
-            Info<< "\nSolving for fluid region "
-                << fluidRegions[i].name() << endl;
-#           include "setRegionFluidFields.H"
-#           include "readFluidMultiRegionPISOControls.H"
-#           include "solveFluid.H"
+            forAll(fluidRegions, i)
+            {
+#               include "setRegionFluidFields.H"
+#               include "storeOldFluidFields.H"
+            }
         }
 
-        forAll(solidRegions, i)
+
+        // --- PIMPLE loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
-            Info<< "\nSolving for solid region "
-                << solidRegions[i].name() << endl;
-#           include "setRegionSolidFields.H"
-#           include "readSolidMultiRegionPISOControls.H"
-#           include "solveSolid.H"
+            forAll(fluidRegions, i)
+            {
+                Info<< "\nSolving for fluid region "
+                    << fluidRegions[i].name() << endl;
+#               include "setRegionFluidFields.H"
+#               include "readFluidMultiRegionPIMPLEControls.H"
+#               include "solveFluid.H"
+            }
+
+            forAll(solidRegions, i)
+            {
+                Info<< "\nSolving for solid region "
+                    << solidRegions[i].name() << endl;
+#               include "setRegionSolidFields.H"
+#               include "readSolidMultiRegionPIMPLEControls.H"
+#               include "solveSolid.H"
+            }
         }
 
         runTime.write();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C
deleted file mode 100644
index 1668e1144c46f13024e37813296b96642fbcbc92..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C
+++ /dev/null
@@ -1,186 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "coupleManager.H"
-#include "OFstream.H"
-#include "regionProperties.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::coupleManager::coupleManager(const fvPatch& patch)
-:
-    patch_(patch),
-    neighbourRegionName_("undefined-neighbourRegionName"),
-    neighbourPatchName_("undefined-neighbourPatchName"),
-    neighbourFieldName_("undefined-neighbourFieldName"),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-Foam::coupleManager::coupleManager
-(
-    const fvPatch& patch,
-    const dictionary& dict
-)
-:
-    patch_(patch),
-    neighbourRegionName_(dict.lookup("neighbourRegionName")),
-    neighbourPatchName_(dict.lookup("neighbourPatchName")),
-    neighbourFieldName_(dict.lookup("neighbourFieldName")),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-Foam::coupleManager::coupleManager
-(
-    const coupleManager& cm
-)
-:
-    patch_(cm.patch()),
-    neighbourRegionName_(cm.neighbourRegionName()),
-    neighbourPatchName_(cm.neighbourPatchName()),
-    neighbourFieldName_(cm.neighbourFieldName()),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::coupleManager::~coupleManager()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::coupleManager::regionOwner() const
-{
-    const fvMesh& nbrRegion = neighbourRegion();
-
-    const regionProperties& props =
-        localRegion_.objectRegistry::parent().lookupObject<regionProperties>
-        (
-            "regionProperties"
-        );
-
-    label myIndex = findIndex(props.fluidRegionNames(), localRegion_.name());
-    if (myIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), localRegion_.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn("coupleManager::regionOwner() const")
-                << "Cannot find region " << localRegion_.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        myIndex = props.fluidRegionNames().size() + i;
-    }
-    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
-    if (nbrIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn("coupleManager::regionOwner() const")
-                << "Cannot find region " << nbrRegion.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        nbrIndex = props.fluidRegionNames().size() + i;
-    }
-
-    return myIndex < nbrIndex;
-}
-
-
-void Foam::coupleManager::checkCouple() const
-{
-    Info<< "neighbourRegionName_ = " << neighbourRegionName_ << endl;
-    Info<< "neighbourPatchName_ = " << neighbourPatchName_ << endl;
-    Info<< "neighbourFieldName_ = " << neighbourFieldName_ << endl;
-
-    const fvPatch& nPatch = neighbourPatch();
-
-    if (patch_.size() != nPatch.size())
-    {
-        FatalErrorIn("Foam::coupleManager::checkCouple()")
-            << "Unequal patch sizes:" << nl
-            << "    patch name (size) = " << patch_.name()
-            << "(" << patch_.size() << ")" << nl
-            << "    neighbour patch name (size) = "
-            << nPatch.name() << "(" << patch_.size() << ")" << nl
-            << abort(FatalError);
-    }
-}
-
-
-void Foam::coupleManager::coupleToObj() const
-{
-    const fvPatch& nPatch = neighbourPatch();
-
-    OFstream obj
-        (
-             patch_.name() + "_to_" + nPatch.name() + "_couple.obj"
-        );
-    const vectorField& c1 = patch_.Cf();
-    const vectorField& c2 = neighbourPatch().Cf();
-
-    if (c1.size() != c2.size())
-    {
-        FatalErrorIn("coupleManager::coupleToObj() const")
-            << "Coupled patches are of unequal size:" << nl
-            << "    patch0 = " << patch_.name()
-            << "(" << patch_.size() <<  ")" << nl
-            << "    patch1 = " << nPatch.name()
-            << "(" << nPatch.size() <<  ")" << nl
-            << abort(FatalError);
-    }
-
-    forAll(c1, i)
-    {
-        obj << "v " << c1[i].x() << " " << c1[i].y() << " " << c1[i].z() << nl
-            << "v " << c2[i].x() << " " << c2[i].y() << " " << c2[i].z() << nl
-            << "l " << (2*i + 1) << " " << (2*i + 2) << endl;
-    }
-}
-
-
-void Foam::coupleManager::writeEntries(Ostream& os) const
-{
-    os.writeKeyword("neighbourRegionName");
-    os << neighbourRegionName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("neighbourPatchName");
-    os << neighbourPatchName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("neighbourFieldName");
-    os << neighbourFieldName_ << token::END_STATEMENT << nl;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H
deleted file mode 100644
index f1f73d23efe8505664cb5886dd3b37a5dc6727d2..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H
+++ /dev/null
@@ -1,170 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    coupleManager
-
-Description
-    Object to handle the coupling of region patches. It can be queried to
-    return the neighbour information.
-
-SourceFiles
-    coupleManager.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef coupleManager_H
-#define coupleManager_H
-
-#include "Ostream.H"
-#include "dictionary.H"
-#include "fvPatch.H"
-#include "fvMesh.H"
-#include "Time.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class coupleManager Declaration
-\*---------------------------------------------------------------------------*/
-
-class coupleManager
-{
-    // Private data
-
-        //- Reference to the local fvPatch
-        const fvPatch& patch_;
-
-        //- Name of neighbour region
-        word neighbourRegionName_;
-
-        //- Name of patch on the neighbour region
-        word neighbourPatchName_;
-
-        //- Name of field on the neighbour region
-        word neighbourFieldName_;
-
-        //- Reference to the local region
-        const fvMesh& localRegion_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise assignment
-        void operator=(const coupleManager&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct from fvPatch
-        coupleManager(const fvPatch& patch);
-
-        //- Construct from fvPatch and dictionary
-        coupleManager(const fvPatch& patch, const dictionary& dict);
-
-        //- Copy constructor
-        coupleManager(const coupleManager& cm);
-
-
-    // Destructor
-
-        ~coupleManager();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return a reference to the local patch
-            inline const fvPatch& patch() const;
-
-            //- Return the name of the neighbour region
-            inline const word& neighbourRegionName() const;
-
-            //- Return the name of the patch on the neighbour region
-            inline const word& neighbourPatchName() const;
-
-            //- Return the name of the field on the neighbour region
-            inline const word& neighbourFieldName() const;
-
-            //- Return a reference to the neighbour mesh
-            inline const fvMesh& neighbourRegion() const;
-
-            //- Return the neighbour patch ID
-            inline label neighbourPatchID() const;
-
-            //- Return a reference to the neighbour patch
-            inline const fvPatch& neighbourPatch() const;
-
-            //- Return a reference to the neighbour patch field
-            template<class Type>
-            inline const fvPatchField<Type>& neighbourPatchField() const;
-
-            //- Am I owner (= first to evaluate) of this region interface? 
-            bool regionOwner() const;
-
-            //- Check that the couple is valid
-            void checkCouple() const;
-
-
-        // Edit
-
-            //- Return the name of the neighbour region
-            word& neighbourRegionName();
-
-            //- Return the name of the patch on the neighbour region
-            word& neighbourPatchName();
-
-            //- Return the name of the field on the neighbour region
-            word& neighbourFieldName();
-
-
-        // Write
-
-           //- Write couple to obj file
-           void coupleToObj() const;
-
-           //- Write entries for patches
-           void writeEntries(Ostream& os) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "coupleManagerI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManagerI.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManagerI.H
deleted file mode 100644
index 45e059ebae27d425e90403927352f7f36c5e8cde..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManagerI.H
+++ /dev/null
@@ -1,98 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-inline const Foam::fvPatch& Foam::coupleManager::patch() const
-{
-    return patch_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourRegionName() const
-{
-    return neighbourRegionName_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourPatchName() const
-{
-    return neighbourPatchName_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourFieldName() const
-{
-     return neighbourFieldName_;
-};
-
-
-inline const Foam::fvMesh& Foam::coupleManager::neighbourRegion() const
-{
-    return localRegion_.objectRegistry::parent()
-        .lookupObject<fvMesh>(neighbourRegionName_);
-}
-
-
-inline Foam::label Foam::coupleManager::neighbourPatchID() const
-{
-    return neighbourRegion().boundaryMesh().findPatchID(neighbourPatchName_);
-}
-
-
-inline const Foam::fvPatch& Foam::coupleManager::neighbourPatch() const
-{
-    return neighbourRegion().boundary()[neighbourPatchID()];
-}
-
-
-template<class Type>
-inline const Foam::fvPatchField<Type>&
-Foam::coupleManager::neighbourPatchField() const
-{
-    return neighbourPatch().lookupPatchField
-        <GeometricField<Type, fvPatchField, volMesh>, Type>
-            (neighbourFieldName_);
-}
-
-
-inline Foam::word& Foam::coupleManager::neighbourRegionName()
-{
-    return neighbourRegionName_;
-};
-
-
-inline Foam::word& Foam::coupleManager::neighbourPatchName()
-{
-    return neighbourPatchName_;
-};
-
-
-inline Foam::word& Foam::coupleManager::neighbourFieldName()
-{
-    return neighbourFieldName_;
-};
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
deleted file mode 100644
index 07bcac29a683edc4dcb95805f63df358f7613d77..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
+++ /dev/null
@@ -1,150 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-
-#include "solidWallTemperatureCoupledFvPatchScalarField.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    coupleManager_(p),
-    KName_("undefined-K")
-{}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const solidWallHeatFluxTemperatureCoupledFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
-    KName_(ptf.KName_)
-{}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
-    KName_(dict.lookup("K"))
-{
-    if (dict.found("value"))
-    {
-        fvPatchField<scalar>::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-    }
-    else
-    {
-        evaluate();
-    }
-}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const solidWallHeatFluxTemperatureCoupledFvPatchScalarField& whftcsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(whftcsf, iF),
-    coupleManager_(whftcsf.coupleManager_),
-    KName_(whftcsf.KName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const fvPatchField<scalar>& neighbourField =
-        coupleManager_.neighbourPatchField<scalar>();
-
-    const fvPatchField<scalar>& K =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
-
-    gradient() = -refCast<const solidWallTemperatureCoupledFvPatchScalarField>
-        (neighbourField).flux()/K;
-
-    fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
-    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H
deleted file mode 100644
index b7867ca8ce4a125415906382d43c43cedb916339..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H
+++ /dev/null
@@ -1,165 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-
-Description
-    Fixed heat-flux boundary condition for temperature, to be used by the
-    conjugate heat transfer solver.
-
-    Example usage:
-        myInterfacePatchName
-        {
-            type                solidWallHeatFluxTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
-            neighbourFieldName  T;
-            K                   K;
-            value               uniform 300;
-        }
-
-SourceFiles
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallHeatFluxTemperatureCoupledFvPatchScalarField_H
-#define solidWallHeatFluxTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "fixedGradientFvPatchFields.H"
-#include "coupleManager.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-   Class solidWallHeatFluxTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-:
-    public fixedGradientFvPatchScalarField
-{
-    // Private data
-
-        //- Couple manager object
-        coupleManager coupleManager_;
-
-        //- Name of thermal conductivity field
-        word KName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("solidWallHeatFluxTemperatureCoupled");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        // solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        //  onto a new patch
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallHeatFluxTemperatureCoupledFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-                    (
-                        *this
-                    )
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallHeatFluxTemperatureCoupledFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-                    (
-                        *this,
-                        iF
-                    )
-            );
-        }
-
-
-    // Member functions
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
index b6bcfffa3f63f947ec0a84734b12df7ff39893a4..5bec5673108e4cd0c0684e041b6443f769c0ff55 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
@@ -28,8 +28,62 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
+#include "directMappedPatchBase.H"
 #include "regionProperties.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::interfaceOwner
+(
+    const polyMesh& nbrRegion
+) const
+{
+    const fvMesh& myRegion = patch().boundaryMesh().mesh();
+
+    const regionProperties& props =
+        myRegion.objectRegistry::parent().lookupObject<regionProperties>
+        (
+            "regionProperties"
+        );
+
+    label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
+    if (myIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), myRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn
+            (
+                "solidWallMixedTemperatureCoupledFvPatchScalarField"
+                "::interfaceOwner(const polyMesh&) const"
+            )   << "Cannot find region " << myRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        myIndex = props.fluidRegionNames().size() + i;
+    }
+    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
+    if (nbrIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn("coupleManager::interfaceOwner(const polyMesh&) const")
+                << "Cannot find region " << nbrRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        nbrIndex = props.fluidRegionNames().size() + i;
+    }
+
+    return myIndex < nbrIndex;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
@@ -40,7 +94,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
-    coupleManager_(p),
+    neighbourFieldName_("undefined-neighbourFieldName"),
     KName_("undefined-K")
 {
     this->refValue() = 0.0;
@@ -60,7 +114,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
+    neighbourFieldName_(ptf.neighbourFieldName_),
     KName_(ptf.KName_),
     fixesValue_(ptf.fixesValue_)
 {}
@@ -75,14 +129,46 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
+    neighbourFieldName_(dict.lookup("neighbourFieldName")),
     KName_(dict.lookup("K"))
 {
+    if (!isA<directMappedPatchBase>(this->patch().patch()))
+    {
+        FatalErrorIn
+        (
+            "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+            "solidWallMixedTemperatureCoupledFvPatchScalarField\n"
+            "(\n"
+            "    const fvPatch& p,\n"
+            "    const DimensionedField<scalar, volMesh>& iF,\n"
+            "    const dictionary& dict\n"
+            ")\n"
+        )   << "\n    patch type '" << p.type()
+            << "' not type '" << directMappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << dimensionedInternalField().name()
+            << " in file " << dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
     fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
-    refValue() = static_cast<scalarField>(*this);
-    refGrad() = 0.0;
-    valueFraction() = 1.0;
-    fixesValue_ = true;
+
+    if (dict.found("refValue"))
+    {
+        // Full restart
+        refValue() = scalarField("refValue", dict, p.size());
+        refGrad() = scalarField("refGradient", dict, p.size());
+        valueFraction() = scalarField("valueFraction", dict, p.size());
+        fixesValue_ = readBool(dict.lookup("fixesValue"));
+    }
+    else
+    {
+        // Start from user entered data. Assume fixedValue.
+        refValue() = *this;
+        refGrad() = 0.0;
+        valueFraction() = 1.0;
+        fixesValue_ = true;
+    }
 }
 
 
@@ -94,7 +180,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(wtcsf, iF),
-    coupleManager_(wtcsf.coupleManager_),
+    neighbourFieldName_(wtcsf.neighbourFieldName_),
     KName_(wtcsf.KName_),
     fixesValue_(wtcsf.fixesValue_)
 {}
@@ -116,33 +202,111 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
         return;
     }
 
+    // Get the coupling information from the directMappedPatchBase
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
+    (
+        patch().patch()
+    );
+    const polyMesh& nbrMesh = mpp.sampleMesh();
+
     tmp<scalarField> intFld = patchInternalField();
 
+    if (interfaceOwner(nbrMesh))
+    {
+        // Note: other side information could be cached - it only needs
+        // to be updated the first time round the iteration (i.e. when
+        // switching regions) but unfortunately we don't have this information.
+
+        const mapDistribute& distMap = mpp.map();
+        const fvPatch& nbrPatch = refCast<const fvMesh>
+        (
+            nbrMesh
+        ).boundary()[mpp.samplePolyPatch().index()];
+
+
+        // Calculate the temperature by harmonic averaging
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        const solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
+        refCast<const solidWallMixedTemperatureCoupledFvPatchScalarField>
+        (
+            nbrPatch.lookupPatchField<volScalarField, scalar>
+            (
+                neighbourFieldName_
+            )
+        );
+
+        // Swap to obtain full local values of neighbour internal field
+        scalarField nbrIntFld = nbrField.patchInternalField();
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            distMap.constructSize(),
+            distMap.subMap(),           // what to send
+            distMap.constructMap(),     // what to receive
+            nbrIntFld
+        );
+
+        // Swap to obtain full local values of neighbour K*delta
+        scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            distMap.constructSize(),
+            distMap.subMap(),           // what to send
+            distMap.constructMap(),     // what to receive
+            nbrKDelta
+        );
+
+
+        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
+
+        // Calculate common wall temperature. Reuse *this to store common value.
+        scalarField Twall
+        (
+            (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
+          / (myKDelta() + nbrKDelta)
+        );
+        // Assign to me
+        fvPatchScalarField::operator=(Twall);
+        // Distribute back and assign to neighbour
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            nbrField.size(),
+            distMap.constructMap(),     // reverse : what to send
+            distMap.subMap(),
+            Twall
+        );
+        const_cast<solidWallMixedTemperatureCoupledFvPatchScalarField&>
+        (
+            nbrField
+        ).fvPatchScalarField::operator=(Twall);
+    }
+
+
+    // Switch between fixed value (of harmonic avg) or gradient
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     label nFixed = 0;
 
     // Like snGrad but bypass switching on refValue/refGrad.
-    tmp<scalarField> normalGradient =
-        (*this-intFld())
-      * patch().deltaCoeffs();
+    tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
 
     if (debug)
     {
+        scalar Q = gSum(K()*patch().magSf()*normalGradient());
+
         Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
             << "updateCoeffs() :"
             << " patch:" << patch().name()
+            << " heatFlux:" << Q
             << " walltemperature "
-            << " min:"
-            <<  returnReduce
-                (
-                    (this->size() > 0 ? min(*this) : VGREAT),
-                    minOp<scalar>()
-                )
-            << " max:"
-            <<  returnReduce
-                (
-                    (this->size() > 0 ? max(*this) : -VGREAT),
-                    maxOp<scalar>()
-                )
+            << " min:" << gMin(*this)
+            << " max:" << gMax(*this)
             << " avg:" << gAverage(*this)
             << endl;
     }
@@ -150,7 +314,7 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
     forAll(*this, i)
     {
         // if outgoing flux use fixed value.
-        if (intFld()[i] > operator[](i))
+        if (normalGradient()[i] < 0.0)
         {
             this->refValue()[i] = operator[](i);
             this->refGrad()[i] = 0.0;   // not used
@@ -185,80 +349,16 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
 }
 
 
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::evaluate
-(
-    const Pstream::commsTypes
-)
-{
-    if (!this->updated())
-    {
-        this->updateCoeffs();
-    }
-
-    if (!coupleManager_.regionOwner())
-    {
-        // I am the last one to evaluate.
-
-        tmp<scalarField> intFld = patchInternalField();
-
-        const fvPatch& nbrPatch = coupleManager_.neighbourPatch();
-
-        solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
-        refCast<solidWallMixedTemperatureCoupledFvPatchScalarField>
-        (
-            const_cast<fvPatchField<scalar>&>
-            (
-                coupleManager_.neighbourPatchField<scalar>()
-            )
-        );
-        tmp<scalarField> nbrIntFld = nbrField.patchInternalField();
-        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
-        tmp<scalarField> nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
-
-        // Calculate common wall temperature and assign to both sides
-        scalarField::operator=
-        (
-            (myKDelta()*intFld + nbrKDelta()*nbrIntFld)
-          / (myKDelta() + nbrKDelta())
-        );
-
-        nbrField.scalarField::operator=(*this);
-
-        if (debug)
-        {
-            Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-                << "updateCoeffs() :"
-                << " patch:" << patch().name()
-                << " setting master and slave to wall temperature "
-                << " min:"
-                <<  returnReduce
-                    (
-                        (this->size() > 0 ? min(*this) : VGREAT),
-                        minOp<scalar>()
-                    )
-                << " max:"
-                <<  returnReduce
-                    (
-                        (this->size() > 0 ? max(*this) : -VGREAT),
-                        maxOp<scalar>()
-                    )
-                << " avg:" << gAverage(*this)
-                << endl;
-        }
-    }
-
-    fvPatchScalarField::evaluate();
-}
-
-
 void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::write
 (
     Ostream& os
 ) const
 {
     mixedFvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
+    os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
+        << token::END_STATEMENT << nl;
     os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
 }
 
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
index 34c32b1abf4204864d1cc7268ecf3e5874a25d36..a358815d3ea7e741baa09e37200fc86ea66119e9 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
@@ -37,13 +37,26 @@ Description
         myInterfacePatchName
         {
             type                solidWallMixedTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
             neighbourFieldName  T;
             K                   K;
             value               uniform 300;
         }
 
+    Needs to be on underlying directMapped(Wall)FvPatch.
+
+
+    Note: lags interface data so both sides use same data.
+    Old scheme:
+    - last-to-solve calculates harmonic average
+    - gets used next time on both regions.
+
+    New scheme:
+    - problem: schedule to calculate average would interfere
+    with standard processor swaps.
+    - so: updateCoeffs sets both to same Twall. Only need to do
+    this for last outer iteration but don't have access to this.
+
+
 SourceFiles
     solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
@@ -54,7 +67,6 @@ SourceFiles
 
 #include "fvPatchFields.H"
 #include "mixedFvPatchFields.H"
-#include "coupleManager.H"
 #include "fvPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -72,14 +84,21 @@ class solidWallMixedTemperatureCoupledFvPatchScalarField
 {
     // Private data
 
-        //- Couple manager object
-        coupleManager coupleManager_;
-
+        //- Name of field on the neighbour region
+        const word neighbourFieldName_;
+        
         //- Name of thermal conductivity field
-        word KName_;
+        const word KName_;
 
         bool fixesValue_;
 
+
+    // Private Member Functions
+
+        //- Am I or neighbour owner of interface
+        bool interfaceOwner(const polyMesh& nbrRegion) const;
+
+
 public:
 
     //- Runtime type information
@@ -162,12 +181,6 @@ public:
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 
-        //- Evaluate the patch field
-        virtual void evaluate
-        (
-            const Pstream::commsTypes commsType=Pstream::blocking
-        );
-
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
deleted file mode 100644
index 05bbdd1ba3fb75adce7642f1afd732f8c87d3c3d..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
+++ /dev/null
@@ -1,156 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "solidWallTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(p, iF),
-    coupleManager_(p),
-    KName_("undefined-K")
-{}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const solidWallTemperatureCoupledFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
-    KName_(ptf.KName_)
-{}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedValueFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
-    KName_(dict.lookup("K"))
-{
-    if (dict.found("value"))
-    {
-        fvPatchField<scalar>::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-    }
-    else
-    {
-        evaluate();
-    }
-}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const solidWallTemperatureCoupledFvPatchScalarField& wtcsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(wtcsf, iF),
-    coupleManager_(wtcsf.coupleManager_),
-    KName_(wtcsf.KName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::solidWallTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const fvPatchField<scalar>& neighbourField =
-        coupleManager_.neighbourPatchField<scalar>();
-
-    operator==(neighbourField);
-
-    fixedValueFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallTemperatureCoupledFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
-    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::solidWallTemperatureCoupledFvPatchScalarField::flux() const
-{
-    const fvPatchScalarField& Kw =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
-
-    const fvPatchScalarField& Tw = *this;
-
-    return Tw.snGrad()*Kw;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    solidWallTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H
deleted file mode 100644
index 22fd168175134eb812ac447caac3d397dee21df3..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H
+++ /dev/null
@@ -1,160 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    solidWallHeatFluxCoupledFvPatchScalarField
-
-Description
-    Fixed value boundary condition for temperature, to be used by the
-    conjugate heat transfer solver.
-
-    Example usage:
-        myInterfacePatchName
-        {
-            type                solidWallTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
-            neighbourFieldName  T;
-            K                   K;
-            value               uniform 300;
-        }
-
-SourceFiles
-    solidWallTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallTemperatureCoupledFvPatchScalarField_H
-#define solidWallTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "fixedValueFvPatchFields.H"
-#include "coupleManager.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-        Class solidWallTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallTemperatureCoupledFvPatchScalarField
-:
-    public fixedValueFvPatchScalarField
-{
-    // Private data
-
-        //- Couple manager object
-        coupleManager coupleManager_;
-
-        //- Name of thermal conductivity field
-        word KName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("solidWallTemperatureCoupled");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given solidWallTemperatureCoupledFvPatchScalarField
-        //  onto a new patch
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallTemperatureCoupledFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallTemperatureCoupledFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallTemperatureCoupledFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallTemperatureCoupledFvPatchScalarField(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        //- (intensive) flux
-        tmp<scalarField> flux() const;
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
index e719d92433926797b7850dc356ade232228ba502..65467f80864609f9225e9d058b9568890c500988 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
@@ -14,12 +14,10 @@
          (
             UEqn()
          ==
-           -fvc::reconstruct
+            fvc::reconstruct
             (
-                (
-                    fvc::snGrad(pd)
-                  + ghf*fvc::snGrad(rho)
-                ) * mesh.magSf()
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              - fvc::snGrad(p)*mesh.magSf()
             )
         );
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index 1f6e50de0a3872b76c239e920f054d5e1cb54ea2..fd9624685f8fcc07cc4f82d014676fb469bbd064 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -6,43 +6,15 @@
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
     PtrList<volScalarField> DpDtFluid(fluidRegions.size());
-    PtrList<volScalarField> ghFluid(fluidRegions.size());
-    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
-    PtrList<volScalarField> pdFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
 
-    dimensionedScalar pRef
-    (
-        "pRef",
-        dimensionSet(1, -1, -2, 0, 0),
-        rp.lookup("pRef")
-    );
-
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
     {
         Info<< "*** Reading fluid mesh thermophysical properties for region "
             << fluidRegions[i].name() << nl << endl;
 
-        Info<< "    Adding to pdFluid\n" << endl;
-        pdFluid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "pd",
-                    runTime.timeName(),
-                    fluidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                fluidRegions[i]
-            )
-        );
-
         Info<< "    Adding to thermoFluid\n" << endl;
         thermoFluid.set
         (
@@ -145,6 +117,7 @@
             i,
             new volScalarField
             (
+                "DpDt",
                 fvc::DDt
                 (
                     surfaceScalarField
@@ -157,36 +130,5 @@
             )
         );
 
-        const dictionary& environmentalProperties =
-            fluidRegions[i].lookupObject<IOdictionary>
-                ("environmentalProperties");
-        dimensionedVector g(environmentalProperties.lookup("g"));
-
-        Info<< "    Adding to ghFluid\n" << endl;
-        ghFluid.set
-        (
-            i,
-            new volScalarField
-            (
-                "gh",
-                 g & fluidRegions[i].C()
-            )
-        );
-        ghfFluid.set
-        (
-            i,
-            new surfaceScalarField
-            (
-                "ghf",
-                 g & fluidRegions[i].Cf()
-            )
-        );
-
-        Info<< "    Updating p from pd\n" << endl;
-        thermoFluid[i].p() == pdFluid[i] + rhoFluid[i]*ghFluid[i] + pRef;
-        thermoFluid[i].correct();
-
         initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
     }
-
-
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
index 7f2202d5937d4e17ab988aef292979339269aabb..e070537db2c701e27baafd6cf64ef015e3898100 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
@@ -1,5 +1,5 @@
 {
-    tmp<fvScalarMatrix> hEqn
+    fvScalarMatrix hEqn
     (
         fvm::ddt(rho, h)
       + fvm::div(phi, h)
@@ -7,8 +7,16 @@
      ==
         DpDt
     );
-    hEqn().relax();
-    hEqn().solve();
+    if (oCorr == nOuterCorr-1)
+    {
+        hEqn.relax();
+        hEqn.solve(mesh.solver("hFinal"));
+    }
+    else
+    {
+        hEqn.relax();
+        hEqn.solve();
+    }
 
     thermo.correct();
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index e297989809aaa930b9c0be00e3049842869e44aa..a264b68fe5eab2388d38ca1aaf75a3f416db2330 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -1,5 +1,5 @@
 {
-    bool closedVolume = pd.needReference();
+    bool closedVolume = p.needReference();
 
     rho = thermo.rho();
 
@@ -17,31 +17,34 @@
         )
     );
 
-    phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
+    phi = phiU + fvc::interpolate(rho)*(g & mesh.Sf())*rhorUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, pd)
-          + fvc::ddt(psi)*pRef
-          + fvc::ddt(psi, rho)*gh
+            fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, pd)
+          - fvm::laplacian(rhorUAf, p)
         );
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        if
+        (
+            oCorr == nOuterCorr-1
+         && corr == nCorr-1
+         && nonOrth == nNonOrthCorr
+        )
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pdEqn.flux();
+            phi += pEqn.flux();
         }
     }
 
@@ -49,27 +52,24 @@
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
 
-    // Update pressure field (including bc)
-    p == pd + rho*gh + pRef;
+    // Update pressure substantive derivative
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
     // Solve continuity
-#   include "rhoEqn.H"
+    #include "rhoEqn.H"
 
     // Update continuity errors
-#   include "compressibleContinuityErrors.H"
+    #include "compressibleContinuityErrors.H"
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
     {
-        p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi);
+        p += (massIni - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
         rho = thermo.rho();
     }
 
     // Update thermal conductivity
     K = thermoFluid[i].Cp()*turb.alphaEff();
-
-    // Update pd (including bc)
-    pd == p - (rho*gh + pRef);
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..413c0225f0a34a82b3af1751011d7e8af7b935ba
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
@@ -0,0 +1,9 @@
+    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
+
+    int nCorr(readInt(pimple.lookup("nCorrectors")));
+
+    int nNonOrthCorr =
+        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
+
+    bool momentumPredictor =
+        pimple.lookupOrDefault<Switch>("momentumPredictor", true);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H
deleted file mode 100644
index a6b0ac9fe1c0c27d1f2479f27ab3715c87537cd9..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H
+++ /dev/null
@@ -1 +0,0 @@
-    solve(fvm::ddt(rho) + fvc::div(phi));
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 35a64418a8c6cbd1d488970f1d665420a2f752b6..72c9bbc4faf7e81c0fd759466c210a536929ff2c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -7,12 +7,15 @@
     surfaceScalarField& phi = phiFluid[i];
     compressible::turbulenceModel& turb = turbulence[i];
     volScalarField& DpDt = DpDtFluid[i];
-    const volScalarField& gh = ghFluid[i];
-    const surfaceScalarField& ghf = ghfFluid[i];
-    volScalarField& pd = pdFluid[i];
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
     volScalarField& h = thermo.h();
 
     const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);
+
+    const dictionary& environmentalProperties =
+        fluidRegions[i].lookupObject<IOdictionary>
+        ("environmentalProperties");
+
+    const dimensionedVector g(environmentalProperties.lookup("g"));
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
index 19ec50cac253986322e9155df1c2bcc0d7632c4a..86dd4344c15310f219bfb9bda406b22f2ffe7585 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
@@ -1,15 +1,18 @@
-#   include "rhoEqn.H"
-    for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
-    {
-    #   include "UEqn.H"
+if (oCorr == 0)
+{
+    #include "rhoEqn.H"
+}
 
-    #   include "hEqn.H"
+#include "UEqn.H"
 
-     // --- PISO loop
+#include "hEqn.H"
 
-        for (int corr=0; corr<nCorr; corr++)
-        {
-    #       include "pEqn.H"
-        }
-    }
-    turb.correct();
+// --- PISO loop
+for (int corr=0; corr<nCorr; corr++)
+{
+    #include "pEqn.H"
+}
+
+turb.correct();
+
+rho = thermo.rho();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..f63e85458e25253209094f92995e0c3c91858fc5
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H
@@ -0,0 +1,2 @@
+    p.storePrevIter();
+    rho.storePrevIter();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..42793d9b9fdb82a890b5994673f07ade7745a659
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H
@@ -0,0 +1,7 @@
+    // We do not have a top-level mesh. Construct the fvSolution for
+    // the runTime instead.
+    fvSolution solutionDict(runTime);
+
+    const dictionary& pimple = solutionDict.subDict("PIMPLE");
+
+    int nOuterCorr(readInt(pimple.lookup("nOuterCorrectors")));
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..e23883c5fae11f298e0459feaacacc4c73af8dbe
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H
@@ -0,0 +1,4 @@
+    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
+
+    int nNonOrthCorr =
+        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H
index 6373e79d49c0869238a30a138aef81a7af905d05..ce6a1c5bb2626f82eb411bc20b4f57b64e3271ff 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H
@@ -1,4 +1,4 @@
-    dictionary piso = solidRegions[i].solutionDict().subDict("PISO");
+    const dictionary& piso = solidRegions[i].solutionDict().subDict("PISO");
 
     int nNonOrthCorr = 0;
     if (piso.found("nNonOrthogonalCorrectors"))
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
index 04c90d2c4c108f479145e86e64c980f7d6296cba..f9e80e3d72385e7f6ae80107cc3b65921904c35f 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
@@ -1,4 +1,4 @@
-//    fvMesh& mesh = solidRegions[i];
+    fvMesh& mesh = solidRegions[i];
 
     volScalarField& rho = rhos[i];
     volScalarField& cp = cps[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
index 5fa731824f375adfeb3e2175f280de55da5cd4b2..ce8b1d0f408d4ec033f88c37846a8c7ab80e61b0 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
@@ -1,10 +1,13 @@
 {
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        solve
+        tmp<fvScalarMatrix> TEqn
         (
-            fvm::ddt(rho*cp, T) - fvm::laplacian(K, T)
+            fvm::ddt(rho*cp, T)
+          - fvm::laplacian(K, T)
         );
+        TEqn().relax();
+        TEqn().solve();
     }
 
     Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;