CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Fourth Order Runge Kutta time integration

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

Like Tree22Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   October 26, 2005, 09:07
Default Hello, Has someone implemen
  #1
Senior Member
 
Frank Bos
Join Date: Mar 2009
Location: The Netherlands
Posts: 340
Rep Power: 10
lr103476 is on a distinguished road
Hello,

Has someone implemented 4th order Runge Kutta time integration. I think this method will be more efficient than the 2nd order CrankNickolson.

Maybe someone tried it and could give me some hints in order to develop it myself.

Regards, Frank
kiddmax and songwukong like this.
__________________
Frank Bos
lr103476 is offline   Reply With Quote

Old   December 16, 2009, 09:19
Default
  #2
New Member
 
Join Date: Nov 2009
Posts: 17
Rep Power: 8
misakagan is on a distinguished road
I'm also wondering if someone has runge-kutta implementation for incompressible turbulent flows (LES or DNS). Apparently the only application dealing such flows is Pisofoam and it uses implicit time stepping, which is very slow.
misakagan is offline   Reply With Quote

Old   June 14, 2010, 12:07
Default
  #3
New Member
 
Michael B Martell Jr
Join Date: Feb 2010
Location: Amherst, MA
Posts: 18
Rep Power: 8
theory37 is on a distinguished road
I too am wondering about this. I am attempting to implement a 3rd order (low storage) RK scheme for an RSTM I am working on. Any ideas?
theory37 is offline   Reply With Quote

Old   September 19, 2010, 16:56
Default Runge-Kutta 4 to top-level OpenFOAM
  #4
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
Hi! I am currently working on compressible flows and writing a RK4 method to OpenFOAM. I am currently learning many things about top-level programming
(mostly from the codes written by other people) but would like to share a simple example of one way to program
RK4 in the top-level code. The following code (that can be put inside the main for-loop of any existing solver to test it) solves the simple advection equation for a variable rho (that can represent of course anything).
Here, we btw can assume for the moment being that
U is a constant field i.e. the velocity.

rhoOld = rho;
phiv = fvc::interpolate(U)& mesh.Sf();
k1 = -runTime.deltaT()*fvc::div(phiv, rhoOld);
k2 = -runTime.deltaT()*fvc::div(phiv, rhoOld + 0.5*k1);
k3 = -runTime.deltaT()*fvc::div(phiv, rhoOld + 0.5*k2);
k4 = -runTime.deltaT()*fvc::div(phiv, rhoOld + k3);

rho = rhoOld + a1*k1 + a2*k2 + a3*k3 + a4*k4;

// ai are the RK4 coefficients

Of the following I would like to hear some further comments about and hopefully the more experienced people could further comment on
these issues (or point out a proper link to a discussion).

When programming explicit code as above the correctBoundaryConditions() function should be used after each update because otherwise there might be
inconsistencies in BC's and also in the processor
BC's. I guess the reason for this is that field operations
such as the ones above have no influence on what
is happening on the boundary; right ?
One can also explicitly update the BC's for a certain
quantity (say e.g. rhoE that is often solved for in
compressible computations) by typing

rhoE.boundaryField() =
rho.boundaryField()*
(
e.boundaryField() + 0.5*magSqr(U.boundaryField())
);

Of course, it remains as user's responsibility
that everything stays consistent when doing
top-level OF solvers.

Regarding the previous question about an incompressible RK4 solver I do not see any problem
of why the above-presented approach for advection
equation would not work also for the
incompressible NS-equations .

Best regards,
Ville
Xinze, kiddmax, eRzBeNgEl and 6 others like this.
ville is offline   Reply With Quote

Old   September 20, 2010, 02:47
Default
  #5
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,904
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by ville View Post
When programming explicit code as above the correctBoundaryConditions() function should be used after each update because otherwise there might be
inconsistencies in BC's and also in the processor
BC's. I guess the reason for this is that field operations
such as the ones above have no influence on what
is happening on the boundary; right ?
Yes and no
Yes, you have to do something to update the boundaries, if you need that. No, what you have to do is not necessarily an explicit call to correctBoundaryConditions().

If you update the value of a field, and you also want to update the corresponding boundaryField, all you have to do is to replace

=

with

==

in the assignment. For example:

k1 == ...;

This should work in OF 1.6 and following. I am not sure about the previous versions, since I noticed this syntax for the first time in 1.6.

Of course, if you do not have an assignment, but a sum with +=, like in the velocity corrector step, you have to call correctBoundaryConditions().

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 20, 2010, 04:54
Default
  #6
Member
 
Juho Peltola
Join Date: Mar 2009
Location: Finland
Posts: 82
Rep Power: 9
juho is on a distinguished road
Quote:
Originally Posted by alberto View Post
Yes and no
This should work in OF 1.6 and following. I am not sure about the previous versions, since I noticed this syntax for the first time in 1.6.
I've used it in 1.4.1, so I would assume it also works in all the following versions.
juho is offline   Reply With Quote

Old   September 20, 2010, 10:26
Default
  #7
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,904
Rep Power: 27
alberto will become famous soon enoughalberto will become famous soon enough
OK. Thanks for the info
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 20, 2010, 10:46
Default
  #8
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
Hi and thanks for the comments!
As said, I am currently considering how to nicely implement the boundary conditions for a fully explicit, density based RK4 solver (say the hardest case of subsonic inflow, outflow for the time being).
Currently, in the prototype version, I define the BC's for p, T and U as they are rather convenient to give. The variables that are solved for are rho, rhoU and rhoE. Now, the BC's for rho, rhoU and rhoE would be needed. In e.g. subsonic inflow the BC for rho would need to be determined by the solution. Thus, p and
T may be used for determining the boundary value of rho.
After this the boundary fields of rhoU and rhoE may be constructed. Any ideas of how to conveniently do this?
How would the more experienced OF-people consider simply
updating the boundary field in the top-level code as is done in
e.g. rhoCentralFoam? Another option would be
defining a new BC type for rho, rhoU and rhoE that is constructed from p, T and U.
Best,
Ville
ville is offline   Reply With Quote

Old   October 31, 2012, 12:00
Default Runge-Kutta 4 density based LES solver implemented
  #9
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
Hi,
to get a closure: I have now implemented into OpenFOAM a RK4 based
fully explicit compressible solver. Works as smoothly as it only can
I've also written RK4 solvers for incompressible flows based on the projection method
which allows us to get rid of the PISO solvers if so desired.
Work based on the incompressible solver was published recently in Computers & Fluids
and can be found currently in the "Articles in Press" section of the journal.

Vuorinen V., Schlatter P., Boersma B., Larmi M., and Fuchs L., A Scale-Selective, Low-
Dissipative Discretization Scheme for the Navier-Stokes Equation, (to appear in Computers and Fluids)

Best,
Ville
ashvinc9 likes this.
ville is offline   Reply With Quote

Old   February 27, 2014, 03:29
Default Publication: Runge-Kutta 4 method for compressible and incompressible flows
  #10
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
Hi,
probably the first published paper on the topic including practical instructions
on how to implement, theory, numerical validation
On the implementation of low-dissipative Runge–Kutta projection methods for time dependent flows using OpenFOAM

Vuorinen et al.
http://www.sciencedirect.com/science...45793014000334

Best,
Ville
thg and songwukong like this.
ville is offline   Reply With Quote

Old   November 11, 2014, 09:27
Default Fluid dynamical part of the code shown herein
  #11
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
Large-eddy simulation in a complex hill terrain enabled by a compact fractional step OpenFOAM® solver


http://www.sciencedirect.com/science...65997814001513

Best wishes,
Ville
ville is offline   Reply With Quote

Old   August 19, 2015, 15:46
Default
  #12
Senior Member
 
Ehsan Asgari
Join Date: Apr 2010
Posts: 292
Rep Power: 9
syavash is on a distinguished road
Quote:
Originally Posted by ville View Post
Large-eddy simulation in a complex hill terrain enabled by a compact fractional step OpenFOAM® solver


http://www.sciencedirect.com/science...65997814001513

Best wishes,
Ville
Dear Ville,

Is it possible to share your incompressible solver?! I am sure many people like me seek for an explicit low-dissipation solver for LES-like simulation.

Syavash.
syavash is offline   Reply With Quote

Old   August 20, 2015, 03:27
Default
  #13
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
Hi,
the functional part of the code is given in the above link entirely inside the article.
You just need to copy that text and modify e.g. pisoFoam to get a working solver.
Note that the projection pressure units are a bit different in the rk4projectionFoam solver version than pisoFoam since we apply the projection method. This is just a matter of convention and the way the pressure is introduced to the system. In the end the units on
LHS and RHS of NS eqs are the same.
Best regards, Ville
ville is offline   Reply With Quote

Old   August 20, 2015, 08:24
Default
  #14
Senior Member
 
Ehsan Asgari
Join Date: Apr 2010
Posts: 292
Rep Power: 9
syavash is on a distinguished road
Quote:
Originally Posted by ville View Post
Hi,
the functional part of the code is given in the above link entirely inside the article.
You just need to copy that text and modify e.g. pisoFoam to get a working solver.
Note that the projection pressure units are a bit different in the rk4projectionFoam solver version than pisoFoam since we apply the projection method. This is just a matter of convention and the way the pressure is introduced to the system. In the end the units on
LHS and RHS of NS eqs are the same.
Best regards, Ville
Thanks Ville,

I have proceeded as the steps in your paper have suggested, but I have encountered some problems in creating the new solver:
1-The variables Uold, Uc, and dU are not defined, so I constructed them in createFields.H as volVectorField. Is it OK?!
2-I have renamed pRef in CreatePoissonMatrix.H to pRefValue because the latter was defined in pisoFoam
3-I have difficulty in defining dt. How should I define this variable? I tried : scalar dt, but OF throws me an error. I think I should consider a dimensionedScalar but do not know the right syntax.
4-Where should I define a1,a2,a3, and a4? I have currently defined them simply as scalar at the beginning of the while-loop.

At the moment the above issues come to my mind. I greatly appreciate if you help me compile the new solver.

Thanks,
Syavash
syavash is offline   Reply With Quote

Old   August 20, 2015, 08:32
Default
  #15
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
Hi,

>1-Uold and Uc variables are not defined, so I constructed them in createFields.H. Is it >OK?!

Of course. They are dummy fields which you can define with something like:

volVectorField Uold
(
IOobject
(
"Uold",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
U
);

volVectorField dU
(
IOobject
(
"dU",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
U
);


>2-I have renamed pRef in CreatePoissonMatrix.H to pRefValue because the latter >was defined in pisoFoam

Sure

>3-I have difficutly in defining dt. How should I define this variable? I tried : scalar dt, >but OF throws me an error. I think I should consider a dimensionedScalar but do not >know the right syntax.

You could replace it with runTime.deltaT() or define e.g. a dimensioned scalar
dt which you set to runTime.deltaT() at the beginning of each timestep. I just wrote
dt in the paper to make it more straightforward


>4-Where should I define a1,a2,a3, and a4? I have currently defined them simply as >scalar at the beginning of the while-loop.

For example you could define a file called rk4coeff.H which you "include" with
#include rk4coeff.H before main loop starts. There you could write something like

Info << "\nDefine RK4 coeff." <<endl;

const scalar a1 = 0.166666666667;
const scalar a2 = 0.333333333333;
const scalar a3 = 0.333333333333;
const scalar a4 = 0.166666666667;

Info << "\n a1 = " <<a1<< "\n a2 = " <<a2<<"\n a3 = " <<a3<<"\n a4 = " <<a4<< endl;

Got it working ?
ville is offline   Reply With Quote

Old   August 20, 2015, 08:48
Default
  #16
Senior Member
 
Ehsan Asgari
Join Date: Apr 2010
Posts: 292
Rep Power: 9
syavash is on a distinguished road
Quote:
Originally Posted by ville View Post
Hi,

>1-Uold and Uc variables are not defined, so I constructed them in createFields.H. Is it >OK?!

Of course. They are dummy fields which you can define with something like:

volVectorField Uold
(
IOobject
(
"Uold",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
U
);

volVectorField dU
(
IOobject
(
"dU",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
U
);


>2-I have renamed pRef in CreatePoissonMatrix.H to pRefValue because the latter >was defined in pisoFoam

Sure

>3-I have difficutly in defining dt. How should I define this variable? I tried : scalar dt, >but OF throws me an error. I think I should consider a dimensionedScalar but do not >know the right syntax.

You could replace it with runTime.deltaT() or define e.g. a dimensioned scalar
dt which you set to runTime.deltaT() at the beginning of each timestep. I just wrote
dt in the paper to make it more straightforward


>4-Where should I define a1,a2,a3, and a4? I have currently defined them simply as >scalar at the beginning of the while-loop.

For example you could define a file called rk4coeff.H which you "include" with
#include rk4coeff.H before main loop starts. There you could write something like

Info << "\nDefine RK4 coeff." <<endl;

const scalar a1 = 0.166666666667;
const scalar a2 = 0.333333333333;
const scalar a3 = 0.333333333333;
const scalar a4 = 0.166666666667;

Info << "\n a1 = " <<a1<< "\n a2 = " <<a2<<"\n a3 = " <<a3<<"\n a4 = " <<a4<< endl;

Got it working ?
Thanks for quick reply!!

Now I am getting an error like this:

Code:
syavash@syavash-VPCF11DGX:~/OpenFOAM/OpenFOAM-2.3.1/applications/solvers/incompressible/rk4projectionFoam$ wmake
options:2:66: warning: backslash and newline separated by space [enabled by default]
     -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \    
 ^
Making dependency list for source file rk4projectionFoam.C
SOURCE=rk4projectionFoam.C ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository -ftemplate-depth-100 -I/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/turbulenceModels/incompressible/turbulenceModel -I/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/../applications/solvers/incompressible/pisoFoam -I/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/transportModels -I/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/transportModels/incompressible/singlePhaseTransportModel -I/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude -IlnInclude -I. -I/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude -I/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/rk4projectionFoam.o
In file included from rk4projectionFoam.C:58:0:
/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/setDeltaT.H: In function ‘int main(int, char**)’:
/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/setDeltaT.H:36:35: error: ‘CoNum’ was not declared in this scope
     scalar maxDeltaTFact = maxCo/(CoNum + SMALL);
                                   ^
In file included from rk4projectionFoam.C:46:0:
/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/initContinuityErrs.H:37:8: warning: unused variable ‘cumulativeContErr’ [-Wunused-variable]
 scalar cumulativeContErr = 0;
        ^
make: *** [Make/linux64GccDPOpt/rk4projectionFoam.o] Error 1
Any idea where I migh be wrong??!

Edit: I could compile the code by adding #include "CourantNo.H" just after #include "readTimeControls.H".
But this warning still persists:

In file included from rk4projectionFoam.C:46:0:
/home/syavash/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/initContinuityErrs.H:37:8: warning: unused variable ‘cumulativeContErr’ [-Wunused-variable]
scalar cumulativeContErr = 0;

Another question: Can I adjust time step by giving courant number as in pimpleFoam?!
syavash is offline   Reply With Quote

Old   August 20, 2015, 09:06
Default
  #17
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
The turbulence model warning would be a matter of some normal include statements
that could be copied from pisoFoam. As you can see, you have now created an OpenFOAM code from scratch and this piece of code does not really assume too many
things: there are fields which are updated in time. Thus, the rk4projectionFoam solver
is simply a field update scheme with explicit time integration and finite volume discretization. About time step control: why could you not do it ? Of course one
needs to understand the algorithm: at which point of the main loop you update
it etc but otherwise you would have quite a freedom to do that.
syavash likes this.
ville is offline   Reply With Quote

Old   August 20, 2015, 09:13
Default
  #18
Senior Member
 
Ehsan Asgari
Join Date: Apr 2010
Posts: 292
Rep Power: 9
syavash is on a distinguished road
Quote:
Originally Posted by ville View Post
The turbulence model warning would be a matter of some normal include statements
that could be copied from pisoFoam. As you can see, you have now created an OpenFOAM code from scratch and this piece of code does not really assume too many
things: there are fields which are updated in time. Thus, the rk4projectionFoam solver
is simply a field update scheme with explicit time integration and finite volume discretization. About time step control: why could you not do it ? Of course one
needs to understand the algorithm: at which point of the main loop you update
it etc but otherwise you would have quite a freedom to do that.
Dear Ville,

Now I am really willing to compare runtime of pisoFoam and the new solver together.
Do you mind if I post my observations here?!

Thanks,
Syavash
syavash is offline   Reply With Quote

Old   August 20, 2015, 09:25
Default
  #19
Member
 
ville vuorinen
Join Date: Mar 2009
Posts: 67
Rep Power: 9
ville is on a distinguished road
Sure. Please bare in mind that the conclusions I've made on runtime differences
were mostly for turbulent flows in parallel runs. Full conclusions are probably
depending on the number of processors, the parallel system which you use, the linear solver, the case (e.g. laminar vs turbulent). Good to start with lid driven cavity and check if you can reproduce the Ghia's data.
ville is offline   Reply With Quote

Old   August 20, 2015, 10:14
Default
  #20
Senior Member
 
Ehsan Asgari
Join Date: Apr 2010
Posts: 292
Rep Power: 9
syavash is on a distinguished road
Quote:
Originally Posted by ville View Post
Sure. Please bare in mind that the conclusions I've made on runtime differences
were mostly for turbulent flows in parallel runs. Full conclusions are probably
depending on the number of processors, the parallel system which you use, the linear solver, the case (e.g. laminar vs turbulent). Good to start with lid driven cavity and check if you can reproduce the Ghia's data.
All right,

As something that migh matter, should any modifications be applied in controlDict, fvScheme, or fvSolution?!

Edit: I have encountered the following error during runtime,

Code:
--> FOAM FATAL IO ERROR: 
keyword div(U) is undefined in dictionary "/media/syavash/science/PHD_Thesis/New/system/fvSchemes.divSchemes"

file: /media/syavash/science/PHD_Thesis/New/system/fvSchemes.divSchemes from line 30 to line 36.

    From function dictionary::lookupEntry(const word&, bool, bool) const
    in file db/dictionary/dictionary.C at line 437.

FOAM exiting
I think modifying of fvScheme seems to be necessary. Could you tell me what actions should I make here?!

Thanks
syavash is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Runge Kutta 4th Order Source Code sugu Main CFD Forum 4 October 26, 2012 03:15
Runge-Kutta 4rd Order method help for 6DoF in CFD siw Main CFD Forum 0 August 29, 2008 06:08
runge kutta Shuo Main CFD Forum 0 January 7, 2008 20:29
4th and 5th Order TVD Runge-Kutta Methods saygin Main CFD Forum 2 January 30, 2006 12:45
Runge Kutta vs adams bashforth time marching vasanth Main CFD Forum 5 January 1, 2006 01:17


All times are GMT -4. The time now is 22:04.