multiphase: Ensure that phi at fixedValue U BCs are up to date

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -58,8 +58,6 @@
         dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
     );
 
-    phi = dimensionedScalar("phi", phi.dimensions(), 0);
-
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -100,8 +98,6 @@
         );
         mrfZones.relativeFlux(phiHbyAs[phasei]);
 
-        phi += alphafs[phasei]*phiHbyAs[phasei];
-
         phiHbyAs[phasei] +=
             rAlphaAUfs[phasei]
            *(
@@ -161,6 +157,9 @@
         phasei++;
     }
 
+    // Reset phi BCs
+    phi.boundaryField() = 0;
+
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -172,6 +171,10 @@
         mrfZones.relativeFlux(phase.phi().oldTime());
         mrfZones.relativeFlux(phase.phi());
 
+        // Update phi BCs before pEqn
+        phi.boundaryField() +=
+            alphafs[phasei].boundaryField()*phase.phi().boundaryField();
+
         phasei++;
     }
 

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -8,9 +8,12 @@
     surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1));
     surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rAU2));
 
-    // Update the phi BCs from U before p BCs are updated
-    phi1.boundaryField() == (mesh.Sf().boundaryField() & U1.boundaryField());
-    phi2.boundaryField() == (mesh.Sf().boundaryField() & U2.boundaryField());
+    // Update the phi BCs from Us before p BCs are updated
+    phi.boundaryField() =
+        alpha1f.boundaryField()
+       *(mesh.Sf().boundaryField() & U1.boundaryField())
+      + alpha2f.boundaryField()
+       *(mesh.Sf().boundaryField() & U2.boundaryField());
 
     volVectorField HbyA1
     (
@@ -174,8 +177,14 @@
         {
             surfaceScalarField mSfGradp(pEqnIncomp.flux()/Dp);
 
+            phi1.boundaryField() ==
+                (mesh.Sf().boundaryField() & U1.boundaryField());
             phi1 = phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1);
+
+            phi2.boundaryField() ==
+                (mesh.Sf().boundaryField() & U2.boundaryField());
             phi2 = phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2);
+
             phi = alpha1f*phi1 + alpha2f*phi2;
 
             dgdt =