MRF: Separate MRF from fvOptions

fvOptions does not have the appropriate structure to support MRF as it
is based on option selection by user-specified fields whereas MRF MUST
be applied to all velocity fields in the particular solver.  A
consequence of the particular design choices in fvOptions made it
difficult to support MRF for multiphase and it is easier to support
frame-related and field related options separately.

Currently the MRF functionality provided supports only rotations but
the structure will be generalized to support other frame motions
including linear acceleration, SRF rotation and 6DoF which will be
run-time selectable.

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -80,6 +80,7 @@ int main(int argc, char *argv[])
     pisoControl potentialFlow(mesh, "potentialFlow");
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -90,8 +91,7 @@ int main(int argc, char *argv[])
     // function objects so do it ourselves
     runTime.functionObjects().start();
 
-    fvOptions.makeRelative(phi);
+    MRF.makeRelative(phi);
-
     adjustPhi(phi, U, p);
 
     // Non-orthogonal velocity potential corrector loop
@@ -113,7 +113,7 @@ int main(int argc, char *argv[])
         }
     }
 
-    fvOptions.makeAbsolute(phi);
+    MRF.makeAbsolute(phi);
 
     Info<< "Continuity error = "
         << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -92,6 +92,7 @@ int main(int argc, char *argv[])
     #include "readCombustionProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/combustion/PDRFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        betav*fvm::ddt(rho, U)
+        betav*fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         betav*rho*g

applications/solvers/combustion/XiFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         rho*g

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -72,6 +72,7 @@ int main(int argc, char *argv[])
     #include "readCombustionProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/combustion/XiFoam/pEqn.H

@@ -18,7 +18,7 @@ if (pimple.transonic())
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -50,7 +50,7 @@ else
         )
     );
 
-    fvOptions.makeRelative(phiHbyA);
+    MRF.makeRelative(phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -50,6 +50,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "initContinuityErrs.H"

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -74,6 +74,7 @@ int main(int argc, char *argv[])
 
     #include "readCombustionProperties.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "initContinuityErrs.H"

applications/solvers/combustion/engineFoam/pEqn.H

@@ -21,7 +21,7 @@ if (pimple.transonic())
     );
 
     fvc::makeRelative(phid, psi, U);
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -54,7 +54,7 @@ else
     );
 
     fvc::makeRelative(phiHbyA, rho, U);
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/combustion/fireFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         parcels.SU(U)

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -55,6 +55,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createClouds.H"
     #include "createSurfaceFilmModel.H"

applications/solvers/combustion/fireFoam/pEqn.H

@@ -18,7 +18,7 @@ surfaceScalarField phiHbyA
   + phig
 );
 
-fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 // Update the fixedFluxPressure BCs to ensure flux consistency
 setSnGrad<fixedFluxPressureFvPatchScalarField>
@@ -26,7 +26,7 @@ setSnGrad<fixedFluxPressureFvPatchScalarField>
     p_rgh.boundaryField(),
     (
         phiHbyA.boundaryField()
-      - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+      - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
        *rho.boundaryField()
     )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
 );

applications/solvers/combustion/reactingFoam/LTSReactingFoam/LTSReactingFoam.C

@@ -50,6 +50,7 @@ int main(int argc, char *argv[])
 
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/combustion/reactingFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         rho*g

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -18,7 +18,7 @@ if (pimple.transonic())
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -50,7 +50,7 @@ else
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -48,6 +48,7 @@ int main(int argc, char *argv[])
 
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         fvOptions(rho, U)

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H

@@ -23,7 +23,7 @@
       + phig
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
@@ -31,7 +31,7 @@
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
            *rho.boundaryField()
         )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
     );

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C

@@ -49,6 +49,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H

@@ -22,7 +22,7 @@
             )
         );
 
-        fvOptions.makeRelative(phiHbyA);
+        MRF.makeRelative(phiHbyA);
 
         surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
 
@@ -63,7 +63,7 @@
             )
         );
 
