Merge branch 'master' of github.com:OpenFOAM/OpenFOAM-2.3.x

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

@@ -53,8 +53,6 @@
 
         tpEqn().setReference(pRefCell, pRefValue);
 
-        fvOptions.constrain(tpEqn(), rho.name());
-
         tpEqn().solve();
 
         if (simple.finalNonOrthogonalIter())

applications/solvers/compressible/sonicFoam/EEqn.H

@@ -1,11 +1,21 @@
 {
-    solve
+    fvScalarMatrix EEqn
     (
         fvm::ddt(rho, e) + fvm::div(phi, e)
       + fvc::ddt(rho, K) + fvc::div(phi, K)
       + fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)")
       - fvm::laplacian(turbulence->alphaEff(), e)
+     ==
+        fvOptions(rho, e)
     );
 
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(e);
+
     thermo.correct();
 }

applications/solvers/compressible/sonicFoam/Make/options

@@ -1,6 +1,10 @@
 EXE_INC = \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
@@ -9,5 +13,7 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -lfiniteVolume

applications/solvers/compressible/sonicFoam/UEqn.H

@@ -3,7 +3,18 @@ fvVectorMatrix UEqn
     fvm::ddt(rho, U)
   + fvm::div(phi, U)
   + turbulence->divDevRhoReff(U)
+ ==
+    fvOptions(rho, U)
 );
 
-solve(UEqn == -fvc::grad(p));
-K = 0.5*magSqr(U);
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (pimple.momentumPredictor())
+{
+    solve(UEqn == -fvc::grad(p));
+
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -16,6 +16,8 @@ surfaceScalarField phid
     )
 );
 
+fvOptions.makeRelative(fvc::interpolate(psi), phid);
+
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
 {
@@ -24,8 +26,12 @@ while (pimple.correctNonOrthogonal())
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
       - fvm::laplacian(rhorAUf, p)
+     ==
+        fvOptions(psi, p, rho.name())
     );
 
+    fvOptions.constrain(pEqn);
+
     pEqn.solve();
 
     if (pimple.finalNonOrthogonalIter())
@@ -39,4 +45,5 @@ while (pimple.correctNonOrthogonal())
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+fvOptions.correct(U);
 K = 0.5*magSqr(U);

applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options

@@ -1,10 +1,13 @@
 EXE_INC = \
     -I.. \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -12,6 +15,8 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
     -lfiniteVolume \
-    -ldynamicMesh \
-    -lmeshTools
+    -ldynamicMesh

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

@@ -19,6 +19,7 @@ surfaceScalarField phid
 );
 
 fvc::makeRelative(phid, psi, U);
+fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
@@ -27,8 +28,12 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
       - fvm::laplacian(rhorAUf, p)
+     ==
+        fvOptions(psi, p, rho.name())
     );
 
+    fvOptions.constrain(pEqn);
+
     pEqn.solve();
 
     phi = pEqn.flux();
@@ -38,6 +43,8 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
 
 {
     rhoUf = fvc::interpolate(rho*U);

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ Description
 #include "turbulenceModel.H"
 #include "motionSolver.H"
 #include "pimpleControl.H"
+#include "fvIOoptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "initContinuityErrs.H"
 

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ Description
 #include "psiThermo.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fvIOoptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,6 +44,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     pimpleControl pimple(mesh);

applications/solvers/compressible/sonicFoam/sonicLiquidFoam/compressibleContinuityErrs.H

@@ -1,5 +1,5 @@
 {
-#   include "rhoEqn.H"
+    solve(fvm::ddt(rho) + fvc::div(phi));
 }
 {
     scalar sumLocalContErr =

applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C

@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
         #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
 
-        #include "rhoEqn.H"
+        solve(fvm::ddt(rho) + fvc::div(phi));
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.H

@@ -125,6 +125,12 @@ public:
 
     // Member Functions
 
+        //- Mixture properties
+        const incompressibleTwoPhaseInteractingMixture& mixture() const
+        {
+            return mixture_;
+        }
+
         //- Return the mixture mean density
         tmp<volScalarField> rho() const;
 

applications/solvers/multiphase/potentialFreeSurfaceFoam/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean
+wclean potentialFreeSurfaceDyMFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/potentialFreeSurfaceFoam/Allwmake

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake
+wmake potentialFreeSurfaceDyMFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H

@@ -1,5 +1,5 @@
 volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/Make/files

@@ -0,0 +1,3 @@
+potentialFreeSurfaceDyMFoam.C
+
+EXE = $(FOAM_APPBIN)/potentialFreeSurfaceDyMFoam

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/Make/options

@@ -0,0 +1,26 @@
+EXE_INC = \
+    -I.. \
+    -I../../interFoam/interDyMFoam \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+
+EXE_LIBS = \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/correctPhi.H

@@ -0,0 +1,58 @@
+if (mesh.changing())
+{
+    forAll(U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].fixesValue())
+        {
+            U.boundaryField()[patchI].initEvaluate();
+        }
+    }
+
+    forAll(U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].fixesValue())
+        {
+            U.boundaryField()[patchI].evaluate();
+
+            phi.boundaryField()[patchI] =
+                U.boundaryField()[patchI]
+              & mesh.Sf().boundaryField()[patchI];
+        }
+    }
+}
+
+{
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", p_gh.dimensions(), 0.0),
+        pcorrTypes
+    );
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
+        );
+
+        pcorrEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+    }
+
+    #include "continuityErrs.H"
+}

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/createPghCorrTypes.H

@@ -0,0 +1,13 @@
+    wordList pcorrTypes
+    (
+        p_gh.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    for (label i=0; i<p_gh.boundaryField().size(); i++)
+    {
+        if (p_gh.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H

@@ -0,0 +1,74 @@
+{
+    rAU = 1.0/UEqn().A();
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
+
+    if (pimple.nCorrPISO() <= 1)
+    {
+        UEqn.clear();
+    }
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + rAUf*fvc::ddtCorr(U, Uf)
+    );
+
+    if (p_gh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_gh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_gh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    // Non-orthogonal pressure corrector loop
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_ghEqn
+        (
+            fvm::laplacian(rAUf, p_gh) == fvc::div(phiHbyA)
+        );
+
+        p_ghEqn.setReference(p_ghRefCell, p_ghRefValue);
+
+        p_ghEqn.solve(mesh.solver(p_gh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_ghEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            p_gh.relax();
+
+            U = HbyA - rAU*fvc::grad(p_gh);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    {
+        Uf = fvc::interpolate(U);
+        surfaceVectorField n(mesh.Sf()/mesh.magSf());
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
+    }
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
+    p = p_gh + (g & mesh.C());
+}

applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/potentialFreeSurfaceDyMFoam.C

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    potentialFreeSurfaceDyMFoam
+
+Description
+    Incompressible Navier-Stokes solver with inclusion of a wave height field
+    to enable single-phase free-surface approximations.
+
+    Wave height field, zeta, used by pressure boundary conditions.
+
+    Optional mesh motion and mesh topology changes including adaptive
+    re-meshing.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "createPghCorrTypes.H"
+
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rAUf", dimTime, 1.0)
+    );
+
+    #include "correctPhi.H"
+    #include "createUf.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
+                }
+
+                if (mesh.changing() && correctPhi)
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
+
+                    #include "correctPhi.H"
+
+                    // Make the flux relative to the mesh motion
+                    fvc::makeRelative(phi, U);
+                }
+
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H

@@ -10,10 +10,10 @@
     fvScalarMatrix he1Eqn
     (
         fvm::ddt(alpha1, rho1, he1) + fvm::div(alphaRhoPhi1, he1)
-      - fvm::Sp(fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1), he1)
+      - fvm::Sp(contErr1, he1)
 
       + fvc::ddt(alpha1, rho1, K1) + fvc::div(alphaRhoPhi1, K1)
-      - (fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1))*K1
+      - contErr1*K1
 
       + (
             he1.name() == thermo1.phasePropertyName("e")
@@ -33,15 +33,16 @@
         heatTransferCoeff*(thermo2.T() - thermo1.T())
       + heatTransferCoeff*he1/Cpv1
       - fvm::Sp(heatTransferCoeff/Cpv1, he1)
+      + fvOptions(alpha1, rho1, he1)
     );
 
     fvScalarMatrix he2Eqn
     (
         fvm::ddt(alpha2, rho2, he2) + fvm::div(alphaRhoPhi2, he2)
-      - fvm::Sp(fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2), he2)
+      - fvm::Sp(contErr2, he2)
 
       + fvc::ddt(alpha2, rho2, K2) + fvc::div(alphaRhoPhi2, K2)
-      - (fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2))*K2
+      - contErr2*K2
 
       + (
             he2.name() == thermo2.phasePropertyName("e")
@@ -61,12 +62,15 @@
         heatTransferCoeff*(thermo1.T() - thermo2.T())
       + heatTransferCoeff*he2/Cpv2
       - fvm::Sp(heatTransferCoeff/Cpv2, he2)
+      + fvOptions(alpha2, rho2, he2)
     );
 
     he1Eqn.relax();
+    fvOptions.constrain(he1Eqn);
     he1Eqn.solve();
 
     he2Eqn.relax();
+    fvOptions.constrain(he2Eqn);
     he2Eqn.solve();
 
     thermo1.correct();

applications/solvers/multiphase/twoPhaseEulerFoam/Make/options

@@ -9,7 +9,10 @@ EXE_INC = \
     -IphaseCompressibleTurbulenceModels/lnInclude \
     -IinterfacialModels/lnInclude \
     -ItwoPhaseSystem/lnInclude \
-    -Iaveraging
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -21,4 +24,6 @@ EXE_LIBS = \
     -lcompressibleTwoPhaseSystem \
     -lcompressibleEulerianInterfacialModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling

applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H

