ENH:

Update of overRhoPimpleDyMFoam and overInterDyMFoam solvers.
Adding corresponding tutorials with best possible settings
The main effort was put on reducing pressure spikes as the
stencil change with hole cells on the background mesh.

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/compressibleContinuityErrs.H

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    continuityErrs
+
+Description
+    Calculates and prints the continuity errors.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    dimensionedScalar totalMass = fvc::domainIntegrate(cellMask*rho);
+
+    scalar sumLocalContErr =
+    (
+        fvc::domainIntegrate(mag(cellMask*(rho - thermo.rho())))/totalMass
+    ).value();
+
+    scalar globalContErr =
+    (
+        fvc::domainIntegrate(cellMask*(rho - thermo.rho()))/totalMass
+    ).value();
+
+    cumulativeContErr += globalContErr;
+
+    Info<< "time step continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}
+
+// ************************************************************************* //

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H

@@ -1,11 +1,84 @@
-CorrectPhi
-(
-    U,
-    phi,
-    p,
-    rho,
-    psi,
-    dimensionedScalar("rAUf", dimTime, 1),
-    divrhoU,
-    pimple
-);
+if (mesh.changing())
+{
+    volVectorField::Boundary& bfld = U.boundaryFieldRef();
+    forAll(bfld, patchi)
+    {
+        if (bfld[patchi].fixesValue())
+        {
+            bfld[patchi].initEvaluate();
+        }
+    }
+
+    surfaceScalarField::Boundary& phiBfld = phi.boundaryFieldRef();
+    forAll(bfld, patchi)
+    {
+        if (bfld[patchi].fixesValue())
+        {
+            bfld[patchi].evaluate();
+
+            phiBfld[patchi] =
+                rho.boundaryField()[patchi]
+              * (
+                    bfld[patchi]
+                  & mesh.Sf().boundaryField()[patchi]
+                );
+        }
+    }
+}
+    // Initialize BCs list for pcorr to zero-gradient
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    // Set BCs of pcorr to fixed-value for patches at which p is fixed
+    forAll(p.boundaryField(), patchi)
+    {
+        if (p.boundaryField()[patchi].fixesValue())
+        {
+            pcorrTypes[patchi] = fixedValueFvPatchScalarField::typeName;
+        }
+    }
+
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(p.dimensions(), Zero),
+        pcorrTypes
+    );
+
+    mesh.setFluxRequired(pcorr.name());
+
+{
+    dimensionedScalar rAUf("rAUf", dimTime, 1.0);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::ddt(psi, pcorr)
+          + fvc::div(phi)
+          - fvm::laplacian(rAUf, pcorr)
+         ==
+            divrhoU()
+        );
+
+        pcorrEqn.solve(mesh.solver(pcorr.select(pimple.finalInnerIter())));
+        //Bypass virtual layer
+        //mesh.fvMesh::solve(pcorrEqn, d);
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi += pcorrEqn.flux();
+        }
+    }
+}

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H

@@ -1,11 +1,4 @@
-#include "createTimeControls.H"
-
-bool correctPhi
+bool ddtCorr
 (
-    pimple.dict().lookupOrDefault("correctPhi", true)
-);
-
-bool checkMeshCourantNo
-(
-    pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
+    pimple.dict().lookupOrDefault("ddtCorr", true)
 );

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H

