fvOptions: Correct handling of density and add multiphase support

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/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.H

@@ -132,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
             );
@@ -155,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")),
+    UName_(coeffs_.lookupOrDefault<word>("UName", "U"))
-    rhoName_(coeffs_.lookupOrDefault<word>("rhoName", "rho"))
 {
     initialise();
 }
@@ -104,19 +103,20 @@ void Foam::fv::MRFSource::addSup
     const label fieldI
 )
 {
-    if (eqn.dimensions() == dimForce)
+    // Add to rhs of equation
-    {
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
-
-        // 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,30 +115,11 @@ 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
         (
@@ -150,6 +131,30 @@ void Foam::fv::actuationDiskSource::addSup
         );
     }
 }
+
+
+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>
+            //- Source term to momentum equation
-            virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
+            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>
+            //- Scalar
-            virtual void addSup(fvMatrix<scalar>& eqn, const label fieldI);
+            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,19 +102,24 @@ 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, 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
+            //- Add explicit contribution to momentum equation
-            virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
+            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,30 +68,11 @@ 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
         (
@@ -103,6 +84,30 @@ void Foam::fv::radialActuationDiskSource::addSup
         );
     }
 }
+
+
+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>
+        //- Source term to momentum equation
-        virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
+        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

@@ -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
@@ -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;
@@ -435,7 +434,6 @@ Foam::fv::rotorDiskSource::rotorDiskSource
 )
 :
     option(name, modelType, dict, mesh),
-    rhoName_("none"),
     rhoRef_(1.0),
     omega_(0.0),
     nBlades_(0),
@@ -465,8 +463,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,
@@ -475,8 +475,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;
@@ -545,11 +543,7 @@ void Foam::fv::rotorDiskSource::calculate
             scalar tipFactor = neg(radius/rMax_ - tipEffect_);
 
             // calculate forces perpendicular to blade
-            scalar pDyn = 0.5*magSqr(Uc);
+            scalar pDyn = 0.5*rho[cellI]*magSqr(Uc);
-            if (compressible)
-            {
-                pDyn *= trho()[cellI];
-            }
 
             scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;
             vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl);
@@ -571,7 +565,6 @@ void Foam::fv::rotorDiskSource::calculate
         }
     }
 
-
     if (output)
     {
         reduce(AOAmin, minOp<scalar>());
@@ -594,42 +587,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_,
+            mesh_
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedVector("zero", dims, vector::zero)
+        dimensionedVector
+        (
+            "zero",
+            eqn.dimensions()/dimVolume,
+            vector::zero
+        )
     );
 
-    const volVectorField& U = eqn.psi();
+    // Read the reference density for incompressible flow
-
+    coeffs_.lookup("rhoRef") >> rhoRef_;
-    const vectorField Uin(inflowVelocity(U));
 
+    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())
@@ -653,7 +673,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

@@ -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,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
 
@@ -95,7 +94,6 @@ SourceFiles
 #include "bladeModel.H"
 #include "profileModelList.H"
 #include "volFieldsFwd.H"
-#include "dimensionSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -149,10 +147,7 @@ protected:
 
     // Protected data
 
-        //- Name of density field
+        //- Reference density for incompressible case
-        word rhoName_;
-
-        //- Reference density for rhoName = 'none'
         scalar rhoRef_;
 
         //- Rotational speed [rad/s]
@@ -258,7 +253,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]
@@ -272,18 +267,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,
@@ -294,8 +285,20 @@ public:
 
         // Source term addition
 
-            //- Source term to fvMatrix<vector>
+            //- Source term to momentum equation
-            virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
+            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);
+    rotor_.calculate(rho, U, thetag, force, false, false);
-
-    bool compressible = rotor_.compressible();
-    tmp<volScalarField> trho = rotor_.rho();
 
     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);
+            scalar coeff2 = rho[cellI]*coeff1*pow4(radius);
-            if (compressible)
-            {
-                coeff2 *= trho()[cellI];
-            }
 
             // 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)
+    correctTrim(geometricOneField(), U, force);
-    {
-        word calcType = "forces";
-        if (useCoeffs_)
-        {
-            calcType = "coefficients";
 }
 
-        Info<< type() << ":" << nl
-            << "    solving for target trim " << calcType << nl;
 
-        const scalar rhoRef = rotor_.rhoRef();
+void Foam::targetCoeffTrim::correct
-
+(
-        // iterate to find new pitch angles to achieve target force
+    const volScalarField rho,
-        scalar err = GREAT;
+    const vectorField& U,
-        label iter = 0;
+    vectorField& force
-        tensor J(tensor::zero);
+)
-
-        vector old = vector::zero;
-        while ((err > tol_) && (iter < nIter_))
 {
-            // cache initial theta vector
+    correctTrim(rho, U, force);
-            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;
-    }
 }
 
 

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

@@ -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,8 +170,60 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
 
     fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
 
-    if (eqn.dimensions() == dimForce)
+    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
+        (
+            name_ + ":UNbr",
+            nbrMesh.time().timeName(),
+            nbrMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        nbrMesh,
+        dimensionedVector("zero", U.dimensions(), vector::zero)
+    );
+
+    // map local velocity onto neighbour region
+    meshInterp().mapSrcToTgt
+    (
+        U.internalField(),
+        plusEqOp<vector>(),
+        UNbr.internalField()
+    );
+
+    fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
+
     volScalarField rhoNbr
     (
         IOobject
@@ -201,9 +252,6 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
         dimensionedScalar("zero", dimViscosity, 0.0)
     );
 
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
-
     const volScalarField& mu =
         mesh_.lookupObject<volScalarField>(muName_);
 
@@ -224,11 +272,6 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
     );
 
     porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);
-    }
-    else
-    {
-        porosityPtr_->addResistance(nbrEqn);
-    }
 
     // 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");
@@ -245,11 +245,9 @@ void Foam::fv::interRegionHeatTransferModel::addSup
                     "   const label "
                     ")"
                 )   << " on mesh " << mesh_.name()
-                    << " could not find object fluidThermo."
+                    << " could not find object basicThermo."
-                    << " The available objects : "
+                    << " 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>
+            //- Source term to energy equation
-            virtual void addSup(fvMatrix<scalar>& eqn, const label fieldI);
+            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;