
[Sponsors] 
Crank Nicolson scheme implemented wrong? 

LinkBack  Thread Tools  Search this Thread  Display Modes 
April 17, 2020, 22:20 
Crank Nicolson scheme implemented wrong?

#1 
New Member
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 9 
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 CrankNicolson temporal discretization scheme for any variable is defined as , where is the value two timesteps earlier. There is potential error in OpenFOAM v7 CrankNicolsonDdtScheme implementation. The offcenter coefficient for blending with Euler implicit scheme is defined as Code:
ddtSchemes { default CrankNicolson ocCoeff; // ocCoeff [0, 1] } Code:
cnCoeff = 1.0/(1.0 + ocCoeff); However, within the implementation cnCoeff is not declared/defined, and a separate variable coef defined as Code:
coef_ = 1 + ocCoeff; Let us investigate how this apparent bug affects the discretization scheme (for fvcDdt member function and mesh not moving): Code:
// Offcentering coefficient function // 1 > CN, less than one blends with EI autoPtr<Function1<scalar>> ocCoeff_; // Return the coefficient for Euler scheme for the first timestep // for and CN thereafter coef_ = 1 + ocCoeff(); // Return the reciprocal timestep coefficient for Euler for the // first timestep and CN thereafter rDtCoef_ = coef_/deltaT(); // Return ddt0 multiplied by the offcentreing 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()) Code:
fvcDdt(vf) = (1 + ocCoeff) * (vf  vf.oldTime())/deltaT  ocCoeff * fvcDdt(vf.oldTime()) Approximating the earlier temporal derivative as implicit Euler Code:
fvcDdt(vf.oldTime()) = (vf  vf.oldTime())/deltaT0 For uniform time step, the standard Crank Nicolson scheme (ocCoeff = 1) becomes , which is very different from Crank Nicolson scheme. This could be solved if the cnCoeff is used instead of coeff! With uniform time step, the standard Crank Nicolson scheme (ocCoeff = 1) can be written as 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: Proposed correction: 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) 

April 27, 2020, 08:24 
Verification: temporal discretization coefficients

#2 
New Member
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 9 
To verify the CrankNicolson 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; } Code:
deltaT: 0.2 ocCoeff: 0.5, coef: 1.5, rDtCoef: 7.5, For corrected implementation, the coefficients are: Code:
deltaT: 0.2 cnCoeff: 0.67 coefft: 0.67 coefft0: 0.33 coefft00: 0.33, 

April 27, 2020, 08:43 
use of previous temporal derivative is potentially inaccurate

#3 
New Member
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 9 
The Crank Nicolson implementation in OpenFOAM uses previous temporal derivative instead of . Either way, the scheme need to store an additional field. There is no additional benefit of saving , although it may lead to potential inaccurate temporal derivative!
Here is how (simplified using uniform time step and ocCoef): OF CrankNiconson implementation Similarly, the . Thus, at any time step n, the temporal derivative involves all phi values from 0 to n. which is clearly a deviation from CrankNicolson discrete temporal discretization scheme. 

April 27, 2020, 10:37 

#4 
Member
Rodrigo
Join Date: Mar 2010
Posts: 91
Rep Power: 13 
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 CrankNicolson was accidental or even prosecuted. Either the way, clearing this may be in benefit for the OFcommunity. 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! 

April 27, 2020, 11:58 

#5 
Super Moderator

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:
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.
__________________
Keep foaming, Tobias Holzmann 

April 27, 2020, 13:41 

#6 
New Member
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 9 
Hello Tobi,
I have implemented the CrankNicolson 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/978331916 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 onedimensional wave equation. Your kind suggestion is much appreciated. Thanks again both of you! Regards, Rajib 

April 29, 2020, 23:24 
Comparison of temporal discretization schemes in OpenFOAM

#7 
New Member
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 9 
A 1D scalar transport equation is used to verify the temporal discretization schemes, where a passive scalar (phi) is convected along the Cartesian xaxis using a uniform flow velocity of (1, 0, 0):
The schemes used for temporal discretization are: Domain extends [0,1] in a Cartesian XY plane. A sinusoidal field is applied on the left (x = 0) at time 0, and convected toward the right (along x axis). After 1 characteristic time, scalar field values are plotted along line , 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 CrankNicolson (CNMod) scheme implemented by the author (as described in the first post) is oscillation free. However, the modified CrankNicolson 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. 

April 30, 2020, 12:13 

#8 
Senior Member
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 898
Rep Power: 9 
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.
__________________
The OpenFOAM community is the biggest contributor to OpenFOAM: User guide/Wiki1/Wiki2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESIOpenCFDTrademark) Bugs/Features: FOAMExtend (WikkiFSB) Bugs: OpenFOAM.org How to create a MWE New: v2106 (Community) 

May 1, 2020, 11:09 
Update on CrankNiconson scheme spatial discretization treatment

#9 
New Member
Rajib Roy
Join Date: Jun 2014
Location: Laramie, Wyoming
Posts: 18
Rep Power: 9 
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 CrankNicolson (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: and where is spatial discretization operator. If we introduce, a CN coefficient (offcenter coefficient > ocCoeff) from blending with implicit Euler scheme, the reformulated schemes will be: The correction I proposed in the earlier post, following the OpenFOAM implementation, is neither semiimplicit 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 semiimplicit 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. 

May 4, 2020, 20:00 
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: 9 
The CrankNicolson temporal discretization scheme can be viewed as summation of implicit and explicit Euler discretization schemes.
For any variable , let assume the linear system for the partial differential equation is written as: where is the spatial discretization operator. The temporal derivatives can discretized as : Implicit Euler Explicit Euler Crank Nicolson Adding equations {Eqn:implEuler} and {Eqn:explEuler}, we get With central difference approximation for , the final form of the discretization becomes: 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 timestep. Thus, equation {Eqn:governingEqn} at any time step can be approximated as, where is the final residual of the linear system at timestep . 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 Substituting, value in the CrankNicolson equation {Eqn:implEuler}, we get OpenFOAM applies this form of the CrankNicolson scheme with assumption that the linear system residual are driven to machine zero 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 CrankNicolson 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 CrankNicolson scheme in the format of equation{Eqn:CNOF} is unstable, but original format (equation {Eqn:CN} is stable. Last edited by rajibroy; May 5, 2020 at 14:35. 

May 5, 2020, 09:57 

#11 
Member
Rodrigo
Join Date: Mar 2010
Posts: 91
Rep Power: 13 

Tags 
cranknicolson, openfoam, temporal discretization 
Thread Tools  Search this Thread 
Display Modes  


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 