ENH: Added new hydrostatic pressure function object

Calculates and outputs the pressure fields p_rgh and ph_rgh based on the
option that was previously hard-coded in the fireFoam solver

Usage
    Example of function object specification to calculate hydrostatic pressure:
    \verbatim
    hydrostaticPressure1
    {
        type        hydrostaticPressure;
        libs        ("libinitialisationFunctionObjects.so");
        ...
    }
    \endverbatim

    Where the entries comprise:
    \table
        Property     | Description             | Required    | Default value
        log          | Log to standard output  | no          | yes
        p_rgh        | Name of p_rgh field     | no          | p_rgh
        ph_rgh       | Name of ph_rgh field    | no          | ph_rgh
        pRef         | Name of pressure ref field | no       | pRef
        rho          | Name of density field   | no          | rho
        U            | Name of velocity field  | no          | U
        gh           | Name of gravity*height volume field  | no | gh
        ghf          | Name of gravity*height surface field | no | ghf
        nCorrectors  | Number of correctors when solving ph_rgh | no | 5
    \endtable

Note
    Calculates the hydrostatic pressure on construction/re-initialisation;
    the execute and write functions are not used.

src/functionObjects/Allwmake

@@ -6,6 +6,7 @@ cd ${0%/*} || exit 1                        # Run from this directory
 
 wmake $targetType field
 wmake $targetType forces
+wmake $targetType initialisation
 wmake $targetType lagrangian
 wmake $targetType utilities
 wmake $targetType solvers

src/functionObjects/doc/functionObjects.dox

@@ -41,6 +41,7 @@ The current range of features comprises of:
 - \ref grpFieldFunctionObjects
 - \ref grpForcesFunctionObjects
 - \ref grpGraphicsFunctionObjects
+- \ref grpInitialisationFunctionObjects
 - \ref grpLagrangianFunctionObjects
 - \ref grpSolversFunctionObjects
 - \ref grpUtilitiesFunctionObjects

src/functionObjects/initialisation/Make/files

@@ -0,0 +1,3 @@
+hydrostaticPressure/hydrostaticPressure.C
+
+LIB = $(FOAM_LIBBIN)/libinitialisationFunctionObjects

src/functionObjects/initialisation/Make/options

@@ -0,0 +1,43 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/fvAgglomerationMethods/pairPatchAgglomeration/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -llagrangian \
+    -ldistributionModels \
+    -lsampling \
+    -lsurfMesh \
+    -lfluidThermophysicalModels \
+    -lsolidThermo \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lcompressibleTransportModels \
+    -lincompressibleTurbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lchemistryModel \
+    -lreactionThermophysicalModels \
+    -lpairPatchAgglomeration

src/functionObjects/initialisation/hydrostaticPressure/hydrostaticPressure.C

@@ -0,0 +1,243 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "hydrostaticPressure.H"
+#include "basicThermo.H"
+#include "uniformDimensionedFields.H"
+#include "volFields.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvmLaplacian.H"
+#include "fvcSnGrad.H"
+#include "constrainPressure.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(hydrostaticPressure, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        hydrostaticPressure,
+        dictionary
+    );
+}
+}
+
+
+Foam::dimensionedScalar
+Foam::functionObjects::hydrostaticPressure::pRef() const
+{
+    if (pRefName_ == "none")
+    {
+        return dimensionedScalar("pRef", dimPressure, pRefValue_);
+    }
+    else
+    {
+        return mesh_.lookupObject<uniformDimensionedScalarField>(pRefName_);
+    }
+}
+
+
+void Foam::functionObjects::hydrostaticPressure::calculateAndWrite()
+{
+    const auto& pRef = this->pRef();
+    const auto& U = mesh_.lookupObject<volVectorField>(UName_);
+    const auto& gh = mesh_.lookupObject<volScalarField>(ghName_);
+    const auto& ghf = mesh_.lookupObject<surfaceScalarField>(ghfName_);
+    auto& rho = mesh_.lookupObjectRef<volScalarField>(rhoName_);
+    auto& thermo = mesh_.lookupObjectRef<basicThermo>(basicThermo::dictName);
+    auto& p_rgh = mesh_.lookupObjectRef<volScalarField>(p_rghName_);
+    auto& ph_rgh = mesh_.lookupObjectRef<volScalarField>(ph_rghName_);
+
+    auto& p = thermo.p();
+
+    Info<< "Performing hydrostatic pressure initialisation";
+    if (mesh_.name() != polyMesh::defaultRegion)
+    {
+        Info<< "for region " << mesh_.name();
+    }
+
+
+    if (thermo.incompressible())
+    {
+        Info<< ": incompressible" << endl;
+
+        // Constant density and temperature
+
+        thermo.correct();
+        rho = thermo.rho();
+        p = ph_rgh + rho*gh + pRef;
+        p_rgh = ph_rgh;
+    }
+    else
+    {
+        Info<< ": compressible" << endl;
+
+        p = ph_rgh + rho*gh + pRef;
+        thermo.correct();
+        rho = thermo.rho();
+
+        for (label i=0; i<nCorrectors_; ++i)
+        {
+            surfaceScalarField rhof("rhof", fvc::interpolate(rho));
+
+            surfaceScalarField phig
+            (
+                "phig",
+               -rhof*ghf*fvc::snGrad(rho)*mesh_.magSf()
+            );
+
+            // Update the pressure BCs to ensure flux consistency
+            constrainPressure(ph_rgh, rho, U, phig, rhof);
+
+            fvScalarMatrix ph_rghEqn
+            (
+                fvm::laplacian(rhof, ph_rgh) == fvc::div(phig)
+            );
+
+            ph_rghEqn.relax();
+
+            ph_rghEqn.solve();
+
+            p = ph_rgh + rho*gh + pRef;
+            thermo.correct();
+            rho = thermo.rho();
+
+            Info<< "Hydrostatic pressure variation "
+                << (max(ph_rgh) - min(ph_rgh)).value() << endl;
+        }
+
+        p_rgh = ph_rgh;
+
+        Log << "    writing field " << ph_rgh.name() << nl;
+        ph_rgh.write();
+    }
+
+    Log << "    writing field " << rho.name() << nl;
+    rho.write();
+
+    Log << "    writing field " << p_rgh.name() << nl;
+    p_rgh.write();
+
+    Log << "    writing field " << p.name() << nl;
+    p.write();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::hydrostaticPressure::hydrostaticPressure
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    p_rghName_("p_rgh"),
+    ph_rghName_("ph_rgh"),
+    pRefName_("pRef"),
+    pRefValue_(0),
+    rhoName_("rho"),
+    UName_("U"),
+    ghName_("gh"),
+    ghfName_("ghf"),
+    nCorrectors_(5)
+{
+    if (read(dict))
+    {
+        // Read and store the initial ph_rgh field
+        volScalarField* ph_rghPtr =
+            new volScalarField
+            (
+                IOobject
+                (
+                    ph_rghName_,
+                    runTime.timeName(),
+                    mesh_,
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE // To enable restart
+                ),
+                mesh_
+            );
+
+        mesh_.objectRegistry::store(ph_rghPtr);
+
+        bool reInitialise = dict.lookupOrDefault<bool>("reInitialise", false);
+
+        if (runTime.timeIndex() == 0 || reInitialise)
+        {
+            calculateAndWrite();
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::hydrostaticPressure::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict))
+    {
+        dict.readIfPresent("p_rgh", p_rghName_);
+        dict.readIfPresent("ph_rgh", ph_rghName_);
+        dict.readIfPresent("pRef", pRefName_);
+        dict.readIfPresent("rho", rhoName_);
+        dict.readIfPresent("U", UName_);
+        dict.readIfPresent("gh", ghName_);
+        dict.readIfPresent("ghf", ghfName_);
+        dict.readIfPresent("nCorrectors", nCorrectors_);
+
+        if (pRefName_ == "none")
+        {
+            pRefValue_ = dict.get<scalar>("pRefValue");
+        }
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::hydrostaticPressure::execute()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::hydrostaticPressure::write()
+{
+    return true;
+}
+
+
+// ************************************************************************* //

src/functionObjects/initialisation/hydrostaticPressure/hydrostaticPressure.H

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::hydrostaticPressure
+
+Group
+    grpInitialisationFunctionObjects
+
+Description
+    Calculates and outputs the pressure fields p_rgh and ph_rgh
+
+Usage
+    Example of function object specification to calculate hydrostatic pressure:
+    \verbatim
+    hydrostaticPressure1
+    {
+        type        hydrostaticPressure;
+        libs        ("libinitialisationFunctionObjects.so");
+        ...
+    }
+    \endverbatim
+
+    Where the entries comprise:
+    \table
+        Property     | Description             | Required    | Default value
+        log          | Log to standard output  | no          | yes
+        p_rgh        | Name of p_rgh field     | no          | p_rgh
+        ph_rgh       | Name of ph_rgh field    | no          | ph_rgh
+        pRef         | Name of pressure ref field | no       | pRef
+        rho          | Name of density field   | no          | rho
+        U            | Name of velocity field  | no          | U
+        gh           | Name of gravity*height volume field  | no | gh
+        ghf          | Name of gravity*height surface field | no | ghf
+        nCorrectors  | Number of correctors when solving ph_rgh | no | 5
+    \endtable
+
+Note
+    Calculates the hydrostatic pressure on construction/re-initialisation;
+    the execute and write functions are not used.
+
+SourceFiles
+    hydrostaticPressure.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_hydrostaticPressure_H
+#define functionObjects_hydrostaticPressure_H
+
+#include "fvMeshFunctionObject.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class hydrostaticPressure Declaration
+\*---------------------------------------------------------------------------*/
+
+class hydrostaticPressure
+:
+    public fvMeshFunctionObject
+{
+protected:
+
+    // Protected data
+
+        //- Name of p - rho*g*h field, default is "p_rgh"
+        word p_rghName_;
+
+        //- Name of p_hydrostatic - rho*g*h field, default is "ph_rgh"
+        word ph_rghName_;
+
+        //- Name of uniform pressure reference field, default is "pRef"
+        //  Set to "none" to set to zero
+        word pRefName_;
+
+        //- Reference pressure if pRef is set to "none"
+        scalar pRefValue_;
+
+        //- Name of density field, default is "rho"
+        word rhoName_;
+
+        //- Name of velocity field, default is "ph_rgh"
+        word UName_;
+
+        //- Name of g*h volume field, default is "gh"
+        word ghName_;
+
+        //- Name of g*h surface field, default is "ghf"
+        word ghfName_;
+
+        //- Number of correctors when solving for ph_rgh
+        label nCorrectors_;
+
+
+    // Protected Member Functions
+
+        //- Helper function to return the reference pressure
+        //  Returns zero if pRefName = none; otherwise returns the value
+        //  from the mesh database
+        dimensionedScalar pRef() const;
+
+        //- Calculate the fields and write
+        void calculateAndWrite();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("hydrostaticPressure");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        hydrostaticPressure
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~hydrostaticPressure() = default;
+
+
+    // Member Functions
+
+        //- Read the hydrostaticPressure data
+        virtual bool read(const dictionary& dict);
+
+        //- Calculate the p_rgh field
+        virtual bool execute();
+
+        //- Write the p_rgh and derived fields
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //