
[Sponsors] 
Suggested unsteady, implicit solver stable with arbitrarily large time steps 

LinkBack  Thread Tools  Display Modes 
February 22, 2010, 03:47 

#21 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
I think we are automatically subscribed to all the threads where we post a comment, so everybody should be here :D
P.S. Hi Niklas. I just noticed you got the super powers. Still using twoPhaseEuler?
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image. OpenQBMM  An opensource implementation of quadraturebased moment methods Last edited by alberto; February 22, 2010 at 03:48. Reason: Added PS 

February 22, 2010, 03:51 
ok, I see the link

#22 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
... yes, ok, and now I have seen the link to the new location in the 'bug forum' Niklas created. Sorry for this 'offtopic' post


February 22, 2010, 05:32 

#23  
New Member
Charles Mockett
Join Date: Dec 2009
Location: Berlin, Germany
Posts: 9
Rep Power: 7 
Dear all,
This thread remains very interesting reading! Just a quick remark, admittedly somewhat off topic, but underlining some of our motivation in this matter. Quote:
Best regards, Charlie. 

February 22, 2010, 05:55 

#24  
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 19 
Quote:
In fact, the only time I was involved in that piece of code was when I was helping a new PhD student and added the framework for the kinetic theory, so he had something to start from. Unfortunately he quit for a much better paid job and I decided to release the rather unfinished code to the public, hoping that someone would find it useful, since I knew I wouldn't have time to work with it. ...and the reason I move the topic is because I see this topic as a feature enhancement and code development. 

February 22, 2010, 06:14 
DetachedEddy simulation

#25 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Charlie, dear all,
that is indeed important. It now should be obvious to anyone that we both are aiming at calibrating our C_DESconstants in OpenFOAM for our hybrid RANS LES implementations. We use different background models for the DetachedEddy Simulations we perform and we would like to achieve what we have achieved and both published with ELAN in FLOMANIA or DESider also with OpenFOAM. For this procedure where we compute the decay of isotropic turbulence in LES mode my suggestion now is that we do not use any relaxation in OpenFOAM. I remember this also worked in ELAN  Charlie, am I right? Did we not experience that relaxation was not required for this case? I am not sure about the more complicated cases we are aiming at where above mentioned grid configurations will be applied and where the computations are already running for a long time. There we might have to underrelax our equations. A remedy would be to slightly increase the relaxation factors for momentum only during the next days for the running computations and see whether we reach URF_U=1 without divergence. Apart from that, we should go on looking into that issue and try to find out what we have to do to generally avoid this dependence on relaxation. One idea I had was to relax the U equation only for computing the intermediate velocities, but not manipulating the coefficients of the momentum equation that are used for calculating the pressure correction. An alternative is explicitly relaxing the intermediate velocities before the pressure correction is computed. In this case, the matrix coefficients are not manipulated. However, this did not help, but we have to check the details of what we did an check how the relaxation factors are applied to the matrix coefficients for momentum. The results for the decay should be dominated by modeled and numerical diffusion. Thus, with all settings equal, the relaxation factors should NOT have any influence. Best regards, Ulf. Last edited by djbungee; February 22, 2010 at 11:32. Reason: typos 

February 23, 2010, 19:45 
Solution to avoid dependence on relaxation

#26 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear all,
my initial idea not to apply relaxation on the coefficients of the momentum equation used for computing the pressure correction was correct, but not correctly implemented at first. Below please find a code adjustment suggested to avoid dependence of results on relaxation for the procedure listed above. The basic idea is to use tmp<fvVectorMatrix> MomPredEqn for computing the intermediate velocities with an arbitrary relaxation. With tmp<fvVectorMatrix> UEqn only the coefficients are computed that are used for pressure correction. The momentum equation is not solved for a second time. To avoid relaxation a line with UEqn().relax(1); is added before the coefficients for the pressure equation are assembled. The attached picture demonstrates that the initial dependence on relaxation ('ORIG') that has been shown using relaxation factors of 1 and 0.5 for all equations is now completely gone ('NORELAX_ON_PRESSCORR'). Suggestions, comments and additional tests are welcome! This should also be applied, e.g., in pimpleFoam! Best regards, Ulf. Code:
. . . while( (maxResidual > convergenceCriterion)  (pFinalize == 1) ) { eqnResidual = 1; MaxResidualPrev=maxResidual; maxResidual = 0; IIcount++; tmp<fvVectorMatrix> MomPredEqn ( fvm::ddt(U) + fvm::div(phi, U) + sgsModel>divDevBeff(U) ); if (pFinalize != 1) { MomPredEqn().relax(); } if (pFinalize == 1) { MomPredEqn().relax(1); } eqnResidual = solve ( MomPredEqn() == fvc::reconstruct(fvc::snGrad(p)*mesh.magSf()) ).initialResidual(); maxResidual = max(eqnResidual, maxResidual); tmp<fvVectorMatrix> UEqn ( fvm::ddt(U) + fvm::div(phi, U) + sgsModel>divDevBeff(U) ); UEqn().relax(1); p.boundaryField().updateCoeffs(); volScalarField rUA = 1.0/UEqn().A(); U = rUA*UEqn().H(); UEqn.clear(); phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); adjustPhi(phi, U, p); p.storePrevIter(); scalar StoppPressCorr =1; int CountPressCorr =0; scalar PrevPressEqnResidual=1.E20; while ( (StoppPressCorr > convergenceCriterion) && (CountPressCorr < 10) ) . . . Last edited by djbungee; February 23, 2010 at 19:58. Reason: added remark on pimple etc. 