@@ -9,7 +9,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
 
 {
     volScalarField virtualMassCoeff(fluid.virtualMassCoeff());
-
     volVectorField liftForce(fluid.liftForce());
     volVectorField wallLubricationForce(fluid.wallLubricationForce());
     volVectorField turbulentDispersionForce(fluid.turbulentDispersionForce());
@@ -19,15 +18,14 @@ volScalarField dragCoeff(fluid.dragCoeff());
         (
             fvm::ddt(alpha1, rho1, U1)
           + fvm::div(alphaRhoPhi1, U1)
-          - fvm::Sp(fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1), U1)
+          - fvm::Sp(contErr1, U1)
+          + mrfZones(alpha1*rho1 + virtualMassCoeff, U1)
           + phase1.turbulence().divDevRhoReff(U1)
          ==
           - fvm::Sp(dragCoeff, U1)
-          - (
-                liftForce
-              + wallLubricationForce
-              + turbulentDispersionForce
-            )
+          - liftForce
+          - wallLubricationForce
+          - turbulentDispersionForce
           - virtualMassCoeff
            *(
                 fvm::ddt(U1)
@@ -35,9 +33,10 @@ volScalarField dragCoeff(fluid.dragCoeff());
               - fvm::Sp(fvc::div(phi1), U1)
               - DDtU2
             )
+          + fvOptions(alpha1, rho1, U1)
         );
-        mrfZones.addCoriolis(alpha1*rho1 + virtualMassCoeff, U1Eqn);
         U1Eqn.relax();
+        fvOptions.constrain(U1Eqn);
     }
 
     {
@@ -45,15 +44,14 @@ volScalarField dragCoeff(fluid.dragCoeff());
         (
             fvm::ddt(alpha2, rho2, U2)
           + fvm::div(alphaRhoPhi2, U2)
-          - fvm::Sp(fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2), U2)
+          - fvm::Sp(contErr2, U2)
+          + mrfZones(alpha2*rho2 + virtualMassCoeff, U2)
           + phase2.turbulence().divDevRhoReff(U2)
          ==
           - fvm::Sp(dragCoeff, U2)
-          + (
-                liftForce
-              + wallLubricationForce
-              + turbulentDispersionForce
-            )
+          + liftForce
+          + wallLubricationForce
+          + turbulentDispersionForce
           - virtualMassCoeff
            *(
                 fvm::ddt(U2)
@@ -61,8 +59,9 @@ volScalarField dragCoeff(fluid.dragCoeff());
               - fvm::Sp(fvc::div(phi2), U2)
               - DDtU1
             )
+          + fvOptions(alpha2, rho2, U2)
         );
-        mrfZones.addCoriolis(alpha2*rho2 + virtualMassCoeff, U2Eqn);
         U2Eqn.relax();
+        fvOptions.constrain(U2Eqn);
     }
 }

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -130,6 +130,7 @@
         pEqnComp1 =
             (
                 fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
+              - (fvOptions(alpha1, rho1)&rho1)
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
             )/rho1
           + (alpha1/rho1)*correction
@@ -143,7 +144,9 @@
         pEqnComp2 =
             (
                 fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
+              - (fvOptions(alpha2, rho2)&rho2)
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+
             )/rho2
           + (alpha2/rho2)*correction
             (
@@ -158,6 +161,7 @@
         pEqnComp1 =
             (
                 fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
+              - (fvOptions(alpha1, rho1)&rho1)
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
             )/rho1
           + (alpha1*psi1/rho1)*correction(fvm::ddt(p));
@@ -165,6 +169,7 @@
         pEqnComp2 =
             (
                 fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
+              - (fvOptions(alpha2, rho2)&rho2)
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
             )/rho2
           + (alpha2*psi2/rho2)*correction(fvm::ddt(p));
@@ -227,6 +232,7 @@
                   - phiP1
                 );
             U1.correctBoundaryConditions();
+            fvOptions.correct(U1);
 
             U2 = HbyA2
               + fvc::reconstruct
@@ -239,6 +245,7 @@
                   - phiP2
                 );
             U2.correctBoundaryConditions();
+            fvOptions.correct(U2);
 
             U = fluid.U();
         }

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C

@@ -35,6 +35,7 @@ Description
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
 #include "IOMRFZoneList.H"
+#include "fvIOoptionList.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -48,6 +49,7 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRFZones.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "CourantNos.H"
@@ -74,8 +76,21 @@ int main(int argc, char *argv[])
             fluid.solve();
             fluid.correct();
 
-            #include "EEqns.H"
+            volScalarField contErr1
+            (
+                fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
+              - (fvOptions(alpha1, rho1)&rho1)
+            );
+
+            volScalarField contErr2
+            (
+                fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
+               - (fvOptions(alpha2, rho2)&rho2)
+            );
+
+
             #include "UEqns.H"
+            #include "EEqns.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

src/finiteVolume/cfdTools/compressible/rhoEqn.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,7 +30,19 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    solve(fvm::ddt(rho) + fvc::div(phi));
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        fvOptions(rho)
+    );
+
+    fvOptions.constrain(rhoEqn);
+
+    rhoEqn.solve();
+
+    fvOptions.correct(rho);
 }
 
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/MRF/MRFZone.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -338,7 +338,6 @@ void Foam::MRFZone::addCoriolis
     }
 
     const labelList& cells = mesh_.cellZones()[cellZoneID_];
-    const scalarField& V = mesh_.V();
     vectorField& ddtUc = ddtU.internalField();
     const vectorField& Uc = U.internalField();
 
@@ -347,7 +346,7 @@ void Foam::MRFZone::addCoriolis
     forAll(cells, i)
     {
         label celli = cells[i];
-        ddtUc[celli] += V[celli]*(Omega ^ Uc[celli]);
+        ddtUc[celli] += (Omega ^ Uc[celli]);
     }
 }
 

src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -160,6 +160,46 @@ void Foam::MRFZoneList::addCoriolis
 }
 
 
+Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::operator()
+(
+    const volVectorField& U
+)
+{
+    tmp<volVectorField> tacceleration
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "MRFZoneList:acceleration",
+                U.mesh().time().timeName(),
+                U.mesh()
+            ),
+            U.mesh(),
+            dimensionedVector("0", U.dimensions()/dimTime, vector::zero)
+        )
+    );
+    volVectorField& acceleration = tacceleration();
+
+    forAll(*this, i)
+    {
+        operator[](i).addCoriolis(U, acceleration);
+    }
+
+    return tacceleration;
+}
+
+
+Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::operator()
+(
+    const volScalarField& rho,
+    const volVectorField& U
+)
+{
+    return rho*operator()(U);
+}
+
+
 void Foam::MRFZoneList::makeRelative(volVectorField& U) const
 {
     forAll(*this, i)

src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,6 +102,19 @@ public:
         //- Add the Coriolis force contribution to the momentum equation
         void addCoriolis(const volScalarField& rho, fvVectorMatrix& UEqn) const;
 
+        //- Return the frame acceleration
+        tmp<volVectorField> operator()
+        (
+            const volVectorField& U
+        );
+
+        //- Return the frame acceleration force
+        tmp<volVectorField> operator()
+        (
+            const volScalarField& rho,
+            const volVectorField& U
+        );
+
         //- Make the given absolute velocity relative within the MRF region
         void makeRelative(volVectorField& U) const;
 

src/fvOptions/Make/options

@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \

src/fvOptions/fvOptions/fvIOoptionList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
     IOobject io
     (
         "fvOptions",
-        mesh.time().system(),
+        mesh.time().constant(),
         mesh,
         IOobject::MUST_READ,
         IOobject::NO_WRITE
@@ -45,7 +45,8 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
 
     if (io.headerOk())
     {
-        Info<< "Creating finite volume options from " << io.name() << nl
+        Info<< "Creating finite volume options from "
+            << io.instance()/io.name() << nl
             << endl;
 
         io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
@@ -53,10 +54,25 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
     }
     else
     {
-        Info<< "No finite volume options present" << nl << endl;
+        // Check if the fvOptions file is in system
+        io.instance() = mesh.time().system();
 
-        io.readOpt() = IOobject::NO_READ;
-        return io;
+        if (io.headerOk())
+        {
+            Info<< "Creating finite volume options from "
+                << io.instance()/io.name() << nl
+                << endl;
+
+            io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
+            return io;
+        }
+        else
+        {
+            Info<< "No finite volume options present" << nl << endl;
+
+            io.readOpt() = IOobject::NO_READ;
+            return io;
+        }
     }
 }
 

src/fvOptions/fvOptions/fvOption.C

@@ -403,13 +403,21 @@ void Foam::fv::option::correct(volTensorField& fld)
 }
 
 
-void Foam::fv::option::addSup(fvMatrix<scalar>& eqn, const label fieldI)
+void Foam::fv::option::addSup
+(
+    fvMatrix<scalar>& eqn,
+    const label fieldI
+)
 {
     // do nothing
 }
 
 
-void Foam::fv::option::addSup(fvMatrix<vector>& eqn, const label fieldI)
+void Foam::fv::option::addSup
+(
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
 {
     // do nothing
 }
@@ -425,18 +433,141 @@ void Foam::fv::option::addSup
 }
 
 
-void Foam::fv::option::addSup(fvMatrix<symmTensor>& eqn, const label fieldI)
+void Foam::fv::option::addSup
+(
+    fvMatrix<symmTensor>& eqn,
+    const label fieldI
+)
 {
     // do nothing
 }
 
 
-void Foam::fv::option::addSup(fvMatrix<tensor>& eqn, const label fieldI)
+void Foam::fv::option::addSup
+(
+    fvMatrix<tensor>& eqn,
+    const label fieldI
+)
 {
     // do nothing
 }
 
 
+void Foam::fv::option::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldI
+)
+{
+    // do nothing
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    // do nothing
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<sphericalTensor>& eqn,
+    const label fieldI
+)
+{
+    // do nothing
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<symmTensor>& eqn,
+    const label fieldI
+)
+{
+    // do nothing
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<tensor>& eqn,
+    const label fieldI
+)
+{
+    // do nothing
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldI
+)
+{
+    addSup(alpha*rho, eqn, fieldI);
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    addSup(alpha*rho, eqn, fieldI);
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<sphericalTensor>& eqn,
+    const label fieldI
+)
+{
+    addSup(alpha*rho, eqn, fieldI);
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<symmTensor>& eqn,
+    const label fieldI
+)
+{
+    addSup(alpha*rho, eqn, fieldI);
+}
+
+
+void Foam::fv::option::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<tensor>& eqn,
+    const label fieldI
+)
+{
+    addSup(alpha*rho, eqn, fieldI);
+}
+
+
 void Foam::fv::option::setValue(fvMatrix<scalar>& eqn, const label fieldI)
 {
     // do nothing

src/fvOptions/fvOptions/fvOption.H

@@ -384,6 +384,97 @@ public:
                 );
 
 
+            // Add explicit and implicit contributions to compressible equations
+
+                //- Scalar
+                virtual void addSup
+                (
+                    const volScalarField& rho,
+                    fvMatrix<scalar>& eqn,
+                    const label fieldI
+                );
+
+                //- Vector
+                virtual void addSup
+                (
+                    const volScalarField& rho,
+                    fvMatrix<vector>& eqn,
+                    const label fieldI
+                );
+
+                //- Spherical tensor
+                virtual void addSup
+                (
+                    const volScalarField& rho,
+                    fvMatrix<symmTensor>& eqn,
+                    const label fieldI
+                );
+
+                //- Symmetric tensor
+                virtual void addSup
+                (
+                    const volScalarField& rho,
+                    fvMatrix<sphericalTensor>& eqn,
+                    const label fieldI
+                );
+
+                //- Tensor
+                virtual void addSup
+                (
+                    const volScalarField& rho,
+                    fvMatrix<tensor>& eqn,
+                    const label fieldI
+                );
+
+
+            // Add explicit and implicit contributions to phase equations
+
+                //- Scalar
+                virtual void addSup
+                (
+                    const volScalarField& alpha,
+                    const volScalarField& rho,
+                    fvMatrix<scalar>& eqn,
+                    const label fieldI
+                );
+
+                //- Vector
+                virtual void addSup
+                (
+                    const volScalarField& alpha,
+                    const volScalarField& rho,
+                    fvMatrix<vector>& eqn,
+                    const label fieldI
+                );
+
+                //- Spherical tensor
+                virtual void addSup
+                (
+                    const volScalarField& alpha,
+                    const volScalarField& rho,
+                    fvMatrix<symmTensor>& eqn,
+                    const label fieldI
+                );
+
+                //- Symmetric tensor
+                virtual void addSup
+                (
+                    const volScalarField& alpha,
+                    const volScalarField& rho,
+                    fvMatrix<sphericalTensor>& eqn,
+                    const label fieldI
+                );
+
+                //- Tensor
+                virtual void addSup
+                (
+                    const volScalarField& alpha,
+                    const volScalarField& rho,
+                    fvMatrix<tensor>& eqn,
+                    const label fieldI
+                );
+
+
             // Set values directly
 
                 //- Scalar

src/fvOptions/fvOptions/fvOptionList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,10 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "fvOptionList.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvMesh.H"
 #include "surfaceFields.H"
-#include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -42,6 +39,37 @@ namespace fv
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+const Foam::dictionary& Foam::fv::optionList::optionsDict
+(
+    const dictionary& dict
+) const
+{
+    if (dict.found("options"))
+    {
+        return dict.subDict("options");
+    }
+    else
+    {
+        return dict;
+    }
+}
+
+
+bool Foam::fv::optionList::readOptions(const dictionary& dict)
+{
+    checkTimeIndex_ = mesh_.time().timeIndex() + 2;
+
+    bool allOk = true;
+    forAll(*this, i)
+    {
+        option& bs = this->operator[](i);
+        bool ok = bs.read(dict.subDict(bs.name()));
+        allOk = (allOk && ok);
+    }
+    return allOk;
+}
+
+
 void Foam::fv::optionList::checkApplied() const
 {
     if (mesh_.time().timeIndex() == checkTimeIndex_)
@@ -63,7 +91,7 @@ Foam::fv::optionList::optionList(const fvMesh& mesh, const dictionary& dict)
     mesh_(mesh),
     checkTimeIndex_(mesh_.time().startTimeIndex() + 2)
 {
-    reset(dict);
+    reset(optionsDict(dict));
 }
 
 
@@ -171,16 +199,7 @@ void Foam::fv::optionList::makeAbsolute
 
 bool Foam::fv::optionList::read(const dictionary& dict)
 {
-    checkTimeIndex_ = mesh_.time().timeIndex() + 2;
-
-    bool allOk = true;
-    forAll(*this, i)
-    {
-        option& bs = this->operator[](i);
-        bool ok = bs.read(dict.subDict(bs.name()));
-        allOk = (allOk && ok);
-    }
-    return allOk;
+    return readOptions(optionsDict(dict));
 }
 
 

src/fvOptions/fvOptions/fvOptionList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,12 @@ protected:
 
     // Protected Member Functions
 
+        //- Return the "options" sub-dictionary if present otherwise return dict
+        const dictionary& optionsDict(const dictionary& dict) const;
+
+        //- Read options dictionary
+        bool readOptions(const dictionary& dict);
+
         //- Check that all sources have been applied
         void checkApplied() const;
 
@@ -126,18 +132,37 @@ public:
             );
 
             //- Return source for equation
-            template<class Type, class RhoType>
+            template<class Type>
             tmp<fvMatrix<Type> > operator()
             (
-                const RhoType& rho,
+                const volScalarField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& fld
             );
 
             //- Return source for equation with specified name
-            template<class Type, class RhoType>
+            template<class Type>
             tmp<fvMatrix<Type> > operator()
             (
-                const RhoType& rho,
+                const volScalarField& rho,
+                GeometricField<Type, fvPatchField, volMesh>& fld,
+                const word& fieldName
+            );
+
+            //- Return source for equation
+            template<class Type>
+            tmp<fvMatrix<Type> > operator()
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                GeometricField<Type, fvPatchField, volMesh>& fld
+            );
+
+            //- Return source for equation with specified name
+            template<class Type>
+            tmp<fvMatrix<Type> > operator()
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& fld,
                 const word& fieldName
             );
@@ -149,10 +174,6 @@ public:
             template<class Type>
             void constrain(fvMatrix<Type>& eqn);
 
-            //- Apply constraints to equation with specified name
-            template<class Type>
-            void constrain(fvMatrix<Type>& eqn, const word& fieldName);
-
 
         // Flux manipulations
 

src/fvOptions/fvOptions/fvOptionListTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,10 +80,8 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
     const dimensionSet ds = fld.dimensions()/dimTime*dimVolume;
 
     tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(fld, ds));
-
     fvMatrix<Type>& mtx = tmtx();
 
-
     forAll(*this, i)
     {
         option& source = this->operator[](i);
@@ -111,10 +109,10 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
 }
 
 
-template<class Type, class RhoType>
+template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
 (
-    const RhoType& rho,
+    const volScalarField& rho,
     GeometricField<Type, fvPatchField, volMesh>& fld
 )
 {
@@ -122,10 +120,10 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
 }
 
 
-template<class Type, class RhoType>
+template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
 (
-    const RhoType& rho,
+    const volScalarField& rho,
     GeometricField<Type, fvPatchField, volMesh>& fld,
     const word& fieldName
 )
@@ -135,10 +133,8 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
     const dimensionSet ds = rho.dimensions()*fld.dimensions()/dimTime*dimVolume;
 
     tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(fld, ds));
-
     fvMatrix<Type>& mtx = tmtx();
 
-
     forAll(*this, i)
     {
         option& source = this->operator[](i);
@@ -157,7 +153,63 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
                         << fieldName << endl;
                 }
 
-                source.addSup(mtx, fieldI);
+                source.addSup(rho, mtx, fieldI);
+            }
+        }
+    }
+
+    return tmtx;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    GeometricField<Type, fvPatchField, volMesh>& fld
+)
+{
+    return this->operator()(alpha, rho, fld, fld.name());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    GeometricField<Type, fvPatchField, volMesh>& fld,
+    const word& fieldName
+)
+{
+    checkApplied();
+
+    const dimensionSet ds =
+        alpha.dimensions()*rho.dimensions()*fld.dimensions()/dimTime*dimVolume;
+
+    tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(fld, ds));
+    fvMatrix<Type>& mtx = tmtx();
+
+    forAll(*this, i)
+    {
+        option& source = this->operator[](i);
+
+        label fieldI = source.applyToField(fieldName);
+
+        if (fieldI != -1)
+        {
+            source.setApplied(fieldI);
+
+            if (source.isActive())
+            {
+                if (debug)
+                {
+                    Info<< "Applying source " << source.name() << " to field "
+                        << fieldName << endl;
+                }
+
+                source.addSup(alpha, rho, mtx, fieldI);
             }
         }
     }
@@ -168,17 +220,6 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
 
 template<class Type>
 void Foam::fv::optionList::constrain(fvMatrix<Type>& eqn)
