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

LinkBack  Thread Tools  Display Modes 
February 11, 2010, 10:10 
Suggested unsteady, implicit solver stable with arbitrarily large time steps

#1 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear all,
as already discussed in the forum "Questions & Answers related to OpenFOAM®", in the article thread "pimpleFoam, transientSimpleFoam, simpleFoam, turbFoam, oodles, transientSimpleOodles", in the XING OpenFOAMgroup, the CFLnumber in most of the transient solvers provided with the OpenFOAM library is not required, and to limit SIMPLE to steady simulations is also wrong. Moreover, in these solvers, you do actually not get convergence on an (arbitrarily) large time step and, therefore, experience convergence problems for larger time steps. What I consider important in order to achieve an implicit procedure with arbitrary timestep sizes is that all discretization is done in an implicit manner (fvm::....) and that the sequentially solved equations are solved for _several_ times in a single time step to reach a converged _SYSTEM_ of equations on a time step. Below and attached you will find an early development version of our code with fixed convergence criteria and some features that still need some improvements. However, the basic idea should become apparent. With this code you can achieve large time steps, and if the initial residuals in the inner loop below the timestep loop do not decrease, the timestep is probably physically too 'coarse' or too large and should be decreased. Suggestions and discussions are welcome. We want to discuss this on the Stammtisch in Berlin tomorrow, for anyone interested in the bofsession 'Hybrid RANS/LES and DES'. Best regards, Ulf. Code:
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright held by original author \\/ M anipulation   License This file is part of OpenFOAM. OpenFOAM is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 021101301 USA Application transientSimpleOodles Description Incompressible LES solver using transient SIMPLE. \**/ #include "fvCFD.H" #include "singlePhaseTransportModel.H" #include "incompressible/transportModel/transportModel.H" #include "LESModel.H" #include "IFstream.H" #include "OFstream.H" #include "Random.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H" #include "createTime.H" #include "createMesh.H"//#include "createMeshNoClear.H" #include "createFields.H" #include "initContinuityErrs.H" #include "readTransportProperties.H" Info<< "\nStarting time loop\n" << endl; for (runTime++; !runTime.end(); runTime++) { Info<< "Time = " << runTime.timeName() << nl << endl; #include "readPISOControls.H" #include "CourantNo.H" scalar eqnResidual = 1; scalar convergenceCriterion = 1.E20; scalar maxResidual = convergenceCriterion + 1; scalar MaxResidualChange = 0; scalar MaxResidualPrev = 0; int IIcount=0; // Pressurevelocity SIMPLE corrector //for (int corr=0; corr<nCorr; corr++) while(maxResidual > convergenceCriterion) { // Momentum predictor eqnResidual = 1; MaxResidualPrev=maxResidual; maxResidual = 0; IIcount++; tmp<fvVectorMatrix> UEqn ( fvm::ddt(U) + fvm::div(phi, U) + sgsModel>divDevBeff(U) ); UEqn().relax(); eqnResidual = solve ( UEqn() == fvc::grad(p) ).initialResidual(); // Retain maximum residual for global convergence check maxResidual = max(eqnResidual, maxResidual); 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); // Nonorthogonal pressure corrector loop if (nCorr != 1) //Ulf aus pimple { p.storePrevIter(); } ////// scalar StoppPressCorr =1; int CountPressCorr =0; scalar PrevPressEqnResidual=1.E20; //for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) while ( (StoppPressCorr > convergenceCriterion) && (CountPressCorr < 10) ) { CountPressCorr++; fvScalarMatrix pEqn ( fvm::laplacian(rUA, p) == fvc::div(phi) ); pEqn.setReference(pRefCell, pRefValue); eqnResidual = pEqn.solve().initialResidual(); StoppPressCorr = eqnResidual; // Stop loop if residual increases, i.e. set loop counter to 10 if (eqnResidual > PrevPressEqnResidual) CountPressCorr=10; PrevPressEqnResidual = eqnResidual; if (StoppPressCorr <= convergenceCriterion) { phi = pEqn.flux(); maxResidual = max(eqnResidual, maxResidual); } if (CountPressCorr == 10) { // Do not save maxResidual so that the outer, Ubased maxResidual is not influenced if convergence is not reached phi = pEqn.flux(); } } //// # include "continuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); // Ulf, nach unten geschoben: DAS SOLLTE HINTER MOMENTUM CORRECTOR: sgsModel>correct(); // Momentum corrector U = rUA*fvc::grad(p); U.correctBoundaryConditions(); sgsModel>correct(); // Info << " > maxResidual(U,p) " << maxResidual << endl; // CONVERGENCE CHECK MaxResidualChange=fabs(MaxResidualPrev  maxResidual); if (IIcount > 100 ) MaxResidualChange = MaxResidualChange * 100./IIcount; // If residual does not change anymore, stop loop if ( (MaxResidualChange/maxResidual < 0.001) && (IIcount!=1) ) { Info << " > maxResidual(U,p) " << maxResidual << endl; Info << " > MaxResidualChange (U,p) " <<MaxResidualChange <<endl; Info << " > MaxResidualChange/maxResidual " <<MaxResidualChange/maxResidual <<endl; Info << " > NOT reached convergence criterion " << convergenceCriterion << " on time step " << runTime.timeName() << " after " << IIcount <<" cycles " <<endl; Info << " > maxResidual (U,p) set to 0 to end loop due to low MaxResidualChange" <<endl; maxResidual=0; } // FINALLY REACHED CONVERGENCE if (maxResidual < convergenceCriterion) { Info << " > reached convergence criterion " << convergenceCriterion << " on time step " << runTime.timeName() << " after " << IIcount <<" cycles " <<endl; } } runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return(0); } // ************************************************************************* // Last edited by djbungee; February 15, 2010 at 11:10. Reason: adjusted sourcecode formatting for enhanced readability 

