View Issue Details

IDProjectCategoryView StatusLast Update
0003490OpenFOAMBugpublic2020-05-05 08:07
ReporterrajibroyAssigned Tohenry 
PrioritynormalSeveritymajorReproducibilityalways
Status closedResolutionno change required 
PlatformLinuxOSUbuntuOS Version18.04
Product Version6 
Fixed in Version 
Summary0003490: Crank Nicolson ddt scheme implementation inaccurary
DescriptionThe Crank-Nicolson 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 cfd-online OpenFOAM forum post #1 (https://www.cfd-online.com/Forums/openfoam-programming-development/226125-crank-nicolson-scheme-implemented-wrong.html)

Furthermore, use of the previous temporal derivative instead of old-old time step has a potential flaw. Any temporal derivative at time n includes all time-steps from start-Time to nth time-step. Details are available at same cfd-online forum post #3.

An easy fix would to use cnCoef instead of coef and use two previous time-step fields instead of a temporal derivative.

Steps To ReproduceA 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 semi-implicit 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 InformationI 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
TagsNo tags attached.

Activities

henry

2020-05-03 19:32

manager   ~0011311

> To my understanding, the intended CN scheme is explicit type.

No the implementation is implicit, at least the part relating to the new time-step is implicit.

henry

2020-05-03 19:38

manager   ~0011312

> An easy fix would to use cnCoef instead of coef and use two previous time-step fields instead of a temporal derivative.

Doesn't that make the scheme bacward rather than Crank-Nicolson? Backward is more stable than Crank-Nicolson and does not require any off-centering.

rajibroy

2020-05-03 19:46

reporter   ~0011313

> Doesn't that make the scheme backward rather than Crank-Nicolson?

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.

rajibroy

2020-05-03 19:58

reporter   ~0011314

>> To my understanding, the intended CN scheme is explicit type.

> No the implementation is implicit, at least the part relating to the new time-step is implicit.

You are correct (the new time-step is implicit)!
I was referring that with new time-step 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 time-step and reast of the half at old time-step)

henry

2020-05-03 20:16

manager   ~0011315

Last edited: 2020-05-03 20:17

View 3 revisions

> 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 phase-fraction in the VoF solvers in which Crank Nicolson has a special implementation for the optional semi-implicit approach.

Could you provide a scalarTransportFoam test-case 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?

rajibroy

2020-05-03 20:39

reporter   ~0011316

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 Crank-Nicolson 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

henry

2020-05-03 21:01

manager   ~0011317

> If the previous temporal derivative is assumed as Euler implicit

This is only the case for the second time-step, not for subsequent time-steps.

henry

2020-05-03 21:04

manager   ~0011318

> 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 second-order? Have you found that the current Crank-Nicolson is not?

rajibroy

2020-05-03 21:15

reporter   ~0011319

>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.

henry

2020-05-03 21:15

manager   ~0011320

> ddt(\phi) = [c1*\phi +c2*\phi_0 + c3*\phi_{00} ]/dt

I don't think this is the correct form for Crank-Nicolson, it is the form for backward differencing. In Crank-Nicolson the time-derivative is equated to the average of the old and new time RHS.

henry

2020-05-03 21:17

manager   ~0011321

The results of the test-case demonstrate that Crank-Nicolson behaves similarly to backward which is what I would expect, it does not demonstrate that the implementation of either is incorrect.

rajibroy

2020-05-03 21:30

reporter   ~0011322

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/guide-schemes-time-crank-nicolson.html) has the same implementation, but mentioned different discretization coefficients, which also matches with book that I added in the additional reference.

henry

2020-05-03 21:37

manager   ~0011323

> Could it be that it's a different version of backward scheme instead of a CrankNicolson scheme!

No, they are both 2nd-order schemes and should behave similarly, particularly if a small amount of off-centering is used in Crank-Nicolson

henry

2020-05-03 21:39

manager   ~0011324

Last edited: 2020-05-03 22:15

View 3 revisions

The either the ESI documentation is incorrect or they have re-implemented CrankNicolson differently.

henry

2020-05-03 21:45

manager   ~0011325

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

rajibroy

2020-05-03 23:17

reporter   ~0011326

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)

henry

2020-05-03 23:41

manager   ~0011327

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

rajibroy

2020-05-03 23:58

reporter   ~0011328

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.

henry

2020-05-04 00:37

manager   ~0011329

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

rajibroy

2020-05-05 01:06

reporter   ~0011330

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)

henry

2020-05-05 08:07

manager   ~0011331

The case provided demonstrates that the current Crank-Nicolson implementation is correct and the proposed change is not correct and would reduce the order and accuracy of the scheme.

Issue History

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