February 23, 2010, 21:18 

#27  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hello all,
Quote:
Code:
volScalarField rUA = 1.0/UEqn().A(); U = rUA*UEqn().H(); UEqn.clear(); The relaxation is applied to the matrix as shown here (see comments in the code): http://foam.sourceforge.net/doc/Doxy...ce.html#l00483 Best, 

February 24, 2010, 02:23 

#28 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hello,
I quickly run some test in a cavity case, using what Ferziger and Peric did as example. I used pimpleFoam, with the same settings for the numerical schemes used in the pimpleFoam tutorial, but changing the time integration to backward. Results here: http://dl.dropbox.com/u/659842/OpenF...Relaxation.pdf I added the plots of the difference between the solutions along y. As you can see the difference is low, and then has a peak in a specific point for all the solutions. The pressure has the peak of difference at the centre of the cavity, while the velocity has the peak of difference around the minimum value of the pressure. I want to run the same cases again with stricter convergence requirements to be sure of the results. Of course any comment is welcome Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image. OpenQBMM  An opensource implementation of quadraturebased moment methods Last edited by alberto; February 24, 2010 at 02:57. Reason: Added comment 

February 24, 2010, 04:36 

#29  
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Alberto, dear all,
let me give some details. Quote:
Quote:
I know, the difference might be small and dependent on the time step size. But as was explained earlier, also by Charlie, our aim is to use an arbitrarily large time step that is a result of 'physical' consideration avoiding as many errors as possible. Completely dispensing with relaxation for the U Equation is a solution, but in many of the unsteady cases we run, this would never work. We have to relax the momentum predictor. Best regards, Ulf. Last edited by djbungee; February 24, 2010 at 05:16. Reason: directly addressed Alberto, details in last sentence 

February 24, 2010, 12:17 

#30 
New Member
Ana Cubero
Join Date: Dec 2009
Location: University of Zaragoza
Posts: 2
Rep Power: 0 
Dear all,
I had to cope with some of the issues described in this thread during my PhD, when I developed an inhouse code with a coupled and implicit solver. Maybe, my point of view (although I am not at all an expert on OpenFoam) could be in some way helpful. As already mentioned, the Rhie and Chow formulation (technique also known as Momentum Interpolation), expressed in its original form, produces converged solutions which depend on the linear relaxation factor and the (both pseudo or real) time step. Although the results accuracy is not significantly affected, this technique can even fail in its main objective of avoiding pressure wiggles on colocated grids. I believe that the pressure equation used in OF can be seen as a Momentum Interpolation (I mean, flux are calculated from the Rhie and Chow procedure and introduced on the continuity equation). Furthermore, a remedy for the case of small time steps is provided throughout the "ddtPhiCorr" function (please, correct me if I am wrong). However, no corrections for linear relaxation factors are available (maybe because of its little practical importance) I do not know if you could be interested in our paper about the momentum interpolation technique, where it is described in some detail, showing reasons for its inconsistency when unsteady problems. There, we propose a compact formulation for both unsteady problems and including linear relaxation. I could try to implement it on OF and check it. (The paper is: Ana Cubero and Norberto Fueyo, A compact momentum interpolation procedure for unsteady flows and relaxation, Numerical Heat Transfer, part B, 52, 6, 507529.) On the other hand, I also believe that it is good to include an inner loop bellow the real time step. However, differently from discussed above, I would be more interested in a dual time step, where you have a pseudotime iteration bellow the realtime loop (this is usually done as preconditioning of solvers for compressible flows at low Mach numbers). I do not have a clear reason for that, but I do not feel good sensations with linear relaxation factors... By using a the dual step procedure, we would have two temporal terms, the false time step being limited by convergence and the real time step by physical temporal scales. Thank you for reading until here, Ana 

February 24, 2010, 14:32 

#31 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
@djbungee: OK. Thanks.
You case (turbulence decay) is probably much more sensitive to this problem than what can be observed in other applications, I agree. It would be interesting, at least for me, to see a measure of the difference in the solution, and compare it with the error causes by other assumptions, to have a feeling of the relative order of magnitude. @ana: Thanks for the reference. If you plan to code your work in OpenFOAM, feel free to poke me. I also agree with you about the twolevel time stepping. Best, 

February 25, 2010, 04:48 

#32 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Ana, dear Alberto,
yes, dual time stepping is an alternative and a more flexible and general way, but it will not further help us here. The problem regarding dependence on relaxation caused by using 'relaxed' or 'nonrelaxed' coefficients has been solved and the coupling of equations on a time step works. So the 'local goal' of this small thread has been reached. What might help to overcome the weaknesses of our segregated approach in solving our equations is a blockmatrix solver. http://powerlab.fsb.hr/ped/kturbo/Op...ringer2009.pdf I expect that we would be able to achieve timedependent results of equal quality with larger time steps and higher stability with that. @Alberto, the measure of error for the decay introduced, e.g. by the convection scheme, is visibly much larger. I mentioned earlier: The results for the decay should be dominated by modeled and numerical diffusion. With these computations we want to calibrate our modeled diffusion, i.a. a constant C_DES for a DetachedEddy simulation for a fixed numerical scheme and procedure. This is actually only valid for the applied spatial resolution, or in other words, it is not a constant value. Therefore, we compute many different spatial resolutions. This way we get a good guess.... So for us it was important to avoid this error as there are enough other errors introduced that are not so easily corrected. Best regards, Ulf. 

February 25, 2010, 13:23 

#33  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
An additional reference which I found recently, where many improvements are suggested to the standard RhieChow interpolation. S. Zhang, X. Zhao, General formulations for the RhieChow interpolation, ASME Heat Transfer/Fluid Engineering Summer Conference, 2004. They developed a formulation of RhieChow independent from underrelaxation, time stepping and they give guidelines to deal with discontinuities, multiphase flows and porous zones. Best, 

February 25, 2010, 16:01 

#34 
Senior Member
Fabian Braennstroem
Join Date: Mar 2009
Posts: 407
Rep Power: 10 
Hello Ulf,
I have some trouble understanding your lines (eps. the pFinalize part) from yesterday. How do you set pFinalize and shouldn't it be vice versa, i.e. if (pFinalize == 1) { MomPredEqn().relax(); } if (pFinalize != 1) { MomPredEqn().relax(1); } Or did I misunderstand your explanation? Gruesse! Fabian 

February 26, 2010, 03:37 

#35 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Hi Fabian,
the implementation is meant to work like this, that for the loop before the _final_ loop, when convergence is 'detected', pFinalize is set to 1. After that the final loop starts as the entering condition for the dowhile loop allows for the inner loop on the time step to be entered for either pFinalize==1 or 'no convergence'. Therefore, when pFinalize==1, relaxation should be applied for neither the momentum predictor nor the coefficients for the pressure correction (the latter is the case for every loop, of course). Actually, this might not be necessary anymore after the possibility of relaxation has been removed from the determination of the coefficients for the pressure equation. But in case you do want to avoid relaxation in the final loop, you have to say (IF I UNDERSTAND THE CALL CORRECTLY ....) Code:
if (pFinalize =! 1) { MomPredEqn().relax(); } if (pFinalize == 1) { MomPredEqn().relax(1); } Hope that this answers your question. If not, we should have a telephone call on the weekend  I prefer to write my programs Saturday night, when all our children are in bed  have you got time? Have a great day ahead ... Grüße, Ulf. 

February 26, 2010, 06:23 

#36 
New Member
Ana Cubero
Join Date: Dec 2009
Location: University of Zaragoza
Posts: 2
Rep Power: 0 
Dear Ulf,
I am sorry, but I still think that solutions will depend on the (pseudo or real) time step used (and not only because of the different resolution). I mean, if we solve the evolution of a flow towards the steady state, the final, converged solution will depend on the time step used. If you look at the expression for the (predicted) flux (for an incompressible flow and using the firtorder Euler scheme): phi = interp(U) + ddtPhiCorr = interp(U) + interp(1/A/DeltaT)*(phi.oldTimeU.oldTime), where U = H/A. At convergence, we have phi=phi.oldTime and U=U.oldTime and then phi = interp(H/A). As the A term does include DeltaT ( A = A_steady + coef/DeltaT, right?), our predicted flux will depend on the time step. Regarding the block matrix solver, are you using structures provided by OF for building the matrix or are you creating it by yourselves? Best regards, Ana 

March 4, 2010, 05:57 

#37 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Ana,
this answer I write down here is based on assumptions only, but let me give you my ideas. I am not too deep in the details of the implementation, but I worked with and performed development work in other codes, so this is how I remember it. The flux you have mentioned is only used for the pressure correction equation. After convergence the flux in the pressurecorrection equation is then used as the 'global flux' (phi = pEqn.flux()) and the momentum is corrected by the adjusted pressure. At convergence of the _system_ of equations this correction tends to zero. There are many possibilities on how to define the coefficients for the pressure correction equation. As I remember it there are implementations that neglect terms with neighbor coefficients in order to make the solution of the pressure equation as easy as possible, e.g. solving explicit expressions instead of solving an equation system. However, the pressure correction found goes to zero for any procedure and you get conservation of momentum and mass. Therefore, these simplifications do not play an important role as might be the case for Dt in the center coefficient A. I might be wrong about this in general. However, we did not observe any dependence on time step size beyond the dependence on whether you resolve the physics or not. In these observations might be included what you described. Therefore, your considerations might be rather theoretical and I have the same question Alberto asked us: What is the relative error, can it be quantified? Actually, that is also the same question we were asked initially when starting this thread and the discussion about dependence on relaxation came up. We want to establish a solver with arbitrarily large, physically justified times steps. We have achieved that at least in a way that we can use larger time steps now and there is no dependence on relaxation factors left! However, when you have very small time steps by physical considerations and, e.g. when you do LES and want to resolve all time scale, our corrections to the procedure are not important. However, in this case, I would prefer a completely explicit procedure to save computation time anyway. Of course, the observations and examinations we made tells us nothing about the amount of dependence on the time step _apart_ from resolution as you mentioned. However, we want to compute additional test cases and describe our observations. Fabian has made a very interesting suggestion about simulating the decay of a TaylorGreen vortex where an analytical solution is available to identify the advantages of our implementation. Regarding the block matrix solver: we are not involved, just 'observing' what is going on. Therefore I added the link. Best regards, Ulf. Last edited by djbungee; March 8, 2010 at 06:42. Reason: removed ambigous formulation 

May 11, 2010, 10:47 

#38 
Member
Cedric Van Holsbeke
Join Date: Dec 2009
Location: Belgium
Posts: 81
Rep Power: 7 
Is there any chance that a fully implicit (PIMPLE?) solver will be available in the next release of OpenFOAM? I have a lot of cases where accuracy is not a priority, but speed is.
In Fluent I can use very large time steps (0.005 s) which give good results for my cases with lagrangian particle tracking and are super fast. In OpenFOAM, I need time steps < 1e6 s in order to be stable. 

May 11, 2010, 15:30 

#39  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
Quote:
Best, 

May 12, 2010, 02:51 

#40  
Member
Cedric Van Holsbeke
Join Date: Dec 2009
Location: Belgium
Posts: 81
Rep Power: 7 
Quote:
Sensitivity analysis in Fluent showed that a time step 0.005s is more than accurate enough for my purposes, but pimpleFoam does not seem to manage. That's why a implicit solver would be a nice addition. Last edited by CedricVH; May 12, 2010 at 03:28. 

Tags 
implicit, pimple, relaxation, simple, unsteady 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Implicit versus Explicit  Deepak  Main CFD Forum  16  November 29, 2014 21:54 
Working directory via command line  Luiz  CFX  4  March 6, 2011 21:02 
Doubt on Implicit Methods  analyse In India  Main CFD Forum  10  March 9, 2007 04:01 
Time scales and Fluent's unsteady solver issue  Freeman  FLUENT  6  December 13, 2005 15:30 
Convergence with coupled implicit solver  Henrik StrÃ¶m  FLUENT  1  October 29, 2005 03:57 