Updates and extensions to the MRF and SRF handling.

applications/solvers/incompressible/pisoFoam/pisoFoam.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    turbFoam
+    pisoFoam
 
 Description
     Transient solver for incompressible flow.
@@ -40,15 +40,14 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "createTime.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -56,8 +55,8 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
-#       include "CourantNo.H"
+        #include "CourantNo.H"
 
         // Pressure-velocity PISO corrector
         {
@@ -120,7 +119,7 @@ int main(int argc, char *argv[])
                     }
                 }
 
-#               include "continuityErrs.H"
+                #include "continuityErrs.H"
 
                 U -= rUA*fvc::grad(p);
                 U.correctBoundaryConditions();

bin/foamPackSource

@@ -77,9 +77,9 @@ find -H $packDir               \
  -a ! -name "core.[1-9]*"      \
  -a ! -name "log[0-9]*"        \
  -a ! -name "libccmio*"        \
- -a ! -name "\.ebrowse"        \
 | sed                          \
  -e "\@$packDir/.git/@d"       \
+ -e "\@/.tags/@d"              \
  -e "\@$packDir/lib/@d"        \
  -e "\@libccmio.*/@d"          \
  -e '\@applications/bin/@d'    \

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

@@ -32,53 +32,96 @@ License
 #include "syncTools.H"
 #include "faceSet.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::MRFZone, 0);
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::MRFZone::setMRFFaces
+void Foam::MRFZone::setMRFFaces()
-(
-    labelList& faceType,
-    const labelList& excludedPatchIDs
-)
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    // Knock out coupled patches
+    // Type per face:
-    forAll(patches, patchi)
+    //  0:not in zone
+    //  1:moving with frame
+    //  2:other
+    labelList faceType(mesh_.nFaces(), 0);
+
+    // Determine faces in cell zone
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // (without constructing cells)
+
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+    // Cells in zone
+    boolList zoneCell(mesh_.nCells(), false);
+
+    if (cellZoneID_ != -1)
     {
-        if (patches[patchi].coupled())
+        const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
+        forAll(cellLabels, i)
         {
-            const polyPatch& pp = patches[patchi];
+            zoneCell[cellLabels[i]] = true;
+        }
+    }
 
-            forAll(pp, j)
+
+    label nZoneFaces = 0;
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
+        {
+            faceType[faceI] = 1;
+            nZoneFaces++;
+        }
+    }
+
+
+    labelHashSet excludedPatches(excludedPatchLabels_);
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled() || excludedPatches.found(patchI))
+        {
+            forAll(pp, i)
             {
-                label faceI = pp.start()+j;
+                label faceI = pp.start()+i;
 
-                if (faceType[faceI] == 1)
+                if (zoneCell[own[faceI]])
                 {
                     faceType[faceI] = 2;
+                    nZoneFaces++;
+                }
+            }
+        }
+        else if (!isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+
+                if (zoneCell[own[faceI]])
+                {
+                    faceType[faceI] = 1;
+                    nZoneFaces++;
                 }
             }
         }
     }
 
-    // All explicitly provided exclusions
+    // Now we have for faceType:
-    forAll(excludedPatchIDs, i)
+    //  0   : face not in cellZone
-    {
+    //  1   : internal face or normal patch face
-        const polyPatch& pp = patches[excludedPatchIDs[i]];
+    //  2   : coupled patch face or excluded patch face
 
-        forAll(pp, j)
+    // Sort into lists per patch.
-        {
-            label faceI = pp.start()+j;
 
-            if (faceType[faceI] == 1)
-            {
-                faceType[faceI] = 2;
-            }
-        }
-    }
-
-    // Collect into lists per patch.
     internalFaces_.setSize(mesh_.nFaces());
     label nInternal = 0;
 
@@ -142,148 +185,42 @@ void Foam::MRFZone::setMRFFaces
         }
     }
 
-    //if (debug)
-    //{
-    //    faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
-    //    Pout<< "Writing internal faces in MRF zone to faceSet "
-    //        << internalFaces.name() << endl;
-    //    internalFaces.write();
-    //}
-    //{
-    //    faceSet MRFFaces(mesh_, "includedFaces", 100);
-    //    forAll(includedFaces_, patchi)
-    //    {
-    //        forAll(includedFaces_[patchi], i)
-    //        {
-    //            label patchFacei = includedFaces_[patchi][i];
-    //            MRFFaces.insert(patches[patchi].start()+patchFacei);
-    //        }
-    //    }
-    //    Pout<< "Writing patch faces in MRF zone to faceSet "
-    //        << MRFFaces.name() << endl;
-    //    MRFFaces.write();
-    //}
-    //{
-    //    faceSet excludedFaces(mesh_, "excludedFaces", 100);
-    //    forAll(excludedFaces_, patchi)
-    //    {
-    //        forAll(excludedFaces_[patchi], i)
-    //        {
-    //            label patchFacei = excludedFaces_[patchi][i];
-    //            excludedFaces.insert(patches[patchi].start()+patchFacei);
-    //        }
-    //    }
-    //    Pout<< "Writing faces in MRF zone with special handling to faceSet "
-    //        << excludedFaces.name() << endl;
-    //    excludedFaces.write();
-    //}
-}
 
