From ef44df91f237537491b2c1e002a66103e207b354 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Thu, 23 Nov 2023 16:25:38 +0100
Subject: [PATCH] ENH: support direct lookup of solver controls

    OLD:
        pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
        pEqn.solve(mesh.solver("Yi"));

    NEW:
        pEqn.solve(p.select(piso.finalInnerIter()));
        pEqn.solve("Yi");
---
 applications/solvers/DNS/dnsFoam/dnsFoam.C    |  2 +-
 .../solvers/combustion/PDRFoam/pEqn.H         |  4 +-
 .../combustion/XiFoam/XiDyMFoam/pEqn.H        |  4 +-
 .../combustion/XiFoam/XiEngineFoam/pEqn.H     |  4 +-
 applications/solvers/combustion/XiFoam/pEqn.H |  4 +-
 .../solvers/combustion/chemFoam/YEqn.H        |  2 +-
 .../solvers/combustion/fireFoam/YEEqn.H       |  2 +-
 .../solvers/combustion/fireFoam/pEqn.H        |  2 +-
 .../solvers/combustion/reactingFoam/YEqn.H    |  2 +-
 .../solvers/combustion/reactingFoam/pEqn.H    |  4 +-
 .../solvers/combustion/reactingFoam/pcEqn.H   |  4 +-
 .../rhoReactingBuoyantFoam/pEqn.H             |  2 +-
 .../rhoPimpleAdiabaticFoam/pEqn.H             |  2 +-
 .../overRhoPimpleDyMFoam/correctPhi.H         |  2 +-
 .../rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H |  4 +-
 .../solvers/compressible/rhoPimpleFoam/pEqn.H |  4 +-
 .../compressible/rhoPimpleFoam/pcEqn.H        |  4 +-
 .../electromagnetics/mhdFoam/mhdFoam.C        |  4 +-
 .../buoyantBoussinesqPimpleFoam/pEqn.H        |  2 +-
 .../overBuoyantPimpleDyMFoam/correctPhi.H     |  2 +-
 .../overBuoyantPimpleDyMFoam/pEqn.H           |  2 +-
 .../heatTransfer/buoyantPimpleFoam/pEqn.H     |  2 +-
 .../solid/solveSolid.H                        |  2 +-
 .../chtMultiRegionFoam/fluid/EEqn.H           |  2 +-
 .../chtMultiRegionFoam/fluid/UEqn.H           |  2 +-
 .../chtMultiRegionFoam/fluid/YEqn.H           |  2 +-
 .../chtMultiRegionFoam/fluid/pEqn.H           | 11 ++--
 .../chtMultiRegionFoam/solid/solveSolid.H     |  2 +-
 .../solvers/incompressible/icoFoam/icoFoam.C  |  2 +-
 .../nonNewtonianIcoFoam/nonNewtonianIcoFoam.C |  2 +-
 .../pimpleFoam/SRFPimpleFoam/pEqn.H           |  2 +-
 .../pimpleFoam/overPimpleDyMFoam/pEqn.H       |  2 +-
 .../solvers/incompressible/pimpleFoam/pEqn.H  |  2 +-
 .../solvers/incompressible/pisoFoam/pEqn.H    |  2 +-
 .../shallowWaterFoam/shallowWaterFoam.C       |  2 +-
 .../lagrangian/DPMFoam/DPMDyMFoam/pEqn.H      |  2 +-
 .../solvers/lagrangian/DPMFoam/pEqn.H         |  2 +-
 .../lagrangian/coalChemistryFoam/YEqn.H       |  2 +-
 .../lagrangian/coalChemistryFoam/pEqn.H       |  4 +-
 .../lagrangian/kinematicParcelFoam/pEqn.H     |  2 +-
 .../lagrangian/reactingParcelFoam/YEqn.H      |  2 +-
 .../lagrangian/reactingParcelFoam/pEqn.H      |  2 +-
 .../simpleReactingParcelFoam/YEqn.H           |  2 +-
 .../lagrangian/simpleCoalParcelFoam/YEqn.H    |  2 +-
 .../solvers/lagrangian/sprayFoam/YEqn.H       |  2 +-
 .../solvers/lagrangian/sprayFoam/pEqn.H       |  4 +-
 .../lagrangian/sprayFoam/sprayDyMFoam/pEqn.H  |  4 +-
 .../solvers/multiphase/MPPICInterFoam/pEqn.H  |  2 +-
 .../cavitatingFoam/cavitatingDyMFoam/pEqn.H   |  2 +-
 .../solvers/multiphase/cavitatingFoam/pEqn.H  |  2 +-
 .../compressibleInterDyMFoam/pEqn.H           |  2 +-
 .../compressibleInterFilmFoam/pEqn.H          |  2 +-
 .../compressibleInterIsoFoam/pEqn.H           |  2 +-
 .../overCompressibleInterDyMFoam/pEqn.H       |  2 +-
 .../multiphase/compressibleInterFoam/pEqn.H   |  2 +-
 .../compressibleMultiphaseInterFoam/pEqn.H    |  2 +-
 .../driftFluxFoam/alphaEqnSubCycle.H          |  2 +-
 .../solvers/multiphase/driftFluxFoam/pEqn.H   |  2 +-
 .../icoReactingMultiphaseInterFoam/pEqn.H     |  2 +-
 .../interCondensatingEvaporatingFoam/pEqn.H   |  2 +-
 .../interFoam/overInterDyMFoam/pEqn.H         |  2 +-
 .../solvers/multiphase/interFoam/pEqn.H       |  2 +-
 .../interPhaseChangeDyMFoam/pEqn.H            |  2 +-
 .../overInterPhaseChangeDyMFoam/pEqn.H        |  2 +-
 .../multiphase/interPhaseChangeFoam/pEqn.H    |  2 +-
 .../multiphase/multiphaseEulerFoam/pEqn.H     |  2 +-
 .../potentialFreeSurfaceFoam/pEqn.H           |  2 +-
 .../potentialFreeSurfaceDyMFoam/pEqn.H        |  2 +-
 .../reactingMultiphaseEulerFoam/YEqns.H       |  2 +-
 .../reactingMultiphaseEulerFoam/pU/pEqn.H     |  6 +--
 .../reactingMultiphaseEulerFoam/pUf/pEqn.H    |  6 +--
 .../reactingTwoPhaseEulerFoam/YEqns.H         |  4 +-
 .../reactingTwoPhaseEulerFoam/pUf/pEqn.H      |  6 +--
 .../multiphase/twoLiquidMixingFoam/pEqn.H     |  2 +-
 .../multiphase/twoPhaseEulerFoam/pU/pEqn.H    |  2 +-
 .../multiphase/twoPhaseEulerFoam/pUf/pEqn.H   |  2 +-
 src/finiteArea/faMatrices/faMatrix/faMatrix.C | 54 +++++++++++++++++--
 src/finiteArea/faMatrices/faMatrix/faMatrix.H | 44 ++++++++++++---
 .../faMatrices/faMatrix/faMatrixSolve.C       | 13 +++--
 .../cfdTools/general/CorrectPhi/CorrectPhi.C  |  4 +-
 .../fvMatrices/fvMatrix/fvMatrix.C            | 43 +++++++++++++--
 .../fvMatrices/fvMatrix/fvMatrix.H            | 42 ++++++++++++---
 .../fvMatrices/fvMatrix/fvMatrixSolve.C       | 16 +++++-
 .../solvers/energyTransport/energyTransport.C |  4 +-
 .../solvers/scalarTransport/scalarTransport.C |  6 +--
 ...lacementComponentLaplacianFvMotionSolver.C |  2 +-
 ...velocityComponentLaplacianFvMotionSolver.C |  2 +-
 .../displacementSBRStressFvMotionSolver.C     |  2 +-
 .../displacementLaplacianFvMotionSolver.C     |  2 +-
 ...dBodyDisplacementLaplacianFvMotionSolver.C |  2 +-
 .../surfaceAlignedSBRStressFvMotionSolver.C   |  2 +-
 .../velocityLaplacianFvMotionSolver.C         |  2 +-
 .../MultiComponentPhaseModel.C                |  2 +-
 93 files changed, 292 insertions(+), 151 deletions(-)

diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index 2a00f53b457..b63a2cf7e72 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -118,7 +118,7 @@ int main(int argc, char *argv[])
                 fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
             );
 
-            pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+            pEqn.solve(p.select(piso.finalInnerIter()));
 
             phi = phiHbyA - pEqn.flux();
 
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index c7308f8b122..f37e8f5171b 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -27,7 +27,7 @@ if (pimple.transonic())
             betav*fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -57,7 +57,7 @@ else
             betav*fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H b/applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H
index 34a822a0a28..cd0565c2325 100644
--- a/applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H
@@ -30,7 +30,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -66,7 +66,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/combustion/XiFoam/XiEngineFoam/pEqn.H b/applications/solvers/combustion/XiFoam/XiEngineFoam/pEqn.H
index fedaea71067..a1f5317d39c 100644
--- a/applications/solvers/combustion/XiFoam/XiEngineFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/XiEngineFoam/pEqn.H
@@ -35,7 +35,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -71,7 +71,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index afba129157e..f21ecd8039a 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -29,7 +29,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -64,7 +64,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/combustion/chemFoam/YEqn.H b/applications/solvers/combustion/chemFoam/YEqn.H
index f31ff16496f..8cfb8cfeb5c 100644
--- a/applications/solvers/combustion/chemFoam/YEqn.H
+++ b/applications/solvers/combustion/chemFoam/YEqn.H
@@ -6,7 +6,7 @@
         solve
         (
             fvm::ddt(rho, Yi) - chemistry.RR(specieI),
-            mesh.solver("Yi")
+            "Yi"
         );
     }
 }
diff --git a/applications/solvers/combustion/fireFoam/YEEqn.H b/applications/solvers/combustion/fireFoam/YEEqn.H
index adc849bc2da..e4f6afa17ab 100644
--- a/applications/solvers/combustion/fireFoam/YEEqn.H
+++ b/applications/solvers/combustion/fireFoam/YEEqn.H
@@ -35,7 +35,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YiEqn);
 
-            YiEqn.solve(mesh.solver("Yi"));
+            YiEqn.solve("Yi");
 
             fvOptions.correct(Yi);
 
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 19cafee7cea..9061b514b17 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -36,7 +36,7 @@ while (pimple.correctNonOrthogonal())
       + fvOptions(psi, p_rgh, rho.name())
     );
 
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/combustion/reactingFoam/YEqn.H b/applications/solvers/combustion/reactingFoam/YEqn.H
index ea8681628a1..5d3ed162072 100644
--- a/applications/solvers/combustion/reactingFoam/YEqn.H
+++ b/applications/solvers/combustion/reactingFoam/YEqn.H
@@ -34,7 +34,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YiEqn);
 
-            YiEqn.solve(mesh.solver("Yi"));
+            YiEqn.solve("Yi");
 
             fvOptions.correct(Yi);
 
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index f4ff96931bd..b7240face8a 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -37,7 +37,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -72,7 +72,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/combustion/reactingFoam/pcEqn.H b/applications/solvers/combustion/reactingFoam/pcEqn.H
index 5ca9b4adca1..c9b0171faa6 100644
--- a/applications/solvers/combustion/reactingFoam/pcEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pcEqn.H
@@ -49,7 +49,7 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -89,7 +89,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
index d5feab6e538..03e0623dee8 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@@ -41,7 +41,7 @@ while (pimple.correctNonOrthogonal())
       - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/compressible/rhoPimpleAdiabaticFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleAdiabaticFoam/pEqn.H
index 27398f42bd3..13f2db4ffa0 100644
--- a/applications/solvers/compressible/rhoPimpleAdiabaticFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleAdiabaticFoam/pEqn.H
@@ -63,7 +63,7 @@
                 fvOptions(psi, p, rho.name())
             );
 
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+            pEqn.solve(p.select(pimple.finalInnerIter()));
 
             // Rhie & Chow interpolation (part 2)
             if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H
index a332be82849..190763d2793 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H
@@ -72,7 +72,7 @@ if (mesh.changing())
             divrhoU()
         );
 
-        pcorrEqn.solve(mesh.solver(pcorr.select(pimple.finalInnerIter())));
+        pcorrEqn.solve(pcorr.select(pimple.finalInnerIter()));
         //Bypass virtual layer
         //mesh.fvMesh::solve(pcorrEqn, d);
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
index 227610a01ba..4cb64661b4d 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
@@ -56,7 +56,7 @@ if (pimple.transonic())
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -78,7 +78,7 @@ else
     {
         fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 8e2c3d836d3..9942c8eddd3 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -54,7 +54,7 @@ if (pimple.transonic())
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -76,7 +76,7 @@ else
     {
         fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index 86812465538..f2861b769c3 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -65,7 +65,7 @@ if (pimple.transonic())
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -90,7 +90,7 @@ else
     {
         fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
index eb3ba594823..478d8ebee5a 100644
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@@ -127,7 +127,7 @@ int main(int argc, char *argv[])
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
-                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+                    pEqn.solve(p.select(piso.finalInnerIter()));
 
                     if (piso.finalNonOrthogonalIter())
                     {
@@ -167,7 +167,7 @@ int main(int argc, char *argv[])
                     fvm::laplacian(rABf, pB) == fvc::div(phiB)
                 );
 
-                pBEqn.solve(mesh.solver(pB.select(bpiso.finalInnerIter())));
+                pBEqn.solve(pB.select(bpiso.finalInnerIter()));
 
                 if (bpiso.finalNonOrthogonalIter())
                 {
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index d6f04136dd1..3832ad47608 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@@ -27,7 +27,7 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H
index a206c943cdb..708a7254d2a 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/correctPhi.H
@@ -72,7 +72,7 @@ if (mesh.changing())
             divrhoU()
         );
 
-        //pcorrEqn.solve(mesh.solver(pcorr.select(pimple.finalInnerIter())));
+        //pcorrEqn.solve(pcorr.select(pimple.finalInnerIter()));
         //Bypass virtual layer
         const dictionary& d = mesh.solver
             (
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H
index 06b6f99fc9b..cd5f97790a2 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H
@@ -50,7 +50,7 @@ while (pimple.correctNonOrthogonal())
       - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 733e97028d7..58a84e38b44 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -52,7 +52,7 @@ while (pimple.correctNonOrthogonal())
         compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
     );
 
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/solid/solveSolid.H
index 165345f6ce2..f63e3335ea3 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/solid/solveSolid.H
@@ -22,7 +22,7 @@ if (finalIter)
 
         fvOptions.constrain(hEqn);
 
-        hEqn.solve(mesh.solver(h.select(finalIter)));
+        hEqn.solve(h.select(finalIter));
 
         fvOptions.correct(h);
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
index 39d238abcee..5939182d1f4 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
@@ -33,7 +33,7 @@
     }
     else
     {
-        EEqn.solve(mesh.solver(he.select(finalIter)));
+        EEqn.solve(he.select(finalIter));
         fvOptions.correct(he);
 
         thermo.correct();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
index 65654eeefed..33eb61e69ae 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
@@ -28,7 +28,7 @@
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
             ),
-            mesh.solver(U.select(finalIter))
+            U.select(finalIter)
         );
 
         fvOptions.correct(U);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
index b42c7986abb..1c17b76eeb7 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
@@ -44,7 +44,7 @@ if (Y.size())
 
             fvOptions.constrain(YiEqn);
 
-            YiEqn.solve(mesh.solver("Yi"));
+            YiEqn.solve("Yi");
 
             fvOptions.correct(Yi);
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 206c32eeb2b..8b3f629fa71 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -52,15 +52,12 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
         p_rghEqn.solve
         (
-            mesh.solver
+            p_rgh.select
             (
-                p_rgh.select
                 (
-                    (
-                        oCorr == nOuterCorr-1
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
+                    oCorr == nOuterCorr-1
+                 && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
                 )
             )
         );
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
index 2d11bfd21df..466e3a8a4d3 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
@@ -28,7 +28,7 @@
             mesh.data().setFinalIteration(true);
         }
 
-        hEqn.solve(mesh.solver(h.select(finalIter)));
+        hEqn.solve(h.select(finalIter));
 
         fvOptions.correct(h);
 
diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C
index 72bb4ec51ad..f03ccf80767 100644
--- a/applications/solvers/incompressible/icoFoam/icoFoam.C
+++ b/applications/solvers/incompressible/icoFoam/icoFoam.C
@@ -141,7 +141,7 @@ int main(int argc, char *argv[])
 
                 pEqn.setReference(pRefCell, pRefValue);
 
-                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+                pEqn.solve(p.select(piso.finalInnerIter()));
 
                 if (piso.finalNonOrthogonalIter())
                 {
diff --git a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
index c84a168c2d7..85041130ff6 100644
--- a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
+++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
@@ -114,7 +114,7 @@ int main(int argc, char *argv[])
 
                 pEqn.setReference(pRefCell, pRefValue);
 
-                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+                pEqn.solve(p.select(piso.finalInnerIter()));
 
                 if (piso.finalNonOrthogonalIter())
                 {
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
index 3cb0580ba46..13120b50867 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
@@ -40,7 +40,7 @@ while (pimple.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+    pEqn.solve(p.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H
index 5f8971e2a0b..76381446e00 100644
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H
@@ -32,7 +32,7 @@ while (pimple.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+    pEqn.solve(p.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index cca90fdc9ac..6908b5259ae 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -45,7 +45,7 @@ while (pimple.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+    pEqn.solve(p.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/incompressible/pisoFoam/pEqn.H b/applications/solvers/incompressible/pisoFoam/pEqn.H
index 866d84e9457..0e0e32a2e69 100644
--- a/applications/solvers/incompressible/pisoFoam/pEqn.H
+++ b/applications/solvers/incompressible/pisoFoam/pEqn.H
@@ -26,7 +26,7 @@ while (piso.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+    pEqn.solve(p.select(piso.finalInnerIter()));
 
     if (piso.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
index e303c8a747a..ebd0c4fc2c9 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
@@ -135,7 +135,7 @@ int main(int argc, char *argv[])
                       - fvm::laplacian(ghrAUf, h)
                     );
 
-                    hEqn.solve(mesh.solver(h.select(pimple.finalInnerIter())));
+                    hEqn.solve(h.select(pimple.finalInnerIter()));
 
                     if (pimple.finalNonOrthogonalIter())
                     {
diff --git a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H
index dd4d382befa..aa345ca122a 100644
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H
@@ -34,7 +34,7 @@
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/lagrangian/DPMFoam/pEqn.H b/applications/solvers/lagrangian/DPMFoam/pEqn.H
index ea0bd1c101d..304d69b723a 100644
--- a/applications/solvers/lagrangian/DPMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/DPMFoam/pEqn.H
@@ -32,7 +32,7 @@
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
index f3cbb702a7a..59bb2aae167 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
@@ -36,7 +36,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YiEqn);
 
-            YiEqn.solve(mesh.solver("Yi"));
+            YiEqn.solve("Yi");
 
             fvOptions.correct(Yi);
 
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 0182e453ba1..0267c10adcc 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -33,7 +33,7 @@ if (pimple.transonic())
           + fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -69,7 +69,7 @@ else
           + fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/lagrangian/kinematicParcelFoam/pEqn.H b/applications/solvers/lagrangian/kinematicParcelFoam/pEqn.H
index d32c1dc45c5..79725f23851 100644
--- a/applications/solvers/lagrangian/kinematicParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/pEqn.H
@@ -34,7 +34,7 @@ while (pimple.correctNonOrthogonal())
 
     pEqn.setReference(pRefCell, pRefValue);
 
-    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+    pEqn.solve(p.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
index ecbf1208873..c0179eadb86 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
@@ -37,7 +37,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YEqn);
 
-            YEqn.solve(mesh.solver("Yi"));
+            YEqn.solve("Yi");
 
             fvOptions.correct(Yi);
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 81f5e5447ad..30ad3e4b7c0 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -47,7 +47,7 @@ while (pimple.correctNonOrthogonal())
       - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
index db2d6c9b9c4..925c55df325 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
@@ -34,7 +34,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YEqn);
 
-            YEqn.solve(mesh.solver("Yi"));
+            YEqn.solve("Yi");
 
             fvOptions.correct(Yi);
 
diff --git a/applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H b/applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H
index a6f70bf4f2a..27c1fd46db5 100644
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H
@@ -34,7 +34,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YEqn);
 
-            YEqn.solve(mesh.solver("Yi"));
+            YEqn.solve("Yi");
 
             fvOptions.correct(Yi);
 
diff --git a/applications/solvers/lagrangian/sprayFoam/YEqn.H b/applications/solvers/lagrangian/sprayFoam/YEqn.H
index e80f884202e..f33c64497dd 100644
--- a/applications/solvers/lagrangian/sprayFoam/YEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/YEqn.H
@@ -35,7 +35,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             fvOptions.constrain(YEqn);
 
-            YEqn.solve(mesh.solver("Yi"));
+            YEqn.solve("Yi");
 
             fvOptions.correct(Yi);
 
diff --git a/applications/solvers/lagrangian/sprayFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/pEqn.H
index 838044ba733..198bf00d52c 100644
--- a/applications/solvers/lagrangian/sprayFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/pEqn.H
@@ -41,7 +41,7 @@ if (pimple.transonic())
           + fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -77,7 +77,7 @@ else
           + fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
index 54a8ae0322b..a7850a5e903 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
@@ -42,7 +42,7 @@ if (pimple.transonic())
           + fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -77,7 +77,7 @@ else
           + fvOptions(psi, p, rho.name())
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/MPPICInterFoam/pEqn.H b/applications/solvers/multiphase/MPPICInterFoam/pEqn.H
index 2e1330b3b8b..a455c7ac9a7 100644
--- a/applications/solvers/multiphase/MPPICInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/MPPICInterFoam/pEqn.H
@@ -38,7 +38,7 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index f5b1d8ef01e..568d2ef4189 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -36,7 +36,7 @@
           - fvm::laplacian(rhorAUf, p)
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index 2cb7b0fc5ad..2e1b449a7e0 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -33,7 +33,7 @@
           - fvm::laplacian(rhorAUf, p)
         );
 
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+        pEqn.solve(p.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index d4c0f41aab8..d56a53bf914 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -99,7 +99,7 @@
               + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
             )
           + p_rghEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            p_rgh.select(pimple.finalInnerIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H
index 5cdd9737d0c..0069c53995e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H
@@ -108,7 +108,7 @@
         solve
         (
             p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            p_rgh.select(pimple.finalInnerIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/pEqn.H
index f3fc5b546c4..0aad99d4769 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/pEqn.H
@@ -124,7 +124,7 @@
          solve
          (
              p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
-             mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+             p_rgh.select(pimple.finalInnerIter())
          );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H
index b69ce2e75d0..5b704a0f4bd 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H
@@ -110,7 +110,7 @@
         solve
         (
             p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            p_rgh.select(pimple.finalInnerIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 5854748e99b..381e45e75bc 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -108,7 +108,7 @@
         solve
         (
             p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            p_rgh.select(pimple.finalInnerIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
index 479bc2cd34f..107fe3c0cca 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@@ -80,7 +80,7 @@
         (
             p_rghEqnComp
           + p_rghEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            p_rgh.select(pimple.finalInnerIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
index 115f9c67388..b38144700cc 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
@@ -54,7 +54,7 @@
           - fvm::laplacian(turbulence->nut(), alpha1)
         );
 
-        alpha1Eqn.solve(mesh.solver("alpha1Diffusion"));
+        alpha1Eqn.solve("alpha1Diffusion");
 
         alphaPhi += alpha1Eqn.flux();
         alpha2 = 1.0 - alpha1;
diff --git a/applications/solvers/multiphase/driftFluxFoam/pEqn.H b/applications/solvers/multiphase/driftFluxFoam/pEqn.H
index 7a2be2efd05..b627d519139 100644
--- a/applications/solvers/multiphase/driftFluxFoam/pEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/pEqn.H
@@ -32,7 +32,7 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
index bf2c86525b0..8906a58a616 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
@@ -56,7 +56,7 @@
 
         p_rghEqn.setReference(pRefCell, pRefValue);
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/pEqn.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/pEqn.H
index d291705f495..bdcf35449d2 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/pEqn.H
@@ -39,7 +39,7 @@
 
         p_rghEqn.setReference(pRefCell, pRefValue);
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H
index 3b4df83409e..3a401339454 100644
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H
@@ -40,7 +40,7 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index e443f269848..b4e3658a3ef 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -47,7 +47,7 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
index 8d992e80b0e..152d4f2111f 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
@@ -44,7 +44,7 @@
 
         p_rghEqn.setReference(pRefCell, pRefValue);
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H
index 850888b9f4f..d5904e179dd 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H
@@ -45,7 +45,7 @@
         //p_rghEqn.setReference(pRefCell, pRefValue);
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index c08deda4038..44ad5d8ec07 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -38,7 +38,7 @@
 
         p_rghEqn.setReference(pRefCell, pRefValue);
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index c0790b96989..02ae8778b9f 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -223,7 +223,7 @@
             //   + (alpha2/rho2)*pEqnComp2()
             // ) +
             pEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            p_rgh.select(pimple.finalInnerIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
index 0e259d88092..6e4db195d6f 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
@@ -30,7 +30,7 @@ while (pimple.correctNonOrthogonal())
 
     p_ghEqn.setReference(p_ghRefCell, p_ghRefValue);
 
-    p_ghEqn.solve(mesh.solver(p_gh.select(pimple.finalInnerIter())));
+    p_ghEqn.solve(p_gh.select(pimple.finalInnerIter()));
 
     if (pimple.finalNonOrthogonalIter())
     {
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
index 9f7da03a11c..d89b6462b52 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
@@ -37,7 +37,7 @@
 
         p_ghEqn.setReference(p_ghRefCell, p_ghRefValue);
 
-        p_ghEqn.solve(mesh.solver(p_gh.select(pimple.finalInnerIter())));
+        p_ghEqn.solve(p_gh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/YEqns.H b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/YEqns.H
index 73eb99dd8c7..75e071490fd 100644
--- a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/YEqns.H
+++ b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/YEqns.H
@@ -24,7 +24,7 @@
             );
 
             YiEqn.relax();
-            YiEqn.solve(mesh.solver("Yi"));
+            YiEqn.solve("Yi");
         }
     }
 }
diff --git a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/pU/pEqn.H
index 5f65eb51297..6d212446c63 100644
--- a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -370,11 +370,7 @@ while (pimple.correct())
                 }
             }
 
-            solve
-            (
-                pEqn,
-                mesh.solver(p_rgh.select(pimple.finalInnerIter()))
-            );
+            pEqn.solve(p_rgh.select(pimple.finalInnerIter()));
         }
 
         // Correct fluxes and velocities on last non-orthogonal iteration
diff --git a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/pUf/pEqn.H
index f2e3d510675..08d97031b04 100644
--- a/applications/solvers/multiphase/reactingMultiphaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingMultiphaseEulerFoam/pUf/pEqn.H
@@ -358,11 +358,7 @@ while (pimple.correct())
                 }
             }
 
-            solve
-            (
-                pEqn,
-                mesh.solver(p_rgh.select(pimple.finalInnerIter()))
-            );
+            pEqn.solve(p_rgh.select(pimple.finalInnerIter()));
         }
 
         // Correct fluxes and velocities on last non-orthogonal iteration
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/YEqns.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/YEqns.H
index ba1ee2328d7..6df93827c4e 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/YEqns.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/YEqns.H
@@ -20,7 +20,7 @@
             );
 
             Y1iEqn.relax();
-            Y1iEqn.solve(mesh.solver("Yi"));
+            Y1iEqn.solve("Yi");
         }
     }
 
@@ -39,7 +39,7 @@
             );
 
             Y2iEqn.relax();
-            Y2iEqn.solve(mesh.solver("Yi"));
+            Y2iEqn.solve("Yi");
         }
     }
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 0be323eb88d..929f228a7af 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -341,11 +341,7 @@ while (pimple.correct())
                 pEqn += pEqnComp2();
             }
 
-            solve
-            (
-                pEqn,
-                mesh.solver(p_rgh.select(pimple.finalInnerIter()))
-            );
+            pEqn.solve(p_rgh.select(pimple.finalInnerIter()));
         }
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index d801fb8fa5b..24fa24adfa0 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@@ -29,7 +29,7 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
 
         if (pimple.finalNonOrthogonalIter())
         {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
index 71d75bd26e0..02c401708c1 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
@@ -312,7 +312,7 @@ while (pimple.correct())
         solve
         (
             pEqnComp1() + pEqnComp2() + pEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            p_rgh.select(pimple.finalInnerIter())
         );
 
         // Correct fluxes and velocities on last non-orthogonal iteration
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
index 0e1b3b7ff6f..e4c30664fb4 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
@@ -291,7 +291,7 @@ while (pimple.correct())
         solve
         (
             pEqnComp1() + pEqnComp2() + pEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            p_rgh.select(pimple.finalInnerIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.C b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
index 9dc7a9f6430..55aa863bf71 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrix.C
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
@@ -743,6 +743,26 @@ Foam::faMatrix<Type>::flux() const
 }
 
 
+template<class Type>
+const Foam::dictionary& Foam::faMatrix<Type>::solverDict
+(
+    const word& name
+) const
+{
+    return psi_.mesh().solverDict(name);
+}
+
+
+template<class Type>
+const Foam::dictionary& Foam::faMatrix<Type>::solverDict() const
+{
+    return psi_.mesh().solverDict
+    (
+        psi_.name()
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -1110,11 +1130,11 @@ void Foam::checkMethod
 template<class Type>
 Foam::SolverPerformance<Type> Foam::solve
 (
-    faMatrix<Type>& fam,
+    faMatrix<Type>& mat,
     const dictionary& solverControls
 )
 {
-    return fam.solve(solverControls);
+    return mat.solve(solverControls);
 }
 
 
@@ -1134,9 +1154,35 @@ Foam::SolverPerformance<Type> Foam::solve
 
 
 template<class Type>
-Foam::SolverPerformance<Type> Foam::solve(faMatrix<Type>& fam)
+Foam::SolverPerformance<Type> Foam::solve
+(
+    faMatrix<Type>& mat,
+    const word& name
+)
+{
+    return mat.solve(name);
+}
+
+
+template<class Type>
+Foam::SolverPerformance<Type> Foam::solve
+(
+    const tmp<faMatrix<Type>>& tmat,
+    const word& name
+)
+{
+    SolverPerformance<Type> solverPerf(tmat.constCast().solve(name));
+
+    tmat.clear();
+
+    return solverPerf;
+}
+
+
+template<class Type>
+Foam::SolverPerformance<Type> Foam::solve(faMatrix<Type>& mat)
 {
-    return fam.solve();
+    return mat.solve();
 }
 
 
diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.H b/src/finiteArea/faMatrices/faMatrix/faMatrix.H
index d8bc359f139..47b19d76048 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrix.H
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -448,6 +448,10 @@ public:
             //  Use the given solver controls
             SolverPerformance<Type> solve(const dictionary&);
 
+            //- Solve returning the solution statistics.
+            //  Uses \p name solver controls from faSolution
+            SolverPerformance<Type> solve(const word& name);
+
             //- Solve returning the solution statistics.
             //  Solver controls read from faSolution
             SolverPerformance<Type> solve();
@@ -467,6 +471,12 @@ public:
             //- Return the face-flux field from the matrix
             tmp<GeometricField<Type, faePatchField, edgeMesh>> flux() const;
 
+            //- Return the solver dictionary (from faSolution) for \p name
+            const dictionary& solverDict(const word& name) const;
+
+            //- Return the solver dictionary for psi
+            const dictionary& solverDict() const;
+
 
     // Member Operators
 
@@ -578,27 +588,45 @@ void checkMethod
 
 
 //- Solve returning the solution statistics given convergence tolerance
-// Use the given solver controls
+//  Use the given solver controls
 template<class Type>
-SolverPerformance<Type> solve(faMatrix<Type>&, const dictionary&);
+SolverPerformance<Type> solve
+(
+    faMatrix<Type>&,
+    const dictionary& solverControls
+);
+
+//- Solve returning the solution statistics given convergence tolerance,
+//- deleting temporary matrix after solution.
+//  Use the given solver controls
+template<class Type>
+SolverPerformance<Type> solve
+(
+    const tmp<faMatrix<Type>>&,
+    const dictionary& solverControls
+);
 
 
+//- Solve returning the solution statistics given convergence tolerance
+//  Uses \p name solver controls from faSolution
+template<class Type>
+SolverPerformance<Type> solve(faMatrix<Type>&, const word& name);
+
 //- Solve returning the solution statistics given convergence tolerance,
 //- deleting temporary matrix after solution.
-// Use the given solver controls
+//  Uses \p name solver controls from faSolution
 template<class Type>
-SolverPerformance<Type> solve(const tmp<faMatrix<Type>>&, const dictionary&);
+SolverPerformance<Type> solve(const tmp<faMatrix<Type>>&, const word& name);
 
 
 //- Solve returning the solution statistics given convergence tolerance
-//  Solver controls read faSolution
+//  Uses solver controls from faSolution
 template<class Type>
 SolverPerformance<Type> solve(faMatrix<Type>&);
 
-
 //- Solve returning the solution statistics given convergence tolerance,
 //- deleting temporary matrix after solution.
-//  Solver controls read faSolution
+//  Uses solver controls from faSolution
 template<class Type>
 SolverPerformance<Type> solve(const tmp<faMatrix<Type>>&);
 
diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C b/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C
index 4de12f14cb5..87d5ef0e03d 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -171,14 +171,21 @@ Foam::SolverPerformance<Type> Foam::faMatrix<Type>::solve
 template<class Type>
 Foam::SolverPerformance<Type> Foam::faMatrix<Type>::faSolver::solve()
 {
-    return solve(faMat_.psi().mesh().solverDict(faMat_.psi().name()));
+    return solve(faMat_.solverDict());
+}
+
+
+template<class Type>
+Foam::SolverPerformance<Type> Foam::faMatrix<Type>::solve(const word& name)
+{
+    return this->solve(solverDict(name));
 }
 
 
 template<class Type>
 Foam::SolverPerformance<Type> Foam::faMatrix<Type>::solve()
 {
-    return solve(this->psi().mesh().solverDict(this->psi().name()));
+    return this->solve(solverDict());
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
index c92d10b68c0..06dcfb73077 100644
--- a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
+++ b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
@@ -107,7 +107,7 @@ void Foam::CorrectPhi
 
         pcorrEqn.solve
         (
-            mesh.solver(pcorr.select(pimple.finalNonOrthogonalIter()))
+            pcorr.select(pimple.finalNonOrthogonalIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
@@ -182,7 +182,7 @@ void Foam::CorrectPhi
 
         pcorrEqn.solve
         (
-            mesh.solver(pcorr.select(pimple.finalNonOrthogonalIter()))
+            pcorr.select(pimple.finalNonOrthogonalIter())
         );
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 59a95ccc3d7..6b15f83645a 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -1528,6 +1528,16 @@ flux() const
 }
 
 
+template<class Type>
+const Foam::dictionary& Foam::fvMatrix<Type>::solverDict
+(
+    const word& name
+) const
+{
+    return psi_.mesh().solverDict(name);
+}
+
+
 template<class Type>
 const Foam::dictionary& Foam::fvMatrix<Type>::solverDict() const
 {
@@ -1918,11 +1928,11 @@ void Foam::checkMethod
 template<class Type>
 Foam::SolverPerformance<Type> Foam::solve
 (
-    fvMatrix<Type>& fvm,
+    fvMatrix<Type>& mat,
     const dictionary& solverControls
 )
 {
-    return fvm.solve(solverControls);
+    return mat.solve(solverControls);
 }
 
 template<class Type>
@@ -1941,9 +1951,34 @@ Foam::SolverPerformance<Type> Foam::solve
 
 
 template<class Type>
-Foam::SolverPerformance<Type> Foam::solve(fvMatrix<Type>& fvm)
+Foam::SolverPerformance<Type> Foam::solve
+(
+    fvMatrix<Type>& mat,
+    const word& name
+)
+{
+    return mat.solve(name);
+}
+
+template<class Type>
+Foam::SolverPerformance<Type> Foam::solve
+(
+    const tmp<fvMatrix<Type>>& tmat,
+    const word& name
+)
+{
+    SolverPerformance<Type> solverPerf(tmat.constCast().solve(name));
+
+    tmat.clear();
+
+    return solverPerf;
+}
+
+
+template<class Type>
+Foam::SolverPerformance<Type> Foam::solve(fvMatrix<Type>& mat)
 {
-    return fvm.solve();
+    return mat.solve();
 }
 
 template<class Type>
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
index b6d2528aef5..c00236f7c9f 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -603,6 +603,10 @@ public:
             //  Use the given solver controls
             SolverPerformance<Type> solveSegregatedOrCoupled(const dictionary&);
 
+            //- Solve segregated or coupled returning the solution statistics.
+            //  Solver controls read from fvSolution
+            SolverPerformance<Type> solveSegregatedOrCoupled();
+
             //- Solve segregated returning the solution statistics.
             //  Use the given solver controls
             SolverPerformance<Type> solveSegregated(const dictionary&);
@@ -615,6 +619,10 @@ public:
             //  Use the given solver controls
             SolverPerformance<Type> solve(const dictionary&);
 
+            //- Solve returning the solution statistics.
+            //  Uses \p name solver controls from fvSolution
+            SolverPerformance<Type> solve(const word& name);
+
             //- Solve returning the solution statistics.
             //  Solver controls read from fvSolution
             SolverPerformance<Type> solve();
@@ -641,7 +649,11 @@ public:
             tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
                 flux() const;
 
-            //- Return the solver dictionary taking into account finalIteration
+            //- Return the solver dictionary (from fvSolution) for \p name
+            const dictionary& solverDict(const word& name) const;
+
+            //- Return the solver dictionary for psi,
+            //- taking into account finalIteration
             const dictionary& solverDict() const;
 
 
@@ -757,8 +769,11 @@ void checkMethod
 //- Solve returning the solution statistics given convergence tolerance
 //  Use the given solver controls
 template<class Type>
-SolverPerformance<Type> solve(fvMatrix<Type>&, const dictionary&);
-
+SolverPerformance<Type> solve
+(
+    fvMatrix<Type>&,
+    const dictionary& solverControls
+);
 
 //- Solve returning the solution statistics given convergence tolerance,
 //- deleting temporary matrix after solution.
@@ -767,19 +782,30 @@ template<class Type>
 SolverPerformance<Type> solve
 (
     const tmp<fvMatrix<Type>>&,
-    const dictionary&
+    const dictionary& solverControls
 );
 
 
 //- Solve returning the solution statistics given convergence tolerance
-//  Solver controls read fvSolution
+//  Uses \p name solver controls from fvSolution
 template<class Type>
-SolverPerformance<Type> solve(fvMatrix<Type>&);
+SolverPerformance<Type> solve(fvMatrix<Type>&, const word& name);
+
+//- Solve returning the solution statistics given convergence tolerance,
+//- deleting temporary matrix after solution.
+//  Uses \p name solver controls from fvSolution
+template<class Type>
+SolverPerformance<Type> solve(const tmp<fvMatrix<Type>>&, const word& name);
 
 
+//- Solve returning the solution statistics given convergence tolerance
+//  Uses solver controls from fvSolution
+template<class Type>
+SolverPerformance<Type> solve(fvMatrix<Type>&);
+
 //- Solve returning the solution statistics given convergence tolerance,
 //- deleting temporary matrix after solution.
-//  Solver controls read fvSolution
+//  Uses solver controls from fvSolution
 template<class Type>
 SolverPerformance<Type> solve(const tmp<fvMatrix<Type>>&);
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
index 6522f4d0cf7..8fbd1150031 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -314,6 +314,13 @@ Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solveCoupled
 }
 
 
+template<class Type>
+Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solveSegregatedOrCoupled()
+{
+    return this->solveSegregatedOrCoupled(solverDict());
+}
+
+
 template<class Type>
 Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solve
 (
@@ -339,6 +346,13 @@ Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::fvSolver::solve()
 }
 
 
+template<class Type>
+Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solve(const word& name)
+{
+    return this->solve(solverDict(name));
+}
+
+
 template<class Type>
 Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solve()
 {
diff --git a/src/functionObjects/solvers/energyTransport/energyTransport.C b/src/functionObjects/solvers/energyTransport/energyTransport.C
index 477c09d4768..ea5b022f117 100644
--- a/src/functionObjects/solvers/energyTransport/energyTransport.C
+++ b/src/functionObjects/solvers/energyTransport/energyTransport.C
@@ -395,7 +395,7 @@ bool Foam::functionObjects::energyTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(mesh_.solverDict(schemesField_));
+            sEqn.solve(schemesField_);
         }
     }
     else if (phi.dimensions() == dimVolume/dimTime)
@@ -433,7 +433,7 @@ bool Foam::functionObjects::energyTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(mesh_.solverDict(schemesField_));
+            sEqn.solve(schemesField_);
         }
     }
     else
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.C b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
index f486afeb124..ec369782760 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.C
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
@@ -295,7 +295,7 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             sEqn.relax(relaxCoeff);
             fvOptions_.constrain(sEqn);
-            sEqn.solve(mesh_.solverDict(schemesField_));
+            sEqn.solve(schemesField_);
 
             tTPhiUD = sEqn.flux();
         }
@@ -333,7 +333,7 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(mesh_.solverDict(schemesField_));
+            sEqn.solve(schemesField_);
         }
     }
     else if (phi.dimensions() == dimVolume/dimTime)
@@ -353,7 +353,7 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(mesh_.solverDict(schemesField_));
+            sEqn.solve(schemesField_);
         }
     }
     else
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
index 4d3c735d21e..c0656e0516e 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
@@ -242,7 +242,7 @@ void Foam::displacementComponentLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-    TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    TEqn.solveSegregatedOrCoupled();
     fvOptions.correct(cellDisplacement_);
 }
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
index d9c53c4d7b3..68d95d1f981 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
@@ -146,7 +146,7 @@ void Foam::velocityComponentLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-    TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    TEqn.solveSegregatedOrCoupled();
     fvOptions.correct(cellMotionU_);
 }
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
index a50f0e1b898..c220b6635c8 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
@@ -242,7 +242,7 @@ void Foam::displacementSBRStressFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-    TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    TEqn.solveSegregatedOrCoupled();
     fvOptions.correct(cellDisplacement_);
 }
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
index f7e53f1883e..e29a6e48652 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
@@ -351,7 +351,7 @@ void Foam::displacementLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-    TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    TEqn.solveSegregatedOrCoupled();
     fvOptions.correct(cellDisplacement_);
 }
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/solidBodyDisplacementLaplacian/solidBodyDisplacementLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/solidBodyDisplacementLaplacian/solidBodyDisplacementLaplacianFvMotionSolver.C
index 7861c2215e4..79ddc0d9bed 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/solidBodyDisplacementLaplacian/solidBodyDisplacementLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/solidBodyDisplacementLaplacian/solidBodyDisplacementLaplacianFvMotionSolver.C
@@ -359,7 +359,7 @@ void Foam::solidBodyDisplacementLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-    TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    TEqn.solveSegregatedOrCoupled();
     fvOptions.correct(cellDisplacement_);
 }
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
index ddc2da8c2f2..0d6c19e100c 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
@@ -423,7 +423,7 @@ void Foam::surfaceAlignedSBRStressFvMotionSolver::solve()
         fvOptions.constrain(DEqn);
 
         // Note: solve uncoupled
-        DEqn.solveSegregatedOrCoupled(DEqn.solverDict());
+        DEqn.solveSegregatedOrCoupled();
 
         fvOptions.correct(cellDisp);
     }
diff --git a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
index 35ba1f3a5ec..2459be05e63 100644
--- a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
@@ -146,7 +146,7 @@ void Foam::velocityLaplacianFvMotionSolver::solve()
         );
 
         fvOptions.constrain(UEqn);
-        UEqn.solveSegregatedOrCoupled(UEqn.solverDict());
+        UEqn.solveSegregatedOrCoupled();
         fvOptions.correct(cellMotionU_);
     }
 }
diff --git a/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
index 3fd085f7e86..e6ab27f2680 100644
--- a/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
+++ b/src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
@@ -420,7 +420,7 @@ void Foam::MultiComponentPhaseModel<BasePhaseModel, phaseThermo>::solveYi
                     )
                 );
 
-                YiDiffEqn.solve(mesh.solver("diffusion" + Yi.name()));
+                YiDiffEqn.solve("diffusion" + Yi.name());
             }
 
             Yt += Yi;
-- 
GitLab