CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Crank Nicolson scheme implemented wrong?

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree14Likes
  • 1 Post By rajibroy
  • 1 Post By rajibroy
  • 1 Post By rajibroy
  • 2 Post By guin
  • 1 Post By Tobi
  • 4 Post By rajibroy
  • 1 Post By HPE
  • 2 Post By rajibroy
  • 1 Post By guin

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 17, 2020, 22:20
Default Crank Nicolson scheme implemented wrong?
  #1
New Member
 
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 8
rajibroy is on a distinguished road
Edit 2: 2020 May 04

OpenFOAM CrankNicolson implementation is correct when at each time step the linear system residual is zero. For derivation, please see forum post 10.


Original post:
The Crank-Nicolson temporal discretization scheme for any variable \phi is defined as

\frac{\partial{\phi}}{\partial t} = \frac{\phi - \phi_{00}}{2 \Delta t}, where \phi_{00} is the \phi value two time-steps earlier.

There is potential error in OpenFOAM v7 CrankNicolsonDdtScheme implementation. The off-center coefficient for blending with Euler implicit scheme is defined as

Code:
ddtSchemes
{
     default         CrankNicolson ocCoeff; // ocCoeff [0, 1]
}
The Crank Nicolson coefficient is defined by

Code:
cnCoeff = 1.0/(1.0 + ocCoeff);
ocCoeff = 1 (cnCoeff = 0.5) mean the standard CrankNicolson scheme and ocCoeff = 0 (cnCoeff = 1) means an inplicit Euler scheme. Any value in between 0 and 1 will produce a linearly blended scheme.

However, within the implementation cnCoeff is not declared/defined, and a separate variable coef defined as

Code:
coef_ = 1 + ocCoeff;
, which if used instead of cnCoeff will create a very different temporal discretization scheme.


Let us investigate how this apparent bug affects the discretization scheme (for fvcDdt member function and mesh not moving):

Code:
//- Off-centering coefficient function
//  1 -> CN, less than one blends with EI
autoPtr<Function1<scalar>> ocCoeff_;


//- Return the coefficient for Euler scheme for the first time-step
//  for and CN thereafter
coef_ = 1 + ocCoeff();


//- Return the reciprocal time-step coefficient for Euler for the
//  first time-step and CN thereafter
rDtCoef_ = coef_/deltaT();


//- Return ddt0 multiplied by the off-centreing coefficient
offCentre_(ddt0()) = ocCoeff()*ddt0;


// previous temporal derivative
ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime()) - offCentre_(ddt0());

// explicit temporal derivative calculation
fvcDdt(vf) = rDtCoef*(vf - vf.oldTime()) - offCentre_(ddt0())
Condensing into one equation, we get
Code:
fvcDdt(vf) = (1 + ocCoeff) * (vf - vf.oldTime())/deltaT - ocCoeff * fvcDdt(vf.oldTime())
\frac{\partial{\phi}}{\partial t} =    (1 + ocCoeff) *\frac{\phi - \phi_{0}}{\Delta t} - ocCoeff *  \frac{\partial{\phi_0}}{\partial t}

Approximating the earlier temporal derivative as implicit Euler
Code:
fvcDdt(vf.oldTime()) = (vf - vf.oldTime())/deltaT0
The discretization can be written as:
\frac{\partial{\phi}}{\partial t} =    (1 + ocCoeff) *\frac{\phi - \phi_{0}}{\Delta t}  - ocCoeff  *\frac{\phi_{0} - \phi_{00}}{\Delta t0}

For uniform time step, the standard Crank Nicolson scheme (ocCoeff = 1) becomes

\frac{\partial{\phi}}{\partial t} =    2*\frac{\phi - \phi_{0}}{\Delta t}  - 1 *\frac{\phi_{0} - \phi_{00}}{\Delta t}

\frac{\partial{\phi}}{\partial t} =    \frac{2 \phi - 3 \phi_{0} + \phi_{00}}{\Delta t}, which is very different from Crank Nicolson scheme.

This could be solved if the cnCoeff is used instead of coeff!
\frac{\partial{\phi}}{\partial t} =    \frac{1}{(1 + ocCoeff)} *\frac{\phi - \phi_{0}}{\Delta t}  + \frac{ocCoeff}{(1 + ocCoeff)}  *\frac{\phi_{0} - \phi_{00}}{\Delta t0}

With uniform time step, the standard Crank Nicolson scheme (ocCoeff = 1) can be written as
\frac{\partial{\phi}}{\partial t} =    \frac{\phi - \phi_{0}}{2 \Delta t}  + \frac{\phi_{0} - \phi_{00}}{ 2 \Delta t}

\frac{\partial{\phi}}{\partial t} =    \frac{\phi - \phi_{00}}{2 \Delta t}

I have confidence in my own analysis, but it is hard to believe OpenFOAM developers and users missed such a crucial bug. I will be happy to proven wrong!

Please review, and am looking forward to your feedback!

Edit: For (ocCoeff = 0.5), the scheme becomes:

OpenFOAM implementation:

\frac{\partial{\phi}}{\partial t} =  \frac{3}{2}  \frac{\phi - \phi_{0}}{\Delta t}  - \frac{1}{2}  \frac{\phi_{0} - \phi_{00}}{ \Delta t}

\frac{\partial{\phi}}{\partial t} =    \frac{1.5 \phi - 2 \phi_{0} + 0.5 \phi_{00}}{\Delta t}

Proposed correction:

\frac{\partial{\phi}}{\partial t} =  \frac{2}{3}  \frac{\phi - \phi_{0}}{\Delta t}  + \frac{1}{3}  \frac{\phi_{0} - \phi_{00}}{ \Delta t}


\frac{\partial{\phi}}{\partial t} = \frac{  \frac{2}{3} \phi - \frac{1}{3} \phi_{0} - \frac{1}{3} \phi_{00}  } { \Delta t}
Tobi likes this.

Last edited by rajibroy; May 4, 2020 at 20:03. Reason: Edit 1: Added an example; Edit 2: OF CN scheme is correct (see post 10)
rajibroy is offline   Reply With Quote

Old   April 27, 2020, 08:24
Default Verification: temporal discretization coefficients
  #2
New Member
 
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 8
rajibroy is on a distinguished road
To verify the Crank-Nicolson scheme, we can printout the the discretization coefficients:

Code:
template<class Type>
tmp<fvMatrix<Type>>
CrankNicolsonTestDdtScheme<Type>::fvmDdt
(
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
.
.
.
Info  << "deltaT" << deltaT()
        << "\tocCoeff: "     << ocCoeff() 
        << ",\tcoef: "     << coef_(ddt0)
        << ",\trDtCoef: "  << rDtCoef
         // << ",\trDtCoef0: " << rDtCoef0_(ddt0).value() 
         << endl;
}
For off-center coefficient 0.5 and uniform time-step of\Delta t = 0.2 s, the OF v7 implementation provides:

Code:
deltaT: 0.2  ocCoeff: 0.5,	coef: 1.5,	rDtCoef: 7.5,
which verifies inaccurate implementation as described in the post above (example with off-center coefficient 0.5).

For corrected implementation, the coefficients are:
Code:
deltaT: 0.2 cnCoeff: 0.67 coefft: 0.67 coefft0: -0.33 coefft00: -0.33,
which matches the correct implementation above.
Tobi likes this.
rajibroy is offline   Reply With Quote

Old   April 27, 2020, 08:43
Default use of previous temporal derivative is potentially inaccurate
  #3
New Member
 
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 8
rajibroy is on a distinguished road
The Crank Nicolson implementation in OpenFOAM uses previous temporal derivative \frac{\partial \phi_{0}}{\partial t_{0}} instead of \phi_{00}. Either way, the scheme need to store an additional field. There is no additional benefit of saving \frac{\partial \phi_{0}}{\partial t_{0}} , although it may lead to potential inaccurate temporal derivative!

Here is how (simplified using uniform time step and ocCoef):

OF Crank-Niconson implementation
\frac{\partial{\phi}}{\partial t} =    (1 + ocCoeff) *\frac{\phi - \phi_{0}}{\Delta t} - ocCoeff *  \frac{\partial{\phi_0}}{\partial t}

\frac{\partial{\phi}}{\partial t} = f \left(  \phi, \phi_{0}, \frac{\partial{\phi_0}}{\partial t} \right)

Similarly, the \frac{\partial{\phi_{0}}}{\partial t} = f \left(  \phi_{0}, \phi_{00}, \frac{\partial{\phi_{00}}}{\partial t} \right).

Thus, at any time step n, the temporal derivative involves all phi values from 0 to n.

\frac{\partial{\phi_{n}}}{\partial t} = f \left(  \phi_{n}, \phi_{n-1},  \phi_{n-2}, \phi_{n-3}, \cdot \cdot \cdot , \phi_{2}, \phi_{1}, \phi_{0}  \right),

which is clearly a deviation from Crank-Nicolson discrete temporal discretization scheme.
Tobi likes this.
rajibroy is offline   Reply With Quote

Old   April 27, 2020, 10:37
Default
  #4
Member
 
Rodrigo
Join Date: Mar 2010
Posts: 73
Rep Power: 12
guin is on a distinguished road
Hi Rajib,

Did you think about bringing this up to the development team?
https://bugs.openfoam.org/rules.php

It may turn out that this deviation from original Crank-Nicolson was accidental or even prosecuted. Either the way, clearing this may be in benefit for the OF-community.

Hint: You will be probably asked for a case or a reproducible set of instructions that probe this deviated behaviour.

Best & thanks for your effort!
Tobi and rajibroy like this.
guin is offline   Reply With Quote

Old   April 27, 2020, 11:58
Default
  #5
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,390
Blog Entries: 6
Rep Power: 42
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Dear Rajib,

first of all, thank you for your clear post and the usage of code tags. Your hints are clear and we can easily follow. I checked the code just the last 5 minutes but I cannot dive into that scheme right now for two reasons:
  • I am not sure why we need the old-old time step here?

As far as I remember, the Crank Nic. scheme uses the actual and old value. Regarding the coefficient. I need to find some reference for that.
rajibroy likes this.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   April 27, 2020, 13:41
Default
  #6
New Member
 
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 8
rajibroy is on a distinguished road
Hello Tobi,
I have implemented the Crank-Nicolson scheme following the explanation in the section 13.3..5, page 512 of the book

Moukalled, F., Mangani, L., & Darwish, M. (2015). The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAMŪ and Matlab. https://doi.org/10.1007/978-3-319-16

If you need a copy of the chapter, could you please pm an email address? If you come across references for CrankNicolson scheme for FVM implementation, would you please post in the thread?

---------------------------------------------------------------------------
Hello Rodrigo,
Thank you for your kind suggestion.

Prior to reaching out to the OF developers, I am reaching out to the community users to get feedback so that bug report is comprehensive.

Regarding a test case, I might need to develop a one-dimensional wave equation. Your kind suggestion is much appreciated.

Thanks again both of you!

Regards,
Rajib
rajibroy is offline   Reply With Quote

Old   April 29, 2020, 23:24
Default Comparison of temporal discretization schemes in OpenFOAM
  #7
New Member
 
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 8
rajibroy is on a distinguished road
A 1D scalar transport equation is used to verify the temporal discretization schemes, where a passive scalar (phi) is convected along the Cartesian x-axis using a uniform flow velocity of (1, 0, 0):

The schemes used for temporal discretization are:

\frac{\partial \phi}{\partial t} + U \frac{\partial \phi}{\partial x} = 0.

Domain extends [0,1] in a Cartesian XY plane. A sinusoidal field [ \phi (y) = sin(2 \pi y) ] is applied on the left (x = 0) at time 0, and convected toward the right (along x axis).


\text{ Euler implicit} \quad \frac{\partial \phi}{\partial t} = \frac{\phi -\phi_0}{\Delta t}

\text{ backward (SOU)} \quad  \frac{\partial \phi}{\partial t} = \frac{3\phi -4\phi_0 + \phi_{00}}{2 \Delta t}

\text{ Crank Nicolson (OFv6)} \quad \frac{\partial \phi} {\partial t} = 2* \frac{\phi -\phi_0}{\Delta t} -  \frac{\partial \phi_{0}}{\partial t}

\text{ Crank Nicolson (Modified)} \quad  \frac{\partial \phi}{\partial t} = \frac{\phi -\phi_{00}}{2 \Delta t}

After 1 characteristic time, scalar field values are plotted along line y = 0.25, which is analytically 1. The scheme comparison plot is attached.

Analysis:
The CrankNicolson scheme implemented in OpemFOAM is behaving oscillatory similar to the 2nd order upwind (aka backward) scheme. The modified Crank-Nicolson (CNMod) scheme implemented by the author (as described in the first post) is oscillation free.

However, the modified Crank-Nicolson scheme have the same accuracy as the Euler scheme. I was expecting better accuracy form the CNMod scheme.

The case study is available at GitHub https://github.com/rroyCFD/temporalD...rification.git

Please review and I look forward to your comments.
Attached Images
File Type: jpg temporalSchemeComparison.jpg (76.5 KB, 26 views)
rajibroy is offline   Reply With Quote

Old   April 30, 2020, 12:13
Default
  #8
HPE
Senior Member
 
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 619
Rep Power: 7
HPE is on a distinguished road
is there any chance for you to report your findings in the issue tracker of one of the OpenFOAM variants? links below may guide you.
rajibroy likes this.
HPE is offline   Reply With Quote

Old   May 1, 2020, 11:09
Default Update on Crank-Niconson scheme spatial discretization treatment
  #9
New Member
 
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 8
rajibroy is on a distinguished road
After further study, I realized that any temporal discretization scheme shall be discussed in conjunction with explicit/implicit treatment of the spatial discretization operators.

The Crank-Nicolson (CN) scheme implementation in OpenFOAM doesn't interact with spatial discretization operators. Hence, the CN scheme can't be applied in the default solvers (i.e. the user have to modify the momentum/transport equation of the solvers).

The two types of CN schemes are:

\text{CN semi-implicit: } \frac{\phi - \phi_{0}} {\Delta t} = 0.5 \left( \mathbf{S} + \mathbf{S_{0}} \right)
and
\text{CN explicit: } \frac{\phi - \phi_{00}} {2 \Delta t} = \mathbf{S_{0}},

where \mathbf{S} is spatial discretization operator.

If we introduce, a CN coefficient (off-center coefficient > ocCoeff) from blending with implicit Euler scheme, the reformulated schemes will be:

ocCoef \in [1, 0] \quad 1: \text{semi-Implicit CN} \quad 0: \text{implicit Euler}

cnCoef = \frac{1}{1 + ocCoef}; \quad ocCoef \in [0.5, 1]

\text{CN semi-implicit: } \frac{\phi - \phi_{0}} {\Delta t} = \left( cnCoef* \mathbf{S} + (1 -cnCoef )* \mathbf{S_{0}} \right)

The correction I proposed in the earlier post, following the OpenFOAM implementation, is neither semi-implicit nor explicit. The implementation intended to blend temporall discretization operators of the explicit CN and implicit Euler schemes while all spatial operators are treated as implicit (where possible).

Please find the attached plot of the scalar transport of the sinusoidal wave. The explicit CN scheme is not stable (super oscillatory) in absence of diffusivity (thus not plotted). Both inaccurate OpenFOAM implementation and semi-implicit implementation of the CN scheme give same result as the backward scheme (which might misled the OF developers).

I look forward to your comments, and planning to submit a bug report with OpenFOAM.
Attached Images
File Type: jpg temporalSchemeComparisonNew.jpg (72.6 KB, 16 views)
rajibroy is offline   Reply With Quote

Old   May 4, 2020, 20:00
Default OpenFOAM CrankNicolson implementation is correct when linear system residual is zero
  #10
New Member
 
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 8
rajibroy is on a distinguished road
The Crank-Nicolson temporal discretization scheme can be viewed as summation of implicit and explicit Euler discretization schemes.

For any variable \phi, let assume the linear system for the partial differential equation is written as:

\frac{ \partial{\phi}}{\partial t} = \mathbf{S}(\phi),
where \mathbf{S} is the spatial discretization operator.

The temporal derivatives can discretized as :

Implicit Euler
\frac{ \partial{\phi^{t +  \frac{1}{2}}}} {\partial t} =  \mathbf{S}(\phi^{t +  \frac{1}{2}})

\frac{\phi^{t +  \frac{1}{2}} - \phi^t}{(\Delta t/2)} =  \mathbf{S}(\phi^{t +  \frac{1}{2}})



Explicit Euler
\frac{ \partial{\phi^{t + 1}}}{\partial t} =  \mathbf{S}(\phi^{t +  \frac{1}{2}})

\frac{\phi^{t + 1} - \phi^{t +  \frac{1}{2}}}{(\Delta t/2)} =  \mathbf{S}(\phi^{t +  \frac{1}{2}})


Crank Nicolson
Adding equations {Eqn:implEuler} and {Eqn:explEuler}, we get

\frac{\phi^{t +  \frac{1}{2}} - \phi^t}{(\Delta t/2)} + \frac{\phi^{t + 1} - \phi^{t +  \frac{1}{2}}}{(\Delta t/2)} = 2 * \mathbf{S}(\phi^{t +  \frac{1}{2}})



\frac{\phi^{t + 1} - \phi^t}{(\Delta t)} = \mathbf{S}(\phi^{t +  \frac{1}{2}})


With central difference approximation for \mathbf{S}(\phi^{t +1/2}), the final form of the discretization becomes:

\frac{\phi^{t + 1} - \phi^t}{\Delta t} =  \frac{1}{2} \left[ \mathbf{S}(\phi^{t + 1}) + \mathbf{S}(\phi^{t}) \right]


The temporal term is similar to the implicit Euler scheme, and key difference is in semi implicit treatment of the spatial operators.

Now, let's assume that an unsteady simulation is performed where the linear system is solved to tight tolerance at each time step before proceeding to next time-step. Thus, equation {Eqn:governingEqn} at any time step t can be approximated as,

\frac{ \partial{\phi^t}} {\partial t} = \mathbf{S}(\phi^t) + \mathbf{R}(t),

where \mathbf{R}(t) is the final residual of the linear system at time-step t. We can also express the spatial operators in terms of the temporal derivative and the final residual. For example, the equation {Eqn:approxEqn} can be reformulated as

\mathbf{S}(\phi^t) = \frac{ \partial{\phi^t}}{\partial t} - \mathbf{R}(t).

Substituting, \mathbf{S}(\phi^t) value in the Crank-Nicolson equation {Eqn:implEuler}, we get

\frac{\phi^{t + 1} - \phi^t}{\Delta t} =  \frac{1}{2} \left[ \mathbf{S}(\phi^{t + 1}) + \frac{ \partial{\phi^t}} {\partial t} - \mathbf{R}(t) \right],


2 \frac{\phi^{t + 1} - \phi^t}{\Delta t} - \frac{ \partial{\phi^t}} {\partial t} + \mathbf{R}(t) = \mathbf{S}(\phi^{t + 1}).

OpenFOAM applies this form of the CrankNicolson scheme with assumption that the linear system residual are driven to machine zero \mathbf{R}(t) \rightarrow 0 or atleast insignificantly small.

Usually in PISO like segregated algorithm, driving the linear system residual to machine zero is very expensive and computationally inefficient. It might be the reason the CrankNicolson scheme is recommended to be used with small blending with implicit Euler scheme.

I tested the Crank-Nicolson scheme in a projection algorithm (one outer loop only) where the initial residual are small but higher than PISO type algorithm multiple outer loops. The Crank-Nicolson scheme in the format of equation{Eqn:CNOF} is unstable, but original format (equation {Eqn:CN} is stable.
Attached Files
File Type: pdf CrankNicolsonDdtScheme_deriverdBy_RajibRoy.pdf (92.2 KB, 13 views)
guin and HPE like this.

Last edited by rajibroy; May 5, 2020 at 14:35.
rajibroy is offline   Reply With Quote

Old   May 5, 2020, 09:57
Default
  #11
Member
 
Rodrigo
Join Date: Mar 2010
Posts: 73
Rep Power: 12
guin is on a distinguished road
Here the closed report for future reference:

https://bugs.openfoam.org/view.php?id=3490
HPE likes this.
guin is offline   Reply With Quote

Reply

Tags
crank-nicolson, openfoam, temporal discretization

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
udf error srihari FLUENT 1 October 31, 2016 14:18
Crank Nicolson scheme msrinath80 OpenFOAM Running, Solving & CFD 6 November 14, 2010 13:59
Could someone explain the implementation of the convection scheme liuhuafei OpenFOAM Running, Solving & CFD 0 January 20, 2008 19:18
extrapolation in MUSCL scheme Chandra Main CFD Forum 6 February 14, 2007 11:21
High order difference scheme for KIVA3v2??? LiQiang Main CFD Forum 1 May 9, 2005 13:50


All times are GMT -4. The time now is 11:13.