ENH: relative velocity usage for porosity in MRF regions

- New optional keyword 'MRFRelativeVelocity' (default=true)
usable by the porosity models. For example, in 'constant/porosityProperties':

    porosity1
    {
        type            DarcyForchheimer;
        ...

        MRFRelativeVelocity    yes;
        ...
    }

src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2018 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -314,12 +314,6 @@ Foam::porosityModels::powerLawLopesdaCosta::powerLawLopesdaCosta
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::porosityModels::powerLawLopesdaCosta::~powerLawLopesdaCosta()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::scalarField&
@@ -344,47 +338,47 @@ void Foam::porosityModels::powerLawLopesdaCosta::calcForce
     scalarField Udiag(U.size(), Zero);
     const scalarField& V = mesh_.V();
 
-    apply(Udiag, V, rho, U);
+    apply(V, rho, U, Udiag);
 
-    force = Udiag*U;
+    force += Udiag*U;
 }
 
 
 void Foam::porosityModels::powerLawLopesdaCosta::correct
 (
-    fvVectorMatrix& UEqn
+    const volVectorField& U,
+    scalarField& Udiag,
+    vectorField& Usource,
+    bool compressible
 ) const
 {
-    const vectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
 
-    if (UEqn.dimensions() == dimForce)
+    if (compressible)
     {
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
+        const auto& rho = mesh_.lookupObject<volScalarField>(rhoName_);
 
-        apply(Udiag, V, rho, U);
+        apply(V, rho, U, Udiag);
     }
     else
     {
-        apply(Udiag, V, geometricOneField(), U);
+        apply(V, geometricOneField(), U, Udiag);
     }
 }
 
 
 void Foam::porosityModels::powerLawLopesdaCosta::correct
 (
-    fvVectorMatrix& UEqn,
+    const volVectorField& U,
     const volScalarField& rho,
-    const volScalarField& mu
+    const volScalarField& mu,
+    scalarField& Udiag,
+    vectorField& Usource
 ) const
 {
-    const vectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
 
-    apply(Udiag, V, rho, U);
+    apply(V, rho, U, Udiag);
 }
 
 
@@ -394,18 +388,17 @@ void Foam::porosityModels::powerLawLopesdaCosta::correct
     volTensorField& AU
 ) const
 {
-    const vectorField& U = UEqn.psi();
+    const volVectorField& U = UEqn.psi();
 
     if (UEqn.dimensions() == dimForce)
     {
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
+        const auto& rho = mesh_.lookupObject<volScalarField>(rhoName_);
 
-        apply(AU, rho, U);
+        apply(rho, U, AU);
     }
     else
     {
-        apply(AU, geometricOneField(), U);
+        apply(geometricOneField(), U, AU);
     }
 }
 

src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -133,19 +134,19 @@ class powerLawLopesdaCosta
         template<class RhoFieldType>
         void apply
         (
-            scalarField& Udiag,
             const scalarField& V,
             const RhoFieldType& rho,
-            const vectorField& U
+            const vectorField& U,
+            scalarField& Udiag
         ) const;
 
         //- Apply resistance
         template<class RhoFieldType>
         void apply
         (
-            tensorField& AU,
             const RhoFieldType& rho,
-            const vectorField& U
+            const vectorField& U,
+            tensorField& AU
         ) const;
 
 
@@ -172,7 +173,7 @@ public:
     );
 
     //- Destructor
-    virtual ~powerLawLopesdaCosta();
+    virtual ~powerLawLopesdaCosta() = default;
 
 
     // Member Functions
@@ -190,14 +191,22 @@ public:
         ) const;
 
         //- Add resistance
-        virtual void correct(fvVectorMatrix& UEqn) const;
+        virtual void correct
+        (
+            const volVectorField& U,
+            scalarField& Udiag,
+            vectorField& Usource,
+            bool compressible
+        ) const;
 
         //- Add resistance
         virtual void correct
         (
-            fvVectorMatrix& UEqn,
+            const volVectorField& U,
             const volScalarField& rho,
-            const volScalarField& mu
+            const volScalarField& mu,
+            scalarField& Udiag,
+            vectorField& Usource
         ) const;
 
         //- Add resistance