-
+    if (debug)
-void Foam::MRFZone::setMRFFaces()
-{
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    // Type per face:
-    //  0:not in zone
-    //  1:moving with frame
-    //  2:other
-    labelList faceType(mesh_.nFaces(), 0);
-
-    bool faceZoneFound = (faceZoneID_ != -1);
-    reduce(faceZoneFound, orOp<bool>());
-
-    if (faceZoneFound)
     {
-        // Explicitly provided faces.
+        faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
-        if (faceZoneID_ != -1)
+        Pout<< "Writing " << internalFaces.size()
+            << " internal faces in MRF zone to faceSet "
+            << internalFaces.name() << endl;
+        internalFaces.write();
+
+        faceSet MRFFaces(mesh_, "includedFaces", 100);
+        forAll(includedFaces_, patchi)
         {
-            const labelList& zoneFaces = mesh_.faceZones()[faceZoneID_];
+            forAll(includedFaces_[patchi], i)
-
-            forAll(zoneFaces, i)
             {
-                faceType[zoneFaces[i]] = 1;
+                label patchFacei = includedFaces_[patchi][i];
-            }
+                MRFFaces.insert(patches[patchi].start()+patchFacei);
-
-            if (allPatchesMove_)
-            {
-                // Explicitly provided patches that do not move
-                setMRFFaces(faceType, patchLabels_);
-            }
-            else
-            {
-                setMRFFaces(faceType, labelList(0));
             }
         }
-    }
+        Pout<< "Writing " << MRFFaces.size()
-    else
+            << " patch faces in MRF zone to faceSet "
-    {
+            << MRFFaces.name() << endl;
-        // Determine faces in cell zone
+        MRFFaces.write();
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        // (without constructing cells)
 
-        const labelList& own = mesh_.faceOwner();
+        faceSet excludedFaces(mesh_, "excludedFaces", 100);
-        const labelList& nei = mesh_.faceNeighbour();
+        forAll(excludedFaces_, patchi)
-
-        // Cells in zone
-        boolList zoneCell(mesh_.nCells(), false);
-
-        if (cellZoneID_ != -1)
         {
-            const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
+            forAll(excludedFaces_[patchi], i)
-            forAll(cellLabels, i)
             {
-                zoneCell[cellLabels[i]] = true;
+                label patchFacei = excludedFaces_[patchi][i];
+                excludedFaces.insert(patches[patchi].start()+patchFacei);
             }
         }
-
+        Pout<< "Writing " << excludedFaces.size()
-
+            << " faces in MRF zone with special handling to faceSet "
-        label nZoneFaces = 0;
+            << excludedFaces.name() << endl;
-
+        excludedFaces.write();
-        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-        {
-            if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
-            {
-                faceType[faceI] = 1;
-                nZoneFaces++;
-            }
-        }
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if (!isA<emptyPolyPatch>(pp))
-            {
-                forAll(pp, i)
-                {
-                    label faceI = pp.start()+i;
-
-                    if (zoneCell[own[faceI]])
-                    {
-                        faceType[faceI] = 1;
-                        nZoneFaces++;
-                    }
-                }
-            }
-        }
-        syncTools::syncFaceList(mesh_, faceType, maxEqOp<label>(), false);
-
-
-        Info<< nl
-            << "MRFZone " << name_ << " : did not find a faceZone; using "
-            << returnReduce(nZoneFaces, sumOp<label>())
-            << " faces from the cellZone instead." << endl;
-
-
-        if (allPatchesMove_)
-        {
-            // Explicitly provided excluded patches
-            setMRFFaces(faceType, patchLabels_);
-        }
-        else
-        {
-            setMRFFaces(faceType, labelList(0));
-        }
     }
 }
 
@@ -296,36 +233,41 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
     name_(is),
     dict_(is),
     cellZoneID_(mesh_.cellZones().findZoneID(name_)),
