PDRFoam: create HbyA and phiHbyA fields rather than reset U and phi

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = invA & UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = invA & UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -34,19 +35,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)*
         (
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             betav*fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*invA, p)
         );
 
@@ -54,7 +58,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -62,7 +66,7 @@ else
 #include "rhoEqn.H"
 #include "continuityErrs.H"
 
-U -= invA & (betav*fvc::grad(p));
+U = HbyA - (invA & (betav*fvc::grad(p)));
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);