diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index 2a00f53b457575c28f97579cb1a04b4c9de49302..b63a2cf7e72384950cfdadb133d5e1a768e02108 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 c7308f8b122d90a98222f20ea6b7c9403c9cac0a..f37e8f5171bc770bdeb94fb7d7ab3a7904045e69 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 34a822a0a282e540fbd908758d50b8bafc8fddfb..cd0565c23252ec7517e560c34e0bae1cbcb4b18a 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 fedaea7106779c6da777cc654fb3fe3770d2f547..a1f5317d39cce7bcbf88396b6cf687b075dd1ac1 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 afba129157eabffbd7a24e257dda287d533721a3..f21ecd8039a7bd7221a941d64026abe3605326f0 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 f31ff16496f56f646fda61532712b9ce33762496..8cfb8cfeb5cb8b5b2095855efd8c8d22a17d11bf 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 adc849bc2da6e79e0545edd5293edbb8d01eae99..e4f6afa17ab1ac609bea4ac0415f1d9eaddf5d53 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 19cafee7cea7b22e998dc92130bbd48f287864f1..9061b514b17c451cf2899502f66b75362c116cea 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 ea8681628a12ca529766c6c1177d5bd2f6c50072..5d3ed162072576dd37574a9d754d78ecc68b5255 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 f4ff96931bdd8df2e8fd2553c03bfe8213c9059d..b7240face8ad48122ae4b0ce9ac7201c978dc475 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 5ca9b4adca1f940f1ef1167a2336e5117a62d229..c9b0171faa678b9f040a7fd9b240b9d78a1379ea 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 d5feab6e538b07d0f5b487346004604390c8dd76..03e0623dee846cd300fc8705715b5682879676f4 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 27398f42bd335f3bdcfe78076b66f076ab42df49..13f2db4ffa0de236336d95c72c495d41faa42de1 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 a332be82849e4e75318b51ac9d0dcd52603fd70c..190763d27936ba6a728d14fca90b8b24dfe40690 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 227610a01ba41bc67fc22ece219b07c1fc5fe579..4cb64661b4dc618a667551905b154f2b101b96c2 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 8e2c3d836d3672284b2d30a908c1a7dc93a6bfcc..9942c8eddd3893627536b12eb16d51c09161deef 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 868124655382da67cd0b9841ff355a6f76ad9e84..f2861b769c3616644ccdfd4cf46c543a154ba76f 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 eb3ba594823c5f414defde7c0c77210eecc6e556..478d8ebee5a73c8fa726b09837685edf89e7133f 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 d6f04136dd146fffd4c4e36f1302834600140403..3832ad47608b34597a3ac588400f05abf0766f6d 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 a206c943cdbaa77320753b3112b8ee8567d0e508..708a7254d2a4d75041ed825959f26f282bd6bb6b 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 06b6f99fc9b234bce809512633d2b183a3b8ed24..cd5f97790a287a4e26a3793f5287216180907fe6 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 733e97028d781aba3e1079bbdfea8fb3fbfd5c92..58a84e38b44ce53fbc439d040da23a4c34f68626 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 165345f6ce258f0fb6c607ea6b8abc75ec622fdc..f63e3335ea3c02afe3ad7ffcfced66f6364832e3 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 39d238abcee325aefa71263f94efcf4b4d44033b..5939182d1f488b10edea539f3bfbc63d1aa82dea 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 65654eeefedf427e874c9083019361757544c6c1..33eb61e69ae50ff87040caadf848eade8302a990 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 b42c7986abb53c8f174851a1c3f83af774974c19..1c17b76eeb7aec7890d53a6fd9869cc325ed59aa 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 206c32eeb2b2ad0054b847e1cd1d10b404e60c72..8b3f629fa71ee79c9cd50821d230c6481c9f616e 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 2d11bfd21df350f2657e368dfe470d577fcce6f8..466e3a8a4d39e2c9c79f522f4c5cc3ad9e69f4e7 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 72bb4ec51adf51e3f2539769166f68e59e9ee669..f03ccf80767f84cd61eb608e1b098b3c580fd54a 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 c84a168c2d76b562ce15b4f4b7eea6c17ac7a1ac..85041130ff65c23b76088512b780b96e7b495ae9 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 3cb0580ba465828c513bfb1fb07a4b6af7e9335b..13120b50867c2b5dd4141ab01fa186b0294eb4df 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 5f8971e2a0b8ed87dda5dd1515de716a4d3d2b65..76381446e00920c58454bcdeb8b8dc91553b1476 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 cca90fdc9ac36d1373627284680d714600c661a6..6908b5259ae5cfe1844731c9b4298fbe1d770f53 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 866d84e94571f60e8a1621e7676775d08908da99..0e0e32a2e694a61970a3c36735f1252f0bd09ca9 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 e303c8a747a3fd3a79eb2d024969f64a178051a2..ebd0c4fc2c9bb2ef50b94cdec244d2dc4a7c5467 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 dd4d382befaed467343166d527cb864dc9b6d72e..aa345ca122a27bfd978ef37400acb6106ae576ab 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 ea0bd1c101d9e27a89eea1db2f37bfb66e27d92f..304d69b723a592e50b1d981ecd396090384e45b1 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 f3cbb702a7a7b27ee64d4a03777ab30f60b23b7e..59bb2aae16725529d397a9b630c85513bcc7aa8f 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 0182e453ba1324cfc1ab28b5ec9322c62e837208..0267c10adcc87a808c3c5fc17da4bec63a99076c 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 d32c1dc45c59481cc462eec3b3ccfef317883df5..79725f238514167cbcd1871bbb106cc5f8752fe7 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 ecbf12088739bbb67cc90d7c66c3c681af81becd..c0179eadb860d572bbfa457878efce1061e2595e 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 81f5e5447adcadc497f7dc868e9ee2c146a344a5..30ad3e4b7c0fe3091527e1af02e98463009abaef 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 db2d6c9b9c4f547162a8fd655a75876b3ee6886c..925c55df3253f646f5c3e4a29f2c662b9adfd471 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 a6f70bf4f2a47452790a906c4819cf5ffcedf2bd..27c1fd46db57e274f10bfcb60eb42032d7a4c501 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 e80f884202e797601981900eab113fb7ca12f90a..f33c64497dd1becff376c826a07b42333010d350 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 838044ba73313746825cb37c3aa96d8886ba560f..198bf00d52caab6aae55d8e2075ddacae20bfa94 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 54a8ae0322b02b2d2798da6399375c605e486058..a7850a5e903257f9a921d475fc91edb2148ca694 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 2e1330b3b8bbf0364a67f0fbbc3f1fcb4a0f1fdc..a455c7ac9a7674737888773a04d287e2cf5cff5d 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 f5b1d8ef01e945d3c5c93be98bb169b71e2aeafa..568d2ef4189b5a4659e31c142eb73add76479f55 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 2cb7b0fc5adc35f1b10ec4089baac32f47650989..2e1b449a7e025271e3ce3f9a382368b9b4affea3 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 d4c0f41aab897174f9f69176c9aa8c6bff0ae5c2..d56a53bf9142fc0f070449bc7e92c5d4edc52c24 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 5cdd9737d0cdc4caf1158a2fdb48455e3108b676..0069c53995ec152618ecc00a686fbed9197d3f73 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 f3fc5b546c46cb659cac43eddbc7601613202cb4..0aad99d4769bfce6d0601235a97a0f7abb9734a6 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 b69ce2e75d08d89ddc17f542202bc9007b6c6885..5b704a0f4bd162da8c7aae1eb47f502708b44dc2 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 5854748e99b83398c49cc52a2fce4c0666c80eb7..381e45e75bc95a241bcb5c5cd1e8c0ceffa807a9 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 479bc2cd34f4f57439c1dddc4cad338b3a830545..107fe3c0ccafc5b9c8de8334dad97c47fee5b8bd 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 115f9c67388e54828d4727dca36a6f72fe250065..b38144700cc079bdfd9156771ad102d7d09a8bd6 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 7a2be2efd05b733fd97bf3624c139a2933a8b739..b627d519139d8a4afb2d87150c57bce15667c709 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 bf2c86525b0e9013136afd911bb9a7372972c0c9..8906a58a616aebb81e5bec7681c3900a300ac74d 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 d291705f49575da742a73908f32676caddf4d88f..bdcf35449d243569f4d819069890e77088979eae 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 3b4df83409e38152de696ffe223c47dad18ed34e..3a4013394544e4829a74613e55a7978692fe3a00 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 e443f269848af443a24bff674a839ac86e245ad3..b4e3658a3ef2212562d3fd86c63e79a2d38a97fe 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 8d992e80b0ef9d90e2d490b97add5d2e5538b56a..152d4f2111f92bc025c773787a8bacbccdd6b529 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 850888b9f4fab372e9b5630f41c35b2f232e2e60..d5904e179dd49941462db886fcec9456d445a6ff 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 c08deda40381ded552f40e65218610dcda6d0408..44ad5d8ec0792d593e48203aeb92b611ba2e98e3 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 c0790b9698967aee8a984a1c635f68406c58c58d..02ae8778b9f5e762b9e77f226f053a825010154e 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 0e259d880925f241db007f841548d2e2ef9ed1c0..6e4db195d6f314a346cb9fa154cd0324347c9f9b 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 9f7da03a11c28a750ec91642f0460e8f5e07a94f..d89b6462b529fe25403694c950ab2331bd80e58f 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 73eb99dd8c7edaf5bd676711977c5c6227547fd4..75e071490fdc0b2dd666dc6b2c59dff30c045a6e 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 5f65eb5129764c325bd80fe1f5dcd13da1493cd9..6d212446c6362cca0357373a6e60c6b2102f1fc1 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 f2e3d510675c8d1ba5e31e43b5af10b76de91839..08d97031b044fd0379c63b333fde4638e4cdc742 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 ba1ee2328d7f3fb2c5e3f8e8949d649a15a7bc82..6df93827c4ee02585def86a345603587a1b4aa0e 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 0be323eb88dcfb76bd557efcc4df0224279c43a8..929f228a7afef7d725e1dca5008fe8290869a64c 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 d801fb8fa5b066b33f9748e2b5521a5dbd4cb076..24fa24adfa04927714958ac8d3c35584495c195c 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 71d75bd26e04d1f4b80e4294203595fe4b7b6737..02c401708c1dee8ba74f33069090a84c6b5c8235 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 0e1b3b7ff6fdc5bbc3b194fc9ace32adb2b4e700..e4c30664fb482ffd96cb436760c7d66f4e960515 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 9dc7a9f6430f19430733ed84c724b12fcb5b4252..55aa863bf71724c4fcecbb41c0c21b080b183645 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 d8bc359f1399197f266930bdd329ea70517ec70b..47b19d76048d6b0e454716c00176f390d5794413 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 4de12f14cb5b3fac9825264581899a4004ae48af..87d5ef0e03d90717476d6eb160d52ac195b25dc0 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 c92d10b68c0651deb5fbeb984bd49df5515d97db..06dcfb73077c30d985ee177facc5a52c839a1187 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 59a95ccc3d78c85abeb033a81f58404ec6f38783..6b15f83645ad712ad3820fe17e074193b5a0d217 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 b6d2528aef50abee3adbf79638a50cafaeb4abb6..c00236f7c9f8a5da6539e4fefb304872df54a5f9 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 6522f4d0cf7aa06b1a10b50833bc321ccc1178f6..8fbd1150031b1d231bb0fdb1743087646e9cdce2 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 477c09d47687ad38e71f9aa7c72de1e014b38f26..ea5b022f117853cb267d42d18c7ca1c25b41a7cf 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 f486afeb12403844dc4d950b12015a09a96373db..ec36978276013c0150671fadddad2263deb44362 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 4d3c735d21e5499c04a4958b29b80e5d2cfae9fe..c0656e0516e149060fc79d80ef0ad39e9b6c04bc 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 d9c53c4d7b34948546f56fbc593fdef04d64c0be..68d95d1f981a185bde6a83352ceb587c6aaed84c 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 a50f0e1b898c0045dd78afa11457e5bcf5336991..c220b6635c808c85e052c7c461e0485b10d10b3d 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 f7e53f1883ef8b15a8e03658eb3273409125a274..e29a6e486520e1a535267a8e67cd3808b641001a 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 7861c2215e4e95ea8f0682bd28b223112ca0a188..79ddc0d9beddc825555178727e8b9d4713a4c75f 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 ddc2da8c2f2c51294339fb84750753fcd736adaf..0d6c19e100ce7088224aa1e595f860bc377f091a 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 35ba1f3a5ece45bf21546928c9266f37c7baa289..2459be05e6319e46df44558706f572e6b33720dd 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 3fd085f7e86c0627af82af1237ac2d630964f3f7..e6ab27f26804d5013cac6f39d6281afa4072679e 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;