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-2016 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 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  With a coefficient of 1 the scheme is fully centred and second-order,
43  with a coefficient of 0 the scheme is equivalent to Euler-implicit.
44  A coefficient of 0.9 has been found to be suitable for a range of cases for
45  which higher-order accuracy in time is needed and provides similar accuracy
46  and stability to the backward scheme. However, the backward scheme has
47  been found to be more robust and provides formal second-order accuracy in
48  time.
49 
50  The advantage of the Crank-Nicolson scheme over backward is that only the
51  new and old-time values are needed, the additional terms relating to the
52  fluxes and sources are evaluated at the mid-point of the time-step which
53  provides the opportunity to limit the fluxes in such a way as to ensure
54  boundedness while maintaining greater accuracy in time compared to the
55  Euler-implicit scheme. This approach is now used with MULES in the
56  interFoam family of solvers. Boundedness cannot be guaranteed with the
57  backward scheme.
58 
59 Note
60  The Crank-Nicolson coefficient for the implicit part of the RHS is related
61  to the off-centering coefficient by
62  \verbatim
63  cnCoeff = 1.0/(1.0 + ocCoeff);
64  \endverbatim
65 
66 See also
67  Foam::Euler
68  Foam::backward
69 
70 SourceFiles
71  CrankNicolsonDdtScheme.C
72 
73 \*---------------------------------------------------------------------------*/
74 
75 #ifndef CrankNicolsonDdtScheme_H
76 #define CrankNicolsonDdtScheme_H
77 
78 #include "ddtScheme.H"
79 
80 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
81 
82 namespace Foam
83 {
84 
85 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
86 
87 namespace fv
88 {
89 
90 /*---------------------------------------------------------------------------*\
91  Class CrankNicolsonDdtScheme Declaration
92 \*---------------------------------------------------------------------------*/
93 
94 template<class Type>
96 :
97  public fv::ddtScheme<Type>
98 {
99  // Private Data
100 
101  //- Class to store the ddt0 fields on the objectRegistry for use in the
102  // next time-step. The start-time index of the CN scheme is also
103  // stored to help handle the transition from Euler to CN
104  template<class GeoField>
105  class DDt0Field
106  :
107  public GeoField
108  {
109  label startTimeIndex_;
110 
111  public:
112 
113  //- Constructor from file for restart.
114  DDt0Field
115  (
116  const IOobject& io,
117  const fvMesh& mesh
118  );
119 
120  //- Constructor from components, initisalised to zero with given
121  // dimensions.
122  DDt0Field
123  (
124  const IOobject& io,
125  const fvMesh& mesh,
127  );
128 
129  //- Return the start-time index
130  label startTimeIndex() const;
131 
132  //- Cast to the underlying GeoField
133  GeoField& operator()();
134 
135  //- Assignment to a GeoField
136  void operator=(const GeoField& gf);
137  };
138 
139 
140  //- Off-centering coefficient, 1 -> CN, less than one blends with EI
141  scalar ocCoeff_;
142 
143 
144  // Private Member Functions
145 
146  //- Disallow default bitwise copy construct
148 
149  //- Disallow default bitwise assignment
150  void operator=(const CrankNicolsonDdtScheme&);
151 
152  template<class GeoField>
153  DDt0Field<GeoField>& ddt0_
154  (
155  const word& name,
156  const dimensionSet& dims
157  );
158 
159  //- Check if the ddt0 needs to be evaluated for this time-step
160  template<class GeoField>
161  bool evaluate(const DDt0Field<GeoField>& ddt0) const;
162 
163  //- Return the coefficient for Euler scheme for the first time-step
164  // for and CN thereafter
165  template<class GeoField>
166  scalar coef_(const DDt0Field<GeoField>&) const;
167 
168  //- Return the old time-step coefficient for Euler scheme for the
169  // second time-step and for CN thereafter
170  template<class GeoField>
171  scalar coef0_(const DDt0Field<GeoField>&) const;
172 
173  //- Return the reciprocal time-step coefficient for Euler for the
174  // first time-step and CN thereafter
175  template<class GeoField>
176  dimensionedScalar rDtCoef_(const DDt0Field<GeoField>&) const;
177 
178  //- Return the reciprocal old time-step coefficient for Euler for the
179  // second time-step and CN thereafter
180  template<class GeoField>
181  dimensionedScalar rDtCoef0_(const DDt0Field<GeoField>&) const;
182 
183  //- Return ddt0 multiplied by the off-centreing coefficient
184  template<class GeoField>
185  tmp<GeoField> offCentre_(const GeoField& ddt0) const;
186 
187 
188 public:
189 
190  //- Runtime type information
191  TypeName("CrankNicolson");
192 
193 
194  // Constructors
195 
196  //- Construct from mesh
198  :
199  ddtScheme<Type>(mesh),
200  ocCoeff_(1.0)
201  {}
202 
203  //- Construct from mesh and Istream
205  :
206  ddtScheme<Type>(mesh, is),
207  ocCoeff_(readScalar(is))
208  {
209  if (ocCoeff_ < 0 || ocCoeff_ > 1)
210  {
212  (
213  is
214  ) << "Off-centreing coefficient = " << ocCoeff_
215  << " should be >= 0 and <= 1"
216  << exit(FatalIOError);
217  }
218  }
219 
220 
221  // Member Functions
222 
223  //- Return mesh reference
224  const fvMesh& mesh() const
225  {
226  return fv::ddtScheme<Type>::mesh();
227  }
228 
229  //- Return the off-centreing coefficient
230  scalar ocCoeff() const
231  {
232  return ocCoeff_;
233  }
234 
236  (
237  const dimensioned<Type>&
238  );
239 
241  (
243  );
244 
246  (
247  const dimensionedScalar&,
249  );
250 
252  (
253  const volScalarField&,
255  );
256 
258  (
259  const volScalarField& alpha,
260  const volScalarField& rho,
262  );
263 
265  (
267  );
268 
270  (
271  const dimensionedScalar&,
273  );
274 
276  (
277  const volScalarField&,
279  );
280 
282  (
283  const volScalarField& alpha,
284  const volScalarField& rho,
286  );
289 
291  (
294  );
295 
297  (
299  const fluxFieldType& phi
300  );
301 
303  (
304  const volScalarField& rho,
307  );
308 
310  (
311  const volScalarField& rho,
313  const fluxFieldType& phi
314  );
315 
316 
318  (
320  );
321 };
322 
323 
324 template<>
326 (
329 );
330 
331 template<>
333 (
334  const volScalarField& U,
335  const surfaceScalarField& phi
336 );
337 
338 template<>
340 (
341  const volScalarField& rho,
342  const volScalarField& U,
343  const surfaceScalarField& Uf
344 );
345 
346 template<>
348 (
349  const volScalarField& rho,
350  const volScalarField& U,
351  const surfaceScalarField& phi
352 );
353 
354 
355 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
356 
357 } // End namespace fv
358 
359 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
360 
361 } // End namespace Foam
362 
363 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
364 
365 #ifdef NoRepository
366  #include "CrankNicolsonDdtScheme.C"
367 #endif
368 
369 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
370 
371 #endif
372 
373 // ************************************************************************* //
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
surfaceScalarField & phi
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const fvMesh & mesh() const
Return mesh reference.
Definition: ddtScheme.H:135
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
Uf
Definition: pEqn.H:67
Generic GeometricField class.
Generic dimensioned Type class.
Abstract base class for ddt schemes.
Definition: ddtScheme.H:64
scalar ocCoeff() const
Return the off-centreing coefficient.
Dimension set for the base types.
Definition: dimensionSet.H:118
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)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
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 > &)
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:331
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
const fvMesh & mesh() const
Return mesh reference.
const volScalarField & psi
A class for managing temporary objects.
Definition: PtrList.H:54
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:91
TypeName("CrankNicolson")
Runtime type information.
Namespace for OpenFOAM.
IOerror FatalIOError