@@ -1,10 +1,10 @@
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiThermo> pThermo
+autoPtr<fluidThermo> pThermo
 (
-    psiThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-psiThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField& p = thermo.p();
@@ -39,6 +39,8 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
 dimensionedScalar rhoMax
 (
     dimensionedScalar::lookupOrDefault
@@ -63,42 +65,25 @@ dimensionedScalar rhoMin
 
 mesh.setFluxRequired(p.name());
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar(p.dimensions()/dimTime, Zero)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
+#include "createDpdt.H"
 
+#include "createK.H"
 
 //- Overset specific
 
 // Add solver-specific interpolations
 {
-    dictionary oversetDict;
-    oversetDict.add("U", true);
-    oversetDict.add("p", true);
-    oversetDict.add("HbyA", true);
-    oversetDict.add("grad(p)", true);
+    wordHashSet& nonInt =
+        const_cast<wordHashSet&>(Stencil::New(mesh).nonInterpolatedFields());
 
-    const_cast<dictionary&>
-    (
-        mesh.schemesDict()
-    ).add
-    (
-        "oversetInterpolationRequired",
-        oversetDict,
-        true
-    );
+    nonInt.insert("HbyA");
+    nonInt.insert("grad(p)");
+    nonInt.insert("surfaceIntegrate(phi)");
+    nonInt.insert("surfaceIntegrate(phiHbyA)");
+    nonInt.insert("cellMask");
+    nonInt.insert("cellDisplacement");
+    nonInt.insert("interpolatedCells");
+    nonInt.insert("cellInterpolationWeight");
 }
 
 // Mask field for zeroing out contributions on hole cells

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/overRhoPimpleDyMFoam.C

@@ -38,10 +38,11 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "CorrectPhi.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
@@ -56,13 +57,13 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "createControl.H"
+    #include "createDyMControls.H"
     #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createMRF.H"
     #include "createFvOptions.H"
-    #include "createRhoUf.H"
+    #include "createRhoUfIfPresent.H"
     #include "createControls.H"
 
     turbulence->validate();
@@ -80,66 +81,107 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readControls.H"
+        #include "readDyMControls.H"
 
+
+        // Store divrhoU from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
+        autoPtr<volScalarField> divrhoU;
+        if (correctPhi)
         {
-            // Store divrhoU from the previous mesh so that it can be mapped
-            // and used in correctPhi to ensure the corrected phi has the
-            // same divergence
-            volScalarField divrhoU
+            divrhoU.reset
             (
-                "divrhoU",
-                fvc::div(fvc::absolute(phi, rho, U))
+                new volScalarField
+                (
+                    "divrhoU",
+                    fvc::div(fvc::absolute(phi, rho, U))
+                )
             );
-
-            if (LTS)
-            {
-                #include "setRDeltaT.H"
-            }
-            else
-            {
-                #include "compressibleCourantNo.H"
-                #include "setDeltaT.H"
-            }
-
-            ++runTime;
-
-            Info<< "Time = " << runTime.timeName() << nl << endl;
-
-            // Store momentum to set rhoUf for introduced faces.
-            volVectorField rhoU("rhoU", rho*U);
-
-            // Do any mesh changes
-            mesh.update();
-
-            if (mesh.changing())
-            {
-                #include "setCellMask.H"
-            }
-
-            if (mesh.changing() && correctPhi)
-            {
-                // Calculate absolute flux from the mapped surface velocity
-                phi = mesh.Sf() & rhoUf;
-
-                #include "correctPhi.H"
-
-                // Make the fluxes relative to the mesh-motion
-                fvc::makeRelative(phi, rho, U);
-            }
         }
 
-        if (mesh.changing() && checkMeshCourantNo)
+        if (LTS)
         {
-            #include "meshCourantNo.H"
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "compressibleCourantNo.H"
+            #include "setDeltaT.H"
         }
 
-        #include "rhoEqn.H"
-        Info<< "rhoEqn max/min : " << max(rho).value()
-            << " " << min(rho).value() << endl;
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+
+                // Do any mesh changes
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    MRF.update();
+
+                    #include "setCellMask.H"
+
+                    const surfaceScalarField faceMaskOld
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+                    );
+
+                    // Zero Uf on old faceMask (H-I)
+                    rhoUf() *= faceMaskOld;
+
+                    surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
+
+                    // Update Uf and phi on new C-I faces
+                    rhoUf() += (1-faceMaskOld)*rhoUfint;
+
+                    // Update Uf boundary
+                    forAll(rhoUf().boundaryField(), patchI)
+                    {
+                        rhoUf().boundaryFieldRef()[patchI] =
+                            rhoUfint.boundaryField()[patchI];
+                    }
+
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & rhoUf();
+
+                    if (correctPhi)
+                    {
+                        #include "correctPhi.H"
+                    }
+
+                    // Zero phi on current H-I
+                    const surfaceScalarField faceMask
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask)
+                    );
+
+                    phi *= faceMask;
+                    U   *= cellMask;
+
+                     // Make the fluxes relative to the mesh-motion
+                    fvc::makeRelative(phi, rho, U);
+
+                }
+
+                if (checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            if (pimple.firstIter() && !pimple.SIMPLErho())
+            {
+                #include "rhoEqn.H"
+            }
+
             #include "UEqn.H"
             #include "EEqn.H"
 
@@ -155,6 +197,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         runTime.printExecutionTime(Info);

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H

@@ -1,80 +1,93 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
-surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+volScalarField rAU("rAU", 1.0/UEqn.A());
+mesh.interpolate(rAU);
 
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA("HbyA", U);
-HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
+
+HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
 
 if (pimple.nCorrPISO() <= 1)
 {
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::interpolate(rho)*fvc::flux(HbyA)
+);
+
+if (ddtCorr)
+{
+    surfaceScalarField faceMaskOld
+    (
+        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+    );
+
+    phiHbyA +=
+        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf));
+}
+
+fvc::makeRelative(phiHbyA, rho, U);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
-       *(
-            fvc::flux(HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
-        )
+       (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
 
-    fvc::makeRelative(phid, psi, U);
-    MRF.makeRelative(fvc::interpolate(psi), phid);
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
 
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
+
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
+    fvScalarMatrix pDDtEqn
     (
-        "phiHbyA",
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
     );
 
-    fvc::makeRelative(phiHbyA, rho, U);
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
-        // Pressure corrector
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -91,25 +104,24 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 volVectorField gradP(fvc::grad(p));
 //mesh.interpolate(gradP);
-U = HbyA - rAU*cellMask*gradP;
+U = cellMask*(HbyA - rAU*gradP);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+if (pressureControl.limit(p))
 {
-    rhoUf = fvc::interpolate(rho*U);
-    surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
+    p.correctBoundaryConditions();
+}
+
+thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+rho = thermo.rho();
+
+{
+    // Correct rhoUf if the mesh is moving
+    fvc::correctRhoUf(rhoUf, rho, U, phi);
 }
 
 if (thermo.dpdt())
@@ -121,3 +133,9 @@ if (thermo.dpdt())
         dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
     }
 }
+
+surfaceScalarField faceMask
+(
+    localMin<scalar>(mesh).interpolate(cellMask)
+);
+phi *= faceMask;

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H

@@ -1,6 +1,9 @@
 #include "readTimeControls.H"
 
-correctPhi = pimple.dict().lookupOrDefault("correctPhi", true);
+correctPhi = pimple.dict().lookupOrDefault("correctPhi", false);
 
 checkMeshCourantNo =
     pimple.dict().lookupOrDefault("checkMeshCourantNo", false);
+
+
+ddtCorr = pimple.dict().lookupOrDefault("ddtCorr", true);

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H

@@ -107,7 +107,6 @@ fvOptions.correct(U);
     Uf += n*(phi/mesh.magSf() - (n & Uf));
 }
 
-
 // Make the fluxes relative to the mesh motion
 fvc::makeRelative(phi, U);
 

applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaSuSp.H

@@ -0,0 +1,9 @@
+zeroField Su;
+zeroField Sp;
+
+volScalarField::Internal divU
+(
+    mesh.moving()
+  ? fvc::div(phiCN() + mesh.phi())
+  : fvc::div(phiCN())
+);

applications/solvers/multiphase/interFoam/overInterDyMFoam/correctPhi.H

@@ -132,10 +132,4 @@
             phi -= pcorrEqn.flux();
         }
     }
-
-    //if (runTime.writeTime())
-    //{
-    //    volScalarField("contPhiPcorr", fvc::div(phi)).write();
-    //    pcorr.write();
-    //}
 }

applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H

@@ -30,37 +30,31 @@ volVectorField U
 
 #include "createPhi.H"
 
+//- Overset specific
+
+// Add solver-specific interpolations
+{
+    wordHashSet& nonInt =
+        const_cast<wordHashSet&>(Stencil::New(mesh).nonInterpolatedFields());
+
+    nonInt.insert("HbyA");
+    nonInt.insert("grad(p_rgh)");
+    nonInt.insert("nHat");
+    nonInt.insert("surfaceIntegrate(phi)");
+    nonInt.insert("surfaceIntegrate(phiHbyA)");
+    nonInt.insert("cellMask");
+    nonInt.insert("cellDisplacement");
+    nonInt.insert("interpolatedCells");
+    nonInt.insert("cellInterpolationWeight");
+    nonInt.insert("pcorr");
+}
 
 
-    //- Overset specific
-
-    // Add solver-specific interpolations
-    {
-        dictionary oversetDict;
-        oversetDict.add("U", true);
-        oversetDict.add("p", true);
-        oversetDict.add("HbyA", true);
-        oversetDict.add("p_rgh", true);
-        oversetDict.add("alpha1", true);
-        oversetDict.add("minGradP", true);
-
-        const_cast<dictionary&>
-        (
-            mesh.schemesDict()
-        ).add
-        (
-            "oversetInterpolationRequired",
-            oversetDict,
-            true
-        );
-    }
-
-    // Mask field for zeroing out contributions on hole cells
-    #include "createCellMask.H"
-
-    // Create bool field with interpolated cells
-    #include "createInterpolatedCells.H"
+// Mask field for zeroing out contributions on hole cells
+#include "createCellMask.H"
 
