CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Implement backward time scheme manually

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 29, 2022, 14:43
Post Implement backward time scheme manually
  #1
Senior Member
 
Join Date: Jan 2019
Posts: 112
Blog Entries: 1
Rep Power: 5
Michael@UW is on a distinguished road
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
Michael@UW is online now   Reply With Quote

Old   September 30, 2022, 02:11
Default
  #2
Member
 
Join Date: Jan 2022
Location: Germany
Posts: 51
Rep Power: 2
überschwupper is on a distinguished road
Quote:
Originally Posted by Michael@UW View Post
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)
überschwupper is offline   Reply With Quote

Old   September 30, 2022, 12:37
Default
  #3
Senior Member
 
Join Date: Jan 2019
Posts: 112
Blog Entries: 1
Rep Power: 5
Michael@UW is on a distinguished road
ü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?
Michael@UW is online now   Reply With Quote

Old   October 12, 2022, 06:16
Default
  #4
Member
 
Join Date: Jan 2022
Location: Germany
Posts: 51
Rep Power: 2
überschwupper is on a distinguished road
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)
überschwupper is offline   Reply With Quote

Old   October 12, 2022, 11:51
Default
  #5
Senior Member
 
Join Date: Jan 2019
Posts: 112
Blog Entries: 1
Rep Power: 5
Michael@UW is on a distinguished road
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 is online now   Reply With Quote

Old   October 12, 2022, 11:55
Default
  #6
Senior Member
 
Join Date: Jan 2019
Posts: 112
Blog Entries: 1
Rep Power: 5
Michael@UW is on a distinguished road
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.
Michael@UW is online now   Reply With Quote

Old   October 13, 2022, 05:40
Default
  #7
Member
 
Join Date: Jan 2022
Location: Germany
Posts: 51
Rep Power: 2
überschwupper is on a distinguished road
Quote:
Originally Posted by Michael@UW View Post
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
überschwupper is offline   Reply With Quote

Old   October 13, 2022, 10:52
Default
  #8
Senior Member
 
Join Date: Jan 2019
Posts: 112
Blog Entries: 1
Rep Power: 5
Michael@UW is on a distinguished road
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 is online now   Reply With Quote

Old   October 13, 2022, 11:20
Default
  #9
Senior Member
 
Join Date: Jan 2019
Posts: 112
Blog Entries: 1
Rep Power: 5
Michael@UW is on a distinguished road
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_();
    }
}
Michael@UW is online now   Reply With Quote

Old   October 13, 2022, 11:51
Default
  #10
Member
 
Join Date: Jan 2022
Location: Germany
Posts: 51
Rep Power: 2
überschwupper is on a distinguished road
Quote:
Originally Posted by Michael@UW View Post
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.
überschwupper is offline   Reply With Quote

Old   October 17, 2022, 18:03
Smile
  #11
Senior Member
 
Join Date: Jan 2019
Posts: 112
Blog Entries: 1
Rep Power: 5
Michael@UW is on a distinguished road
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 =0

This means my explicit implementation of fvm::ddt(backward) is correct but remedy should be introduced for the first time step as OpenFoam does.
Michael@UW is online now   Reply With Quote

Reply

Tags
backward time scheme

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 12:27.