src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCostaTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,10 +33,10 @@ License
 template<class RhoFieldType>
 void Foam::porosityModels::powerLawLopesdaCosta::apply
 (
-    scalarField& Udiag,
     const scalarField& V,
     const RhoFieldType& rho,
-    const vectorField& U
+    const vectorField& U,
+    scalarField& Udiag
 ) const
 {
     const scalar C1m1b2 = (C1_ - 1.0)/2.0;
@@ -59,9 +60,9 @@ void Foam::porosityModels::powerLawLopesdaCosta::apply
 template<class RhoFieldType>
 void Foam::porosityModels::powerLawLopesdaCosta::apply
 (
-    tensorField& AU,
     const RhoFieldType& rho,
-    const vectorField& U
+    const vectorField& U,
+    tensorField& AU
 ) const
 {
     const scalar C1m1b2 = (C1_ - 1.0)/2.0;

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -175,79 +175,79 @@ void Foam::porosityModels::DarcyForchheimer::calcForce
     vectorField& force
 ) const
 {
-    scalarField Udiag(U.size(), Zero);
-    vectorField Usource(U.size(), Zero);
     const scalarField& V = mesh_.V();
 
-    apply(Udiag, Usource, V, rho, mu, U);
+    scalarField Udiag(U.size(), Zero);
+    vectorField Usource(U.size(), Zero);
+    apply(V, rho, mu, U, Udiag, Usource);
 
-    force = Udiag*U - Usource;
+    force += Udiag*U - Usource;
 }
 
 
 void Foam::porosityModels::DarcyForchheimer::correct
 (
-    fvVectorMatrix& UEqn
+    const volVectorField& U,
+    scalarField& Udiag,
+    vectorField& Usource,
+    bool compressible
 ) const
 {
-    const volVectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-    vectorField& Usource = UEqn.source();
 
-    word rhoName(IOobject::groupName(rhoName_, U.group()));
-    word muName(IOobject::groupName(muName_, U.group()));
-    word nuName(IOobject::groupName(nuName_, U.group()));
+    const word rhoName(IOobject::groupName(rhoName_, U.group()));
+    const word muName(IOobject::groupName(muName_, U.group()));
+    const word nuName(IOobject::groupName(nuName_, U.group()));
 
-    if (UEqn.dimensions() == dimForce)
+    if (compressible)
     {
         const auto& rho = mesh_.lookupObject<volScalarField>(rhoName);
 
+        tmp<volScalarField> tmu;
         if (mesh_.foundObject<volScalarField>(muName))
         {
-            const auto& mu = mesh_.lookupObject<volScalarField>(muName);
-
-            apply(Udiag, Usource, V, rho, mu, U);
+            tmu = mesh_.lookupObject<volScalarField>(muName);
         }
         else
         {
             const auto& nu = mesh_.lookupObject<volScalarField>(nuName);
-
-            apply(Udiag, Usource, V, rho, rho*nu, U);
+            tmu = rho*nu;
         }
+
+        apply(V, rho, tmu(), U, Udiag, Usource);
     }
     else
     {
+        tmp<volScalarField> tnu;
+
         if (mesh_.foundObject<volScalarField>(nuName))
         {
-            const auto& nu = mesh_.lookupObject<volScalarField>(nuName);
-
-            apply(Udiag, Usource, V, geometricOneField(), nu, U);
+            tnu = mesh_.lookupObject<volScalarField>(nuName);
         }
         else
         {
             const auto& rho = mesh_.lookupObject<volScalarField>(rhoName);
             const auto& mu = mesh_.lookupObject<volScalarField>(muName);
-
-            apply(Udiag, Usource, V, geometricOneField(), mu/rho, U);
+            tnu = mu/rho;
         }
+
+        apply(V, geometricOneField(), tnu(), U, Udiag, Usource);
     }
 }
 
 
 void Foam::porosityModels::DarcyForchheimer::correct
 (
-    fvVectorMatrix& UEqn,
+    const volVectorField& U,
     const volScalarField& rho,
-    const volScalarField& mu
+    const volScalarField& mu,
+    scalarField& Udiag,
+    vectorField& Usource
 ) const
 {
-    const vectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-    vectorField& Usource = UEqn.source();
 
-    apply(Udiag, Usource, V, rho, mu, U);
+    apply(V, rho, mu, U, Udiag, Usource);
 }
 
 
@@ -259,16 +259,16 @@ void Foam::porosityModels::DarcyForchheimer::correct
 {
     const volVectorField& U = UEqn.psi();
 
-    word rhoName(IOobject::groupName(rhoName_, U.group()));
-    word muName(IOobject::groupName(muName_, U.group()));
-    word nuName(IOobject::groupName(nuName_, U.group()));
+    const word rhoName(IOobject::groupName(rhoName_, U.group()));
+    const word muName(IOobject::groupName(muName_, U.group()));
+    const word nuName(IOobject::groupName(nuName_, U.group()));
 
     if (UEqn.dimensions() == dimForce)
     {
         const auto& rho = mesh_.lookupObject<volScalarField>(rhoName);
         const auto& mu = mesh_.lookupObject<volScalarField>(muName);
 
-        apply(AU, rho, mu, U);
+        apply(rho, mu, U, AU);
     }
     else
     {
@@ -276,14 +276,14 @@ void Foam::porosityModels::DarcyForchheimer::correct
         {
             const auto& nu = mesh_.lookupObject<volScalarField>(nuName);
 
-            apply(AU, geometricOneField(), nu, U);
+            apply(geometricOneField(), nu, U, AU);
         }
         else
         {
             const auto& rho = mesh_.lookupObject<volScalarField>(rhoName);
             const auto& mu = mesh_.lookupObject<volScalarField>(muName);
 
-            apply(AU, geometricOneField(), mu/rho, U);
+            apply(geometricOneField(), mu/rho, U, AU);
         }
     }
 }

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -102,22 +103,22 @@ class DarcyForchheimer
         template<class RhoFieldType>
         void apply
         (
-            scalarField& Udiag,
-            vectorField& Usource,
             const scalarField& V,
             const RhoFieldType& rho,
             const scalarField& mu,
-            const vectorField& U
+            const vectorField& U,
+            scalarField& Udiag,
+            vectorField& Usource
         ) const;
 
         //- Apply
         template<class RhoFieldType>
         void apply
         (
-            tensorField& AU,
             const RhoFieldType& rho,
             const scalarField& mu,
-            const vectorField& U
+            const vectorField& U,
+            tensorField& AU
         ) const;
 
         //- No copy construct
@@ -161,14 +162,22 @@ public:
         ) const;
 
         //- Add resistance
-        virtual void correct(fvVectorMatrix& UEqn) const;
+        virtual void correct
+        (
+            const volVectorField& U,
+            scalarField& Udiag,
+            vectorField& Usource,
+            bool compressible
+        ) const;
 
         //- Add resistance
         virtual void correct
         (
-            fvVectorMatrix& UEqn,
+            const volVectorField& U,
             const volScalarField& rho,
-            const volScalarField& mu
+            const volScalarField& mu,
+            scalarField& Udiag,
+            vectorField& Usource
         ) const;
 
         //- Add resistance

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,12 +31,12 @@ License
 template<class RhoFieldType>
 void Foam::porosityModels::DarcyForchheimer::apply
 (
-    scalarField& Udiag,
-    vectorField& Usource,
     const scalarField& V,
     const RhoFieldType& rho,
     const scalarField& mu,
-    const vectorField& U
+    const vectorField& U,
+    scalarField& Udiag,
+    vectorField& Usource
 ) const
 {
     forAll(cellZoneIDs_, zoneI)
@@ -64,10 +65,10 @@ void Foam::porosityModels::DarcyForchheimer::apply
 template<class RhoFieldType>
 void Foam::porosityModels::DarcyForchheimer::apply
 (
-    tensorField& AU,
     const RhoFieldType& rho,
     const scalarField& mu,
-    const vectorField& U
+    const vectorField& U,
+    tensorField& AU
 ) const
 {
     forAll(cellZoneIDs_, zoneI)

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,11 +47,11 @@ namespace Foam
 
 void Foam::porosityModels::fixedCoeff::apply
 (
-    scalarField& Udiag,
-    vectorField& Usource,
     const scalarField& V,
     const vectorField& U,
-    const scalar rho
+    const scalar rho,
+    scalarField& Udiag,
+    vectorField& Usource
 ) const
 {
     forAll(cellZoneIDs_, zoneI)
@@ -77,9 +77,9 @@ void Foam::porosityModels::fixedCoeff::apply
 
 void Foam::porosityModels::fixedCoeff::apply
 (
-    tensorField& AU,
     const vectorField& U,
-    const scalar rho
+    const scalar rho,
+    tensorField& AU
 ) const
 {
     forAll(cellZoneIDs_, zoneI)
@@ -183,51 +183,46 @@ void Foam::porosityModels::fixedCoeff::calcForce
     const scalarField& V = mesh_.V();
     const scalar rhoRef = coeffs_.get<scalar>("rhoRef");
 
-    apply(Udiag, Usource, V, U, rhoRef);
+    apply(V, U, rhoRef, Udiag, Usource);
 
-    force = Udiag*U - Usource;
+    force += Udiag*U - Usource;
 }
 
 
 void Foam::porosityModels::fixedCoeff::correct
 (
-    fvVectorMatrix& UEqn
+    const volVectorField& U,
+    scalarField& Udiag,
+    vectorField& Usource,
+    bool compressible
 ) const
 {
-    const vectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-    vectorField& Usource = UEqn.source();
 
     scalar rho = 1.0;
-    if (UEqn.dimensions() == dimForce)
+    if (compressible)
     {
         coeffs_.readEntry("rhoRef", rho);
     }
 
-    apply(Udiag, Usource, V, U, rho);
+    apply(V, U, rho, Udiag, Usource);
 }
 
 
 void Foam::porosityModels::fixedCoeff::correct
 (
-    fvVectorMatrix& UEqn,
+    const volVectorField& U,
     const volScalarField&,
-    const volScalarField&
+    const volScalarField&,
+    scalarField& Udiag,
+    vectorField& Usource
 ) const
 {
-    const vectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-    vectorField& Usource = UEqn.source();
 
-    scalar rho = 1.0;
-    if (UEqn.dimensions() == dimForce)
-    {
-        coeffs_.readEntry("rhoRef", rho);
-    }
+    const scalar rho = coeffs_.get<scalar>("rhoRef");
 
-    apply(Udiag, Usource, V, U, rho);
+    apply(V, U, rho, Udiag, Usource);
 }
 
 
@@ -237,15 +232,13 @@ void Foam::porosityModels::fixedCoeff::correct
     volTensorField& AU
 ) const
 {
-    const vectorField& U = UEqn.psi();
-
     scalar rho = 1.0;
     if (UEqn.dimensions() == dimForce)
     {
         coeffs_.readEntry("rhoRef", rho);
     }
 
-    apply(AU, U, rho);
+    apply(UEqn.psi(), rho, AU);
 }
 
 

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -82,19 +83,19 @@ class fixedCoeff
         //- Apply
         void apply
         (
-            scalarField& Udiag,
-            vectorField& Usource,
             const scalarField& V,
             const vectorField& U,
-            const scalar rho
+            const scalar rho,
+            scalarField& Udiag,
+            vectorField& Usource
         ) const;
 
         //- Apply
         void apply
         (
-            tensorField& AU,
             const vectorField& U,
-            const scalar rho
+            const scalar rho,
+            tensorField& AU
         ) const;
 
         //- No copy construct
@@ -138,14 +139,22 @@ public:
         ) const;
 
         //- Add resistance
-        virtual void correct(fvVectorMatrix& UEqn) const;
+        virtual void correct
+        (
+            const volVectorField& U,
+            scalarField& Udiag,
+            vectorField& Usource,
+            bool compressible
+        ) const;
 
         //- Add resistance
         virtual void correct
         (
-            fvVectorMatrix& UEqn,
+            const volVectorField& U,
             const volScalarField& rho,
-            const volScalarField& mu
+            const volScalarField& mu,
+            scalarField& Udiag,
+            vectorField& Usource
         ) const;
 
         //- Add resistance

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C

@@ -28,6 +28,7 @@ License
 
 #include "porosityModel.H"
 #include "volFields.H"
+#include "IOMRFZoneList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -97,7 +98,8 @@ Foam::porosityModel::porosityModel
     csysPtr_
     (
         coordinateSystem::New(mesh, coeffs_, coordinateSystem::typeName_())
-    )
+    ),
+    MRFRelativeVelocity_(true)
 {
     if (zoneName_.empty())
     {
@@ -146,12 +148,17 @@ Foam::porosityModel::porosityModel
 
         Info<< "    local bounds: " << bb.span() << nl << endl;
     }
+
+    if (dict.readIfPresent("MRFRelativeVelocity", MRFRelativeVelocity_))
+    {
+        Info<< "    MRFRelativeVelocity " << MRFRelativeVelocity_ << endl;
+    }
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::porosityModel::transformModelData()
+Foam::tmp<Foam::volVectorField> Foam::porosityModel::transformModelData()
 {
     if (!mesh_.upToDatePoints(*this))
     {
@@ -160,6 +167,32 @@ void Foam::porosityModel::transformModelData()
         // set model up-to-date wrt points
         mesh_.setUpToDatePoints(*this);
     }
+
+
+    auto* MRF = mesh_.findObject<IOMRFZoneList>("MRFProperties");
+    if (MRF && MRFRelativeVelocity_)
+    {
+        auto tUref =
+            tmp<volVectorField>::New
+            (
+                IOobject
+                (
+                    IOobject::scopedName(word("Uporous"), name_),
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedVector(dimVelocity, Zero)
+            );
+
+        MRF->makeRelative(tUref.ref());
+
+        return tUref;
+    }
+
+    return nullptr;
 }
 
 
@@ -170,13 +203,18 @@ Foam::tmp<Foam::vectorField> Foam::porosityModel::porosityModel::force
     const volScalarField& mu
 )
 {
-    transformModelData();
+    tmp<volVectorField> tUref = transformModelData();
 
-    tmp<vectorField> tforce(new vectorField(U.size(), Zero));
+    auto tforce = tmp<vectorField>::New(U.size(), Zero);
 
     if (!cellZoneIDs_.empty())
     {
         this->calcForce(U, rho, mu, tforce.ref());
+
+        if (tUref.valid())
+        {
+            this->calcForce(tUref(), rho, mu, tforce.ref());
+        }
     }
 
     return tforce;
@@ -190,16 +228,30 @@ void Foam::porosityModel::addResistance(fvVectorMatrix& UEqn)
         return;
     }
 
-    transformModelData();
-    this->correct(UEqn);
+    bool compressible = UEqn.dimensions() == dimForce;
+
+    tmp<volVectorField> tUref = transformModelData();
+
+    this->correct(UEqn.psi(), UEqn.diag(), UEqn.source(), compressible);
+
+    if (tUref.valid())
+    {
+        const auto& Uref = tUref();
+        scalarField Udiag(Uref.size(), Zero);
+        vectorField Usource(Uref.size(), Zero);
+
+        this->correct(Uref, Udiag, Usource, compressible);
+
+        UEqn.source() -= Udiag*Uref - Usource;
+    }
 }
 
 
 void Foam::porosityModel::addResistance
 (
-    fvVectorMatrix& UEqn,
     const volScalarField& rho,
-    const volScalarField& mu
+    const volScalarField& mu,
+    fvVectorMatrix& UEqn
 )
 {
     if (cellZoneIDs_.empty())
@@ -207,8 +259,20 @@ void Foam::porosityModel::addResistance
         return;
     }
 
-    transformModelData();
-    this->correct(UEqn, rho, mu);
+    tmp<volVectorField> tUref = transformModelData();
+
+    this->correct(UEqn.psi(), rho, mu, UEqn.diag(), UEqn.source());
+
+    if (tUref.valid())
+    {
+        const auto& Uref = tUref();
+        scalarField Udiag(Uref.size(), Zero);
+        vectorField Usource(Uref.size(), Zero);
+
+        this->correct(Uref, rho, mu, Udiag, Usource);
+
+        UEqn.source() -= Udiag*Uref - Usource;
+    }
 }
 
 
@@ -224,7 +288,10 @@ void Foam::porosityModel::addResistance
         return;
     }
 
-    transformModelData();
+    // Note: reference velocity contribution not handled!
+    // tmp<volVectorField> tUref = transformModelData();
+    (void) transformModelData();
+
     this->correct(UEqn, AU);
 
     if (correctAUprocBC)

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H

@@ -96,6 +96,9 @@ protected:
         //- Local coordinate system
         autoPtr<coordinateSystem> csysPtr_;
 
+        //- Flag to use relative velocity for MRF regions; default = yes
+        Switch MRFRelativeVelocity_;
+
 
     // Protected Member Functions
 
@@ -114,13 +117,21 @@ protected:
             vectorField& force
         ) const = 0;
 
