View Issue Details
| ID | Project | Category | View Status | Date Submitted | Last Update |
|---|---|---|---|---|---|
| 0003490 | OpenFOAM | Bug | public | 2020-05-03 19:19 | 2020-05-05 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 |
| Summary | 0003490: Crank Nicolson ddt scheme implementation inaccurary | ||||
| Description | The 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 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 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 Information | 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 | ||||
| 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 time-step is implicit. |
|
|
> 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. |
|
|
> 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. |
|
|
>> 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) |
|
|
> 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? |
|
|
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 |
|
|
> 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. |
|
|
> 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? |
|
|
>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 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. |
|
|
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. |
|
|
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. |
|
|
> 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 |
|
|
The either the ESI documentation is incorrect or they have re-implemented CrankNicolson differently. |
|
|
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. |
|
|
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. |
|
|
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. |
|
|
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 Crank-Nicolson 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. |
|
|
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. |
| 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 | |
| 2020-05-03 20:17 | henry | Note Edited: 0011315 | |
| 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 | |
| 2020-05-03 22:15 | henry | Note Edited: 0011324 | |
| 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 |