February 16, 2010, 05:50 

#2 
Senior Member
Matthias Voß
Join Date: Mar 2009
Location: Berlin, Germany
Posts: 446
Rep Power: 11 
thanks for posting these additional (but VERY necessary) snippet. i will try to set them up in the baratropicCavitatingFoam and report back.
neewbie P.S.: nice workshop in Berlin 

February 16, 2010, 09:12 
addon: Suggested unsteady, implicit solver stable with arbitrarily large time steps

#3 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear all,
as we were not able to present the energyspectra (Engergy vs. wavenumber) for the computations of decay of isotropic turbulence / box turbulence we made on the OpenFOAMStammtisch on Friday for different settings they are given here. The Fouriertransformation binary did not work on our laptop and we did not get a fortran compiler in a few minutes due to a lack of online access that day. Instead, we presented curves for the temporal evolution of point probes for instantaneous velocity at a given point. Significant differences became apparent. These results are confirmed here and now. Please take a look at the picture attached to this post. What do we see in the attached diagram? Energy vs. wave number at a given point of time / after simulating 2 time units. The basics of the results, how the Fourier transformation was done, and the comparability to available experimental data is of no interest for the following conclusions. The results itself are not very good, admittedly, as a very coarse grid was used. The simulations are done with the source code attached to the first post in this thread. We applied different underrelaxation (no underrelaxation: RELAX=1, and 50% underrelaxation: RELAX=0.5). Then we computed one loop per time step ('one loop') and a variable number of loops on a time step until a certain residual was reached ('more loops, MaxRes=1e20). And we varied the timestep size using DT=0.2 and DT=0.01. Thus, we get 8 simulations. What do we conclude?  For 'large' time steps with relaxation (red curves), the results significantly depend on whether you couple the equations by more loops on a time step or not.  For 'small' time steps (blue curves), this is also true.  Without relaxation (green and magenta) the influence of not coupling the equations decreases with decreasing timestep size. The magentacolored curves are even congruent  With coupling of equations on a timestep the results get 'improved' (are closer to each other), see the dashed lines. > Coupling on a time step (inner loops, controlled by residuals) is required. Another effect, but a totally different issue, is the dependence of the results on relaxation. Does this relate to the RhieandChow interpolation? Has anyone tried to evaluate that? Comments are welcome, but, however, please not on the curves themselves. We know that we need much finder grids for decent results especially when applying full central differencing. Anyway, the effects are the same as described here. We have checked that in other cases. Best regards, Ulf. 

February 16, 2010, 17:27 

#4  
New Member
Charles Mockett
Join Date: Dec 2009
Location: Berlin, Germany
Posts: 9
Rep Power: 7 
Dear Ulf,
Thanks for posting your results and conclusions in detail, especially since I couldn't make the workshop and find the topic very important. First of all, I agree with the procedure you have followed for the new transient SIMPLE solver (multiple iterations of the equation system within each time step). To save us having to compare the source code, could you possibly outline the differences between your solver and the transientSimpleFoam and pimpleFoam solvers? Your investigation with isotropic turbulence is very interesting. I find the results partly in line with my intuition, and partly perplexing: Quote:
I wouldn't necessarily expect this. My feeling is that if the time step is fine enough, explicit and implicit time integration should become equivalent... but I'm really not sure about this. Quote:
Quote:
Quote:
Quote:
So, being sufficiently rattled, I decided to repeat your experiment with ELAN (our inhouse code, also a colocated, FV, "transient SIMPLE" type). My setup was identical except for two things:
spectraTimeStepRelax.png So, what are my conclusions?
One thing's for sure, this URF dependency is highly undesirable and it needs to be fixed. Has anybody seen this with one of the standard solvers? All the best, Charlie. 

February 17, 2010, 06:55 

#5 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Charlie, dear all,
so we do or did have the same expectations, perfect. The dependence of results on under relaxation is nothing new, I guess. Angelo mentioned that in Berlin and searching for 'under relaxation Rhie and Chow' e.g. in google gives many links to related publications. There the effect is described and workarounds are given. However, I am not that deep in the OpenFOAM source code. As you mentioned, we could do that in ELAN, where my roots are  this is now your task . However, there you did not observe this effect  lucky you! Regarding your request to give details about the source code changes I will come back to that later. Basically it is only a dowhile loop on the time step that checks for the initial residuals, and, of course, that is very similar as it is done in ELAN. I remember when I implemented the grid deformation algorithm in ELAN together with the ArbitraryEulerLagrange () formulation for all equations that many adjustments were necessary regarding residuals and convergence check  but that is too long ago for me to compare exactly . Best regards, Ulf. Remark: For this simulation we used OpenFOAM, version 1.6, July 2009 ('unpatched'), employing our modified, LLkepsilonbased turbulence model for hybrid RANS/LES Last edited by djbungee; February 17, 2010 at 08:10. Reason: added remark on code version 

February 17, 2010, 12:58 

#6 
New Member
Charles Mockett
Join Date: Dec 2009
Location: Berlin, Germany
Posts: 9
Rep Power: 7 
Hi Ulf,
Thanks for the reply. Nice that we agree on how the solver should behave (of course!). The problem that the results depend on the underrelaxation factor should definitely be pursued further, as it is a serious limitation. Let's stay on the case and invite contributions / advice from the OpenFOAM community! Best, Charlie. 

February 17, 2010, 13:28 

#7  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Hello,
I have some comments, questions and a suggestion Quote:
Quote:
Quote:
Quote:
Quote:
My question is: where the solutions in the different cases actually converged? If so, what criterion was adopted to define the convergence? The suggestion is on the implementation of a different kind of timestep adaption procedure, which has been already successfully used in MFIX (www.mfix.org), with a SIMPLElike procedure. The approach they use is very effective at speeding up the simulation, and the idea is the following:
This algorithm was shown to be very successful in multiphase simulations with very stiff problems, like in granular flows, and can be quite easily implemented changing time class and with minor changes to the solvers Best, 

February 17, 2010, 13:59 

#8  
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Alberto,
thanks for your confirmation! As I see: we succeeded, as we wanted to show something obvious and something that one would expect . A success like that is a decade ago when teaching CFD  in engineering life this is different and all the engineers in the teams in my department and myself never get results we expect, either in crash, statics, dynamics, solid or fluid mechanics or whatever . Quote:
Quote:
The residual in for all equations in the simulations was down to a very low level (except for the simulations with a fixed number of inner iterations, of course), therefore the systems of equations can definitively considered as converged! The number of iterations on each time step naturally increases with increased under relaxation, but the level of convergence was equal for all the simulations we did (we have only shown a small excerpt). Nonetheless, the results differed dependent on the relaxation factor, although we did have equally converged solutions. This should not depend on the way you implement your pressurebased procedure. We rather think that this comes from the way pressurevelocity decoupling on nonstaggered grids is avoided by employing a RhieandChow interpolation. Best regards, Ulf. Last edited by djbungee; February 17, 2010 at 14:05. Reason: improved smileys and remarks 

February 17, 2010, 16:59 

#9  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Hi,
sorry, I didn't take a deep look to the code before... my bad! The forum makes it quite hard to read with the scrolling box My question about convergence was on this point. Thanks for the clarification. Quote:
Quote:
Quote:
Quote:
To stabilize the interpolation in those cases, a modified RhieChow interpolation is used, including the effect of the bodyforce term, pressure gradient and other source terms in the interpolation. You can do something similar in OpenFOAM (and they do that for some compressible solver like buoyantPisoFoam). You might want to give it a try to check if it plays a role also in your case. In your code you do not have any body force, so you simply need to modify the momentum predictor as follows UEqn() == fvc::reconstruct(fvc::snGrad(p)*mesh.magSf()) Some additional change is made to the correction step too, since also the velocity correction is reconstructed from the correction to the flux (see pEqn.H in buoyantPisoFoam) There is an additional term coming from the unsteady term that brings a dependency on the relaxation factors (see for example Chio, Numerical Heat Transfer, Part A, 36:545 ± 550, 1999), which, to my knowledge is ignored in OpenFOAM (and in other codes). In the paper they show that the error that you make when you neglect the term dependent on the relaxation factors in the interpolation is important when dt = 0.001 and URF_u = 0.7. I hope this helps a bit, and I'm very interested in the topic (Very happy you pointed it out and shared your results with us. Thanks!), since I've been fighting with a related interpolation problem for a work on multiphase flows I'm doing 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 17, 2010 at 17:04. Reason: Added explanation 

February 17, 2010, 18:37 

#10 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Maybe you know it already, but some interesting considerations are in
W. Z. Shen, J. A. Michelsen, N. N. Sorensen, J. N. Sorensen, An Improved SIMPLEC method on collocated grids for steady and unsteady flow computations, Numerical Heat Transfer, Part B, 43: 221239, 2003 In summary


February 18, 2010, 08:08 

#11 
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Alberto,
thank you very much for all the information and suggestions! Wow, I am impressed and see that you are very deep in these issues. We will try what you suggested and have to handle all the information you have given. Please give us some time. We will either come back with some more questions or show the improvements we have achieved. Again, thank you very much and best regards, Ulf. 

February 18, 2010, 13:26 

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

February 19, 2010, 04:48 

#13 
New Member
Charles Mockett
Join Date: Dec 2009
Location: Berlin, Germany
Posts: 9
Rep Power: 7 
Also from me thanks to Alberto for this very comprehensive advice  this will surely help us to overcome this problem!
Best regards, Charlie. 

February 19, 2010, 08:55 

#14 
New Member
Angelo Carnarius
Join Date: Feb 2010
Location: Berlin, Germany
Posts: 3
Rep Power: 7 
Dear all,
this is a very interesting topic as I was really surprised to see, that OpenFOAM seems to have problems with such fundamental numerical schemes. Concerning the discussions above, I have some comments: 1. In agreement to Ulf, I am am very sure that all the pressure correction schemes we are talking about (SIMPLE, SIMPLER, PISO) can be applied to steady as well as unsteady flows. There is no reason to restrict them to only one type of time integration. 2. Independent of how the pressure is calculated, there are two reasons why we have to iterate in an inner loop over the system of equations on each time step: a: We solve the coupled equations in a segregated way, which means that the coefficients of one equation are depending on the other equations. b: There is a nonlinear convective term in the momentum equations. This has to be linearised, which requires an iterative procedure even for time explicit schemes. If we do not apply the inner loop, we get an error, which somehow scales with the time step. I am not sure, whether this error tends to zero with decreasing time step. But this might be an academic problem, which is not of major interest, as long as the error is small for small time steps. However, in the above problem (decay of isotropic turbulence) , a Courant number (CFL) of 0.2, which is quite small, still gives a significant error. Note, that this behaviour is very different from explicit time integration, where CFL<1 is required for numerical stability, even if CFL>>1 would be sufficient for numerical accuracy!!! One should also keep in mind that the main advantage of time implicite schemes is that you don't have to pay attention to the Courant number. You simply choose the time step which is appropriate for the physical representation of your problem. For many industrial applications this leads to a very high Courant number. Thus you definitly need the inner loops. 3. Relaxation for the momentum equation means that we take only a part of the new solution of u_i. This is then used in all other equations. Relaxation for the momentum equations can ALWAYS be applied independent of the pressure correction scheme. 4. Relaxation for the SIMPLE algorithm means that the pressure of the last iteration is only corrected by a part of the pressure correction p'. For the correction of the velocities u_i, the FULL pressure correction p' HAS to be applied to achieve mass conservation. For PISO or SIMPLER, you never use relaxation for the pressure, as these schemes are constructed such that you can use the full pressure correction. It is very important not to mix up the above points in order to separate the influence of all the numerical parameters. From my point of view, inner iteration loops are essential for general applicability, no matter which pressure correction scheme is used. This allows the use of much larger time steps and the solution should be independent of the relaxation factor. Some comments on the Rhie&Chow Interpolation: From my knowledge, there are two flaws of the standard formulation. First, the calculation of the face velocities depends on the relaxation factor of the momentum equation. This holds for steady as well as unsteady problems. There exist different cures for this, simple ones as well as more complicated ones. Second, the damping propertie of the Rhie&Chow Interpolation decreases with small time steps. Modifications to avoid this are quite complex. I will post some results of test calculations on this topics today or on monday for further discussion. Best regards, Angelo 

February 19, 2010, 09:31 
... not a problem of OpenFOAM

#15  
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Angelo, dear all,
just one comment on your statement: Quote:
Thanks to all of you for this very fruitful discussion. Best regards, Ulf. 

February 19, 2010, 13:53 

#16  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Quote:
Of course, not having the dependency would be better. However I don't think "moving" the dependency from the URF to an arbitrary parameter, as done in one of the paper I read is any better, especially if this makes the method griddependent. Quote:
Quote:
For RANS, I agree with you, but the PIMPLE implementation is OK for that. If you are familiar with commercial codes, there is a famous one that iterates for a given number of times, exiting if residuals fall below a certain value, but it does not do absolutely anything if that criterion is not respected, which is not very different from what OpenFOAM does in PIMPLE loops, excluding the fixed number of iterations. Quote:
What test case are you using? Best, Alberto 

February 21, 2010, 20:01 

#17  
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear all,
what you mentioned, Alberto, is right: Quote:
Implementing reconstruct for the pressure gradient in the momentum equation (fvc::reconstruct(fvc::snGrad(p)*mesh.magSf())) as suggested earlier in this thread and alternatively or additionally when updating the momentum after pressure iteration did not change the results at all. Some additional test and implementations yielded the following:
Actually, I am not satisfied with these results and suspect that something is overlooked. Now I am really interested, what you all find out. Has anyone checked this for pimpleFoam? Best regards, Ulf. 

February 22, 2010, 01:51 

#18  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Hello and thank you for the update
Quote:
Quote:
Quote:
What do you think? Best, 

February 22, 2010, 03:35 

#19 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 20 
Moved the thread to this forum instead.
It seems more appropriate. 

February 22, 2010, 03:46 

#20  
Member
Ulf Bunge
Join Date: Mar 2009
Location: Wolfsburg, Germany
Posts: 34
Rep Power: 8 
Dear Niklas,
Quote:
Hopefully, other contributers will find the thread again easily  if not, that would be most inappropriate . Therefore, taking the history of a thread it does not really matter where it resides, I think. It is more important that contributers find it again. For forums or groups that I moderate I suggest and ask before I move a post. Of course, if a thread is totally or obviously off topic or violating rules it is moved or deleted. Well, never mind. Thanks and best regards, Ulf. 

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  17  November 7, 2015 14:14 
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 