-        fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
         fvScalarMatrix pDDtEqn
         (

applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C

@@ -49,6 +49,7 @@ int main(int argc, char *argv[])
 
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/compressible/rhoPimpleFoam/UEqn.H

@@ -1,9 +1,11 @@
 // Solve the Momentum equation
 
+MRF.correctBoundaryVelocity(U);
+
 tmp<fvVectorMatrix> UEqn
 (
-    fvm::ddt(rho, U)
+    fvm::ddt(rho, U) + fvm::div(phi, U)
-  + fvm::div(phi, U)
+  + MRF.DDt(rho, U)
   + turbulence->divDevRhoReff(U)
  ==
     fvOptions(rho, U)

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -26,7 +26,7 @@ if (pimple.transonic())
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -58,7 +58,7 @@ else
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/rhoLTSPimpleFoam.C

@@ -55,6 +55,7 @@ int main(int argc, char *argv[])
 
     #include "setInitialrDeltaT.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

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

@@ -27,7 +27,7 @@ if (pimple.transonic())
     );
 
     fvc::makeRelative(phid, psi, U);
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -58,7 +58,7 @@ else
     );
 
     fvc::makeRelative(phiHbyA, rho, U);
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C

@@ -58,6 +58,7 @@ int main(int argc, char *argv[])
 
     #include "readControls.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "CourantNo.H"

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -51,6 +51,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

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

@@ -26,7 +26,7 @@ if (pimple.transonic())
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     surfaceScalarField phic
     (
@@ -69,7 +69,7 @@ else
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);

applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C

@@ -51,6 +51,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/compressible/rhoSimpleFoam/UEqn.H

@@ -1,8 +1,11 @@
     // Solve the Momentum equation
 
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         fvOptions(rho, U)

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -16,7 +16,7 @@
            *(fvc::interpolate(HbyA) & mesh.Sf())
         );
 
-        fvOptions.makeRelative(fvc::interpolate(psi), phid);
+        MRF.makeRelative(fvc::interpolate(psi), phid);
 
         while (simple.correctNonOrthogonal())
         {
@@ -49,7 +49,7 @@
             fvc::interpolate(rho*HbyA) & mesh.Sf()
         );
 
-        fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
         closedVolume = adjustPhi(phiHbyA, U, p);
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H

@@ -1,8 +1,11 @@
     // Construct the Momentum equation
 
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         fvOptions(rho, U)

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H

@@ -22,7 +22,7 @@
         fvc::interpolate(rho*HbyA) & mesh.Sf()
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     closedVolume = adjustPhi(phiHbyA, U, p);
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C

@@ -49,6 +49,7 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createZones.H"
     #include "initContinuityErrs.H"

applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C

@@ -47,6 +47,7 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -17,7 +17,7 @@ if (simple.transonic())
        *(fvc::interpolate(HbyA) & mesh.Sf())
     );
 
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     surfaceScalarField phic
     (
@@ -61,7 +61,7 @@ else
         fvc::interpolate(rho*HbyA) & mesh.Sf()
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     closedVolume = adjustPhi(phiHbyA, U, p);
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C

@@ -49,6 +49,7 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/compressible/sonicFoam/UEqn.H

@@ -1,9 +1,11 @@
 // Solve the Momentum equation
 
+MRF.correctBoundaryVelocity(U);
+
 fvVectorMatrix UEqn
 (
-    fvm::ddt(rho, U)
+    fvm::ddt(rho, U) + fvm::div(phi, U)
-  + fvm::div(phi, U)
+  + MRF.DDt(rho, U)
   + turbulence->divDevRhoReff(U)
  ==
     fvOptions(rho, U)

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -16,7 +16,7 @@ surfaceScalarField phid
     )
 );
 
-fvOptions.makeRelative(fvc::interpolate(psi), phid);
+MRF.makeRelative(fvc::interpolate(psi), phid);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())

applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H

@@ -17,7 +17,7 @@ surfaceScalarField phid
 );
 
 fvc::makeRelative(phid, psi, U);
-fvOptions.makeRelative(fvc::interpolate(psi), phid);
+MRF.makeRelative(fvc::interpolate(psi), phid);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())

applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C