+// Create bool field with interpolated cells
+#include "createInterpolatedCells.H"
 
 
 Info<< "Reading transportProperties\n" << endl;

applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C

@@ -150,23 +150,47 @@ int main(int argc, char *argv[])
                     // Update cellMask field for blocking out hole cells
                     #include "setCellMask.H"
                     #include "setInterpolatedCells.H"
-                }
 
-                if ((mesh.changing() && correctPhi) || mesh.topoChanging())
-                {
-                    // Calculate absolute flux from the mapped surface velocity
-                    // Note: temporary fix until mapped Uf is assessed
-                    Uf = fvc::interpolate(U);
+                    const surfaceScalarField faceMaskOld
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+                    );
+
+                    // Zero Uf on old faceMask (H-I)
+                    Uf *= faceMaskOld;
+
+                    const surfaceVectorField Uint(fvc::interpolate(U));
+                    // Update Uf and phi on new C-I faces
+                    Uf += (1-faceMaskOld)*Uint;
+
+                    // Update Uf boundary
+                    forAll(Uf.boundaryField(), patchI)
+                    {
+                        Uf.boundaryFieldRef()[patchI] =
+                            Uint.boundaryField()[patchI];
+                    }
 
-                    // Calculate absolute flux from the mapped surface velocity
                     phi = mesh.Sf() & Uf;
 
-                    #include "correctPhi.H"
+                    // Correct phi on individual regions
+                    if (correctPhi)
+                    {
+                         #include "correctPhi.H"
+                    }
+
+                    mixture.correct();
+
+                    // Zero phi on current H-I
+                    const surfaceScalarField faceMask
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask)
+                    );
+                    phi *= faceMask;
+                    U   *= cellMask;
 
                     // Make the flux relative to the mesh motion
                     fvc::makeRelative(phi, U);
 
-                    mixture.correct();
                 }
 
                 if (mesh.changing() && checkMeshCourantNo)
@@ -175,9 +199,16 @@ int main(int argc, char *argv[])
                 }
             }
 
+
             #include "alphaControls.H"
             #include "alphaEqnSubCycle.H"
 
+            const surfaceScalarField faceMask
+            (
+                localMin<scalar>(mesh).interpolate(cellMask)
+            );
+            rhoPhi *= faceMask;
+
             mixture.correct();
 
             #include "UEqn.H"

applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H

@@ -1,23 +1,39 @@
 {
     rAU = 1.0/UEqn.A();
+    //mesh.interpolate(rAU);
+
     surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
+    volVectorField H("H", UEqn.H());
+
     volVectorField HbyA("HbyA", U);
     //HbyA = rAU*UEqn.H();
-    HbyA = constrainHbyA(rAU*UEqn.H(), U, p_rgh);
+    HbyA = constrainHbyA(rAU*H, U, p_rgh);
 
     if (massFluxInterpolation)
     {
         #include "interpolatedFaces.H"
     }
 
+    if (runTime.outputTime())
+    {
+        H.write();
+        rAU.write();
+        HbyA.write();
+    }
+
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 
     if (ddtCorr)
     {
-        phiHbyA += fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf);
+        surfaceScalarField faceMaskOld
+        (
+            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+        );
+        phiHbyA +=
+            fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf);
     }
     MRF.makeRelative(phiHbyA);
 
@@ -35,7 +51,6 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
-
     surfaceScalarField phig
     (
         (
@@ -60,7 +75,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(faceMask*rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -73,14 +88,12 @@
 
             p_rgh.relax();
 
-            // Reconstruct body forces (-grad(p) and gh etc)
-            volVectorField minGradP
-            (
-                "minGradP",
-                fvc::reconstruct((phig - p_rghEqn.flux())/rAUf)
-            );
-            //U = HbyA + rAU*cellMask*minGradP;
-            U = fvc::reconstruct(phi);
+            U =
+                cellMask*
+                (
+                   HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf)
+                );
+
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
@@ -97,16 +110,19 @@
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phi, U);
 
+    // Zero faces H-I for transport Eq after pEq
+    phi *= faceMask;
+
     p == p_rgh + rho*gh;
 
     if (p_rgh.needReference())
     {
-        p += dimensionedScalar
+        p_rgh += dimensionedScalar
         (
-            "p",
+            "p_rgh",
             p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
+            pRefValue - getRefCellValue(p_rgh, pRefCell)
         );
-        p_rgh = p - rho*gh;
+        p == p_rgh + rho*gh;
     }
 }