Hello,
Has someone implemen
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 |
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.
|
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?
|
Runge-Kutta 4 to top-level OpenFOAM
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 |
Quote:
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, |
Quote:
|
OK. Thanks for the info :)
|
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 |
Runge-Kutta 4 density based LES solver implemented
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 |
Publication: Runge-Kutta 4 method for compressible and incompressible flows
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 |
Fluid dynamical part of the code shown herein
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 |
Quote:
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. |
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 |
Quote:
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 |
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 ? |
Quote:
Now I am getting an error like this: Code:
syavash@syavash-VPCF11DGX:~/OpenFOAM/OpenFOAM-2.3.1/applications/solvers/incompressible/rk4projectionFoam$ wmake 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?! |
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. |
Quote:
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 |
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. |
Quote:
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: Thanks |
Code:
--> FOAM FATAL IO ERROR: I think div(U) is missing from fvSchemes. I have added it like: Code:
div(U) Gauss linear; Code:
--> FOAM FATAL ERROR: Syavash |
Quote:
But I think a clarification on fvSchemes helps me a lot! Another point is that I can only see pressure residuals on the screen. How should I include velocity residuals either?! Syavash Edit: Can you post an example of fvSchemes & fvSolution files?! |
Quote:
Any comments?! |
Dear All,
Il read this topic and succefully do the implementation of rk4 time scheme. Is someone have some return to do about it ? I have the same error than you, but i dislike the new pressure dimension, may we have not well do some operation. I try the simulation of blasius boundary flow, and the time step down to 1.e-150 any idea. Thank you for answers. |
Quote:
Dear Ville, I have read your paper and your code. you have done a wonderful work in this and help me a lot. but i still have some questions, i cannot figure it out. 1>as we know, the rk4 is a fully explicit method, in openfoam, the explicit equation should use the fvc::ddt. but you use the fvm::ddt in code, such as solve(fvm::ddt(rho) + rhokcycle0 * RK4values[cycle]); (solve rho) why? and what's the meaning? 2> can we implement the time schemes directly in folder? in this way, maybe we can use the RK4 in all of the solvers thank you !!! Best, Frank |
The code lacks a deltaT parameter in the pressure step in "PressureCorrection.H"
Quote:
Obviously, it is a problem of time dimension, the reason is a lack of deltaT parameter on the construction of pressure correction, plz refer to Projection method, see the "Chorin's projection method" part. To help the one who is interested, I paste the corresponding code down here Code:
// PressureCorrection.H |
1 Attachment(s)
Quote:
Kind Regards, Syavash |
Test case
Thank you for sharing the code, Syavash. Could you kindly attach a minimum test case for the code, for example maybe a small case changing the OpenFOAM cavity. I believe this will help the other foamers.
Best regards, Yu CHENG |
2 Attachment(s)
Quote:
Very good idea. I have modified the solver to be able to take a constant source term into account. I have also provided a lid-driven cavity example with the new RK4Foam solver. I hope the other foamers would find it useful and hopefully do some experiments and comparisons. Kind Regards, Syavash |
Good work
Hey Syavash,
You've done a great work! I have compiled the code and run the cavity case, the preliminary results seem to be reasonable. I believe this work will definitely help other foamers. Thank you! Best wishes! Yu. |
Quote:
I have implemented this solver and I'm comparing it to pimpleFoam. It's compiling fine and also seems to be solving. However, I have two questions: (1) What is the justification for multiplying or dividing by runTime.deltaT()? I can see that we need this for dimensional purposes, but how will this change the solver? Will this term factorise out? (2) How will residuals be affected? I seem to have very high initial pressure residuals in the rk4projection method (compared to my pimpleFoam case). I'm using the same fvSolution specs for both cases, so I didn't expect this. Thanks, RH. |
Hey Syavash,
Thank you for the solver and the test case. I was able to compile the solver and run the test case with out a glitch. I tried comparing the solver running time between the RK4Foam and pisoFoam and found out that while the results are almost identical, the RK4Foam took more than 5 times the solver time of the pisoFOAM. The case I run was the cavity case with RE = 2500 and implementing RAS modeling. I was under the impression that RK4Foam results in saving of solver running time. Is that wrong or is there something I am doing wrong? Thanks |
Quote:
|
Quote:
The runtime depends on the number of processes as well. I have compared both pisoFoam and RK4 in LES simulations with cell number in order of 1e+8 and over 100 processes. I found them quite equally fast. However, it should be noted that in RK4 solver, the pressure equation (elliptic) is solved 4 times while in pisoFoam it is solved only 2 times in default! This is the main reason why RK4 could be slower than pisoFoam in some cases. In the end, I stick with pisoFoam as I found it more robust and more accurate in some canonical flow simulations! Best Regards, Syavash |
Quote:
The deltaT term is multiplied by because simply it is there in the original N-S equations! Have a look at some numerical textbook on Runge-Kutta algorithm. Best Regards, Syavash |
Quote:
|
Quote:
Best Regards, Syavash |
Quote:
|
Quote:
|
Quote:
Is it possible to modified the 'RK4' be an convenient 'option' in fvschemes file for ddtSchemes such as 'Euler' or 'backward'? then used for most compressible or incompressible solvers? |
All times are GMT -4. The time now is 10:14. |