Crank Nicolson scheme implemented wrong?

rajibroy · April 17, 2020, 22:20

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}$

rajibroy · April 27, 2020, 08:24

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.

rajibroy · April 27, 2020, 08:43

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.

guin · April 27, 2020, 10:37

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 · April 27, 2020, 11:58

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 · April 27, 2020, 13:41

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 · April 29, 2020, 23:24

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

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

HPE · April 30, 2020, 12:13

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 · May 1, 2020, 11:09

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.

rajibroy · May 4, 2020, 20:00

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

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

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

guin · May 5, 2020, 09:57

Here the closed report for future reference:

https://bugs.openfoam.org/view.php?id=3490

April 17, 2020, 22:20	Crank Nicolson scheme implemented wrong?	#1
rajibroy New Member Rajib Roy Join Date: Jun 2014 Location: Laramie, Wyoming Posts: 18 Rep Power: 11	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)

April 27, 2020, 08:24	Verification: temporal discretization coefficients	#2
rajibroy New Member Rajib Roy Join Date: Jun 2014 Location: Laramie, Wyoming Posts: 18 Rep Power: 11	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.

April 27, 2020, 08:43	use of previous temporal derivative is potentially inaccurate	#3
rajibroy New Member Rajib Roy Join Date: Jun 2014 Location: Laramie, Wyoming Posts: 18 Rep Power: 11	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.

April 27, 2020, 10:37		#4
guin Member Rodrigo Join Date: Mar 2010 Posts: 98 Rep Power: 16	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.

April 27, 2020, 11:58		#5
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51	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

April 27, 2020, 13:41		#6
rajibroy New Member Rajib Roy Join Date: Jun 2014 Location: Laramie, Wyoming Posts: 18 Rep Power: 11	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

April 30, 2020, 12:13		#8
HPE Senior Member Herpes Free Engineer Join Date: Sep 2019 Location: The Home Under The Ground with the Lost Boys Posts: 932 Rep Power: 12	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. __________________ The OpenFOAM community is the biggest contributor to OpenFOAM: User guide/Wiki-1/Wiki-2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESI-OpenCFD-Trademark) Bugs/Features: FOAM-Extend (Wikki-FSB) Bugs: OpenFOAM.org How to create a MWE New: Forkable OpenFOAM mirror

May 5, 2020, 09:57		#11
guin Member Rodrigo Join Date: Mar 2010 Posts: 98 Rep Power: 16	Here the closed report for future reference: https://bugs.openfoam.org/view.php?id=3490 HPE likes this.

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