CrankNicolsonDdtScheme.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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 stabilise 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  template<class GeoField>
164  DDt0Field<GeoField>& ddt0_
165  (
166  const word& name,
167  const dimensionSet& dims
168  );
169 
170  //- Check if the ddt0 needs to be evaluated for this time-step
171  template<class GeoField>
172  bool evaluate(const DDt0Field<GeoField>& ddt0) const;
173 
174  //- Return the coefficient for Euler scheme for the first time-step
175  // for and CN thereafter
176  template<class GeoField>
177  scalar coef_(const DDt0Field<GeoField>&) const;
178 
179  //- Return the old time-step coefficient for Euler scheme for the
180  // second time-step and for CN thereafter
181  template<class GeoField>
182  scalar coef0_(const DDt0Field<GeoField>&) const;
183 
184  //- Return the reciprocal time-step coefficient for Euler for the
185  // first time-step and CN thereafter
186  template<class GeoField>
187  dimensionedScalar rDtCoef_(const DDt0Field<GeoField>&) const;
188 
189  //- Return the reciprocal old time-step coefficient for Euler for the
190  // second time-step and CN thereafter
191  template<class GeoField>
192  dimensionedScalar rDtCoef0_(const DDt0Field<GeoField>&) const;
193 
194  //- Return ddt0 multiplied by the off-centreing coefficient
195  template<class GeoField>
196  tmp<GeoField> offCentre_(const GeoField& ddt0) const;
197 
198 
199 public:
200 
201  //- Runtime type information
202  TypeName("CrankNicolson");
203 
204 
205  // Constructors
206 
207  //- Construct from mesh
209 
210  //- Construct from mesh and Istream
211  CrankNicolsonDdtScheme(const fvMesh& mesh, Istream& is);
212 
213  //- Disallow default bitwise copy construction
215 
216 
217  // Member Functions
218 
219  //- Return mesh reference
220  const fvMesh& mesh() const
221  {
222  return fv::ddtScheme<Type>::mesh();
223  }
224 
225  //- Return the current off-centreing coefficient
226  scalar ocCoeff() const
227  {
228  return ocCoeff_->value(mesh().time().value());
229  }
230 
232  (
233  const dimensioned<Type>&
234  );
235 
237  (
239  );
240 
242  (
243  const dimensionedScalar&,
245  );
246 
248  (
249  const volScalarField&,
251  );
252 
254  (
255  const volScalarField& alpha,
256  const volScalarField& rho,
258  );
259 
260  virtual tmp<fvMatrix<Type>> fvmDdt
261  (
263  );
264 
265  virtual tmp<fvMatrix<Type>> fvmDdt
266  (
267  const dimensionedScalar&,
269  );
270 
271  virtual tmp<fvMatrix<Type>> fvmDdt
272  (
273  const volScalarField&,
275  );
276 
277  virtual tmp<fvMatrix<Type>> fvmDdt
278  (
279  const volScalarField& alpha,
280  const volScalarField& rho,
282  );
285 
287  (
290  );
291 
293  (
295  const fluxFieldType& phi
296  );
297 
299  (
300  const volScalarField& rho,
303  );
304 
306  (
307  const volScalarField& rho,
309  const fluxFieldType& phi
310  );
311 
312 
314  (
316  );
317 
318 
319  // Member Operators
320 
321  //- Disallow default bitwise assignment
322  void operator=(const CrankNicolsonDdtScheme&) = delete;
323 };
324 
325 
326 template<>
328 (
331 );
332 
333 template<>
335 (
336  const volScalarField& U,
337  const surfaceScalarField& phi
338 );
339 
340 template<>
342 (
343  const volScalarField& rho,
344  const volScalarField& U,
345  const surfaceScalarField& Uf
346 );
347 
348 template<>
350 (
351  const volScalarField& rho,
352  const volScalarField& U,
353  const surfaceScalarField& phi
354 );
355 
356 
357 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
358 
359 } // End namespace fv
360 
361 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
362 
363 } // End namespace Foam
364 
365 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
366 
367 #ifdef NoRepository
368  #include "CrankNicolsonDdtScheme.C"
369 #endif
370 
371 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
372 
373 #endif
374 
375 // ************************************************************************* //
virtual tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
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...
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
void evaluate(GeometricField< Type, PatchField, GeoMesh > &result, const Function1< Type > &func, const GeometricField< Type, PatchField, GeoMesh > &x)
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.
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
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:129
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
virtual tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
labelList fv(nPoints)
phi
Definition: correctPhi.H:3
autoPtr< surfaceVectorField > Uf
virtual tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
U
Definition: pEqn.H:72
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
CrankNicolsonDdtScheme(const fvMesh &mesh)
Construct from mesh.
virtual 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
void operator=(const CrankNicolsonDdtScheme &)=delete
Disallow default bitwise assignment.
A class for managing temporary objects.
Definition: PtrList.H:53
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.