diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
index c0b4b73afa..cade87c352 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
@@ -61,74 +61,86 @@ int main(int argc, char *argv[])
 
         surfaceScalarField rho_pos
         (
+            "rho_pos",
             fvc::interpolate(rho, pos, "reconstruct(rho)")
         );
         surfaceScalarField rho_neg
         (
+            "rho_neg",
             fvc::interpolate(rho, neg, "reconstruct(rho)")
         );
 
         surfaceVectorField rhoU_pos
         (
+            "rhoU_pos",
             fvc::interpolate(rhoU, pos, "reconstruct(U)")
         );
         surfaceVectorField rhoU_neg
         (
+            "rhoU_neg",
             fvc::interpolate(rhoU, neg, "reconstruct(U)")
         );
 
         volScalarField rPsi(1.0/psi);
         surfaceScalarField rPsi_pos
         (
+            "rPsi_pos",
             fvc::interpolate(rPsi, pos, "reconstruct(T)")
         );
         surfaceScalarField rPsi_neg
         (
+            "rPsi_neg",
             fvc::interpolate(rPsi, neg, "reconstruct(T)")
         );
 
         surfaceScalarField e_pos
         (
+            "e_pos",
             fvc::interpolate(e, pos, "reconstruct(T)")
         );
         surfaceScalarField e_neg
         (
+            "e_neg",
             fvc::interpolate(e, neg, "reconstruct(T)")
         );
 
-        surfaceVectorField U_pos(rhoU_pos/rho_pos);
-        surfaceVectorField U_neg(rhoU_neg/rho_neg);
+        surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
+        surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
 
-        surfaceScalarField p_pos(rho_pos*rPsi_pos);
-        surfaceScalarField p_neg(rho_neg*rPsi_neg);
+        surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
+        surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
 
-        surfaceScalarField phiv_pos(U_pos & mesh.Sf());
-        surfaceScalarField phiv_neg(U_neg & mesh.Sf());
+        surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
+        surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
 
         volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
         surfaceScalarField cSf_pos
         (
+            "cSf_pos",
             fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf()
         );
         surfaceScalarField cSf_neg
         (
+            "cSf_neg",
             fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf()
         );
 
         surfaceScalarField ap
         (
+            "ap",
             max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
         );
         surfaceScalarField am
         (
+            "am",
             min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
         );
 
-        surfaceScalarField a_pos(ap/(ap - am));
+        surfaceScalarField a_pos("a_pos", ap/(ap - am));
 
         surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
 
-        surfaceScalarField aSf(am*a_pos);
+        surfaceScalarField aSf("aSf", am*a_pos);
 
         if (fluxScheme == "Tadmor")
         {
@@ -136,13 +148,13 @@ int main(int argc, char *argv[])
             a_pos = 0.5;
         }
 
-        surfaceScalarField a_neg(1.0 - a_pos);
+        surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
 
         phiv_pos *= a_pos;
         phiv_neg *= a_neg;
 
-        surfaceScalarField aphiv_pos(phiv_pos - aSf);
-        surfaceScalarField aphiv_neg(phiv_neg + aSf);
+        surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
+        surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
 
         // Reuse amaxSf for the maximum positive and negative fluxes
         // estimated by the central scheme
@@ -166,6 +178,7 @@ int main(int argc, char *argv[])
 
         surfaceScalarField phiEp
         (
+            "phiEp",
             aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
           + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
           + aSf*p_pos - aSf*p_neg
@@ -202,6 +215,7 @@ int main(int argc, char *argv[])
         // --- Solve energy
         surfaceScalarField sigmaDotU
         (
+            "sigmaDotU",
             (
                 fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
               + (mesh.Sf() & fvc::interpolate(tauMC))
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 0a276c6e6f..c66440d5bc 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -77,6 +77,7 @@ void Foam::multiphaseSystem::solveAlphas()
             phasei,
             new surfaceScalarField
             (
+                "phi" + alpha1.name() + "Corr",
                 fvc::flux
                 (
                     phi_,
@@ -106,7 +107,7 @@ void Foam::multiphaseSystem::solveAlphas()
             {
                 surfaceScalarField phic
                 (
-                    (mag(phi_) + mag(phase1.phi() - phase2.phi()))/mesh_.magSf()
+                    (mag(phi_) + mag(phir))/mesh_.magSf()
                 );
 
                 phir += min(cAlpha()*phic, max(phic))*nHatf(phase1, phase2);
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 0f32807132..e8112b300a 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -65,6 +65,7 @@
         const volScalarField& alpha = phase;
 
         alphafs.set(phasei, fvc::interpolate(alpha).ptr());
+        alphafs[phasei].rename("hmm" + alpha.name());
 
         volScalarField dragCoeffi
         (
@@ -227,7 +228,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            surfaceScalarField mSfGradp(pEqnIncomp.flux()/rAUf);
+            surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
 
             phasei = 0;
             phi = dimensionedScalar("phi", phi.dimensions(), 0);
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index e41b012860..1261843677 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -569,6 +569,7 @@ void Foam::multiphaseMixture::solveAlphas
             phasei,
             new surfaceScalarField
             (
+                "phi" + alpha.name() + "Corr",
                 fvc::flux
                 (
                     phi_,