# Implement backward time scheme manually

 Register Blogs Members List Search Today's Posts Mark Forums Read

 September 29, 2022, 14:43 Implement backward time scheme manually #1 Senior Member   Join Date: Jan 2019 Posts: 112 Blog Entries: 1 Rep Power: 5 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

September 30, 2022, 02:11
#2
Member

Join Date: Jan 2022
Location: Germany
Posts: 51
Rep Power: 2
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)
==
- 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)

 September 30, 2022, 12:37 #3 Senior Member   Join Date: Jan 2019 Posts: 112 Blog Entries: 1 Rep Power: 5 ü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?

 October 12, 2022, 06:16 #4 Member   Join Date: Jan 2022 Location: Germany Posts: 51 Rep Power: 2 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 06:25. Reason: deleted (havent read precisely)

 October 12, 2022, 11:51 #5 Senior Member   Join Date: Jan 2019 Posts: 112 Blog Entries: 1 Rep Power: 5 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 tmp> backwardDdtScheme::fvmDdt ( const GeometricField& vf ) { tmp> tfvm ( new fvMatrix ( vf, vf.dimensions()*dimVol/dimTime ) ); fvMatrix& 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 , 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() = However, in the manually implemented version fvm::Sp(1.5/dt, U) will be set the diagonal entry to be fvm.diag() = Here is a difference , 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 12, 2022, 11:55 #6 Senior Member   Join Date: Jan 2019 Posts: 112 Blog Entries: 1 Rep Power: 5 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.

October 13, 2022, 05:40
#7
Member

Join Date: Jan 2022
Location: Germany
Posts: 51
Rep Power: 2
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 tmp> backwardDdtScheme::fvmDdt ( const GeometricField& vf ) { tmp> tfvm ( new fvMatrix ( vf, vf.dimensions()*dimVol/dimTime ) ); fvMatrix& 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 , 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() = Here is a difference , 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

 October 13, 2022, 10:52 #8 Senior Member   Join Date: Jan 2019 Posts: 112 Blog Entries: 1 Rep Power: 5 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 which is not expected. Still need to figure out why there is the difference between OpenFOAM and the manual implementation.

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

October 13, 2022, 11:51
#10
Member

Join Date: Jan 2022
Location: Germany
Posts: 51
Rep Power: 2
Quote:
 Originally Posted by Michael@UW I just noticed Code: scalar deltaT0 = deltaT0_(vf); is 0 which was not what I thought . This means Code: scalar coefft = 1 + deltaT/(deltaT + deltaT0) is zero, hence backward returns to Euler. Question 1 answered! Code: template template scalar backwardDdtScheme::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.

 October 17, 2022, 18:03 #11 Senior Member   Join Date: Jan 2019 Posts: 112 Blog Entries: 1 Rep Power: 5 I just confirmed the difference between the two implementations for the first time step: backward in OpenFOAM turns to be Euler, so ap1 = source1 = manually backward implementation ap2 = source2 = The following is confirmed. ap2 - ap1 = 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.

 Tags backward time scheme