@@ -54,6 +54,7 @@ int main(int argc, char *argv[])
 
     #include "readControls.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "CourantNo.H"

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -50,6 +50,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H

@@ -1,9 +1,11 @@
     // Solve the momentum equation
 
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(U)
+        fvm::ddt(U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(U)
       + turbulence->divDevReff(U)
      ==
         fvOptions(U)

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C

@@ -65,6 +65,7 @@ int main(int argc, char *argv[])
 
     #include "createFields.H"
     #include "createIncompressibleRadiationModel.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -15,7 +15,7 @@
       + phig
     );
 
-    fvOptions.makeRelative(phiHbyA);
+    MRF.makeRelative(phiHbyA);
 
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
@@ -23,7 +23,7 @@
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );
 

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H

@@ -1,8 +1,11 @@
     // Solve the momentum equation
 
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
+      + MRF.DDt(U)
       + turbulence->divDevReff(U)
      ==
         fvOptions(U)

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C

@@ -63,6 +63,7 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -14,7 +14,7 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
     );
 
-    fvOptions.makeRelative(phiHbyA);
+    MRF.makeRelative(phiHbyA);
 
     adjustPhi(phiHbyA, U, p_rgh);
 
@@ -26,7 +26,7 @@
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );
 

applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H

@@ -1,9 +1,11 @@
     // Solve the Momentum equation
 
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         fvOptions(rho, U)

applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C

@@ -53,6 +53,7 @@ int main(int argc, char *argv[])
 
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRadiationModel.H"
     #include "initContinuityErrs.H"

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -23,7 +23,7 @@
       + phig
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
@@ -31,7 +31,7 @@
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
            *rho.boundaryField()
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );

applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H

@@ -1,8 +1,11 @@
     // Solve the Momentum equation
 
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         fvOptions(rho, U)

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -50,6 +50,7 @@ int main(int argc, char *argv[])
 
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRadiationModel.H"
     #include "initContinuityErrs.H"

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -17,7 +17,7 @@
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
 
@@ -29,7 +29,7 @@
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
            *rho.boundaryField()
         )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
     );

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H

@@ -1,7 +1,11 @@
     // Solve the Momentum equation
+
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turb.divDevRhoReff(U)
      ==
         fvOptions(rho, U)

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H

@@ -18,6 +18,7 @@
     PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
     PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
 
+    PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
     PtrList<fv::IOoptionList> fluidFvOptions(fluidRegions.size());
 
     // Populate fluid field pointer lists
@@ -205,6 +206,13 @@
             )
         );
 
+        Info<< "    Adding MRF\n" << endl;
+        MRFfluid.set
+        (
+            i,
+            new IOMRFZoneList(fluidRegions[i])
+        );
+
         Info<< "    Adding fvOptions\n" << endl;
         fluidFvOptions.set
         (

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -19,7 +19,7 @@
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
 
@@ -31,7 +31,7 @@
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
            *rho.boundaryField()
         )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
     );

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -12,6 +12,7 @@
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
 
+    IOMRFZoneList& MRF = MRFfluid[i];
     fv::IOoptionList& fvOptions = fluidFvOptions[i];
 
     const dimensionedScalar initialMass

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H

@@ -1,8 +1,11 @@
     // Solve the Momentum equation
+
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turb.divDevRhoReff(U)
      ==
         fvOptions(rho, U)

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -15,6 +15,7 @@
     List<scalar> initialMassFluid(fluidRegions.size());
     List<bool> frozenFlowFluid(fluidRegions.size(), false);
 
+    PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
     PtrList<fv::IOoptionList> fluidFvOptions(fluidRegions.size());
 
     // Populate fluid field pointer lists
@@ -195,6 +196,13 @@
             fluidRegions[i].solutionDict().subDict("PIMPLE");
         pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
 
+        Info<< "    Adding MRF\n" << endl;
+        MRFfluid.set
+        (
+            i,
+            new IOMRFZoneList(fluidRegions[i])
+        );
+
         Info<< "    Adding fvOptions\n" << endl;
         fluidFvOptions.set
         (

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -23,7 +23,7 @@
       + phig
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
@@ -31,7 +31,7 @@
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
            *rho.boundaryField()
         )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
     );

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -20,6 +20,7 @@
 
     radiation::radiationModel& rad = radiation[i];
 