-        virtual void correct(fvVectorMatrix& UEqn) const = 0;
+        virtual void correct
+        (
+            const volVectorField& U,
+            scalarField& Udiag,
+            vectorField& Usource,
+            bool compressible
+        ) const = 0;
 
         virtual void correct
         (
-            fvVectorMatrix& UEqn,
+            const volVectorField& U,
             const volScalarField& rho,
-            const volScalarField& mu
+            const volScalarField& mu,
+            scalarField& Udiag,
+            vectorField& Usource
         ) const = 0;
 
         virtual void correct
@@ -232,7 +243,8 @@ public:
         const dictionary& dict() const;
 
         //- Transform the model data wrt mesh changes
-        virtual void transformModelData();
+        //  Returns the reference velocity (if applicable)
+        virtual tmp<volVectorField> transformModelData();
 
         //- Return the force over the cell zone(s)
         virtual tmp<vectorField> force
@@ -248,9 +260,9 @@ public:
         //- Add resistance
         virtual void addResistance
         (
-            fvVectorMatrix& UEqn,
             const volScalarField& rho,
-            const volScalarField& mu
+            const volScalarField& mu,
+            fvVectorMatrix& UEqn
         );
 
         //- Add resistance

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2014 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -97,9 +98,8 @@ void Foam::porosityModelList::reset(const dictionary& dict)
 bool Foam::porosityModelList::read(const dictionary& dict)
 {
     bool allOk = true;
-    forAll(*this, i)
+    for (auto& pm : *this)
     {
-        porosityModel& pm = this->operator[](i);
         bool ok = pm.read(dict.subDict(pm.name()));
         allOk = (allOk && ok);
     }
@@ -109,24 +109,21 @@ bool Foam::porosityModelList::read(const dictionary& dict)
 
 bool Foam::porosityModelList::writeData(Ostream& os) const
 {
-    forAll(*this, i)
+    for (const auto& pm : *this)
     {
         os  << nl;
-        this->operator[](i).writeData(os);
+        pm.writeData(os);
     }
 
     return os.good();
 }
 
 
-void Foam::porosityModelList::addResistance
-(
-    fvVectorMatrix& UEqn
-)
+void Foam::porosityModelList::addResistance(fvVectorMatrix& UEqn)
 {
-    forAll(*this, i)
+    for (auto& pm : *this)
     {
-        this->operator[](i).addResistance(UEqn);
+        pm.addResistance(UEqn);
     }
 }
 
@@ -138,9 +135,9 @@ void Foam::porosityModelList::addResistance
     const volScalarField& mu
 )
 {
-    forAll(*this, i)
+    for (auto& pm : *this)
     {
-        this->operator[](i).addResistance(UEqn, rho, mu);
+        pm.addResistance(rho, mu, UEqn);
     }
 }
 
@@ -152,9 +149,9 @@ void Foam::porosityModelList::addResistance
     bool correctAUprocBC
 )
 {
-    forAll(*this, i)
+    for (auto& pm : *this)
     {
-        this->operator[](i).addResistance(UEqn, AU, correctAUprocBC);
+        pm.addResistance(UEqn, AU, correctAUprocBC);
     }
 }
 

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -80,49 +80,50 @@ void Foam::porosityModels::powerLaw::calcForce
     scalarField Udiag(U.size(), Zero);
     const scalarField& V = mesh_.V();
 
-    apply(Udiag, V, rho, U);
+    apply(V, rho, U, Udiag);
 
-    force = Udiag*U;
+    force += Udiag*U;
 }
 
 
 void Foam::porosityModels::powerLaw::correct
 (
-    fvVectorMatrix& UEqn
+    const volVectorField& U,
+    scalarField& Udiag,
+    vectorField& Usource,
+    bool compressible
 ) const
 {
-    const volVectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
 
-    if (UEqn.dimensions() == dimForce)
+    if (compressible)
     {
         const auto& rho = mesh_.lookupObject<volScalarField>
         (
             IOobject::groupName(rhoName_, U.group())
         );
 
-        apply(Udiag, V, rho, U);
+        apply(V, rho, U, Udiag);
     }
     else
     {
-        apply(Udiag, V, geometricOneField(), U);
+        apply(V, geometricOneField(), U, Udiag);
     }
 }
 
 
 void Foam::porosityModels::powerLaw::correct
 (
-    fvVectorMatrix& UEqn,
+    const volVectorField& U,
     const volScalarField& rho,
-    const volScalarField& mu
+    const volScalarField& mu,
+    scalarField& Udiag,
+    vectorField& Usource
 ) const
 {
-    const vectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
 
-    apply(Udiag, V, rho, U);
+    apply(V, rho, U, Udiag);
 }
 
 
@@ -141,11 +142,11 @@ void Foam::porosityModels::powerLaw::correct
             IOobject::groupName(rhoName_, U.group())
         );
 
-        apply(AU, rho, U);
+        apply(rho, U, AU);
     }
     else
     {
-        apply(AU, geometricOneField(), U);
+        apply(geometricOneField(), U, AU);
     }
 }
 

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,19 +84,19 @@ class powerLaw
         template<class RhoFieldType>
         void apply
         (
-            scalarField& Udiag,
             const scalarField& V,
             const RhoFieldType& rho,
-            const vectorField& U
+            const vectorField& U,
+            scalarField& Udiag
         ) const;
 
         //- Apply resistance
         template<class RhoFieldType>
         void apply
         (
-            tensorField& AU,
             const RhoFieldType& rho,
-            const vectorField& U
+            const vectorField& U,
+            tensorField& AU
         ) const;
 
         //- No copy construct
@@ -139,14 +140,22 @@ public:
         ) const;
 
         //- Add resistance
-        virtual void correct(fvVectorMatrix& UEqn) const;
+        virtual void correct
+        (
+            const volVectorField& U,
+            scalarField& Udiag,
+            vectorField& Usource,
+            bool compressible
+        ) const;
 
         //- Add resistance
         virtual void correct
         (
-            fvVectorMatrix& UEqn,
+            const volVectorField& U,
             const volScalarField& rho,
-            const volScalarField& mu
+            const volScalarField& mu,
+            scalarField& Udiag,
+            vectorField& Usource
         ) const;
 
         //- Add resistance

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLawTemplates.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,10 +31,10 @@ License
 template<class RhoFieldType>
 void Foam::porosityModels::powerLaw::apply
 (
-    scalarField& Udiag,
     const scalarField& V,
     const RhoFieldType& rho,
-    const vectorField& U
+    const vectorField& U,
+    scalarField& Udiag
 ) const
 {
     const scalar C0 = C0_;
@@ -56,9 +56,9 @@ void Foam::porosityModels::powerLaw::apply
 template<class RhoFieldType>
 void Foam::porosityModels::powerLaw::apply
 (
-    tensorField& AU,
     const RhoFieldType& rho,
-    const vectorField& U
+    const vectorField& U,
+    tensorField& AU
 ) const
 {
     const scalar C0 = C0_;

src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -62,12 +62,6 @@ Foam::porosityModels::solidification::solidification
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::porosityModels::solidification::~solidification()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::porosityModels::solidification::calcTransformModelData()
@@ -85,49 +79,50 @@ void Foam::porosityModels::solidification::calcForce
     scalarField Udiag(U.size(), Zero);
     const scalarField& V = mesh_.V();
 
-    apply(Udiag, V, rho, U);
+    apply(V, rho, U, Udiag);
 
-    force = Udiag*U;
+    force += Udiag*U;
 }
 
 
 void Foam::porosityModels::solidification::correct
 (
-    fvVectorMatrix& UEqn
+    const volVectorField& U,
+    scalarField& Udiag,
+    vectorField& Usource,
+    bool compressible
 ) const
 {
-    const volVectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
 
-    if (UEqn.dimensions() == dimForce)
+    if (compressible)
     {
         const auto& rho = mesh_.lookupObject<volScalarField>
         (
             IOobject::groupName(rhoName_, U.group())
         );
 
-        apply(Udiag, V, rho, U);
+        apply(V, rho, U, Udiag);
     }
     else
     {
-        apply(Udiag, V, geometricOneField(), U);
+        apply(V, geometricOneField(), U, Udiag);
     }
 }
 
 
 void Foam::porosityModels::solidification::correct
 (
-    fvVectorMatrix& UEqn,
+    const volVectorField& U,
     const volScalarField& rho,
-    const volScalarField& mu
+    const volScalarField& mu,
+    scalarField& Udiag,
+    vectorField& Usource
 ) const
 {
-    const volVectorField& U = UEqn.psi();
     const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
 
-    apply(Udiag, V, rho, U);
+    apply(V, rho, U, Udiag);
 }
 
 
@@ -146,11 +141,11 @@ void Foam::porosityModels::solidification::correct
             IOobject::groupName(rhoName_, U.group())
         );
 
-        apply(AU, rho, U);
+        apply(rho, U, AU);
     }
     else
     {
-        apply(AU, geometricOneField(), U);
+        apply(geometricOneField(), U, AU);
     }
 }
 

src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -120,40 +121,40 @@ class solidification
         template<class AlphaFieldType, class RhoFieldType>
         void apply
         (
-            scalarField& Udiag,
             const scalarField& V,
             const AlphaFieldType& alpha,
             const RhoFieldType& rho,
-            const volVectorField& U
+            const volVectorField& U,
+            scalarField& Udiag
         ) const;
 
         //- Apply resistance
         template<class AlphaFieldType, class RhoFieldType>
         void apply
         (
-            tensorField& AU,
             const AlphaFieldType& alpha,
             const RhoFieldType& rho,
-            const volVectorField& U
+            const volVectorField& U,
+            tensorField& AU
         ) const;
 
         //- Apply resistance
         template<class RhoFieldType>
         void apply
         (
-            scalarField& Udiag,
             const scalarField& V,
             const RhoFieldType& rho,
-            const volVectorField& U
+            const volVectorField& U,
+            scalarField& Udiag
         ) const;
 
         //- Apply resistance
         template<class RhoFieldType>
         void apply
         (
-            tensorField& AU,
             const RhoFieldType& rho,
-            const volVectorField& U
+            const volVectorField& U,
+            tensorField& AU
         ) const;
 
         //- No copy construct
@@ -179,7 +180,7 @@ public:
     );
 
     //- Destructor
-    virtual ~solidification();
+    virtual ~solidification() = default;
 
 
     // Member Functions
@@ -197,14 +198,22 @@ public:
         ) const;
 
         //- Add resistance
-        virtual void correct(fvVectorMatrix& UEqn) const;
+        virtual void correct
+        (
+            const volVectorField& U,
+            scalarField& Udiag,
+            vectorField& Usource,
+            bool compressible
+        ) const;
 
         //- Add resistance
         virtual void correct
         (
-            fvVectorMatrix& UEqn,
+            const volVectorField& U,
             const volScalarField& rho,
-            const volScalarField& mu
+            const volScalarField& mu,
+            scalarField& Udiag,
+            vectorField& Usource
         ) const;
 
         //- Add resistance

src/finiteVolume/cfdTools/general/porosityModel/solidification/solidificationTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,11 +34,11 @@ License
 template<class AlphaFieldType, class RhoFieldType>
 void Foam::porosityModels::solidification::apply
 (
-    scalarField& Udiag,
     const scalarField& V,
     const AlphaFieldType& alpha,
     const RhoFieldType& rho,
-    const volVectorField& U
+    const volVectorField& U,
+    scalarField& Udiag
 ) const
 {
     const auto& T = mesh_.lookupObject<volScalarField>
@@ -61,10 +62,10 @@ void Foam::porosityModels::solidification::apply
 template<class AlphaFieldType, class RhoFieldType>
 void Foam::porosityModels::solidification::apply
 (
-    tensorField& AU,
     const AlphaFieldType& alpha,
     const RhoFieldType& rho,
-    const volVectorField& U
+    const volVectorField& U,
+    tensorField& AU
 ) const
 {
     const auto& T = mesh_.lookupObject<volScalarField>
@@ -88,15 +89,15 @@ void Foam::porosityModels::solidification::apply
 template<class RhoFieldType>
 void Foam::porosityModels::solidification::apply
 (
-    scalarField& Udiag,
     const scalarField& V,
     const RhoFieldType& rho,
-    const volVectorField& U
+    const volVectorField& U,
+    scalarField& Udiag
 ) const
 {
     if (alphaName_ == "none")
     {
-        return apply(Udiag, V, geometricOneField(), rho, U);
+        return apply(V, geometricOneField(), rho, U, Udiag);
     }
     else
     {
@@ -105,7 +106,7 @@ void Foam::porosityModels::solidification::apply
             IOobject::groupName(alphaName_, U.group())
         );
 
-        return apply(Udiag, V, alpha, rho, U);
+        return apply(V, alpha, rho, U, Udiag);
     }
 }
 
@@ -113,14 +114,14 @@ void Foam::porosityModels::solidification::apply
 template<class RhoFieldType>
 void Foam::porosityModels::solidification::apply
 (
-    tensorField& AU,
     const RhoFieldType& rho,
-    const volVectorField& U
+    const volVectorField& U,
+    tensorField& AU
 ) const
 {
     if (alphaName_ == "none")
     {
-        return apply(AU, geometricOneField(), rho, U);
+        return apply(geometricOneField(), rho, U, AU);
     }
     else
     {
@@ -129,7 +130,7 @@ void Foam::porosityModels::solidification::apply
             IOobject::groupName(alphaName_, U.group())
         );
 
-        return apply(AU, alpha, rho, U);
+        return apply(alpha, rho, U, AU);
     }
 }
 

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -270,7 +270,7 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
         muNbr.primitiveFieldRef()
     );
 
-    porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);
+    porosityPtr_->addResistance(rhoNbr, muNbr, nbrEqn);
 
     // Convert source from neighbour to local region
     fvMatrix<vector> porosityEqn(U, eqn.dimensions());