-    faceZoneID_(mesh_.faceZones().findZoneID(name_)),
+    excludedPatchNames_
-    allPatchesMove_(dict_.found("nonRotatingPatches")),
-    patchNames_
     (
-        allPatchesMove_
+        dict_.lookupOrDefault("nonRotatingPatches", wordList(0))
-      ? dict_.lookup("nonRotatingPatches")
-      : dict_.lookup("patches")
     ),
     origin_(dict_.lookup("origin")),
     axis_(dict_.lookup("axis")),
     omega_(dict_.lookup("omega")),
     Omega_("Omega", omega_*axis_)
 {
+    if (dict_.found("patches"))
+    {
+        WarningIn("MRFZone(const fvMesh&, Istream&)")
+            << "Ignoring entry 'patches'\n"
+            << "    By default all patches within the rotating region rotate.\n"
+            << "    Optionally supply excluded patches using 'nonRotatingPatches'."
+            << endl;
+    }
+
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
     axis_ = axis_/mag(axis_);
     Omega_ = omega_*axis_;
 
-    patchLabels_.setSize(patchNames_.size());
+    excludedPatchLabels_.setSize(excludedPatchNames_.size());
 
-    forAll(patchNames_, i)
+    forAll(excludedPatchNames_, i)
     {
-        patchLabels_[i] = patches.findPatchID(patchNames_[i]);
+        excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
 
-        if (patchLabels_[i] == -1)
+        if (excludedPatchLabels_[i] == -1)
         {
             FatalErrorIn
             (
                 "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
-            )   << "cannot find MRF patch " << patchNames_[i]
+            )   << "cannot find MRF patch " << excludedPatchNames_[i]
                 << exit(FatalError);
         }
     }
@@ -364,7 +306,46 @@ void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
 
     forAll(cells, i)
     {
-        Usource[cells[i]] -= V[cells[i]]*(Omega ^ U[cells[i]]);
+        label celli = cells[i];
+        Usource[celli] -= V[celli]*(Omega ^ U[celli]);
+    }
+}
+
+
+void Foam::MRFZone::relativeVelocity(volVectorField& U) const
+{
+    const volVectorField& C = mesh_.C();
+
+    const vector& origin = origin_.value();
+    const vector& Omega = Omega_.value();
+
+    const labelList& cells = mesh_.cellZones()[cellZoneID_];
+
+    forAll(cells, i)
+    {
+        label celli = cells[i];
+        U[celli] -= (Omega ^ (C[celli] - origin));
+    }
+
+    // Included patches
+    forAll(includedFaces_, patchi)
+    {
+        forAll(includedFaces_[patchi], i)
+        {
+            label patchFacei = includedFaces_[patchi][i];
+            U.boundaryField()[patchi][patchFacei] = vector::zero;
+        }
+    }
+
+    // Excluded patches
+    forAll(excludedFaces_, patchi)
+    {
+        forAll(excludedFaces_[patchi], i)
+        {
+            label patchFacei = excludedFaces_[patchi][i];
+            U.boundaryField()[patchi][patchFacei] -=
+                (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
+        }
     }
 }
 
@@ -410,6 +391,49 @@ void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
 }
 
 
+void Foam::MRFZone::absoluteFlux(surfaceScalarField& phi) const
+{
+    const surfaceVectorField& Cf = mesh_.Cf();
+    const surfaceVectorField& Sf = mesh_.Sf();
+
+    const vector& origin = origin_.value();
+    const vector& Omega = Omega_.value();
+
+    // Internal faces
+    forAll(internalFaces_, i)
+    {
+        label facei = internalFaces_[i];
+        phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei];
+    }
+
+    // Included patches
+    forAll(includedFaces_, patchi)
+    {
+        forAll(includedFaces_[patchi], i)
+        {
+            label patchFacei = includedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] +=
+                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+              & Sf.boundaryField()[patchi][patchFacei];
+        }
+    }
+
+    // Excluded patches
+    forAll(excludedFaces_, patchi)
+    {
+        forAll(excludedFaces_[patchi], i)
+        {
+            label patchFacei = excludedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] +=
+                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+              & Sf.boundaryField()[patchi][patchFacei];
+        }
+    }
+}
+
+
 void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
 {
     const vector& origin = origin_.value();

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

@@ -26,7 +26,7 @@ Class
     Foam::MRFZone
 
 Description
-    MRF zone definition based on cell zone and optional face zone and parameters
+    MRF zone definition based on cell zone and parameters
     obtained from a control dictionary constructed from the given stream.
 
     The rotation of the MRF region is defined by an origin and axis of
@@ -74,15 +74,8 @@ class MRFZone
 
         label cellZoneID_;
 
-        //- label of face zone with faces on outside of cell zone.
+        const wordList excludedPatchNames_;
-        label faceZoneID_;
+        labelList excludedPatchLabels_;
-
-        //- Do patches move with frame (true) or are explicitly provided (false,
-        //  old behaviour)
-        Switch allPatchesMove_;
-
-        const wordList patchNames_;
-        labelList patchLabels_;
 
         //- Internal faces that are part of MRF
         labelList internalFaces_;
@@ -101,17 +94,9 @@ class MRFZone
 
     // Private Member Functions
 
-        //- Divide faces in frame according to patch
-        void setMRFFaces
-        (
-            labelList& faceType,
-            const labelList& excludedPatchIDs
-        );
-
         //- Divide faces in frame according to patch
         void setMRFFaces();
 
-
         //- Disallow default bitwise copy construct
         MRFZone(const MRFZone&);
 
@@ -121,6 +106,10 @@ class MRFZone
 
 public:
 
+    // Declare name of the class and its debug switch
+    ClassName("MRFZone");
+
+
     // Constructors
 
         //- Construct from fvMesh and Istream
@@ -165,9 +154,15 @@ public:
         //- Add the Coriolis force contribution to the momentum equation
         void addCoriolis(fvVectorMatrix& UEqn) const;
 
+        //- Make the given absolute velocity relative within the MRF region
+        void relativeVelocity(volVectorField& U) const;
+
         //- Make the given absolute flux relative within the MRF region
         void relativeFlux(surfaceScalarField& phi) const;
 
+        //- Make the given relative flux absolute within the MRF region
+        void absoluteFlux(surfaceScalarField& phi) const;
+
         //- Correct the boundary velocity for the roation of the MRF region
         void correctBoundaryVelocity(volVectorField& U) const;
 

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

@@ -65,6 +65,15 @@ void Foam::MRFZones::addCoriolis(fvVectorMatrix& UEqn) const
 }
 
 
