ENH: limitHeight: new faOption to limit film height

src/faOptions/Make/files

@@ -17,5 +17,6 @@ $(derivedSources)/externalFileSource/externalFileSource.C
 corrections=corrections
 
 $(corrections)/limitVelocity/limitVelocity.C
+$(corrections)/limitHeight/limitHeight.C
 
 LIB = $(FOAM_LIBBIN)/libfaOptions

src/faOptions/corrections/limitHeight/limitHeight.C

@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "limitHeight.H"
+#include "areaFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fa
+{
+    defineTypeNameAndDebug(limitHeight, 0);
+    addToRunTimeSelectionTable
+    (
+        option,
+        limitHeight,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fa::limitHeight::limitHeight
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    faceSetOption(name, modelType, dict, mesh),
+    hName_("h"),
+    max_(0)  // overwritten later
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fa::limitHeight::read(const dictionary& dict)
+{
+    if (!faceSetOption::read(dict))
+    {
+        return false;
+    }
+
+    coeffs_.readIfPresent("h", hName_);
+    coeffs_.readEntry("max", max_);
+
+    fieldNames_.resize(1, hName_);
+
+    applied_.resize(1, false);
+
+    return true;
+}
+
+
+void Foam::fa::limitHeight::correct(areaScalarField& h)
+{
+    // Count nTotFaces ourselves
+    // (maybe only applying on a subset)
+    label nFacesAbove = 0;
+    const label nTotFaces = returnReduce(faces_.size(), sumOp<label>());
+
+    scalarField& hif = h.primitiveFieldRef();
+
+    for (const label facei : faces_)
+    {
+        auto& hval = hif[facei];
+
+        if (hval > max_)
+        {
+            hval *= max_/max(hval, SMALL);
+            ++nFacesAbove;
+        }
+    }
+
+    // Handle boundaries in the case of 'all'
+    label nEdgesAbove = 0;
+
+    if (!faceSetOption::useSubMesh())
+    {
+        for (faPatchScalarField& hp : h.boundaryFieldRef())
+        {
+            if (!hp.fixesValue())
+            {
+                for (auto& hval : hp)
+                {
+                    if (hval > max_)
+                    {
+                        hval *= max_/max(hval, SMALL);
+                        ++nEdgesAbove;
+                    }
+                }
+            }
+        }
+    }
+
+    // Percent, max 2 decimal places
+    const auto percent = [](scalar num, label denom) -> scalar
+    {
+        return (denom ? 1e-2*round(1e4*num/denom) : 0);
+    };
+
+
+    reduce(nFacesAbove, sumOp<label>());
+    reduce(nEdgesAbove, sumOp<label>());
+
+    Info<< type() << ' ' << name_ << " Limited "
+        << nFacesAbove << " ("
+        << percent(nFacesAbove, nTotFaces)
+        << "%) of faces, with max limit " << max_ << endl;
+
+    if (nFacesAbove || nEdgesAbove)
+    {
+        // We've changed internal values so give
+        // boundary conditions opportunity to correct
+        h.correctBoundaryConditions();
+    }
+}
+
+
+// ************************************************************************* //

src/faOptions/corrections/limitHeight/limitHeight.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fa::limitHeight
+
+Description
+    Limits the film height to a specified \c max value.
+
+    The \c limitHeight applies an explicit correction to limit the max film
+    height to the user-specified value. This will ensure that the max value
+    is adhered to, but may become unstable if the limiting instigates an
+    oscillatory/destabilising effect into neighbouring faces that then spreads.
+
+Usage
+    Minimal example by using \c constant/faOptions:
+    \verbatim
+    <faOption>
+    {
+        // Mandatory entries
+        type            limitHeight;
+        active          yes;
+        selectionMode   all;
+        max             <scalar>;
+
+        // Optional entries
+        h               <word>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                   | Type   | Reqd  | Deflt
+      type      | Type name: limitHeight        | word   | yes   | -
+      max       | Maximum height limit [m]      | scalar | yes   | -
+      h         | Name of operand height field  | word   | no    | h
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link faceSetOption.H \endlink
+
+SourceFiles
+    limitHeight.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fa_limitHeight_H
+#define fa_limitHeight_H
+
+#include "faceSetOption.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fa
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class limitHeight Declaration
+\*---------------------------------------------------------------------------*/
+
+class limitHeight
+:
+    public faceSetOption
+{
+protected:
+
+    // Protected Data
+
+        //- Name of operand height field
+        word hName_;
+
+        //- Maximum height [m]
+        scalar max_;
+
+
+    // Protected Member Functions
+
+        //- No copy construct
+        limitHeight(const limitHeight&) = delete;
+
+        //- No copy assignment
+        void operator=(const limitHeight&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("limitHeight");
+
+
+    // Constructors
+
+        //- Construct from components
+        limitHeight
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~limitHeight() = default;
+
+
+    // Member Functions
+
+        //- Read dictionary
+        virtual bool read(const dictionary& dict);
+
+        //- Correct the height field
+        virtual void correct(areaScalarField& h);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //