View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0003490  OpenFOAM  Bug  public  20200503 19:19  20200505 08:07 
Reporter  rajibroy  Assigned To  henry  
Priority  normal  Severity  major  Reproducibility  always 
Status  closed  Resolution  no change required  
Platform  Linux  OS  Ubuntu  OS Version  18.04 
Product Version  6  
Fixed in Version  
Summary  0003490: Crank Nicolson ddt scheme implementation inaccurary  
Description  The CrankNicolson ddt (CN) scheme implemented in OpenFOAM (v6 and v7) is not accurate! To my understanding, the intended CN scheme is explicit type. Without any blending, (i.e. CrankNicolson 1), the temporal derivative is very different than original CN scheme. Although the the blending coefficient cnCoeff = 1.0/(1.0 + ocCoeff) is introduced in code documentation, instead a different coefficient coef_ = 1 + ocCoeff; is used which led to a very different CN scheme. Please see the details outlined in cfdonline OpenFOAM forum post #1 (https://www.cfdonline.com/Forums/openfoamprogrammingdevelopment/226125cranknicolsonschemeimplementedwrong.html) Furthermore, use of the previous temporal derivative instead of oldold time step has a potential flaw. Any temporal derivative at time n includes all timesteps from startTime to nth timestep. Details are available at same cfdonline forum post #3. An easy fix would to use cnCoef instead of coef and use two previous timestep fields instead of a temporal derivative.  
Steps To Reproduce  A simple 1D scalar transport case made to test multiple temporal discretization schemes. The case study is available at GitHub https://github.com/rroyCFD/temporalDiscretizationVerification The modified explicit CrankNicolson scheme is available at https://github.com/rroyCFD/CrankNicolsonModDdtScheme (The Please rename the Function1Mesh class with OF default Function1 class; The Function1Mesh class had added features that I intend to submit as separate feature . the Function1Mesh class can be obtained from https://github.com/rroyCFD/Function1Mesh repository) Please note using CrankNicolson explicit scheme (where the spatial derivative are all explicit) may lead to instability. Introducing blending or treating spatial discrete operators as implicit helps to stabilize the solver. In addition, the semiimplicit CrankNicolson scheme has same discretization as implicit Euler and spatial terms are discretized as half implicit and rest of the half explicit. (see forum post #9)  
Additional Information  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  
Tags  No tags attached.  

> To my understanding, the intended CN scheme is explicit type. No the implementation is implicit, at least the part relating to the new timestep is implicit. 

> An easy fix would to use cnCoef instead of coef and use two previous timestep fields instead of a temporal derivative. Doesn't that make the scheme bacward rather than CrankNicolson? Backward is more stable than CrankNicolson and does not require any offcentering. 

> Doesn't that make the scheme backward rather than CrankNicolson? It doesn't! The coefficients are different. I wish I could write in Latex here! Please refer to forum post #1 where I discussed the temporal coefficients. 

>> To my understanding, the intended CN scheme is explicit type. > No the implementation is implicit, at least the part relating to the new timestep is implicit. You are correct (the new timestep is implicit)! I was referring that with new timestep implicit, the the linear system with CrankNicolson ddt scheme can be explicit (spatial terms are all explicit) or semi implicit (spatial terms are half evaluated at new timestep and reast of the half at old timestep) 

> CrankNicolson ddt scheme can be explicit The implementation in OpenFOAM is for implicit operation, it is not appropriate for explicit solution; all scalars in OpenFOAM are solved implicitly except the phasefraction in the VoF solvers in which Crank Nicolson has a special implementation for the optional semiimplicit approach. Could you provide a scalarTransportFoam testcase which demonstrates the problem with the current Crank Nicolson implementation? Also which coefficient(s) in the current implementation do you think is/are incorrect and what do you think they should be changed to? 

For any variable \phi, let's describe temporal derivative as ddt(\phi) = [c1*\phi +c2*\phi_0 + c3*\phi_{00} ]/dt or ddt(\phi) = [c1*\phi +c2*\phi_0 ]/dt + c4*ddt(\phi_0) The pure CrankNicolson the coefficients shall be (assuming uniform time step): c1 1/2 c2 0 c3 1/2 The current OpenFOAM implementation produces: c1 2 c2 2 c4 1 If the previous temporal derivative is assumed as Euler implicit, the coefficients become: \phi 2 \phi_0 3 \phi_{00} 1 A scalarTransportFoam test case can be download from the Github repository: https://github.com/rroyCFD/temporalDiscretizationVerification 

> If the previous temporal derivative is assumed as Euler implicit This is only the case for the second timestep, not for subsequent timesteps. 

> A scalarTransportFoam test case can be download from the Github repository: https://github.com/rroyCFD/temporalDiscretizationVerification How is the order verified in this case setup? Have you confirmed that backward is secondorder? Have you found that the current CrankNicolson is not? 

>How is the order verified in this case setup? The case doesn't verify the scheme order, only highlight difference between multiple schemes. I don't how to verify the order of a temporal scheme. Shall compare local error convergence order in a mesh refinement study? Please, advise. 

> ddt(\phi) = [c1*\phi +c2*\phi_0 + c3*\phi_{00} ]/dt I don't think this is the correct form for CrankNicolson, it is the form for backward differencing. In CrankNicolson the timederivative is equated to the average of the old and new time RHS. 

The results of the testcase demonstrate that CrankNicolson behaves similarly to backward which is what I would expect, it does not demonstrate that the implementation of either is incorrect. 

Yes, the implemented CN scheme behave like backward scheme! Could it be that it's a different version of backward scheme instead of a CrankNicolson scheme! I was wondering about the reference publication that was used to implement the scheme. For your kind reference, ESI OpenFOAM.com documentation (https://www.openfoam.com/documentation/guides/latest/doc/guideschemestimecranknicolson.html) has the same implementation, but mentioned different discretization coefficients, which also matches with book that I added in the additional reference. 

> Could it be that it's a different version of backward scheme instead of a CrankNicolson scheme! No, they are both 2ndorder schemes and should behave similarly, particularly if a small amount of offcentering is used in CrankNicolson 

The either the ESI documentation is incorrect or they have reimplemented CrankNicolson differently. 

If I run your case with a very small timestep all three schemes give the same results and more or less the same as backward and CrankNicolson with the original larger timestep, demonstrating that the current implementation of backward and CrankNicolson are correct. 

Please see the attached derivation for coefficients of the CrankNicolson scheme and spatial discretization operators. Note, the temporal discretization operators is similar to Euler and spatial discretization operators are semi explicit. To implement the scheme, I had modify the scalarTransportFoam solver UEqn.H files advection and diffusion terms. The scheme result in similar behavior as OpenFOAM CrankNicolson and backward scheme, with very different coefficients for discrete operators. Please advise on any inaccuracy in my derivation. OFbug_3490_rajibroy_CrankNicolsonDdtScheme.pdf (78,816 bytes) OFbug_3490_rajibroy_temporalSchemeComparison.png (314,976 bytes) 

I have found that for the case you provided the CrankNicolson and backward schemes provide results more accurate than Euler for the specified timestep, consistent with them being 2nd order schemes. The Euler scheme can generate the same results but only with a timestep at least an order of magnitude smaller as expected. From the results of your modified CrankNicolson scheme are more similar to those of Euler rather than CrankNicolson or backward indicating that it is a 1storder scheme. From these results I cannot see anything wrong with the current CrankNicolson implementation. 

Agreed that the CN scheme that I proposed earlier behaves like Euler. Please have a look at the 2nd CrankNicolson scheme (red line in the OFbug_3490_rajibroy_temporalSchemeComparison.png image). It behaves like 2nd order scheme. I am not saying OpenFOAM CN scheme is incorrect! I am not able to match the OpenFOAM CN scheme coefficients with the literature. Hence, my 2nd CN scheme following the literature is also 2nd order accurate and very different from OpenFOAM implementation. 

It is no longer clear what you are proposing. Can you provide any case which demonstrates a problem with the current CrankNicolson implementation? 

OpenFOAM CrankNicolson implementation is correct, when at each time step the linear system residual is zero. I have attached my derivation for your kind review. My hypothesis on inaccurate CrankNicolson ddt scheme implementation was wrong. Please consider this issue as resolved. CrankNicolsonDdtScheme_deriverdBy_RajibRoy.pdf (94,386 bytes) 

The case provided demonstrates that the current CrankNicolson implementation is correct and the proposed change is not correct and would reduce the order and accuracy of the scheme. 
Date Modified  Username  Field  Change 

20200503 19:19  rajibroy  New Issue  
20200503 19:32  henry  Note Added: 0011311  
20200503 19:38  henry  Note Added: 0011312  
20200503 19:46  rajibroy  Note Added: 0011313  
20200503 19:58  rajibroy  Note Added: 0011314  
20200503 20:16  henry  Note Added: 0011315  
20200503 20:16  henry  Note Edited: 0011315  View Revisions 
20200503 20:17  henry  Note Edited: 0011315  View Revisions 
20200503 20:39  rajibroy  Note Added: 0011316  
20200503 21:01  henry  Note Added: 0011317  
20200503 21:04  henry  Note Added: 0011318  
20200503 21:15  rajibroy  Note Added: 0011319  
20200503 21:15  henry  Note Added: 0011320  
20200503 21:17  henry  Note Added: 0011321  
20200503 21:30  rajibroy  Note Added: 0011322  
20200503 21:37  henry  Note Added: 0011323  
20200503 21:39  henry  Note Added: 0011324  
20200503 21:45  henry  Note Added: 0011325  
20200503 21:45  henry  Note Edited: 0011324  View Revisions 
20200503 22:15  henry  Note Edited: 0011324  View Revisions 
20200503 23:17  rajibroy  File Added: OFbug_3490_rajibroy_CrankNicolsonDdtScheme.pdf  
20200503 23:17  rajibroy  File Added: OFbug_3490_rajibroy_temporalSchemeComparison.png  
20200503 23:17  rajibroy  Note Added: 0011326  
20200503 23:41  henry  Note Added: 0011327  
20200503 23:58  rajibroy  Note Added: 0011328  
20200504 00:37  henry  Note Added: 0011329  
20200505 01:06  rajibroy  File Added: CrankNicolsonDdtScheme_deriverdBy_RajibRoy.pdf  
20200505 01:06  rajibroy  Note Added: 0011330  
20200505 08:07  henry  Assigned To  => henry 
20200505 08:07  henry  Status  new => closed 
20200505 08:07  henry  Resolution  open => no change required 
20200505 08:07  henry  Note Added: 0011331 