+void Foam::MRFZones::relativeVelocity(volVectorField& U) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).relativeVelocity(U);
+    }
+}
+
+
 void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
 {
     forAll(*this, i)
@@ -74,6 +83,15 @@ void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
 }
 
 
+void Foam::MRFZones::absoluteFlux(surfaceScalarField& phi) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).absoluteFlux(phi);
+    }
+}
+
+
 void Foam::MRFZones::correctBoundaryVelocity(volVectorField& U) const
 {
     forAll(*this, i)

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

@@ -76,9 +76,15 @@ public:
         //- Add the Coriolis force contribution to the momentum equation
         void addCoriolis(fvVectorMatrix& UEqn) const;
 
+        //- Make the given absolute velocity relative within the MRF region
+        void relativeVelocity(volVectorField& U) const;
+
         //- Make the given absolute flux relative within the MRF region
         void relativeFlux(surfaceScalarField& phi) const;
 
+        //- Make the given relative flux absolute within the MRF region
+        void absoluteFlux(surfaceScalarField& phi) const;
+
         //- Correct the boundary velocity for the roation of the MRF region
         void correctBoundaryVelocity(volVectorField& U) const;
 };

src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C

@@ -22,10 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
-    Formulation based on relative velocities
-
 \*---------------------------------------------------------------------------*/
 
 #include "SRFModel.H"

tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C

@@ -63,10 +63,10 @@ int main(int argc, char *argv[])
         // Pressure-velocity SIMPLE corrector
         {
             // Momentum predictor
-
             tmp<fvVectorMatrix> UEqn
             (
                 fvm::div(phi, U)
+              - fvm::Sp(fvc::div(phi), U)
               + turbulence->divDevReff(U)
             );
             mrfZones.addCoriolis(UEqn());

tutorials/incompressible/simpleSRFFoam/simpleSRFFoam/UrelEqn.H

@@ -1,9 +1,9 @@
-    // Momentum predictor
+    // Relative momentum predictor
     tmp<fvVectorMatrix> UrelEqn
     (
         fvm::div(phi, Urel)
-        + turbulence->divDevReff(Urel)
+      + turbulence->divDevReff(Urel)
-        + SRF->Su()
+      + SRF->Su()
     );
 
     UrelEqn().relax();

tutorials/incompressible/simpleSRFFoam/simpleSRFFoam/simpleSRFFoam.C

@@ -27,7 +27,7 @@ Application
 
 Description
     Steady-state solver for incompressible, turbulent flow of non-Newtonian
-    fluids with single rotating frame.
+    fluids in a single rotating frame.
 
 \*---------------------------------------------------------------------------*/
 
@@ -40,16 +40,13 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "setRootCase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-    //mesh.clearPrimitives();
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -57,20 +54,19 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
+        #include "readSIMPLEControls.H"
-#       include "initConvergenceCheck.H"
+        #include "initConvergenceCheck.H"
 
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
         {
-#           include "UEqn.H"
+            #include "UrelEqn.H"
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
 
-
         if (runTime.outputTime())
         {
             volVectorField Uabs
@@ -93,7 +89,7 @@ int main(int argc, char *argv[])
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
 
-#       include "convergenceCheck.H"
+        #include "convergenceCheck.H"
     }
 
     Info<< "End\n" << endl;