ENH: Updated waveSurfacePressure BC to allow for different time integration schemes

src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,34 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "uniformDimensionedFields.H"
+#include "EulerDdtScheme.H"
+#include "CrankNicholsonDdtScheme.H"
+#include "backwardDdtScheme.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::waveSurfacePressureFvPatchScalarField::timeSchemeType,
+        3
+    >::names[] =
+    {
+        fv::EulerDdtScheme<scalar>::typeName.c_str(),
+        fv::CrankNicholsonDdtScheme<scalar>::typeName.c_str(),
+        fv::backwardDdtScheme<scalar>::typeName.c_str()
+    };
+}
+
+
+const Foam::NamedEnum
+<
+    Foam::waveSurfacePressureFvPatchScalarField::timeSchemeType,
+    3
+>   Foam::waveSurfacePressureFvPatchScalarField::timeSchemeTypeNames_;
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -41,10 +69,8 @@ waveSurfacePressureFvPatchScalarField
 :
     fixedValueFvPatchScalarField(p, iF),
     phiName_("phi"),
-    rhoName_("rho"),
     zetaName_("zeta"),
-    zeta0_(p.size(), vector::zero),
+    rhoName_("rho")
-    curTimeIndex_(-1)
 {}
 
 
@@ -58,10 +84,8 @@ waveSurfacePressureFvPatchScalarField
 :
     fixedValueFvPatchScalarField(p, iF),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
     zetaName_(dict.lookupOrDefault<word>("zeta", "zeta")),
-    zeta0_(p.size(), vector::zero),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
-    curTimeIndex_(-1)
 {
     fvPatchField<scalar>::operator=
     (
@@ -81,10 +105,8 @@ waveSurfacePressureFvPatchScalarField
 :
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
     zetaName_(ptf.zetaName_),
-    zeta0_(ptf.zeta0_),
+    rhoName_(ptf.rhoName_)
-    curTimeIndex_(-1)
 {}
 
 
@@ -96,10 +118,8 @@ waveSurfacePressureFvPatchScalarField
 :
     fixedValueFvPatchScalarField(wspsf),
     phiName_(wspsf.phiName_),
-    rhoName_(wspsf.rhoName_),
     zetaName_(wspsf.zetaName_),
-    zeta0_(wspsf.zeta0_),
+    rhoName_(wspsf.rhoName_)
-    curTimeIndex_(wspsf.curTimeIndex_)
 {}
 
 
@@ -112,40 +132,13 @@ waveSurfacePressureFvPatchScalarField
 :
     fixedValueFvPatchScalarField(wspsf, iF),
     phiName_(wspsf.phiName_),
-    rhoName_(wspsf.rhoName_),
     zetaName_(wspsf.zetaName_),
-    zeta0_(wspsf.zeta0_),
+    rhoName_(wspsf.rhoName_)
-    curTimeIndex_(wspsf.curTimeIndex_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::waveSurfacePressureFvPatchScalarField::autoMap
-(
-    const fvPatchFieldMapper& m
-)
-{
-    fixedValueFvPatchScalarField::autoMap(m);
-    zeta0_.autoMap(m);
-}
-
-
-void Foam::waveSurfacePressureFvPatchScalarField::rmap
-(
-    const fvPatchScalarField& ptf,
-    const labelList& addr
-)
-{
-    fixedValueFvPatchScalarField::rmap(ptf, addr);
-
-    const waveSurfacePressureFvPatchScalarField& wspsf =
-        refCast<const waveSurfacePressureFvPatchScalarField>(ptf);
-
-    zeta0_.rmap(wspsf.zeta0_, addr);
-}
-
-
 void Foam::waveSurfacePressureFvPatchScalarField::updateCoeffs()
 {
     if (updated())
@@ -153,65 +146,90 @@ void Foam::waveSurfacePressureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const scalar dt = db().time().deltaTValue();
+    const label patchI = patch().index();
-    const scalar timeI = db().time().timeIndex();
-    const scalar patchI = patch().index();
 
+    const scalar dt = db().time().deltaTValue();
+
+    // retrieve non-const access to zeta field from the database
     volVectorField& zeta =
         const_cast<volVectorField&>
         (
             db().lookupObject<volVectorField>(zetaName_)
         );
-
     vectorField& zetap = zeta.boundaryField()[patchI];
 
-    if (curTimeIndex_ != timeI)
+    // lookup d/dt scheme from database for zeta
-    {
+    const word ddtSchemeName(zeta.mesh().ddtScheme(zeta.name()));
-        zeta0_ = zetap;
+    timeSchemeType timeScheme(timeSchemeTypeNames_[ddtSchemeName]);
-        curTimeIndex_ = timeI;
-    }
 
-    const surfaceScalarField& phi =
+    // retrieve the flux field from the database
-        db().lookupObject<surfaceScalarField>(phiName_);
+    const scalarField& phip =
-
+        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
-    const scalarField& phip = phi.boundaryField()[patchI];
-
-    const uniformDimensionedVectorField& g =
-        db().lookupObject<uniformDimensionedVectorField>("g");
 
+    // cache the patch face-normal vectors
     tmp<vectorField> nf(patch().nf());
 
-    if (phi.dimensions() == dimVelocity*dimArea)
+    // change in zeta due to flux
-    {
+    vectorField dZetap(dt*nf()*phip/patch().magSf());
-        zetap = zeta0_ + nf()*dt*phip/patch().magSf();
 
-        operator==(-g.value() & zetap);
+    if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
-    }
-    else if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
         const scalarField& rhop =
             patch().lookupPatchField<volScalarField, scalar>(rhoName_);
 
-        zetap = zeta0_ + nf()*dt*phip/rhop/patch().magSf();
+        dZetap /= rhop;
-
-        operator==(-rhop*(g.value() & zetap));
     }
-    else
+
+    switch (timeScheme)
     {
-        FatalErrorIn
+        case tsEuler:
-        (
+        case tsCrankNicholson:
-            "waveSurfacePressureFvPatchScalarField::updateCoeffs()"
+        {
-        )
+            zetap = zeta.oldTime().boundaryField()[patchI] + dZetap;
-            << "dimensions of phi are incorrect" << nl
+
-            << "    on patch " << this->patch().name()
+            break;
-            << " of field " << this->dimensionedInternalField().name()
+        }
-            << " in file " << this->dimensionedInternalField().objectPath()
+        case tsBackward:
-            << exit(FatalError);
+        {
+            scalar dt0 = db().time().deltaT0Value();
+
+            scalar c = 1.0 + dt/(dt + dt0);
+            scalar c00 = dt*dt/(dt0*(dt + dt0));
+            scalar c0 = c + c00;
+
+            zetap =         
+                (
+                    c0*zeta.oldTime().boundaryField()[patchI]
+                  - c00*zeta.oldTime().oldTime().boundaryField()[patchI]
+                  + dZetap
+                )/c;
+
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "waveSurfacePressureFvPatchScalarField<Type>::updateCoeffs()"
+            )   << "    Unsupported temporal differencing scheme : "
+                << ddtSchemeName << nl
+                << "    on patch " << this->patch().name()
+                << " of field " << this->dimensionedInternalField().name()
+                << " in file " << this->dimensionedInternalField().objectPath()
+                << abort(FatalError);
+        }
     }
 
-    Info<< "min/max(zetap) = " << gMin(zetap & nf()) << ", "
+
+    Info<< "min/max zetap = " << gMin(zetap & nf()) << ", "
         << gMax(zetap & nf()) << endl;
 
+    // update the surface pressure
+    const uniformDimensionedVectorField& g =
+        db().lookupObject<uniformDimensionedVectorField>("g");
+
+    operator==(-g.value() & zetap);
+
     fixedValueFvPatchScalarField::updateCoeffs();
 }
 
@@ -220,8 +238,8 @@ void Foam::waveSurfacePressureFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
     writeEntryIfDifferent<word>(os, "zeta", "zeta", zetaName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", zetaName_);
     writeEntry("value", os);
 }
 

src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,22 +25,34 @@ Class
     Foam::waveSurfacePressureFvPatchScalarField
 
 Description
-    Applies the surface wave pressure, based on the wave height
+    This is a pressure boundary condition, whose value is calculated as
+    the hydrostatic pressure based on a given displacement:
 
-        p = -rho.g.zeta
+        \f[
+            p = -rho*g*zeta
+        \f]
 
-    Where
+    where
+        \var g    = acceleration due to gravity [m/s2]
+        \var zeta = wave amplitude [m]
 
-        p = pressure - kinematic of dynamic depending on the flux units
+    The wave amplitude is updated as part of the calculation, derived from the
-        zeta = wave height vector [m]
+    local volumetric flux.
-        g = acceleration due to gravity [m/s2]
+    
-
+    Example of the boundary condition specification:
-
+    \verbatim
-    Note:
+        myPatch
-
+        {
-        This boundary also updates the wave height boundary field, which
+            type            waveSurfacePressure;
-        is accessed via lookup from the database
+            phi             phi;        // name of flux field (default = phi)
+            rho             rho;        // name of density field (default = rho)
+            zeta            zeta;       // name amplitude field (default = zeta)
+            value           uniform 0;  // place holder
+        }
+    \endverbatim
 
+    The density field is only required if the flux is mass-based as opposed to
+    volumetric-based.
 
 SourceFiles
     waveSurfacePressureFvPatchScalarField.C
@@ -51,6 +63,7 @@ SourceFiles
 #define waveSurfacePressureFvPatchScalarField_H
 
 #include "fixedValueFvPatchFields.H"
+#include "NamedEnum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,22 +78,34 @@ class waveSurfacePressureFvPatchScalarField
 :
     public fixedValueFvPatchScalarField
 {
+public:
+
+    // Public data
+
+        //- Enumeration defining the available time schemes
+        enum timeSchemeType
+        {
+            tsEuler,
+            tsCrankNicholson,
+            tsBackward
+        };
+
+
+private:
+
     // Private data
 
         //- Flux field name
         word phiName_;
 
-        //- Density field name (if compressible)
-        word rhoName_;
-
         //- Wave height field name
         word zetaName_;
 
-        //- Old-time zeta patch value
+        //- Density field for mass-based flux evaluations
-        vectorField zeta0_;
+        word rhoName_;
 
-        //- Current time index used to store ms0_
+        //- Time scheme type names
-        label curTimeIndex_;
+        static const NamedEnum<timeSchemeType, 3> timeSchemeTypeNames_;
 
 
 public:
@@ -153,22 +178,6 @@ public:
 
     // Member functions
 
-        // Mapping functions
-
-            //- Map (and resize as needed) from self given a mapping object
-            virtual void autoMap
-            (
-                const fvPatchFieldMapper&
-            );
-
-            //- Reverse map the given fvPatchField onto this fvPatchField
-            virtual void rmap
-            (
-                const fvPatchScalarField&,
-                const labelList&
-            );
-
-
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field