CrankNicolsonDdtScheme.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Class
25  Foam::fv::CrankNicolsonDdtScheme
26 
27 Description
28  Second-oder Crank-Nicolson implicit ddt using the current and
29  previous time-step fields as well as the previous time-step ddt.
30 
31  The Crank-Nicolson scheme is often unstable for complex flows in complex
32  geometries and it is necessary to "off-centre" the scheme to stabilize it
33  while retaining greater temporal accuracy than the first-order
34  Euler-implicit scheme. Off-centering is specified via the mandatory
35  coefficient \c ocCoeff in the range [0,1] following the scheme name e.g.
36  \verbatim
37  ddtSchemes
38  {
39  default CrankNicolson 0.9;
40  }
41  \endverbatim
42  or with an optional "ramp" function to transition from the Euler scheme to
43  Crank-Nicolson over a initial period to avoid start-up problems, e.g.
44  \verbatim
45  ddtSchemes
46  {
47  default CrankNicolson
48  ocCoeff
49  {
50  type scale;
51  scale linearRamp;
52  duration 0.01;
53  value 0.9;
54  };
55  }
56  \endverbatim
57  With a coefficient of 1 the scheme is fully centred and second-order,
58  with a coefficient of 0 the scheme is equivalent to Euler-implicit.
59  A coefficient of 0.9 has been found to be suitable for a range of cases for
60  which higher-order accuracy in time is needed and provides similar accuracy
61  and stability to the backward scheme. However, the backward scheme has
62  been found to be more robust and provides formal second-order accuracy in
63  time.
64 
65  The advantage of the Crank-Nicolson scheme over backward is that only the
66  new and old-time values are needed, the additional terms relating to the
67  fluxes and sources are evaluated at the mid-point of the time-step which
68  provides the opportunity to limit the fluxes in such a way as to ensure
69  boundedness while maintaining greater accuracy in time compared to the
70  Euler-implicit scheme. This approach is now used with MULES in the
71  interFoam family of solvers. Boundedness cannot be guaranteed with the
72  backward scheme.
73 
74 Note
75  The Crank-Nicolson coefficient for the implicit part of the RHS is related
76  to the off-centering coefficient by
77  \verbatim
78  cnCoeff = 1.0/(1.0 + ocCoeff);
79  \endverbatim
80 
81 See also
82  Foam::Euler
83  Foam::backward
84 
85 SourceFiles
86  CrankNicolsonDdtScheme.C
87 
88 \*---------------------------------------------------------------------------*/
89 
90 #ifndef CrankNicolsonDdtScheme_H
91 #define CrankNicolsonDdtScheme_H
92 
93 #include "ddtScheme.H"
94 #include "Function1.H"
95 
96 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
97 
98 namespace Foam
99 {
100 
101 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
102 
103 namespace fv
104 {
105 
106 /*---------------------------------------------------------------------------*\
107  Class CrankNicolsonDdtScheme Declaration
108 \*---------------------------------------------------------------------------*/
109 
110 template<class Type>
112 :
113  public fv::ddtScheme<Type>
114 {
115  // Private Data
116 
117  //- Class to store the ddt0 fields on the objectRegistry for use in the
118  // next time-step. The start-time index of the CN scheme is also
119  // stored to help handle the transition from Euler to CN
120  template<class GeoField>
121  class DDt0Field
122  :
123  public GeoField
124  {
125  label startTimeIndex_;
126 
127  public:
128 
129  //- Constructor from file for restart.
130  DDt0Field
131  (
132  const IOobject& io,
133  const fvMesh& mesh
134  );
135 
136  //- Constructor from components, initisalised to zero with given
137  // dimensions.
138  DDt0Field
139  (
140  const IOobject& io,
141  const fvMesh& mesh,
143  );
144 
145  //- Return the start-time index
146  label startTimeIndex() const;
147 
148  //- Cast to the underlying GeoField
149  GeoField& operator()();
150 
151  //- Assignment to a GeoField
152  void operator=(const GeoField& gf);
153  };
154 
155 
156  //- Off-centering coefficient function
157  // 1 -> CN, less than one blends with EI
158  autoPtr<Function1<scalar>> ocCoeff_;
159 
160 
161  // Private Member Functions
162 
163  //- Disallow default bitwise copy construct
165 
166  //- Disallow default bitwise assignment
167  void operator=(const CrankNicolsonDdtScheme&);
168 
169  template<class GeoField>
170  DDt0Field<GeoField>& ddt0_
171  (
172  const word& name,
173  const dimensionSet& dims
174  );
175 
176  //- Check if the ddt0 needs to be evaluated for this time-step
177  template<class GeoField>
178  bool evaluate(const DDt0Field<GeoField>& ddt0) const;
179 
180  //- Return the coefficient for Euler scheme for the first time-step
181  // for and CN thereafter
182  template<class GeoField>
183  scalar coef_(const DDt0Field<GeoField>&) const;
184 
185  //- Return the old time-step coefficient for Euler scheme for the
186  // second time-step and for CN thereafter
187  template<class GeoField>
188  scalar coef0_(const DDt0Field<GeoField>&) const;
189 
190  //- Return the reciprocal time-step coefficient for Euler for the
191  // first time-step and CN thereafter
192  template<class GeoField>
193  dimensionedScalar rDtCoef_(const DDt0Field<GeoField>&) const;
194 
195  //- Return the reciprocal old time-step coefficient for Euler for the
196  // second time-step and CN thereafter
197  template<class GeoField>
198  dimensionedScalar rDtCoef0_(const DDt0Field<GeoField>&) const;
199 
200  //- Return ddt0 multiplied by the off-centreing coefficient
201  template<class GeoField>
202  tmp<GeoField> offCentre_(const GeoField& ddt0) const;
203 
204 
205 public:
206 
207  //- Runtime type information
208  TypeName("CrankNicolson");
209 
210 
211  // Constructors
212 
213  //- Construct from mesh
215 
216  //- Construct from mesh and Istream
217  CrankNicolsonDdtScheme(const fvMesh& mesh, Istream& is);
218 
219 
220  // Member Functions
221 
222  //- Return mesh reference
223  const fvMesh& mesh() const
224  {
225  return fv::ddtScheme<Type>::mesh();
226  }
227 
228  //- Return the current off-centreing coefficient
229  scalar ocCoeff() const
230  {
231  return ocCoeff_->value(mesh().time().value());
232  }
233 
235  (
236  const dimensioned<Type>&
237  );
238 
240  (
242  );
243 
245  (
246  const dimensionedScalar&,
248  );
249 
251  (
252  const volScalarField&,
254  );
255 
257  (
258  const volScalarField& alpha,
259  const volScalarField& rho,
261  );
262 
264  (
266  );
267 
269  (
270  const dimensionedScalar&,
272  );
273 
275  (
276  const volScalarField&,
278  );
279 
281  (
282  const volScalarField& alpha,
283  const volScalarField& rho,
285  );
288 
290  (
293  );
294 
296  (
298  const fluxFieldType& phi
299  );
300 
302  (
303  const volScalarField& rho,
306  );
307 
309  (
310  const volScalarField& rho,
312  const fluxFieldType& phi
313  );
314 
315 
317  (
319  );
320 };
321 
322 
323 template<>
325 (
328 );
329 
330 template<>
332 (
333  const volScalarField& U,
334  const surfaceScalarField& phi
335 );
336 
337 template<>
339 (
340  const volScalarField& rho,
341  const volScalarField& U,
342  const surfaceScalarField& Uf
343 );
344 
345 template<>
347 (
348  const volScalarField& rho,
349  const volScalarField& U,
350  const surfaceScalarField& phi
351 );
352 
353 
354 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
355 
356 } // End namespace fv
357 
358 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
359 
360 } // End namespace Foam
361 
362 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
363 
364 #ifdef NoRepository
365  #include "CrankNicolsonDdtScheme.C"
366 #endif
367 
368 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
369 
370 #endif
371 
372 // ************************************************************************* //
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
surfaceScalarField & phi
const fvMesh & mesh() const
Return mesh reference.
Second-oder Crank-Nicolson implicit ddt using the current and previous time-step fields as well as th...
U
Definition: pEqn.H:83
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
scalar ocCoeff() const
Return the current off-centreing coefficient.
Uf
Definition: pEqn.H:67
Generic GeometricField class.
Generic dimensioned Type class.
Abstract base class for ddt schemes.
Definition: ddtScheme.H:64
const fvMesh & mesh() const
Return mesh reference.
Definition: ddtScheme.H:135
Dimension set for the base types.
Definition: dimensionSet.H:120
ddtScheme< Type >::fluxFieldType fluxFieldType
A class for handling words, derived from string.
Definition: word.H:59
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
labelList fv(nPoints)
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const volScalarField & psi
A class for managing temporary objects.
Definition: PtrList.H:53
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
TypeName("CrankNicolson")
Runtime type information.
Namespace for OpenFOAM.