-{
-    constrain(eqn, eqn.psi().name());
-}
-
-
-template<class Type>
-void Foam::fv::optionList::constrain
-(
-    fvMatrix<Type>& eqn,
-    const word& fieldName
-)
 {
     checkApplied();
 
@@ -186,7 +227,7 @@ void Foam::fv::optionList::constrain
     {
         option& source = this->operator[](i);
 
-        label fieldI = source.applyToField(fieldName);
+        label fieldI = source.applyToField(eqn.psi().name());
 
         if (fieldI != -1)
         {
@@ -197,7 +238,7 @@ void Foam::fv::optionList::constrain
                 if (debug)
                 {
                     Info<< "Applying constraint " << source.name()
-                        << " to field " << fieldName << endl;
+                        << " to field " << eqn.psi().name() << endl;
                 }
 
                 source.setValue(eqn, fieldI);

src/fvOptions/sources/derived/MRFSource/MRFSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,8 +89,7 @@ Foam::fv::MRFSource::MRFSource
 :
     option(name, modelType, dict, mesh),
     mrfPtr_(NULL),
-    UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
-    rhoName_(coeffs_.lookupOrDefault<word>("rhoName", "rho"))
+    UName_(coeffs_.lookupOrDefault<word>("UName", "U"))
 {
     initialise();
 }
@@ -104,19 +103,20 @@ void Foam::fv::MRFSource::addSup
     const label fieldI
 )
 {
-    if (eqn.dimensions() == dimForce)
-    {
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
+    // Add to rhs of equation
+    mrfPtr_->addCoriolis(eqn, true);
+}
 
-        // use 'true' flag to add to rhs of equation
-        mrfPtr_->addCoriolis(rho, eqn, true);
-    }
-    else
-    {
-        // use 'true' flag to add to rhs of equation
-        mrfPtr_->addCoriolis(eqn, true);
-    }
+
+void Foam::fv::MRFSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    // Add to rhs of equation
+    mrfPtr_->addCoriolis(rho, eqn, true);
 }
 
 
@@ -173,7 +173,6 @@ bool Foam::fv::MRFSource::read(const dictionary& dict)
     if (option::read(dict))
     {
         coeffs_.readIfPresent("UName", UName_);
-        coeffs_.readIfPresent("rhoName", rhoName_);
 
         initialise();
 

src/fvOptions/sources/derived/MRFSource/MRFSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,9 +79,6 @@ protected:
         //- Velocity field name, default = U
         word UName_;
 
-        //- Density field name, default = rho
-        word rhoName_;
-
 
     // Protected Member Functions
 
@@ -135,6 +132,17 @@ public:
             );
 
 
+        // Add explicit and implicit contributions to compressible equation
+
+            //- Vector
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
+
+
         // Flux manipulations
 
             //- Make the given absolute flux relative

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,40 +115,45 @@ void Foam::fv::actuationDiskSource::addSup
     const label fieldI
 )
 {
-    bool compressible = false;
-    if (eqn.dimensions() == dimForce)
-    {
-        compressible = true;
-    }
-
     const scalarField& cellsV = mesh_.V();
     vectorField& Usource = eqn.source();
     const vectorField& U = eqn.psi();
 
     if (V() > VSMALL)
     {
-        if (compressible)
-        {
-            addActuationDiskAxialInertialResistance
-            (
-                Usource,
-                cells_,
-                cellsV,
-                mesh_.lookupObject<volScalarField>("rho"),
-                U
-            );
-        }
-        else
-        {
-            addActuationDiskAxialInertialResistance
-            (
-                Usource,
-                cells_,
-                cellsV,
-                geometricOneField(),
-                U
-            );
-        }
+        addActuationDiskAxialInertialResistance
+        (
+            Usource,
+            cells_,
+            cellsV,
+            geometricOneField(),
+            U
+        );
+    }
+}
+
+
+void Foam::fv::actuationDiskSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    const scalarField& cellsV = mesh_.V();
+    vectorField& Usource = eqn.source();
+    const vectorField& U = eqn.psi();
+
+    if (V() > VSMALL)
+    {
+        addActuationDiskAxialInertialResistance
+        (
+            Usource,
+            cells_,
+            cellsV,
+            rho,
+            U
+        );
     }
 }
 

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -187,10 +187,22 @@ public:
             }
 
 
-        // Public Functions
+         // Add explicit and implicit contributions
 
-            //- Source term to fvMatrix<vector>
-            virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
+            //- Source term to momentum equation
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
+
+            //- Source term to compressible momentum equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
 
 
         // I-O

src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -204,6 +204,7 @@ bool Foam::fv::effectivenessHeatExchangerSource::alwaysApply() const
 
 void Foam::fv::effectivenessHeatExchangerSource::addSup
 (
+    const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label
 )

src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -247,10 +247,32 @@ public:
         virtual bool alwaysApply() const;
 
 
-        // Public Functions
+        // Add explicit and implicit contributions
 
-            //- Source term to fvMatrix<scalar>
-            virtual void addSup(fvMatrix<scalar>& eqn, const label fieldI);
+            //- Scalar
+            virtual void addSup
+            (
+                fvMatrix<scalar>& eqn,
+                const label fieldI
+            )
+            {
+                notImplemented
+                (
+                    "effectivenessHeatExchangerSource::addSup(eqn, fieldI): "
+                    "only compressible solvers supported."
+                );
+            }
+
+
+        // Add explicit and implicit contributions to compressible equation
+
+            //- Scalar
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn,
+                const label fieldI
+            );
 
 
         // I-O

src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -86,7 +86,6 @@ Foam::fv::explicitPorositySource::explicitPorositySource
     option(name, modelType, dict, mesh),
     porosityPtr_(NULL),
     UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
-    rhoName_(coeffs_.lookupOrDefault<word>("rhoName", "rho")),
     muName_(coeffs_.lookupOrDefault<word>("muName", "thermo:mu"))
 {
     initialise();
@@ -103,18 +102,23 @@ void Foam::fv::explicitPorositySource::addSup
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
 
-    if (eqn.dimensions() == dimForce)
-    {
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
-        const volScalarField& mu = mesh_.lookupObject<volScalarField>(muName_);
+    porosityPtr_->addResistance(porosityEqn);
 
-        porosityPtr_->addResistance(porosityEqn, rho, mu);
-    }
-    else
-    {
-        porosityPtr_->addResistance(porosityEqn);
-    }
+    eqn -= porosityEqn;
+}
+
+
+void Foam::fv::explicitPorositySource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
+
+    const volScalarField& mu = mesh_.lookupObject<volScalarField>(muName_);
+    porosityPtr_->addResistance(porosityEqn, rho, mu);
 
     eqn -= porosityEqn;
 }
@@ -132,7 +136,6 @@ bool Foam::fv::explicitPorositySource::read(const dictionary& dict)
     if (option::read(dict))
     {
         coeffs_.readIfPresent("UName", UName_);
-        coeffs_.readIfPresent("rhoName", rhoName_);
         coeffs_.readIfPresent("muName", muName_);
 
         return true;

src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -143,13 +143,21 @@ public:
 
         // Add explicit and implicit contributions
 
-            //- Vector
+            //- Add implicit contribution to momentum equation
             virtual void addSup
             (
                 fvMatrix<vector>& eqn,
                 const label fieldI
             );
 
+            //- Add implicit contribution to compressible momentum equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
+
 
         // I-O
 

src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -203,6 +203,17 @@ void Foam::fv::pressureGradientExplicitSource::addSup
 }
 
 
+void Foam::fv::pressureGradientExplicitSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    this->addSup(eqn, fieldI);
+}
+
+
 void Foam::fv::pressureGradientExplicitSource::setValue
 (
     fvMatrix<vector>& eqn,

src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -130,8 +130,20 @@ public:
             //- Correct the pressure gradient
             virtual void correct(volVectorField& U);
 
-            //- Add explicit contribution to equation
-            virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
+            //- Add explicit contribution to momentum equation
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
+
+            //- Add explicit contribution to compressible momentum equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
 
             //- Set 1/A coefficient
             virtual void setValue

src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,40 +68,45 @@ void Foam::fv::radialActuationDiskSource::addSup
     const label fieldI
 )
 {
-    bool compressible = false;
-    if (eqn.dimensions() == dimForce)
-    {
-        compressible = true;
-    }
-
     const scalarField& cellsV = mesh_.V();
     vectorField& Usource = eqn.source();
     const vectorField& U = eqn.psi();
 
     if (V_ > VSMALL)
     {
-        if (compressible)
-        {
-            addRadialActuationDiskAxialInertialResistance
-            (
-                Usource,
-                cells_,
-                cellsV,
-                mesh_.lookupObject<volScalarField>("rho"),
-                U
-            );
-        }
-        else
-        {
-            addRadialActuationDiskAxialInertialResistance
-            (
-                Usource,
-                cells_,
-                cellsV,
-                geometricOneField(),
-                U
-            );
-        }
+        addRadialActuationDiskAxialInertialResistance
+        (
+            Usource,
+            cells_,
+            cellsV,
+            geometricOneField(),
+            U
+        );
+    }
+}
+
+
+void Foam::fv::radialActuationDiskSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    const scalarField& cellsV = mesh_.V();
+    vectorField& Usource = eqn.source();
+    const vectorField& U = eqn.psi();
+
+    if (V_ > VSMALL)
+    {
+        addRadialActuationDiskAxialInertialResistance
+        (
+            Usource,
+            cells_,
+            cellsV,
+            rho,
+            U
+        );
     }
 }
 

src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -145,8 +145,20 @@ public:
 
     // Member Functions
 
-        //- Source term to fvMatrix<vector>
-        virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
+        //- Source term to momentum equation
+        virtual void addSup
+        (
+            fvMatrix<vector>& eqn,
+            const label fieldI
+        );
+
+        //- Source term to compressible momentum equation
+        virtual void addSup
+        (
+            const volScalarField& rho,
+            fvMatrix<vector>& eqn,
+            const label fieldI
+        );
 
 
     // I-O

src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C

@@ -25,10 +25,9 @@ License
 
 #include "rotorDiskSource.H"
 #include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
 #include "trimModel.H"
-#include "unitConversion.H"
 #include "fvMatrices.H"
+#include "geometricOneField.H"
 #include "syncTools.H"
 
 using namespace Foam::constant;
@@ -455,7 +454,6 @@ Foam::fv::rotorDiskSource::rotorDiskSource
 )
 :
     option(name, modelType, dict, mesh),
-    rhoName_("none"),
     rhoRef_(1.0),
     omega_(0.0),
     nBlades_(0),
@@ -486,8 +484,10 @@ Foam::fv::rotorDiskSource::~rotorDiskSource()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class RhoFieldType>
 void Foam::fv::rotorDiskSource::calculate
 (
+    const RhoFieldType& rho,
     const vectorField& U,
     const scalarField& thetag,
     vectorField& force,
@@ -496,8 +496,6 @@ void Foam::fv::rotorDiskSource::calculate
 ) const
 {
     const scalarField& V = mesh_.V();
-    const bool compressible = this->compressible();
-    tmp<volScalarField> trho(rho());
 
     // logging info
     scalar dragEff = 0.0;
@@ -570,11 +568,7 @@ void Foam::fv::rotorDiskSource::calculate
             scalar tipFactor = neg(radius/rMax_ - tipEffect_);
 
             // calculate forces perpendicular to blade
-            scalar pDyn = 0.5*magSqr(Uc);
-            if (compressible)
-            {
-                pDyn *= trho()[cellI];
-            }
+            scalar pDyn = 0.5*rho[cellI]*magSqr(Uc);
 
             scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;
             vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl);
@@ -596,7 +590,6 @@ void Foam::fv::rotorDiskSource::calculate
         }
     }
 
-
     if (output)
     {
         reduce(AOAmin, minOp<scalar>());
@@ -619,42 +612,69 @@ void Foam::fv::rotorDiskSource::addSup
     const label fieldI
 )
 {
-    dimensionSet dims = dimless;
-    if (eqn.dimensions() == dimForce)
-    {
-        coeffs_.lookup("rhoName") >> rhoName_;
-        dims.reset(dimForce/dimVolume);
-    }
-    else
-    {
-        coeffs_.lookup("rhoRef") >> rhoRef_;
-        dims.reset(dimForce/dimVolume/dimDensity);
-    }
-
     volVectorField force
     (
         IOobject
         (
             name_ + ":rotorForce",
             mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
+            mesh_
         ),
         mesh_,
-        dimensionedVector("zero", dims, vector::zero)
+        dimensionedVector
+        (
+            "zero",
+            eqn.dimensions()/dimVolume,
+            vector::zero
+        )
     );
 
-    const volVectorField& U = eqn.psi();
-
-    const vectorField Uin(inflowVelocity(U));
+    // Read the reference density for incompressible flow
+    coeffs_.lookup("rhoRef") >> rhoRef_;
 
+    const vectorField Uin(inflowVelocity(eqn.psi()));
     trim_->correct(Uin, force);
+    calculate(geometricOneField(), Uin, trim_->thetag(), force);
 
-    calculate(Uin, trim_->thetag(), force);
+    // Add source to rhs of eqn
+    eqn -= force;
+
+    if (mesh_.time().outputTime())
+    {
+        force.write();
+    }
+}
 
 
-    // add source to rhs of eqn
+void Foam::fv::rotorDiskSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    volVectorField force
+    (
+        IOobject
+        (
+            name_ + ":rotorForce",
+            mesh_.time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensionedVector
+        (
+            "zero",
+            eqn.dimensions()/dimVolume,
+            vector::zero
+        )
+    );
+
+    const vectorField Uin(inflowVelocity(eqn.psi()));
+    trim_->correct(rho, Uin, force);
+    calculate(rho, Uin, trim_->thetag(), force);
+
+    // Add source to rhs of eqn
     eqn -= force;
 
     if (mesh_.time().outputTime())
@@ -678,7 +698,6 @@ bool Foam::fv::rotorDiskSource::read(const dictionary& dict)
         coeffs_.lookup("fieldNames") >> fieldNames_;
         applied_.setSize(fieldNames_.size(), false);
 
-
         // read co-ordinate system/geometry invariant properties
         scalar rpm(readScalar(coeffs_.lookup("rpm")));
         omega_ = rpm/60.0*mathematical::twoPi;

src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.H

@@ -37,7 +37,6 @@ Description
     rotorDiskSourceCoeffs
     {
         fieldNames      (U);    // names of fields on which to apply source
-        rhoName         rho;    // density field if compressible case
         nBlades         3;      // number of blades
         tipEffect       0.96;   // normalised radius above which lift = 0
 
@@ -96,7 +95,6 @@ SourceFiles
 #include "bladeModel.H"
 #include "profileModelList.H"
 #include "volFieldsFwd.H"
-#include "dimensionSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -150,10 +148,7 @@ protected:
 
     // Protected data
 
-        //- Name of density field
-        word rhoName_;
-
-        //- Reference density for rhoName = 'none'
+        //- Reference density for incompressible case
         scalar rhoRef_;
 
         //- Rotational speed [rad/s]
@@ -262,7 +257,7 @@ public:
 
         // Access
 
-            //- Return the reference density for rhoName = 'none'
+            //- Return the reference density for incompressible case
             inline scalar rhoRef() const;
 
             //- Return the rotational speed [rad/s]
@@ -276,18 +271,14 @@ public:
             //- Return the rotor co-ordinate system (r, theta, z)
             inline const cylindricalCS& coordSys() const;
 
-            //- Return true if solving a compressible case
-            inline bool compressible() const;
-
-            //- Return the density field [kg/m3]
-            inline tmp<volScalarField> rho() const;
-
 
         // Evaluation
 
             //- Calculate forces
+            template<class RhoFieldType>
             void calculate
             (
+                const RhoFieldType& rho,
                 const vectorField& U,
                 const scalarField& thetag,
                 vectorField& force,
@@ -298,8 +289,20 @@ public:
 
         // Source term addition
 
-            //- Source term to fvMatrix<vector>
-            virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
+            //- Source term to momentum equation
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
+
+            //- Source term to compressible momentum equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
 
 
         // I-O

src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,23 +51,4 @@ const Foam::cylindricalCS& Foam::fv::rotorDiskSource::coordSys() const
 }
 
 
-bool Foam::fv::rotorDiskSource::compressible() const
-{
-    return rhoName_ != "none";
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::fv::rotorDiskSource::rho() const
-{
-    if (compressible())
-    {
-        return mesh_.lookupObject<volScalarField>(rhoName_);
-    }
-    else
-    {
-        return volScalarField::null();
-    }
-}
-
-
 // ************************************************************************* //

src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -86,7 +86,21 @@ Foam::tmp<Foam::scalarField> Foam::fixedTrim::thetag() const
 }
 
 
-void Foam::fixedTrim::correct(const vectorField& U, vectorField& force)
+void Foam::fixedTrim::correct
+(
+    const vectorField& U,
+    vectorField& force
+)
+{
+    // do nothing
+}
+
+
+void Foam::fixedTrim::correct
+(
+    const volScalarField rho,
+    const vectorField& U,
+    vectorField& force)
 {
     // do nothing
 }

src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,7 +80,19 @@ public:
         virtual tmp<scalarField> thetag() const;
 
         //- Correct the model
-        virtual void correct(const vectorField& U, vectorField& force);
+        virtual void correct
+        (
+            const vectorField& U,
+            vectorField& force
+        );
+
+        //- Correct the model for compressible flow
+        virtual void correct
+        (
+            const volScalarField rho,
+            const vectorField& U,
+            vectorField& force
+        );
 };
 
 

src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,9 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "targetCoeffTrim.H"
+#include "geometricOneField.H"
 #include "addToRunTimeSelectionTable.H"
-#include "unitConversion.H"
-#include "mathematicalConstants.H"
 
 using namespace Foam::constant;
 
@@ -42,17 +41,16 @@ namespace Foam
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+template<class RhoFieldType>
 Foam::vector Foam::targetCoeffTrim::calcCoeffs
 (
+    const RhoFieldType& rho,
     const vectorField& U,
     const scalarField& thetag,
     vectorField& force
 ) const
 {
-    rotor_.calculate(U, thetag, force, false, false);
-
-    bool compressible = rotor_.compressible();
-    tmp<volScalarField> trho = rotor_.rho();
+    rotor_.calculate(rho, U, thetag, force, false, false);
 
     const labelList& cells = rotor_.cells();
     const vectorField& C = rotor_.mesh().C();
@@ -76,11 +74,7 @@ Foam::vector Foam::targetCoeffTrim::calcCoeffs
         if (useCoeffs_)
         {
             scalar radius = x[i].x();
-            scalar coeff2 = coeff1*pow4(radius);
-            if (compressible)
-            {
-                coeff2 *= trho()[cellI];
-            }
+            scalar coeff2 = rho[cellI]*coeff1*pow4(radius);
 
             // add to coefficient vector
             cf[0] += (fc & yawAxis)/(coeff2 + ROOTVSMALL);
@@ -101,6 +95,99 @@ Foam::vector Foam::targetCoeffTrim::calcCoeffs
 }
 
 
+template<class RhoFieldType>
+void Foam::targetCoeffTrim::correctTrim
+(
+    const RhoFieldType& rho,
+    const vectorField& U,
+    vectorField& force
+)
+{
+    if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
+    {
+        word calcType = "forces";
+        if (useCoeffs_)
+        {
+            calcType = "coefficients";
+        }
+
+        Info<< type() << ":" << nl
+            << "    solving for target trim " << calcType << nl;
+
+        const scalar rhoRef = rotor_.rhoRef();
+
+        // iterate to find new pitch angles to achieve target force
+        scalar err = GREAT;
+        label iter = 0;
+        tensor J(tensor::zero);
+
+        vector old = vector::zero;
+        while ((err > tol_) && (iter < nIter_))
+        {
+            // cache initial theta vector
+            vector theta0(theta_);
+
+            // set initial values
+            old = calcCoeffs(rho, U, thetag(), force);
+
+            // construct Jacobian by perturbing the pitch angles
+            // by +/-(dTheta_/2)
+            for (label pitchI = 0; pitchI < 3; pitchI++)
+            {
+                theta_[pitchI] -= dTheta_/2.0;
+                vector cf0 = calcCoeffs(rho, U, thetag(), force);
+
+                theta_[pitchI] += dTheta_;
+                vector cf1 = calcCoeffs(rho, U, thetag(), force);
+
+                vector ddTheta = (cf1 - cf0)/dTheta_;
+
+                J[pitchI + 0] = ddTheta[0];
+                J[pitchI + 3] = ddTheta[1];
+                J[pitchI + 6] = ddTheta[2];
+
+                theta_ = theta0;
+            }
+
+            // calculate the change in pitch angle vector
+            vector dt = inv(J) & (target_/rhoRef - old);
+
+            // update pitch angles
+            vector thetaNew = theta_ + relax_*dt;
+
+            // update error
+            err = mag(thetaNew - theta_);
+
+            // update for next iteration
+            theta_ = thetaNew;
+            iter++;
+        }
+
+        if (iter == nIter_)
+        {
+            Info<< "    solution not converged in " << iter
+                << " iterations, final residual = " << err
+                << "(" << tol_ << ")" << endl;
+        }
+        else
+        {
+            Info<< "    final residual = " << err << "(" << tol_
+                << "), iterations = " << iter << endl;
+        }
+
+        Info<< "    current and target " << calcType << nl
+            << "        thrust  = " << old[0]*rhoRef << ", " << target_[0] << nl
+            << "        pitch   = " << old[1]*rhoRef << ", " << target_[1] << nl
+            << "        roll    = " << old[2]*rhoRef << ", " << target_[2] << nl
+            << "    new pitch angles [deg]:" << nl
+            << "        theta0  = " << radToDeg(theta_[0]) << nl
+            << "        theta1c = " << radToDeg(theta_[1]) << nl
+            << "        theta1s = " << radToDeg(theta_[2]) << nl
+            << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::targetCoeffTrim::targetCoeffTrim
@@ -185,90 +272,24 @@ Foam::tmp<Foam::scalarField> Foam::targetCoeffTrim::thetag() const
 }
 
 
-void Foam::targetCoeffTrim::correct(const vectorField& U, vectorField& force)
+void Foam::targetCoeffTrim::correct
+(
+    const vectorField& U,
+    vectorField& force
+)
 {
-    if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
-    {
-        word calcType = "forces";
-        if (useCoeffs_)
-        {
-            calcType = "coefficients";
-        }
+    correctTrim(geometricOneField(), U, force);
+}
 
-        Info<< type() << ":" << nl
-            << "    solving for target trim " << calcType << nl;
 
-        const scalar rhoRef = rotor_.rhoRef();
-
-        // iterate to find new pitch angles to achieve target force
-        scalar err = GREAT;
-        label iter = 0;
-        tensor J(tensor::zero);
-
-        vector old = vector::zero;
-        while ((err > tol_) && (iter < nIter_))
-        {
-            // cache initial theta vector
-            vector theta0(theta_);
-
-            // set initial values
-            old = calcCoeffs(U, thetag(), force);
-
-            // construct Jacobian by perturbing the pitch angles
-            // by +/-(dTheta_/2)
-            for (label pitchI = 0; pitchI < 3; pitchI++)
-            {
-                theta_[pitchI] -= dTheta_/2.0;
-                vector cf0 = calcCoeffs(U, thetag(), force);
-
-                theta_[pitchI] += dTheta_;
-                vector cf1 = calcCoeffs(U, thetag(), force);
-
-                vector ddTheta = (cf1 - cf0)/dTheta_;
-
-                J[pitchI + 0] = ddTheta[0];
-                J[pitchI + 3] = ddTheta[1];
-                J[pitchI + 6] = ddTheta[2];
-
-                theta_ = theta0;
-            }
-
-            // calculate the change in pitch angle vector
-            vector dt = inv(J) & (target_/rhoRef - old);
-
-            // update pitch angles
-            vector thetaNew = theta_ + relax_*dt;
-
-            // update error
-            err = mag(thetaNew - theta_);
-
-            // update for next iteration
-            theta_ = thetaNew;
-            iter++;
-        }
-
-        if (iter == nIter_)
-        {
-            Info<< "    solution not converged in " << iter
-                << " iterations, final residual = " << err
-                << "(" << tol_ << ")" << endl;
-        }
-        else
-        {
-            Info<< "    final residual = " << err << "(" << tol_
-                << "), iterations = " << iter << endl;
-        }
-
-        Info<< "    current and target " << calcType << nl
-            << "        thrust  = " << old[0]*rhoRef << ", " << target_[0] << nl
-            << "        pitch   = " << old[1]*rhoRef << ", " << target_[1] << nl
-            << "        roll    = " << old[2]*rhoRef << ", " << target_[2] << nl
-            << "    new pitch angles [deg]:" << nl
-            << "        theta0  = " << radToDeg(theta_[0]) << nl
-            << "        theta1c = " << radToDeg(theta_[1]) << nl
-            << "        theta1s = " << radToDeg(theta_[2]) << nl
-            << endl;
-    }
+void Foam::targetCoeffTrim::correct
+(
+    const volScalarField rho,
+    const vectorField& U,
+    vectorField& force
+)
+{
+    correctTrim(rho, U, force);
 }
 
 

src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -125,13 +125,24 @@ protected:
     // Protected member functions
 
         //- Calculate the rotor force and moment coefficients vector
+        template<class RhoFieldType>
         vector calcCoeffs
         (
+            const RhoFieldType& rho,
             const vectorField& U,
             const scalarField& alphag,
             vectorField& force
         ) const;
 
+        //- Correct the model
+        template<class RhoFieldType>
+        void correctTrim
+        (
+            const RhoFieldType& rho,
+            const vectorField& U,
+            vectorField& force
+        );
+
 
 public:
 
@@ -154,7 +165,19 @@ public:
         virtual tmp<scalarField> thetag() const;
 
         //- Correct the model
-        virtual void correct(const vectorField& U, vectorField& force);
+        virtual void correct
+        (
+            const vectorField& U,
+            vectorField& force
+        );
+
+        //- Correct the model for compressible flow
+        virtual void correct
+        (
+            const volScalarField rho,
+            const vectorField& U,
+            vectorField& force
+        );
 };
 
 

src/fvOptions/sources/derived/rotorDiskSource/trimModel/trimModel/trimModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -120,7 +120,19 @@ public:
         virtual tmp<scalarField> thetag() const = 0;
 
         //- Correct the model
-        virtual void correct(const vectorField& U, vectorField& force) = 0;
+        virtual void correct
+        (
+            const vectorField& U,
+            vectorField& force
+        ) = 0;
+
+        //- Correct the model for compressible flow
+        virtual void correct
+        (
+            const volScalarField rho,
+            const vectorField& U,
+            vectorField& force
+        ) = 0;
 };
 
 

src/fvOptions/sources/general/codedSource/CodedSource.C

@@ -181,6 +181,25 @@ void Foam::fv::CodedSource<Type>::addSup
 }
 
 
+template<class Type>
+void Foam::fv::CodedSource<Type>::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<Type>& eqn,
+    const label fieldI
+)
+{
+    if (debug)
+    {
+        Info<< "CodedSource<"<< pTraits<Type>::typeName
+            << ">::addSup for source " << name_ << endl;
+    }
+
+    updateLibrary(redirectType_);
+    redirectFvOption().addSup(rho, eqn, fieldI);
+}
+
+
 template<class Type>
 void Foam::fv::CodedSource<Type>::setValue
 (

src/fvOptions/sources/general/codedSource/CodedSource.H

@@ -42,7 +42,7 @@ Description
 
     setValue
     (
-         fvMatrix<Type}>& eqn,
+        fvMatrix<Type}>& eqn,
         const label fieldI
     )
 
@@ -206,6 +206,15 @@ public:
                 const label fieldI
             );
 
+            //- Explicit and implicit matrix contributions
+            //  to compressible equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<Type>& eqn,
+                const label fieldI
+            );
+
             //- Set value
             virtual void setValue
             (

src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -193,4 +193,22 @@ void Foam::fv::SemiImplicitSource<Type>::addSup
 }
 
 
+template<class Type>
+void Foam::fv::SemiImplicitSource<Type>::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<Type>& eqn,
+    const label fieldI
+)
+{
+    if (debug)
+    {
+        Info<< "SemiImplicitSource<" << pTraits<Type>::typeName
+            << ">::addSup for source " << name_ << endl;
+    }
+
+    return this->addSup(eqn, fieldI);
+}
+
+
 // ************************************************************************* //

src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,10 +80,10 @@ namespace fv
 
 // Forward declaration of classes
 
-
 template<class Type>
 class SemiImplicitSource;
 
+
 // Forward declaration of friend functions
 
 template<class Type>
@@ -93,6 +93,7 @@ Ostream& operator<<
     const SemiImplicitSource<Type>&
 );
 
+
 /*---------------------------------------------------------------------------*\
                      Class SemiImplicitSource Declaration
 \*---------------------------------------------------------------------------*/
@@ -184,7 +185,19 @@ public:
         // Evaluation
 
             //- Add explicit contribution to equation
-            virtual void addSup(fvMatrix<Type>& eqn, const label fieldI);
+            virtual void addSup
+            (
+                fvMatrix<Type>& eqn,
+                const label fieldI
+            );
+
+            //- Add explicit contribution to compressible equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<Type>& eqn,
+                const label fieldI
+            );
 
 
         // I-O

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -122,7 +122,6 @@ Foam::fv::interRegionExplicitPorositySource::interRegionExplicitPorositySource
     porosityPtr_(NULL),
     firstIter_(-1),
     UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
-    rhoName_(coeffs_.lookupOrDefault<word>("rhoName", "rho")),
     muName_(coeffs_.lookupOrDefault<word>("muName", "thermo:mu"))
 {
     if (active_)
@@ -171,64 +170,108 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
 
     fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
 
-    if (eqn.dimensions() == dimForce)
-    {
-        volScalarField rhoNbr
+    porosityPtr_->addResistance(nbrEqn);
+
+    // convert source from neighbour to local region
+    fvMatrix<vector> porosityEqn(U, eqn.dimensions());
+    scalarField& Udiag = porosityEqn.diag();
+    vectorField& Usource = porosityEqn.source();
+
+    Udiag.setSize(eqn.diag().size(), 0.0);
+    Usource.setSize(eqn.source().size(), vector::zero);
+
+    meshInterp().mapTgtToSrc(nbrEqn.diag(), plusEqOp<scalar>(), Udiag);
+    meshInterp().mapTgtToSrc(nbrEqn.source(), plusEqOp<vector>(), Usource);
+
+    eqn -= porosityEqn;
+}
+
+
+void Foam::fv::interRegionExplicitPorositySource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    initialise();
+
+    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
+
+    const volVectorField& U = eqn.psi();
+
+    volVectorField UNbr
+    (
+        IOobject
         (
-            IOobject
-            (
-                "rho:UNbr",
-                nbrMesh.time().timeName(),
-                nbrMesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            name_ + ":UNbr",
+            nbrMesh.time().timeName(),
             nbrMesh,
-            dimensionedScalar("zero", dimDensity, 0.0)
-        );
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        nbrMesh,
+        dimensionedVector("zero", U.dimensions(), vector::zero)
+    );
 
-        volScalarField muNbr
+    // map local velocity onto neighbour region
+    meshInterp().mapSrcToTgt
+    (
+        U.internalField(),
+        plusEqOp<vector>(),
+        UNbr.internalField()
+    );
+
+    fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
+
+    volScalarField rhoNbr
+    (
+        IOobject
         (
-            IOobject
-            (
-                "mu:UNbr",
-                nbrMesh.time().timeName(),
-                nbrMesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            "rho:UNbr",
+            nbrMesh.time().timeName(),
             nbrMesh,
-            dimensionedScalar("zero", dimViscosity, 0.0)
-        );
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        nbrMesh,
+        dimensionedScalar("zero", dimDensity, 0.0)
+    );
 
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
-
-        const volScalarField& mu =
-            mesh_.lookupObject<volScalarField>(muName_);
-
-        // map local rho onto neighbour region
-        meshInterp().mapSrcToTgt
+    volScalarField muNbr
+    (
+        IOobject
         (
-            rho.internalField(),
-            plusEqOp<scalar>(),
-            rhoNbr.internalField()
-        );
+            "mu:UNbr",
+            nbrMesh.time().timeName(),
+            nbrMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        nbrMesh,
+        dimensionedScalar("zero", dimViscosity, 0.0)
+    );
 
-        // map local mu onto neighbour region
-        meshInterp().mapSrcToTgt
-        (
-            mu.internalField(),
-            plusEqOp<scalar>(),
-            muNbr.internalField()
-        );
+    const volScalarField& mu =
+        mesh_.lookupObject<volScalarField>(muName_);
 
-        porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);
-    }
-    else
-    {
-        porosityPtr_->addResistance(nbrEqn);
-    }
+    // map local rho onto neighbour region
+    meshInterp().mapSrcToTgt
+    (
+        rho.internalField(),
+        plusEqOp<scalar>(),
+        rhoNbr.internalField()
+    );
+
+    // map local mu onto neighbour region
+    meshInterp().mapSrcToTgt
+    (
+        mu.internalField(),
+        plusEqOp<scalar>(),
+        muNbr.internalField()
+    );
+
+    porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);
 
     // convert source from neighbour to local region
     fvMatrix<vector> porosityEqn(U, eqn.dimensions());
@@ -257,7 +300,6 @@ bool Foam::fv::interRegionExplicitPorositySource::read(const dictionary& dict)
     if (option::read(dict))
     {
         coeffs_.readIfPresent("UName", UName_);
-        coeffs_.readIfPresent("rhoName", rhoName_);
         coeffs_.readIfPresent("muName", muName_);
 
         // reset the porosity model?

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,9 +91,6 @@ protected:
         //- Velocity field name, default = U
         word UName_;
 
-        //- Density field name (compressible case only), default = rho
-        word rhoName_;
-
         //- Dynamic viscosity field name (compressible case only)
         //  default = thermo:mu
         word muName_;
@@ -154,6 +151,17 @@ public:
             );
 
 
+        // Add explicit and implicit contributions to compressible equation
+
+            //- Vector
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
+
+
         // I-O
 
             //- Write data

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C

@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "interRegionHeatTransferModel.H"
-#include "fluidThermo.H"
+#include "basicThermo.H"
 #include "fvmSup.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fvcVolumeIntegrate.H"
@@ -216,7 +216,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup
     {
         if (he.dimensions() == dimEnergy/dimMass)
         {
-            if (mesh_.foundObject<fluidThermo>("thermophysicalProperties"))
+            if (mesh_.foundObject<basicThermo>("thermophysicalProperties"))
             {
                 const basicThermo& thermo =
                    mesh_.lookupObject<basicThermo>("thermophysicalProperties");
@@ -237,7 +237,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup
             }
             else
             {
-                 FatalErrorIn
+                FatalErrorIn
                 (
                     "void Foam::fv::interRegionHeatTransferModel::addSup"
                     "("
@@ -245,11 +245,9 @@ void Foam::fv::interRegionHeatTransferModel::addSup
                     "   const label "
                     ")"
                 )   << " on mesh " << mesh_.name()
-                    << " could not find object fluidThermo."
-                    << " The available objects : "
+                    << " could not find object basicThermo."
+                    << " The available objects are: "
                     << mesh_.names()
-                    << " The semi implicit option can only be used for "
-                    << "fluid-fluid inter region heat transfer models "
                     << exit(FatalError);
             }
         }
@@ -265,12 +263,23 @@ void Foam::fv::interRegionHeatTransferModel::addSup
 }
 
 
+void Foam::fv::interRegionHeatTransferModel::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldI
+)
+{
+    addSup(eqn, fieldI);
+}
+
+
 void Foam::fv::interRegionHeatTransferModel::writeData(Ostream& os) const
 {
     os.writeKeyword("name") << this->name() << token::END_STATEMENT << nl;
     os.writeKeyword("nbrRegionName") << nbrRegionName_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("nbrModeleName") << nbrModelName_
+    os.writeKeyword("nbrModelName") << nbrModelName_
         << token::END_STATEMENT << nl;
     os.writeKeyword("master") << master_ << token::END_STATEMENT << nl;
     os.writeKeyword("semiImplicit") << semiImplicit_ << token::END_STATEMENT

src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -157,7 +157,6 @@ public:
 
     // Member Functions
 
-
         // Access
 
             //- Return the heat transfer coefficient
@@ -169,8 +168,20 @@ public:
             //- Return access to the neighbour model
             inline interRegionHeatTransferModel& nbrModel();
 
-            //-Source term to fvMatrix<scalar>
-            virtual void addSup(fvMatrix<scalar>& eqn, const label fieldI);
+            //- Source term to energy equation
+            virtual void addSup
+            (
+                fvMatrix<scalar>& eqn,
+                const label fieldI
+            );
+
+            //- Source term to compressible energy equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn,
+                const label fieldI
+            );
 
             //- Calculate heat transfer coefficient
             virtual void calculateHtc() = 0;

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -285,7 +285,7 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
         SprayParcel<ParcelType>* child = new SprayParcel<ParcelType>(*this);
         child->mass0() = massChild;
         child->d() = dChild;
-        child->nParticle() = massChild/this->rho()*this->volume(dChild);
+        child->nParticle() = massChild/(this->rho()*this->volume(dChild));
 
         const forceSuSp Fcp =
             forces.calcCoupled(*child, dt, massChild, Re, muAv);

src/mesh/extrudeModel/linearNormal/linearNormal.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,7 +101,14 @@ point linearNormal::operator()
     const label layer
 ) const
 {
-    return surfacePoint + layerPoints_[layer - 1]*surfaceNormal;
+    if (layer == 0)
+    {
+        return surfacePoint;
+    }
+    else
+    {
+        return surfacePoint + layerPoints_[layer - 1]*surfaceNormal;
+    }
 }
 
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C

@@ -127,6 +127,10 @@ void Foam::cyclicACMIPolyPatch::setNeighbourFaceAreas() const
 
 void Foam::cyclicACMIPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
+    // initialise the AMI so that base geometry (e.g. cell volumes) are
+    // correctly evaluated
+    resetAMI();
+
     cyclicAMIPolyPatch::initGeometry(pBufs);
 }
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchTemplates.C

@@ -32,13 +32,16 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIPolyPatch::interpolate
     const Field<Type>& fldNonOverlap
 ) const
 {
+    // note: do not scale AMI field as face areas have already been taken
+    // into account
+
     if (owner())
     {
         const scalarField& w = srcMask_;
 
         tmp<Field<Type> > interpField(AMI().interpolateToSource(fldCouple));
 
-        return w*interpField + (1.0 - w)*fldNonOverlap;
+        return interpField + (1.0 - w)*fldNonOverlap;
     }
     else
     {
@@ -49,7 +52,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIPolyPatch::interpolate
             neighbPatch().AMI().interpolateToTarget(fldCouple)
         );
 
-        return w*interpField + (1.0 - w)*fldNonOverlap;
+        return interpField + (1.0 - w)*fldNonOverlap;
     }
 }
 
@@ -74,19 +77,22 @@ void Foam::cyclicACMIPolyPatch::interpolate
     List<Type>& result
 ) const
 {
+    // note: do not scale AMI field as face areas have already been taken
+    // into account
+
     if (owner())
     {
         const scalarField& w = srcMask_;
 
         AMI().interpolateToSource(fldCouple, cop, result);
-        result = w*result + (1.0 - w)*fldNonOverlap;
+        result = result + (1.0 - w)*fldNonOverlap;
     }
     else
     {
         const scalarField& w = neighbPatch().tgtMask();
 
         neighbPatch().AMI().interpolateToTarget(fldCouple, cop, result);
-        result = w*result + (1.0 - w)*fldNonOverlap;
+        result = result + (1.0 - w)*fldNonOverlap;
     }
 }
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -581,11 +581,11 @@ kinematicSingleLayer::kinematicSingleLayer
             "phi",
             time().timeName(),
             regionMesh(),
-            IOobject::READ_IF_PRESENT,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         regionMesh(),
-        dimLength*dimMass/dimTime
+        dimensionedScalar("0", dimLength*dimMass/dimTime, 0.0)
     ),
 
     primaryMassTrans_
@@ -795,7 +795,22 @@ kinematicSingleLayer::kinematicSingleLayer
         correctThermoFields();
 
         deltaRho_ == delta_*rho_;
-        phi_ = fvc::interpolate(deltaRho_*U_) & regionMesh().Sf();
+
+        surfaceScalarField phi0
+        (
+            IOobject
+            (
+                "phi",
+                time().timeName(),
+                regionMesh(),
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            fvc::interpolate(deltaRho_*U_) & regionMesh().Sf()
+        );
+
+        phi_ == phi0;
     }
 }
 

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -554,8 +554,24 @@ thermoSingleLayer::thermoSingleLayer
 
         // Update derived fields
         hs_ == hs(T_);
+
         deltaRho_ == delta_*rho_;
-        phi_ = fvc::interpolate(deltaRho_*U_) & regionMesh().Sf();
+
+        surfaceScalarField phi0
+        (
+            IOobject
+            (
+                "phi",
+                time().timeName(),
+                regionMesh(),
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            fvc::interpolate(deltaRho_*U_) & regionMesh().Sf()
+        );
+
+        phi_ == phi0;
 
         // evaluate viscosity from user-model
         viscosity_->correct(pPrimary_, T_);

tutorials/multiphase/LTSInterFoam/DTCHull/README

@@ -0,0 +1,6 @@
+Reference:
+Duisburg Test Case (DTC)
+el Moctar, O., Shigunov, V., Zorn, T.,
+Duisburg Test Case: Post-Panamax Container Ship for Benchmarking,
+Journal of Ship Technology Research,
+Vol.59, No.3, pp. 50-65, 2012

tutorials/multiphase/interDyMFoam/ras/DTCHull/README

@@ -0,0 +1,6 @@
+Reference:
+Duisburg Test Case (DTC)
+el Moctar, O., Shigunov, V., Zorn, T.,
+Duisburg Test Case: Post-Panamax Container Ship for Benchmarking,
+Journal of Ship Technology Research,
+Vol.59, No.3, pp. 50-65, 2012

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/0.org/U

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    freeSurface
+    {
+        type            pressureInletOutletParSlipVelocity;
+        value           uniform (0 0 0);
+    }
+
+    ".*Wall"
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    "floatingObject.*"
+    {
+        type            movingWallVelocity;
+        value           uniform (0 0 0);
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/0.org/p

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    freeSurface
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    ".*Wall"
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    "floatingObject.*"
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/0.org/p_gh

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p_gh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    freeSurface
+    {
+        type            waveSurfacePressure;
+        value           uniform 0;
+    }
+    ".*Wall"
+    {
+        type            zeroGradient;
+    }
+    "floatingObject.*"
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/0.org/pointDisplacement

@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       pointVectorField;
+    location    "0";
+    object      pointMotionU;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    freeSurface
+    {
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
+    }
+
+    leftWall
+    {
+        type            fixedNormalSlip;
+        n               (-1 0 0);
+    }
+
+    rightWall
+    {
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
+    }
+
+    bottomWall
+    {
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
+    }
+
+    floatingObject
+    {
+        type            fixedNormalSlip;
+        n               (1 0 0);
+    }
+
+    floatingObjectBottom
+    {
+        type            oscillatingDisplacement;
+        amplitude       (0 -0.05 0);
+        omega           6.2831853;
+        value           uniform (0 0 0);
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/Allclean

@@ -0,0 +1,10 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+rm -rf 0
+rm -rf poolHeight poolHeight_vs_time

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/Allrun

@@ -0,0 +1,21 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Set application name
+application=`getApplication`
+
+runApplication blockMesh
+runApplication topoSet
+mv log.topoSet log.topoSet.1
+runApplication subsetMesh -overwrite c0 -patch floatingObject
+runApplication topoSet -dict system/topoSetDict-selectBottom
+runApplication createPatch -overwrite
+
+cp -r 0.org 0 > /dev/null 2>&1
+
+runApplication $application
+
+# ----------------------------------------------------------------- end-of-file

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/constant/dynamicMeshDict

@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      dynamicMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh dynamicMotionSolverFvMesh;
+
+motionSolverLibs ("libfvMotionSolvers.so");
+
+solver displacementLaplacian;
+//solver velocityComponentLaplacian z;
+
+displacementLaplacianCoeffs
+{
+    // diffusivity  uniform;
+    // diffusivity  directional (1 200 0);
+    // diffusivity  motionDirectional (1 1000 0);
+    // diffusivity  file motionDiffusivity;
+    diffusivity  inverseDistance 1(floatingObjectBottom);
+}
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0 );
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/constant/polyMesh/blockMeshDict

@@ -0,0 +1,96 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    ( 0 0 0)
+    (10 0 0)
+    (10 1 0)
+    ( 0 1 0)
+    ( 0 0 0.1)
+    (10 0 0.1)
+    (10 1 0.1)
+    ( 0 1 0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (200 20 1) simpleGrading (10 0.1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    freeSurface
+    {
+        type wall;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    leftWall
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    rightWall
+    {
+        type wall;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    bottomWall
+    {
+        type wall;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    floatingObject
+    {
+        type wall;
+        faces
+        (
+        );
+    }
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/constant/polyMesh/boundary

@@ -0,0 +1,71 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+7
+(
+    freeSurface
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          181;
+        startFace       7429;
+    }
+    leftWall
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          20;
+        startFace       7610;
+    }
+    rightWall
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          20;
+        startFace       7630;
+    }
+    bottomWall
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          200;
+        startFace       7650;
+    }
+    floatingObject
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          18;
+        startFace       7850;
+    }
+    frontAndBack
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          7658;
+        startFace       7868;
+    }
+    floatingObjectBottom
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          19;
+        startFace       15526;
+    }
+)
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/constant/transportProperties

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [0 2 -1 0 0 0 0] 1e-06;
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/constant/turbulenceProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/extractHeightData

@@ -0,0 +1,2 @@
+#!/bin/sh
+awk '{print $1 " " $4}' postProcessing/poolHeight/0/faceSource.dat > poolHeight_vs_time

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/system/controlDict

@@ -0,0 +1,89 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     potentialFreeSurfaceDyMFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         20;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.02;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.4;
+
+maxDeltaT       1;
+
+functions
+{
+    forces
+    {
+        type            forces;
+        functionObjectLibs ("libforces.so");
+        outputControl   outputTime;
+        patches         (floatingObject);
+        pName           p;
+        UName           U;
+        rhoName         rhoInf;
+        log             yes;
+        rhoInf          1000;
+        CofR            (0 0 0);
+    }
+
+    poolHeight
+    {
+        type            faceSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         yes;
+        outputControl   timeStep;
+        outputInterval  1;
+        log             yes;
+        valueOutput     no;
+        source          faceZone;
+        sourceName      f0;
+        operation       areaAverage;
+        fields
+        (
+            zeta
+        );
+    }
+};
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/system/createPatchDict

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createPatchDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Do a synchronisation of coupled points after creation of any patches.
+// Note: this does not work with points that are on multiple coupled patches
+//       with transformations (i.e. cyclics).
+pointSync false;
+
+// Patches to create.
+patches
+(
+    {
+        // Name of new patch
+        name floatingObjectBottom;
+
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type wall;
+        }
+
+        // How to construct: either from 'patches' or 'set'
+        constructFrom set;
+
+        // If constructFrom = set : name of faceSet
+        set floatingObjectBottom;
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/system/fvSchemes

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss upwind;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p_gh;
+    pcorr;
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/system/fvSolution

@@ -0,0 +1,85 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "cellDisplacement.*"
+    {
+        solver          GAMG;
+        tolerance       1e-5;
+        relTol          0;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    "pcorr.*"
+    {
+        solver          GAMG;
+        tolerance       1e-2;
+        relTol          0;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+        maxIter         100;
+    }
+
+    p_gh
+    {
+        $pcorr;
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+
+    p_ghFinal
+    {
+        $p_gh;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+
+    UFinal
+    {
+        $U;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    correctPhi          yes;
+    momentumPredictor   no;
+    nOuterCorrectors    1;
+    nCorrectors         2;
+    nNonOrthogonalCorrectors 0;
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/system/topoSetDict

@@ -0,0 +1,71 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    c0;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (0.1 0.8 -100) (0.4 100 100);
+        }
+    }
+
+    {
+        name    c0;
+        type    cellSet;
+        action  invert;
+    }
+
+    {
+        name    f0;
+        type    faceSet;
+        action  new;
+        source  patchToFace;
+        sourceInfo
+        {
+            name   freeSurface;
+        }
+    }
+
+    {
+        name    f0;
+        type    faceSet;
+        action  subset;
+        source  boxToFace;
+        sourceInfo
+        {
+            box (-100 0.9 -100) (0.2 100 100);
+        }
+    }
+
+    {
+        name    f0;
+        type    faceZoneSet;
+        action  new;
+        source  setToFaceZone;
+        sourceInfo
+        {
+            faceSet    f0;
+        }
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/system/topoSetDict-selectBottom

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    floatingObjectBottom;
+        type    faceSet;
+        action  new;
+        source  patchToFace;
+        sourceInfo
+        {
+            name   floatingObject;
+        }
+    }
+
+    {
+        name    floatingObjectBottom;
+        type    faceSet;
+        action  subset;
+        source  normalToFace;
+        sourceInfo
+        {
+            normal (0 1 0);     // Vector
+            cos     0.01;       // Tolerance (max cos of angle)
+        }
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/T.air

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Tair;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.air;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    walls
+    {
+        type               zeroGradient;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/T.water

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Twater;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.water;
+        inletValue         uniform 300;
+        value              $internalField;
+    }
+    walls
+    {
+        type               zeroGradient;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/U.air

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.air;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/U.water

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.water;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //