Implement backward time scheme manually

Michael@UW · September 29, 2022, 13:43

Suppose the backward ddtSchme is used in fvSchemes file, then we can assemble the momentum equation as follows:

Code:

        fvVectorMatrix UEqn
        (
            fvm::ddt(U) +
            fvm::div(phi, U)
            + turbulence->divDevReff(U)
        ==
          - fvc::grad(p)
        );

I want to implement the backward time scheme manually in the code. According to the formula at https://www.openfoam.com/documentati...-backward.html, the ddt(U) can be replaced by the implicit and explicit sources terms as follows:

Code:

        fvVectorMatrix UEqn
        (
              fvm::div(phi, U)
            + turbulence->divDevReff(U)
        ==
            - fvc::grad(p)
            - fvm::Sp(1.5/dt, U)
           + (2*U.oldTime()  - U.oldTime().oldTime())/dt

        );

However, the second one gives a wrong solution. Can anyone tell me what's wrong with the implementation?

Regards,
Michael

überschwupper · September 30, 2022, 01:11

Quote:

Originally Posted by Michael@UW

Suppose the backward ddtSchme is used in fvSchemes file, then we can assemble the momentum equation as follows:

Code:

        fvVectorMatrix UEqn
        (
            fvm::ddt(U) +
            fvm::div(phi, U)
            + turbulence->divDevReff(U)
        ==
          - fvc::grad(p)
        );

I want to implement the backward time scheme manually in the code. According to the formula at https://www.openfoam.com/documentati...-backward.html, the ddt(U) can be replaced by the implicit and explicit sources terms as follows:
However, the second one gives a wrong solution. Can anyone tell me what's wrong with the implementation?

Regards,
Michael

Code:

        fvVectorMatrix UEqn
        (
              fvm::div(phi, U)
            + turbulence->divDevReff(U)
        ==
            - fvc::grad(p)
            - fvm::Sp(1.5/dt, U)
           + (2*U.oldTime()  - 0.5*U.oldTime(2))/dt

        );

I believe you forgot factor 0.5 for U^n-2. Further I replaced the function with another version of that function (should be equivalent according to the source code, but cleaner)

Michael@UW · September 30, 2022, 11:37

überschwupper, thank you for your reply! Yes, I missed the factor of 0.5. It's good to know the syntax of U.oldTime(2), but it seems to be not available in OF v2106. What version are you using?

I found there is still a difference between the two approaches.

Code:

        fvVectorMatrix UEqn
        (
           fvm::Sp(1.5/dt, U)            
        ==
           (2*U.oldTime() - 0.5*U.oldTime(2)/dt   // U.oldTime(2) not available though in OFv2106
        );

        fvVectorMatrix UEqnDdt
        (
          fvm::ddt(U)   // backward in controlDict
        );

        Info << UEqnDdt - UEqn << endl;
        Info << UEqnDdt.source() - UEqn.source() << endl;

The difference seems to be a constant which is (-ΔV/Δt), and the official OpenFOAM implementation gives a smaller diagonal entry. And the source difference is 0.

How does it come?

überschwupper · October 12, 2022, 05:16

sorry for the late reply, have you found a solution?

I checked the API Guide, unfortunately the ESI Version has not the overloaded function of oldTimes(). In OFv9 there is such a function.

- deleted some part, I will again think about it

Michael@UW · October 12, 2022, 10:51

I should have mentioned that I only tested for the first time step.

In the source code of ddt(backward) at

https://www.openfoam.com/documentati...ce.html#l00444

Code:

template<class Type> 
tmp<fvMatrix<Type>> 
backwardDdtScheme<Type>::fvmDdt 
( 
     const GeometricField<Type, fvPatchField, volMesh>& vf 
) 
{ 
     tmp<fvMatrix<Type>> tfvm 
     ( 
         new fvMatrix<Type> 
         ( 
             vf, 
             vf.dimensions()*dimVol/dimTime 
         ) 
     ); 
 
     fvMatrix<Type>& fvm = tfvm.ref(); 
 
     scalar rDeltaT = 1.0/deltaT_(); 
     scalar deltaT = deltaT_(); 
     scalar deltaT0 = deltaT0_(vf); 
     scalar coefft   = 1 + deltaT/(deltaT + deltaT0); 
     scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); 
     scalar coefft0  = coefft + coefft00; 
 
     fvm.diag() = (coefft*rDeltaT)*mesh().V(); 
 
     if (mesh().moving()) 
     { 
         fvm.source() = ... 
     } 
     else 
     { 
         fvm.source() = ...
     } 
 
     return tfvm; 
}