+    IOMRFZoneList& MRF = MRFfluid[i];
     fv::IOoptionList& fvOptions = fluidFvOptions[i];
 
     const dimensionedScalar initialMass

applications/solvers/incompressible/pimpleFoam/UEqn.H

@@ -1,9 +1,11 @@
 // Solve the Momentum equation
 
+MRF.correctBoundaryVelocity(U);
+
 tmp<fvVectorMatrix> UEqn
 (
-    fvm::ddt(U)
+    fvm::ddt(U) + fvm::div(phi, U)
-  + fvm::div(phi, U)
+  + MRF.DDt(U)
   + turbulence->divDevReff(U)
  ==
     fvOptions(U)

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -15,7 +15,7 @@ surfaceScalarField phiHbyA
   + rAUf*fvc::ddtCorr(U, phi)
 );
 
-fvOptions.makeRelative(phiHbyA);
+MRF.makeRelative(phiHbyA);
 
 adjustPhi(phiHbyA, U, p);
 
@@ -25,7 +25,7 @@ setSnGrad<fixedFluxPressureFvPatchScalarField>
     p.boundaryField(),
     (
         phiHbyA.boundaryField()
-      - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+      - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
     )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
 );
 

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

@@ -15,6 +15,8 @@ surfaceScalarField phiHbyA
   + rAUf*fvc::ddtCorr(U, Uf)
 );
 
+MRF.makeRelative(phiHbyA);
+
 if (p.needReference())
 {
     fvc::makeRelative(phiHbyA, U);
@@ -22,6 +24,17 @@ if (p.needReference())
     fvc::makeAbsolute(phiHbyA, U);
 }
 
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+    p.boundaryField(),
+    (
+        phiHbyA.boundaryField()
+      - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
+    )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+);
+
+// Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C

@@ -39,6 +39,7 @@ Description
 #include "pimpleControl.H"
 #include "CorrectPhi.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,6 +54,7 @@ int main(int argc, char *argv[])
 
     #include "createFields.H"
     #include "createUf.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "readTimeControls.H"
     #include "CourantNo.H"

applications/solvers/incompressible/pimpleFoam/pimpleFoam.C

@@ -39,8 +39,6 @@ Description
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
-#include "IOporosityModelList.H"
-#include "IOMRFZoneList.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -54,6 +52,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/incompressible/simpleFoam/UEqn.H

@@ -1,8 +1,11 @@
     // Momentum predictor
 
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
+      + MRF.DDt(U)
       + turbulence->divDevReff(U)
       ==
         fvOptions(U)

applications/solvers/incompressible/simpleFoam/pEqn.H

@@ -4,7 +4,7 @@
     HbyA = rAU*UEqn().H();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
-    fvOptions.makeRelative(phiHbyA);
+    MRF.makeRelative(phiHbyA);
     adjustPhi(phiHbyA, U, p);
 
     tmp<volScalarField> rAtU(rAU);

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H

@@ -1,5 +1,7 @@
     // Construct the Momentum equation
 
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H

@@ -1,6 +1,3 @@
-    IOMRFZoneList MRF(mesh);
-    MRF.correctBoundaryVelocity(U);
-
     IOporosityModelList pZones(mesh);
     Switch pressureImplicitPorosity(false);
 

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C

@@ -35,7 +35,6 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "simpleControl.H"
-#include "IOMRFZoneList.H"
 #include "IOporosityModelList.H"
 #include "fvIOoptionList.H"
 
@@ -50,8 +49,9 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
+    #include "createPorousZones.H"
     #include "createFvOptions.H"
-    #include "createZones.H"
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/incompressible/simpleFoam/simpleFoam.C

@@ -46,6 +46,7 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/lagrangian/coalChemistryFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -58,6 +58,7 @@ int main(int argc, char *argv[])
 
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createClouds.H"
     #include "createRadiationModel.H"

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -18,7 +18,7 @@ if (pimple.transonic())
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -51,7 +51,7 @@ else
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         parcels.SU(U)

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -18,7 +18,7 @@ surfaceScalarField phiHbyA
   + phig
 );
 
-fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 // Update the fixedFluxPressure BCs to ensure flux consistency
 setSnGrad<fixedFluxPressureFvPatchScalarField>
@@ -26,7 +26,7 @@ setSnGrad<fixedFluxPressureFvPatchScalarField>
     p_rgh.boundaryField(),
     (
         phiHbyA.boundaryField()
-      - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+      - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
        *rho.boundaryField()
     )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
 );

applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C

@@ -53,6 +53,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createClouds.H"
     #include "createRadiationModel.H"

applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -57,6 +57,7 @@ int main(int argc, char *argv[])
 
     #include "readTimeControls.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRDeltaT.H"
     #include "createRadiationModel.H"

applications/solvers/lagrangian/reactingParcelFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -20,7 +20,7 @@
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     fvScalarMatrix pDDtEqn
     (

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -55,6 +55,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createRadiationModel.H"
     #include "createClouds.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H

@@ -1,6 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H

@@ -15,7 +15,7 @@
         fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (simple.correctNonOrthogonal())
     {

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C

@@ -56,6 +56,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createRadiationModel.H"
     #include "createClouds.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/lagrangian/sprayFoam/UEqn.H

@@ -1,9 +1,11 @@
 // Solve the Momentum equation
 
+MRF.correctBoundaryVelocity(U);
+
 tmp<fvVectorMatrix> UEqn
 (
-    fvm::ddt(rho, U)
+    fvm::ddt(rho, U) + fvm::div(phi, U)
-  + fvm::div(phi, U)
+  + MRF.DDt(rho, U)
   + turbulence->divDevRhoReff(U)
  ==
     rho.dimensionedInternalField()*g

applications/solvers/lagrangian/sprayFoam/pEqn.H

@@ -26,7 +26,7 @@ if (pimple.transonic())
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -59,7 +59,7 @@ else
         )
     );
 
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H

@@ -27,7 +27,7 @@ if (pimple.transonic())
     );
 
     fvc::makeRelative(phid, psi, U);
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -59,7 +59,7 @@ else
     );
 
     fvc::makeRelative(phiHbyA, rho, U);
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {

applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C

@@ -56,6 +56,7 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createRhoUf.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createClouds.H"
     #include "createRadiationModel.H"

applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C

@@ -54,6 +54,7 @@ int main(int argc, char *argv[])
 
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "createClouds.H"

applications/solvers/lagrangian/sprayFoam/sprayFoam.C

@@ -52,6 +52,7 @@ int main(int argc, char *argv[])
 
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createClouds.H"
     #include "createRadiationModel.H"

applications/solvers/multiphase/driftFluxFoam/UEqn.H

@@ -1,9 +1,11 @@
     // Solve the Momentum equation
 
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
       + fvc::div(UdmModel.tauDm())
       + turbulence->divDevRhoReff(U)
      ==

applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C

@@ -58,6 +58,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/multiphase/driftFluxFoam/pEqn.H

@@ -11,7 +11,7 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
-    fvOptions.makeRelative(phiHbyA);
+    MRF.makeRelative(phiHbyA);
     adjustPhi(phiHbyA, U, p_rgh);
 
     surfaceScalarField phig
@@ -29,7 +29,7 @@
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );
 

applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C

@@ -63,6 +63,7 @@ int main(int argc, char *argv[])
 
     #include "initContinuityErrs.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "correctPhi.H"
     #include "CourantNo.H"

applications/solvers/multiphase/interFoam/UEqn.H

@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(rho, U)
      ==
         fvOptions(rho, U)

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -58,6 +58,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "readTimeControls.H"
 

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

@@ -11,6 +11,7 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
     );
+    MRF.makeRelative(phiHbyA);
 
     if (p_rgh.needReference())
     {

applications/solvers/multiphase/interFoam/interFoam.C

@@ -62,6 +62,7 @@ int main(int argc, char *argv[])
 
     #include "initContinuityErrs.H"
     #include "createFields.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "readTimeControls.H"
     #include "correctPhi.H"