All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2024 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
76  is related 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, initialised 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-centering 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
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-centering coefficient
226  scalar ocCoeff() const
227  {
228  return ocCoeff_->value(mesh().time().value());
229  }
230 
231  virtual tmp<VolField<Type>> fvcDdt
232  (
233  const dimensioned<Type>&
234  );
235 
236  virtual tmp<VolField<Type>> fvcDdt
237  (
238  const VolField<Type>&
239  );
240 
241  virtual tmp<VolField<Type>> fvcDdt
242  (
243  const dimensionedScalar&,
244  const VolField<Type>&
245  );
246 
247  virtual tmp<VolField<Type>> fvcDdt
248  (
249  const volScalarField&,
250  const VolField<Type>&
251  );
252 
253  virtual tmp<VolField<Type>> fvcDdt
254  (
255  const volScalarField& alpha,
256  const volScalarField& rho,
257  const VolField<Type>& psi
258  );
259 
260  virtual tmp<fvMatrix<Type>> fvmDdt
261  (
262  const VolField<Type>&
263  );
264 
265  virtual tmp<fvMatrix<Type>> fvmDdt
266  (
267  const dimensionedScalar&,
268  const VolField<Type>&
269  );
270 
271  virtual tmp<fvMatrix<Type>> fvmDdt
272  (
273  const volScalarField&,
274  const VolField<Type>&
275  );
276 
277  virtual tmp<fvMatrix<Type>> fvmDdt
278  (
279  const volScalarField& alpha,
280  const volScalarField& rho,
281  const VolField<Type>& psi
282  );
283 
285 
287  (
288  const VolField<Type>& U,
289  const SurfaceField<Type>& Uf
290  );
291 
293  (
294  const VolField<Type>& U,
295  const fluxFieldType& phi
296  );
297 
299  (
300  const volScalarField& rho,
301  const VolField<Type>& U,
302  const SurfaceField<Type>& rhoUf
303  );
304 
306  (
307  const volScalarField& rho,
308  const VolField<Type>& U,
309  const fluxFieldType& phi
310  );
311 
313  (
314  const volScalarField& alpha,
315  const volScalarField& rho,
316  const VolField<Type>& U,
317  const SurfaceField<Type>& Uf
318  )
319  {
321  return fluxFieldType::null();
322  }
323 
325  (
326  const volScalarField& alpha,
327  const volScalarField& rho,
328  const VolField<Type>& U,
329  const fluxFieldType& phi
330  )
331  {
333  return fluxFieldType::null();
334  }
335 
337  (
338  const VolField<Type>&
339  );
340 
341  virtual tmp<scalarField> meshPhi
342  (
343  const VolField<Type>&,
344  const label patchi
345  );
346 
347 
348  // Member Operators
349 
350  //- Disallow default bitwise assignment
351  void operator=(const CrankNicolsonDdtScheme&) = delete;
352 };
353 
354 
355 template<>
357 (
358  const VolField<scalar>& U,
359  const SurfaceField<scalar>& Uf
360 );
361 
362 template<>
364 (
365  const volScalarField& U,
366  const surfaceScalarField& phi
367 );
368 
369 template<>
371 (
372  const volScalarField& rho,
373  const volScalarField& U,
374  const surfaceScalarField& Uf
375 );
376 
377 template<>
379 (
380  const volScalarField& rho,
381  const volScalarField& U,
382  const surfaceScalarField& phi
383 );
384 
385 template<>
387 (
388  const volScalarField& alpha,
389  const volScalarField& rho,
390  const volScalarField& U,
391  const surfaceScalarField& Uf
392 );
393 
394 template<>
396 (
397  const volScalarField& alpha,
398  const volScalarField& rho,
399  const volScalarField& U,
400  const surfaceScalarField& phi
401 );
402 
403 
404 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
405 
406 } // End namespace fv
407 
408 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
409 
410 } // End namespace Foam
411 
412 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
413 
414 #ifdef NoRepository
415  #include "CrankNicolsonDdtScheme.C"
416 #endif
417 
418 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
419 
420 #endif
421 
422 // ************************************************************************* //
Generic GeometricField class.
static const GeometricField< Type, PatchField, GeoMesh > & null()
Return a null geometric field.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
Definition: Istream.H:60
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
Dimension set for the base types.
Definition: dimensionSet.H:125
Generic dimensioned Type class.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:99
Second-oder Crank-Nicolson implicit ddt using the current and previous time-step fields as well as th...
virtual tmp< fvMatrix< Type > > fvmDdt(const VolField< Type > &)
const fvMesh & mesh() const
Return mesh reference.
scalar ocCoeff() const
Return the current off-centering coefficient.
virtual tmp< fluxFieldType > fvcDdtPhiCorr(const VolField< Type > &U, const fluxFieldType &phi)
virtual tmp< surfaceScalarField > meshPhi(const VolField< Type > &)
ddtScheme< Type >::fluxFieldType fluxFieldType
virtual tmp< fluxFieldType > fvcDdtUfCorr(const VolField< Type > &U, const SurfaceField< Type > &Uf)
void operator=(const CrankNicolsonDdtScheme &)=delete
Disallow default bitwise assignment.
virtual tmp< VolField< Type > > fvcDdt(const dimensioned< Type > &)
CrankNicolsonDdtScheme(const fvMesh &mesh)
Construct from mesh.
TypeName("CrankNicolson")
Runtime type information.
Abstract base class for ddt schemes.
Definition: ddtScheme.H:68
const fvMesh & mesh() const
Return mesh reference.
Definition: ddtScheme.H:130
A class for managing temporary objects.
Definition: tmp.H:55
A class for handling words, derived from string.
Definition: word.H:62
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:381
label patchi
const volScalarField & psi
U
Definition: pEqn.H:72
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
Namespace for OpenFOAM.
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
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
labelList fv(nPoints)