So, for $t=0, \Delta t = \Delta t_0$ ,

coefft = 1 + deltaT/(deltaT + deltaT0) = 1.5
coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0)) = 0.5
coefft0 = coefft + coefft0 = 2.0

fvm.diag() = (coefft*rDeltaT)*mesh().V() = $2.0\frac{\Delta V}{\Delta t}$

However, in the manually implemented version fvm::Sp(1.5/dt, U) will be set the diagonal entry to be
fvm.diag() = $1.5\frac{\Delta V}{\Delta t}$

Here is a difference $0.5\frac{\Delta V}{\Delta t}$ , which partially explains why they are different, though I found the difference is twice of this one.

Please let me know if you found some errors in my derivation or have any thoughts.

Michael@UW · October 12, 2022, 10:55

I wonder if it is wrong that OpenFOAM gives fvm.diag() = (coefft*rDeltaT)*mesh().V() = for the first time step when backward temporal scheme is used.

überschwupper · October 13, 2022, 04:40

Quote:

Originally Posted by Michael@UW

I should have mentioned that I only tested for the first time step.

In the source code of ddt(backward) at

https://www.openfoam.com/documentati...ce.html#l00444

Code:

template<class Type> 
tmp<fvMatrix<Type>> 
backwardDdtScheme<Type>::fvmDdt 
( 
     const GeometricField<Type, fvPatchField, volMesh>& vf 
) 
{ 
     tmp<fvMatrix<Type>> tfvm 
     ( 
         new fvMatrix<Type> 
         ( 
             vf, 
             vf.dimensions()*dimVol/dimTime 
         ) 
     ); 
 
     fvMatrix<Type>& fvm = tfvm.ref(); 
 
     scalar rDeltaT = 1.0/deltaT_(); 
     scalar deltaT = deltaT_(); 
     scalar deltaT0 = deltaT0_(vf); 
     scalar coefft   = 1 + deltaT/(deltaT + deltaT0); 
     scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); 
     scalar coefft0  = coefft + coefft00; 
 
     fvm.diag() = (coefft*rDeltaT)*mesh().V(); 
 
     if (mesh().moving()) 
     { 
         fvm.source() = ... 
     } 
     else 
     { 
         fvm.source() = ...
     } 
 
     return tfvm; 
}

So, for $t=0, \Delta t = \Delta t_0$ ,

coefft = 1 + deltaT/(deltaT + deltaT0) = 1.5
coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0)) = 0.5
coefft0 = coefft + coefft0 = 2.0

fvm.diag() = (coefft*rDeltaT)*mesh().V() = 1.5\frac{\Delta V}{\Delta t}

However, in the manually implemented version fvm::Sp(1.5/dt, U) will be set the diagonal entry to be
fvm.diag() = $1.5\frac{\Delta V}{\Delta t}$

Here is a difference $0.5\frac{\Delta V}{\Delta t}$ , which partially explains why they are different, though I found the difference is twice of this one.

Please let me know if you found some errors in my derivation or have any thoughts.

I was also looking inside the code, but had no time yesterday.

I have found an error. You made an mistake while setting in values inside the code equation. The factor is 1.5 not 2.0. It is correct at this point.
there has to be another reason, Maybe treatment of Sp()? I hope finding time to dig into it this evening

Michael@UW · October 13, 2022, 09:52

That's good catch! I did not notice that it is coefft not coefft0 is used to calculate the diagonal entries.
fvm.diag() = (coefft*rDeltaT)*mesh().V();

The questions remain. At the first time step, Backward scheme in OpenFOAM actually is $1.5\frac{dU}{dt}$ which is not expected.
Still need to figure out why there is the difference between OpenFOAM and the manual implementation.

Michael@UW · October 13, 2022, 10:20

I just noticed

Code:

scalar deltaT0 = deltaT0_(vf);

is 0 which was not what I thought $=\Delta t$ . This means

Code:

scalar coefft   = 1 + deltaT/(deltaT + deltaT0)

is zero, hence backward returns to Euler. Question 1 answered!

Code:

template<class Type>
template<class GeoField>
scalar backwardDdtScheme<Type>::deltaT0_(const GeoField& vf) const
{
    if (mesh().time().timeIndex() < 2)
    {
        return GREAT;
    }
    else
    {
        return deltaT0_();
    }
}

überschwupper · October 13, 2022, 10:51

Quote:

Originally Posted by Michael@UW

I just noticed

Code:

scalar deltaT0 = deltaT0_(vf);

is 0 which was not what I thought $=\Delta t$ . This means

Code:

scalar coefft   = 1 + deltaT/(deltaT + deltaT0)

is zero, hence backward returns to Euler. Question 1 answered!

Code:

template<class Type>
template<class GeoField>
scalar backwardDdtScheme<Type>::deltaT0_(const GeoField& vf) const
{
    if (mesh().time().timeIndex() < 2)
    {
        return GREAT;
    }
    else
    {
        return deltaT0_();
    }
}

Good job! i was thinking that this was a definition and not a function call.. But deltaT0 is GREAT and GREAT is 10^6, which would set coefft to 1, coeft00 becomes zero for time indizes below 2, which is logical.

The difference between both is then dependend on what ever oldTime() is storing in the first time steps. Since every time step afterwards is dependend on the previous ones, the error is accumulating.

Michael@UW · October 17, 2022, 17:03

I just confirmed the difference between the two implementations for the first time step:

backward in OpenFOAM turns to be Euler, so
ap1 = $\frac{\Delta V}{\Delta t}$
source1 = $\frac{\Delta V}{\Delta t} U_0$

manually backward implementation
ap2 = $1.5\frac{\Delta V}{\Delta t}$
source2 = $\frac{\Delta V}{\Delta t} U_0$

The following is confirmed.
ap2 - ap1 = $0.5 \frac{\Delta V}{\Delta t}$
source2 - source1 =

This means my explicit implementation of fvm::ddt(backward) is correct but remedy should be introduced for the first time step as OpenFoam does.

September 29, 2022, 13:43	Implement backward time scheme manually	#1
Michael@UW Senior Member Join Date: Jan 2019 Posts: 125 Blog Entries: 1 Rep Power: 0	Suppose the backward ddtSchme is used in fvSchemes file, then we can assemble the momentum equation as follows: Code: fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U) + turbulence->divDevReff(U) == - fvc::grad(p) ); I want to implement the backward time scheme manually in the code. According to the formula at https://www.openfoam.com/documentati...-backward.html, the ddt(U) can be replaced by the implicit and explicit sources terms as follows: Code: fvVectorMatrix UEqn ( fvm::div(phi, U) + turbulence->divDevReff(U) == - fvc::grad(p) *- fvm::Sp(1.5/dt, U) + (2U.oldTime() - U.oldTime().oldTime())/dt** ); However, the second one gives a wrong solution. Can anyone tell me what's wrong with the implementation? Regards, Michael

September 30, 2022, 11:37		#3
Michael@UW Senior Member Join Date: Jan 2019 Posts: 125 Blog Entries: 1 Rep Power: 0	überschwupper, thank you for your reply! Yes, I missed the factor of 0.5. It's good to know the syntax of U.oldTime(2), but it seems to be not available in OF v2106. What version are you using? I found there is still a difference between the two approaches. Code: fvVectorMatrix UEqn ( fvm::Sp(1.5/dt, U) == (2U.oldTime() - 0.5U.oldTime(2)/dt // U.oldTime(2) not available though in OFv2106 ); fvVectorMatrix UEqnDdt ( fvm::ddt(U) // backward in controlDict ); Info << UEqnDdt - UEqn << endl; Info << UEqnDdt.source() - UEqn.source() << endl; The difference seems to be a constant which is (-ΔV/Δt), and the official OpenFOAM implementation gives a smaller diagonal entry. And the source difference is 0. How does it come?

October 12, 2022, 05:16		#4
überschwupper Member Join Date: Jan 2022 Location: Germany Posts: 72 Rep Power: 4	sorry for the late reply, have you found a solution? I checked the API Guide, unfortunately the ESI Version has not the overloaded function of oldTimes(). In OFv9 there is such a function. - deleted some part, I will again think about it Last edited by überschwupper; October 12, 2022 at 05:25. Reason: deleted (havent read precisely)

October 12, 2022, 10:51		#5
Michael@UW Senior Member Join Date: Jan 2019 Posts: 125 Blog Entries: 1 Rep Power: 0	I should have mentioned that I only tested for the first time step. In the source code of ddt(backward) at https://www.openfoam.com/documentati...ce.html#l00444 Code: template<class Type> tmp<fvMatrix<Type>> backwardDdtScheme<Type>::fvmDdt ( const GeometricField<Type, fvPatchField, volMesh>& vf ) { tmp<fvMatrix<Type>> tfvm ( new fvMatrix<Type> ( vf, vf.dimensions()dimVol/dimTime ) ); fvMatrix<Type>& fvm = tfvm.ref(); scalar rDeltaT = 1.0/deltaT_(); scalar deltaT = deltaT_(); scalar deltaT0 = deltaT0_(vf); scalar coefft = 1 + deltaT/(deltaT + deltaT0); scalar coefft00 = deltaTdeltaT/(deltaT0(deltaT + deltaT0)); scalar coefft0 = coefft + coefft00; fvm.diag() = (coefftrDeltaT)mesh().V();* if (mesh().moving()) { fvm.source() = ... } else { fvm.source() = ... } return tfvm; } So, for $t=0, \Delta t = \Delta t_0$ , coefft = 1 + deltaT/(deltaT + deltaT0) = 1.5 coefft00 = deltaTdeltaT/(deltaT0(deltaT + deltaT0)) = 0.5 coefft0 = coefft + coefft0 = 2.0 *fvm.diag() = (coefftrDeltaT)mesh().V() = $2.0\frac{\Delta V}{\Delta t}$* However, in the manually implemented version fvm::Sp(1.5/dt, U) will be set the diagonal entry to be fvm.diag() = $1.5\frac{\Delta V}{\Delta t}$ Here is a difference $0.5\frac{\Delta V}{\Delta t}$ , which partially explains why they are different, though I found the difference is twice of this one. Please let me know if you found some errors in my derivation or have any thoughts.

October 13, 2022, 10:20		#9
Michael@UW Senior Member Join Date: Jan 2019 Posts: 125 Blog Entries: 1 Rep Power: 0	I just noticed Code: scalar deltaT0 = deltaT0_(vf); is 0 which was not what I thought $=\Delta t$ . This means Code: scalar coefft = 1 + deltaT/(deltaT + deltaT0) is zero, hence backward returns to Euler. Question 1 answered! Code: template<class Type> template<class GeoField> scalar backwardDdtScheme<Type>::deltaT0_(const GeoField& vf) const { if (mesh().time().timeIndex() < 2) { return GREAT; } else { return deltaT0_(); } }

October 12, 2022, 10:55		#6
Michael@UW Senior Member Join Date: Jan 2019 Posts: 125 Blog Entries: 1 Rep Power: 0	I wonder if it is wrong that OpenFOAM gives *fvm.diag() = (coefftrDeltaT)mesh().V() =* for the first time step when backward temporal scheme is used.

October 13, 2022, 09:52		#8
Michael@UW Senior Member Join Date: Jan 2019 Posts: 125 Blog Entries: 1 Rep Power: 0	That's good catch! I did not notice that it is coefft not coefft0 is used to calculate the diagonal entries. fvm.diag() = (coefftrDeltaT)mesh().V(); The questions remain. At the first time step, Backward scheme in OpenFOAM actually is $1.5\frac{dU}{dt}$ which is not expected. Still need to figure out why there is the difference between OpenFOAM and the manual implementation.

October 17, 2022, 17:03		#11
Michael@UW Senior Member Join Date: Jan 2019 Posts: 125 Blog Entries: 1 Rep Power: 0	I just confirmed the difference between the two implementations for the first time step: backward in OpenFOAM turns to be Euler, so ap1 = $\frac{\Delta V}{\Delta t}$ source1 = $\frac{\Delta V}{\Delta t} U_0$ manually backward implementation ap2 = $1.5\frac{\Delta V}{\Delta t}$ source2 = $\frac{\Delta V}{\Delta t} U_0$ The following is confirmed. ap2 - ap1 = $0.5 \frac{\Delta V}{\Delta t}$ source2 - source1 = This means my explicit implementation of fvm::ddt(backward) is correct but remedy should be introduced for the first time step as OpenFoam does.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
courant number increases to rather large values	6863523	OpenFOAM Running, Solving & CFD	22	July 5, 2023 23:48
laplacianFoam with source term	Herwig	OpenFOAM Running, Solving & CFD	17	November 19, 2019 13:47
Inconsistencies in reading .dat file during run time in new injection model	Scram_1	OpenFOAM	0	March 23, 2018 22:29
mixerVesselAMI2D's mass is not balancing	sharonyue	OpenFOAM Running, Solving & CFD	6	June 10, 2013 09:34
Upgraded from Karmic Koala 9.10 to Lucid Lynx10.04.3	bookie56	OpenFOAM Installation	8	August 13, 2011 04:03