waves: Added atmospheric boundary layer modelling to the gas

An atmospheric boundary layer velocity can now be added to the gas side
of the wave modelling. The wave superposition class has been given a
run-time selection mechanism, and a derivation added which includes gas
atmospheric boundary layer modelling. This modelling is therefore
available in both the wave boundary conditions, and in setWaves.

This functionality can be selected in the constant/waveProperties file
by supplying a "type" entry and a number of parameters controlling the
boundary layer. For example:

    In constant/waveProperties:

        type            waveAtmBoundaryLayer;

        // properties specifying the wave modelling ...

        UGasRef         (10 0 0);
        hRef            20;
        hWaveMin        -2;
        hWaveMax        3;

UGasRef is the gas velocity relative to the liquid, at the height, hRef,
relative to the wave model origin. hWaveMin and hWaveMax describe the
range of the wave elevation; it is non-trivial to calculate this from
the wave models themselves, so it is required as an input.

The base wave superposition class can be selected with "type wave;", but
also selects by default when the "type" entry is omitted, so the change
is backwards compatible.

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C

@@ -25,14 +25,35 @@ License
 
 #include "atmBoundaryLayer.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-namespace Foam
+const Foam::scalar Foam::atmBoundaryLayer::kappaDefault_ = 0.41;
+
+const Foam::scalar Foam::atmBoundaryLayer::CmuDefault_ = 0.09;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::atmBoundaryLayer::init()
 {
+    if (mag(flowDir_) < small || mag(zDir_) < small)
+    {
+        FatalErrorInFunction
+            << "magnitude of n or z must be greater than zero"
+            << abort(FatalError);
+    }
+
+    // Ensure direction vectors are normalized
+    flowDir_ /= mag(flowDir_);
+    zDir_ /= mag(zDir_);
+
+    Ustar_ = kappa_*Uref_/(log((Zref_ + z0_)/z0_));
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-atmBoundaryLayer::atmBoundaryLayer()
+Foam::atmBoundaryLayer::atmBoundaryLayer()
 :
     flowDir_(Zero),
     zDir_(Zero),
@@ -50,12 +71,49 @@ atmBoundaryLayer::atmBoundaryLayer()
 {}
 
 
-atmBoundaryLayer::atmBoundaryLayer(const vectorField& p, const dictionary& dict)
+Foam::atmBoundaryLayer::atmBoundaryLayer
+(
+    const vector& flowDir,
+    const vector& zDir,
+    const scalar Uref,
+    const scalar Zref,
+    const scalarField& z0,
+    const scalarField& zGround,
+    const scalar kappa,
+    const scalar Cmu,
+    const scalar Ulower,
+    const scalar kLower,
+    const scalar epsilonLower
+)
+:
+    flowDir_(flowDir),
+    zDir_(zDir),
+    kappa_(kappa),
+    Cmu_(Cmu),
+    Uref_(Uref),
+    Zref_(Zref),
+    z0_(z0),
+    zGround_(zGround),
+    Ustar_(z0.size()),
+    offset_(Ulower != 0),
+    Ulower_(Ulower),
+    kLower_(kLower),
+    epsilonLower_(epsilonLower)
+{
+    init();
+}
+
+
+Foam::atmBoundaryLayer::atmBoundaryLayer
+(
+    const vectorField& p,
+    const dictionary& dict
+)
 :
     flowDir_(dict.lookup("flowDir")),
     zDir_(dict.lookup("zDir")),
-    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
-    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", kappaDefault_)),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", CmuDefault_)),
     Uref_(readScalar(dict.lookup("Uref"))),
     Zref_(readScalar(dict.lookup("Zref"))),
     z0_("z0", dict, p.size()),
@@ -66,22 +124,11 @@ atmBoundaryLayer::atmBoundaryLayer(const vectorField& p, const dictionary& dict)
     kLower_(dict.lookupOrDefault<scalar>("kLower", 0)),
     epsilonLower_(dict.lookupOrDefault<scalar>("epsilonLower", 0))
 {
-    if (mag(flowDir_) < small || mag(zDir_) < small)
-    {
-        FatalErrorInFunction
-            << "magnitude of n or z must be greater than zero"
-            << abort(FatalError);
-    }
-
-    // Ensure direction vectors are normalized
-    flowDir_ /= mag(flowDir_);
-    zDir_ /= mag(zDir_);
-
-    Ustar_ = kappa_*Uref_/(log((Zref_ + z0_)/z0_));
+    init();
 }
 
 
-atmBoundaryLayer::atmBoundaryLayer
+Foam::atmBoundaryLayer::atmBoundaryLayer
 (
     const atmBoundaryLayer& abl,
     const fvPatchFieldMapper& mapper
@@ -103,7 +150,7 @@ atmBoundaryLayer::atmBoundaryLayer
 {}
 
 
-atmBoundaryLayer::atmBoundaryLayer(const atmBoundaryLayer& abl)
+Foam::atmBoundaryLayer::atmBoundaryLayer(const atmBoundaryLayer& abl)
 :
     flowDir_(abl.flowDir_),
     zDir_(abl.zDir_),
@@ -123,7 +170,7 @@ atmBoundaryLayer::atmBoundaryLayer(const atmBoundaryLayer& abl)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void atmBoundaryLayer::autoMap(const fvPatchFieldMapper& m)
+void Foam::atmBoundaryLayer::autoMap(const fvPatchFieldMapper& m)
 {
     z0_.autoMap(m);
     zGround_.autoMap(m);
@@ -131,7 +178,7 @@ void atmBoundaryLayer::autoMap(const fvPatchFieldMapper& m)
 }
 
 
-void atmBoundaryLayer::rmap
+void Foam::atmBoundaryLayer::rmap
 (
     const atmBoundaryLayer& blptf,
     const labelList& addr
@@ -143,7 +190,10 @@ void atmBoundaryLayer::rmap
 }
 
 
-tmp<vectorField> atmBoundaryLayer::U(const vectorField& p) const
+Foam::tmp<Foam::vectorField> Foam::atmBoundaryLayer::U
+(
+    const vectorField& p
+) const
 {
     const scalarField Un
     (
@@ -162,7 +212,10 @@ tmp<vectorField> atmBoundaryLayer::U(const vectorField& p) const
 }
 
 
-tmp<scalarField> atmBoundaryLayer::k(const vectorField& p) const
+Foam::tmp<Foam::scalarField> Foam::atmBoundaryLayer::k
+(
+    const vectorField& p
+) const
 {
     tmp<scalarField> tk
     (
@@ -179,7 +232,10 @@ tmp<scalarField> atmBoundaryLayer::k(const vectorField& p) const
 }
 
 
-tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& p) const
+Foam::tmp<Foam::scalarField> Foam::atmBoundaryLayer::epsilon
+(
+    const vectorField& p
+) const
 {
     tmp<scalarField> tepsilon
     (
@@ -196,7 +252,7 @@ tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& p) const
 }
 
 
-void atmBoundaryLayer::write(Ostream& os) const
+void Foam::atmBoundaryLayer::write(Ostream& os) const
 {
     z0_.writeEntry("z0", os) ;
     os.writeKeyword("flowDir")
@@ -226,8 +282,4 @@ void atmBoundaryLayer::write(Ostream& os) const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H

@@ -128,6 +128,15 @@ namespace Foam
 
 class atmBoundaryLayer
 {
+    // Private static data
+
+        //- Default value of the Von Karman constant
+        static const scalar kappaDefault_;
+
+        //- Default value of the turbulent viscosity coefficient
+        static const scalar CmuDefault_;
+
+
     // Private data
 
         //- Flow direction
@@ -171,6 +180,12 @@ class atmBoundaryLayer
         const scalar epsilonLower_;
 
 
+    // Private member functions
+
+        //- Initialisation shared by multiple constructors
+        void init();
+
+
 public:
 
     // Constructors
@@ -178,6 +193,22 @@ public:
         //- Construct null
         atmBoundaryLayer();
 
+        //- Construct from components
+        atmBoundaryLayer
+        (
+            const vector& flowDir,
+            const vector& zDir,
+            const scalar Uref,
+            const scalar Zref,
+            const scalarField& z0,
+            const scalarField& zGround,
+            const scalar kappa = kappaDefault_,
+            const scalar Cmu = CmuDefault_,
+            const scalar ULower = 0,
+            const scalar kLower = 0,
+            const scalar epsilonLower = 0
+        );
+
         //- Construct from the coordinates field and dictionary
         atmBoundaryLayer(const vectorField& p, const dictionary&);
 

src/waves/Make/files

@@ -5,7 +5,9 @@ waveModels/Stokes2/Stokes2.C
 waveModels/Stokes5/Stokes5.C
 waveModels/solitary/solitary.C
 
-waveSuperposition/waveSuperposition.C
+waveSuperpositions/waveSuperposition/waveSuperposition.C
+waveSuperpositions/waveSuperposition/newWaveSuperposition.C
+waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.C
 
 derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C
 derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C

src/waves/Make/options

@@ -1,9 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/atmosphericModels/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
-    -ldynamicMesh
+    -ldynamicMesh \
+    -latmosphericModels

src/waves/waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+     \\/     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 "waveAtmBoundaryLayerSuperposition.H"
+#include "uniformDimensionedFields.H"
+#include "atmBoundaryLayer.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(waveAtmBoundaryLayerSuperposition, 0);
+    addToRunTimeSelectionTable
+    (
+        waveSuperposition,
+        waveAtmBoundaryLayerSuperposition,
+        objectRegistry
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::waveAtmBoundaryLayerSuperposition::waveAtmBoundaryLayerSuperposition
+(
+    const objectRegistry& db
+)
+:
+    waveSuperposition(db),
+    UGasRef_(lookup("UGasRef")),
+    hRef_(readScalar(lookup("hRef"))),
+    hWaveMin_(readScalar(lookup("hWaveMin"))),
+    hWaveMax_(readScalar(lookup("hWaveMax")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::waveAtmBoundaryLayerSuperposition::~waveAtmBoundaryLayerSuperposition()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::vectorField> Foam::waveAtmBoundaryLayerSuperposition::UGas
+(
+    const scalar t,
+    const vectorField& p
+) const
+{
+    const vector gHat =
+        normalised
+        (
+            db().lookupObject<uniformDimensionedVectorField>("g").value()
+        );
+
+    const scalar h0 = - gHat & origin();
+
+    atmBoundaryLayer atm
+    (
+        normalised(UGasRef_),
+      - gHat,
+        mag(UGasRef_),
+        h0 + hRef_,
+        scalarField(p.size(), hWaveMax_ - hWaveMin_),
+        scalarField(p.size(), h0 + hWaveMin_)
+    );
+
+    return waveSuperposition::UGas(t, p) + atm.U(p);
+}
+
+
+void Foam::waveAtmBoundaryLayerSuperposition::write(Ostream& os) const
+{
+    waveSuperposition::write(os);
+
+    os.writeKeyword("UGasRef") << UGasRef_ << token::END_STATEMENT << nl;
+    os.writeKeyword("hRef") << hRef_ << token::END_STATEMENT << nl;
+    os.writeKeyword("hWaveMin") << hWaveMin_ << token::END_STATEMENT << nl;
+    os.writeKeyword("hWaveMax") << hWaveMax_ << token::END_STATEMENT << nl;
+}
+
+
+// ************************************************************************* //

src/waves/waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.H

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+     \\/     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::waveAtmBoundaryLayerSuperposition
+
+Description
+    An extension of waveSuperposition which adds an atmospheric boundary layer
+    model to the gas velocity. The user supplies a gas velocity relative to the
+    mean liquid velocity and a height above the wave coordinate system origin
+    at which that velocity is reached. Also needed are a maximum and minimum
+    wave height which are used to set the surface roughness in the boundary
+    layer model. It is not trivial to determine these from an arbitrary
+    superposition of differently oriented wave models, so they are required as
+    user inputs instead. For a pure sinusoidal wave, the maximum and minimum
+    wave heights can be set to positive and negative amplitude, respectively.
+
+Usage
+    \table
+        Property | Description                             | Req'd? | Default
+        UGasRef  | The gas velocity relative to the \\
+                   liquid at the reference height          | yes    |
+        hRef     | The reference height relative to the \\
+                   origin of the wave coordinate system    | yes    |
+        hWaveMin | The minimum wave elevation              | yes    |
+        hWaveMax | The maximum wave elevation              | yes    |
+    \endtable
+
+    Example specification:
+    \verbatim
+    type        waveAtmBoundaryLayer;
+    // parameters for waveSuperposition ...
+    UGasRef     (10 0 0);
+    hRed        10;
+    hWaveMin    -2;
+    hWaveMax    3;
+    \endverbatim
+
+SourceFiles
+    waveAtmBoundaryLayerSuperposition.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef waveAtmBoundaryLayerSuperposition_H
+#define waveAtmBoundaryLayerSuperposition_H
+
+#include "waveSuperposition.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+             Class waveAtmBoundaryLayerSuperposition Declaration
+\*---------------------------------------------------------------------------*/
+
+class waveAtmBoundaryLayerSuperposition
+:
+    public waveSuperposition
+{
+private:
+
+    // Private member data
+
+        //- The reference velocity
+        const vector UGasRef_;
+
+        //- The reference height
+        const scalar hRef_;
+
+        //- The minimum wave elevation
+        const scalar hWaveMin_;
+
+        //- The maximum wave elevation
+        const scalar hWaveMax_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("waveAtmBoundaryLayer");
+
+
+    // Constructors
+
+        //- Construct from a database
+        waveAtmBoundaryLayerSuperposition(const objectRegistry& db);
+
+
+    //- Destructor
+    ~waveAtmBoundaryLayerSuperposition();
+
+
+    // Member Functions
+
+        //- Get the gas velocity at a given time and global positions
+        virtual tmp<vectorField> UGas
+        (
+            const scalar t,
+            const vectorField& p
+        ) const;
+
+        //- Inherit write from regIOobject
+        using regIOobject::write;
+
+        //- Write
+        void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/waves/waveSuperpositions/waveSuperposition/newWaveSuperposition.C

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+     \\/     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 "waveSuperposition.H"
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+const Foam::waveSuperposition& Foam::waveSuperposition::New
+(
+    const objectRegistry& db
+)
+{
+    if (db.foundObject<waveSuperposition>(dictName))
+    {
+        return db.lookupObject<waveSuperposition>(dictName);
+    }
+
+    const IOdictionary dict
+    (
+        IOobject
+        (
+            dictName,
+            db.time().constant(),
+            db,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    const word type =
+        dict.lookupOrDefault<word>("type", waveSuperposition::typeName);
+
+    objectRegistryConstructorTable::iterator cstrIter =
+        objectRegistryConstructorTablePtr_->find(type);
+
+    if (cstrIter == objectRegistryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown " << waveSuperposition::typeName << " " << type
+            << nl << nl << "Valid types are:" << nl
+            << objectRegistryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    waveSuperposition* ptr = cstrIter()(db).ptr();
+
+    ptr->store();
+
+    return *ptr;
+}
+
+
+// ************************************************************************* //

src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.C

@@ -25,29 +25,22 @@ License
 
 #include "waveSuperposition.H"
 #include "uniformDimensionedFields.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const Foam::word Foam::waveSuperposition::dictName("waveProperties");
 
-
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-
-const Foam::waveSuperposition& Foam::waveSuperposition::New
-(
-    const objectRegistry& db
-)
+namespace Foam
 {
-    if (db.foundObject<waveSuperposition>(dictName))
-    {
-        return db.lookupObject<waveSuperposition>(dictName);
-    }
-    else
-    {
-        waveSuperposition* ptr = new waveSuperposition(db);
-        ptr->store();
-        return *ptr;
-    }
+    defineTypeNameAndDebug(waveSuperposition, 0);
+    defineRunTimeSelectionTable(waveSuperposition, objectRegistry);
+    addToRunTimeSelectionTable
+    (
+        waveSuperposition,
+        waveSuperposition,
+        objectRegistry
+    );
 }
 
 

src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.H

@@ -153,6 +153,21 @@ class waveSuperposition
 
 public:
 
+    //- Runtime type information
+    TypeName("wave");
+
+
+    // Declare runtime construction
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        waveSuperposition,
+        objectRegistry,
+        (const objectRegistry& db),
+        (db)
+    );
+
+
     // Static Data
 
         //- The name of the dictionary
@@ -171,12 +186,6 @@ public:
         //- Construct from a database
         waveSuperposition(const objectRegistry& db);
 
-        //- Construct a copy
-        waveSuperposition(const waveSuperposition& waves);
-
-        //- Construct from a database and gravity
-        waveSuperposition(const objectRegistry& db, const dictionary& dict);
-
 
     //- Destructor
     ~waveSuperposition();
@@ -184,20 +193,54 @@ public:
 
     // Member Functions
 
+        // Access
+
+            //- Return the origin of the wave coordinate system
+            const vector& origin() const
+            {
+                return origin_;
+            }
+
+            //- Return the direction of the wave coordinate system
+            const vector& direction() const
+            {
+                return direction_;
+            }
+
         //- Get the height above the waves at a given time and global positions
-        tmp<scalarField> height(const scalar t, const vectorField& p) const;
+        virtual tmp<scalarField> height
+        (
+            const scalar t,
+            const vectorField& p
+        ) const;
 
         //- Get the liquid velocity at a given time and global positions
-        tmp<vectorField> ULiquid(const scalar t, const vectorField& p) const;
+        virtual tmp<vectorField> ULiquid
+        (
+            const scalar t,
+            const vectorField& p
+        ) const;
 
         //- Get the gas velocity at a given time and global positions
-        tmp<vectorField> UGas(const scalar t, const vectorField& p) const;
+        virtual tmp<vectorField> UGas
+        (
+            const scalar t,
+            const vectorField& p
+        ) const;
 
         //- Get the liquid pressure at a given time and global positions
-        tmp<scalarField> pLiquid(const scalar t, const vectorField& p) const;
+        virtual tmp<scalarField> pLiquid
+        (
+            const scalar t,
+            const vectorField& p
+        ) const;
 
         //- Get the gas pressure at a given time and global positions
-        tmp<scalarField> pGas(const scalar t, const vectorField& p) const;
+        virtual tmp<scalarField> pGas
+        (
+            const scalar t,
+            const vectorField& p
+        ) const;
 
         //- Inherit write from regIOobject
         using regIOobject::write;