Commit 92b6f3dc authored by Henry Weller's avatar Henry Weller
Browse files

CrankNicolsonDdtScheme: Added option "ramp" function to smoothly transition from Euler

    Off-centering is specified via the mandatory coefficient \c ocCoeff in the
    range [0,1] following the scheme name e.g.
    \verbatim
    ddtSchemes
    {
        default         CrankNicolson 0.9;
    }
    \endverbatim
    or with an optional "ramp" function to transition from the Euler scheme to
    Crank-Nicolson over a initial period to avoid start-up problems, e.g.
    \verbatim
    ddtSchemes
    {
        default         CrankNicolson
        ocCoeff
        {
            type scale;
            scale linearRamp;
            duration 0.01;
            value 0.9;
        };
    }
    \endverbatim

Note this functionality is experimental and the specification and implementation
may change if issues arise.
parent aae3705b
......@@ -200,7 +200,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef_
{
if (mesh().time().timeIndex() > ddt0.startTimeIndex())
{
return 1 + ocCoeff_;
return 1 + ocCoeff();
}
else
{
......@@ -218,7 +218,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef0_
{
if (mesh().time().timeIndex() > ddt0.startTimeIndex() + 1)
{
return 1 + ocCoeff_;
return 1 + ocCoeff();
}
else
{
......@@ -256,9 +256,9 @@ tmp<GeoField> CrankNicolsonDdtScheme<Type>::offCentre_
const GeoField& ddt0
) const
{
if (ocCoeff_ < 1)
if (ocCoeff() < 1)
{
return ocCoeff_*ddt0;
return ocCoeff()*ddt0;
}
else
{
......
......@@ -32,13 +32,28 @@ Description
geometries and it is necessary to "off-centre" the scheme to stabilize it
while retaining greater temporal accuracy than the first-order
Euler-implicit scheme. Off-centering is specified via the mandatory
coefficient in the range [0,1] following the scheme name e.g.
coefficient \c ocCoeff in the range [0,1] following the scheme name e.g.
\verbatim
ddtSchemes
{
default CrankNicolson 0.9;
}
\endverbatim
or with an optional "ramp" function to transition from the Euler scheme to
Crank-Nicolson over a initial period to avoid start-up problems, e.g.
\verbatim
ddtSchemes
{
default CrankNicolson
ocCoeff
{
type scale;
scale linearRamp;
duration 0.01;
value 0.9;
};
}
\endverbatim
With a coefficient of 1 the scheme is fully centred and second-order,
with a coefficient of 0 the scheme is equivalent to Euler-implicit.
A coefficient of 0.9 has been found to be suitable for a range of cases for
......@@ -76,6 +91,7 @@ SourceFiles
#define CrankNicolsonDdtScheme_H
#include "ddtScheme.H"
#include "Constant.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -138,7 +154,7 @@ class CrankNicolsonDdtScheme
//- Off-centering coefficient, 1 -> CN, less than one blends with EI
scalar ocCoeff_;
autoPtr<Function1<scalar>> ocCoeff_;
// Private Member Functions
......@@ -197,7 +213,7 @@ public:
CrankNicolsonDdtScheme(const fvMesh& mesh)
:
ddtScheme<Type>(mesh),
ocCoeff_(1.0)
ocCoeff_(new Function1Types::Constant<scalar>("ocCoeff", 1))
{
// Ensure the old-old-time cell volumes are available
// for moving meshes
......@@ -210,17 +226,34 @@ public:
//- Construct from mesh and Istream
CrankNicolsonDdtScheme(const fvMesh& mesh, Istream& is)
:
ddtScheme<Type>(mesh, is),
ocCoeff_(readScalar(is))
ddtScheme<Type>(mesh, is)
{
if (ocCoeff_ < 0 || ocCoeff_ > 1)
token firstToken(is);
if (firstToken.isNumber())
{
FatalIOErrorInFunction
const scalar ocCoeff = firstToken.scalarToken();
if (ocCoeff < 0 || ocCoeff > 1)
{
FatalIOErrorInFunction
(
is
) << "Off-centreing coefficient = " << ocCoeff
<< " should be >= 0 and <= 1"
<< exit(FatalIOError);
}
ocCoeff_ = new Function1Types::Constant<scalar>
(
is
) << "Off-centreing coefficient = " << ocCoeff_
<< " should be >= 0 and <= 1"
<< exit(FatalIOError);
"ocCoeff",
ocCoeff
);
}
else
{
is.putBack(firstToken);
dictionary dict(is);
ocCoeff_ = Function1<scalar>::New("ocCoeff", dict);
}
// Ensure the old-old-time cell volumes are available
......@@ -243,7 +276,7 @@ public:
//- Return the off-centreing coefficient
scalar ocCoeff() const
{
return ocCoeff_;
return ocCoeff_->value(mesh().time().value());
}
tmp<GeometricField<Type, fvPatchField, volMesh>> fvcDdt
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment