What is the physics behind the unconditional stability of implicit schemes?..and though theory says that stability is unconditional,in practice it is not always really unconditional.the maximum permissible time-step value of an implicit scheme is some multiple of the corresponding explicit scheme value.Why is this so?

In order to solve the non-linear PDE's, which govern fluid flow, its normal procedure to use an iterative procedure of which 3 main types exist.

Implicit, Explicit & Semi-Implicit

Pure explicit schemes are conditionally stable, however Semi Implicit schemes are off higher order accuracy than their counterparts and are used as the norm for non-linear coupled equations. In this manner, in order for a scheme to behave like the exact solution, it must behave in a similar manner to the function that it is approximating, and become a better approximation in the next iterate.

Explicit Schemes, allows the calculation of a flow field at an appropriate step from calculations, obtained from previous steps. Since the discretization schemes normally involve some sort of truncation error, the error will accumulate over each step. Hence explicit schemes are known to be stable for a limited step size and due to the accumulation of the error over each time step, normally over estimate the final solution. However in Implicit Schemes the flow field calculations depend on the current time step only. As you can imagine Semi Implicit schemes are obtained by the averaging of implicit and explicit schemes together and in general such schemes are second order accurate and unconditionally stable like implicit schemes Implicit numerical techniques that allow arbitrarily large time-step sizes to be used in calculations are a popular way to reduce CPU time requirements. Unfortunately, these methods are not accurate for convective processes (i.e. Euler portion of the RANS Equations). Implicit methods gain their time-step independence by introducing diffusive effects into the approximating equations. In general for flows of medium Re, the convective time << than the diffusive time and implicit schemes compensate for such by the addition of numerical diffusion. However the addition of such diffusion, normally destroys the phenomena been modelled, from such viable solutions can only be obtained for conservative time steps. (i.e. in the order of the diffusive time)

Hi Radhakrishnan,

I can give you a lead on this. Stability can be looked at from the point of view of dependence of the solution at a future time on the space time domain in the present and past. If the numerical domain of dependence does not include the physical DoD then instability and also maybe ill-posedness results (the precise mathematical constraints are beyond me). It is easy to draw this DoD for different explicit numerical schemes for say some illustrative 1D problems where the solutions are known.

In the case of an implicit method the numerical DoD extends to the whole domain implicitly. This means that the physical DoD is always included in the numerical DoD. This explanation probably does not have the "physical" feel that I think you wanted but this is directly related to determinism and well-posedness. Also, what I say is totally not rigourous from a mathematical point of view so take this just as a pointer.

regards, chidu...

First consider an explicit scheme of the form (for simplicity - a one-dimensional wave equation like system )

du/dt=c*du/dx (I)

where c is the analogy to the sound speed. As one discretizes the equation one gets dx and dt as real finite quantities (delta-x and delta-t). Now the ratio delta-x/delta-t is the speed of the 'numerical' integration. IN order to be able to follow this one dimensional wave equation, which has a speed c, one needs to have that the numerical speed be larger than the physical speed c.

delta-x/delta-t larger than c, this gives the famous Courant Friedrich Levy (CFL) condition that delta-t be smaller than delta-x/c .

The same can be obtained for a second order equation in space of the form:

du/dt = nu*d2u/dx2 (II)

(where nu is for example the coefficient of the viscosity) and one obtains delta-t smaller than nu/delta-x2.

IN order to overcome these stability criteria one needs to look at the discretization scheme. For example for equation (II)

(u(n+1,j)-u(n,j))/delta-t = nu*(u(n,j+1)-u(n,j))/delta-x2, (III)

where index n refers to the time discretization and index j refers to the space discretization. From (III), which is of course an explicit scheme, one can write:

u(n+1,j)=u(n,j)+delta-t*nu*(u(n,j+1)-u(n,j))/delta-x2

when delta-t increases, so does u(n+1,j), and in the limit of very large (going to infinity) delta-t, u(n+1,j) also becomes very large, and the scheme does not converge. In order to remediate to this one uses implicity scheme as follow.

Consider equation (II) with an implicit Crank-Nicholson scheme:

(u(n+1,j)-u(n,j))/delta-t = nu/2*(u(n+1,j+1)-u(n+1,j))/delta-x2 + nu/2*(u(n,j+1)-u(n,j))/delta-x2

this can be re-aranged and leads to

u(n+1,j)=u(n,j) * (1-nu*delta-t/2/delta-x2)/(1+nu*delta-t/2/delta-x2) + ...

Now when delta-t becomes extremely large, this quantity on the right hand side does not diverge any more, but converge. This is what makes the scheme stable.

THe usual way to numerical integrate the implicit equation involves a lot of arithmetic and therefore errors. In addition, one usually uses the inversion of a matrix, which can be made only under certain conditions, these conditions translate into a small time step (but not extrememly small). THis is why, usualy, you still cannot take an infinitely large time step when solving an implicit scheme.

I hope this helps,

Patrick

As Patrick pointed out, implicit schemes are converged using a lot of iterations. These iterations usually involve solving a set of linear system(s) of equations that are either discrete forms of the original linear equations or linearized forms if the original equations are non-linear.

Let the linear system be Ax = b

The convergence depends on the asymptotic condition number (i.e. the spectrum) of A. A term of the order of 1/dt appears in the diagonal of A. If "dt" is small, this term if large and hence we have diagonal dominance. As "dt" is increased, the matrix becomes less and less diagonally dominant and usually (perhaps not always, you may have to ask a linear algebra expert) convergence slows down and eventually iterations may diverge. This is an explanation for why implicit schemes cannot handle large time steps even though theory tells you otherwise. The problem is with our ability in converging the iterations.

There may sometimes be problems with the theories as well. Most stability theories are linear theories which can not predict long term non-linear stability. In stabilities theories (usually), a small disturbance is introduced and if it dies down, the system is said to be stable. Local linearization of the equations is done to predict the growth/decay of the introduced error. The problem with large time steps is that local linearization does not make sense since the solution can change a great deal in just one time step. We can not track the small introduced error using linearized equations.

(1). If you evaluate the spatial terms at the known old time, and express the first order time derivative term by a two point difference between the future time and the old time, then it is explicit method, because the future value of the dependent variable can be computed directly (explicitly) from the known old values. In other words, the future values (and the flow field) depend solely on the old known values (and the old flow field). So, you can march in time this way, easily. But if you take a large time step, you can easily step into a sink hole or man hole. This is exactly the reason why the map comes in different sizes. And that is exactly the reason why it is hard to hit the Mars from the earth by the Mars lander. The gradient in time in the explicit method is evaluated on the right hand side of the equation by the old time spatial variation. This is a real model of the real physics, except that in real world, the time step is continuous and infinitestimal. So, the finite time step in explicit method will deviate from the true path and produce inaccurate answers. And if you keep the steering wheel fixed at certain position when entering the highway ramp, the car will diverge and it will becomes upside down. (no need to buy the Ford SUV)(2). In the implicit method, your are pretending that the spatial derivative term lives in the future time, which is not a right theory either. Because the future flow field is completely unknown, you are in a bigger trouble with this theory. The reason why you have a better chance to get the flow field is that, there exists "boundary conditions" in the future time. The solution theory becomes: how do we guess at the future flow field such that the future spatial derivative will equal to the two point time derivative, bounded by the future B.C's. You adjust (iterate) the future flow field until such condition is satisfied (backward requirement, implicit method, or iterative). Since you are free to adjust the future flow variable (to absorb or re-distribute the errors), you have a better chance to obtain a good solution.(in explicit method, you are not allowed to adjust anything at all) The problem is if the time step is large, and the future boundary condition has changed a lot, then you will be solving a different problem, linking the future and the old time directly by a two point difference, separated by a long time period. Something like traveling into the 22nd century by direct linking. This will fail, because the boundary conditions in the future is no longer related to the old time. If the large time step future boundary condition is such that a train will arrive at the cross road, then there is no need to rush there in one large time step.(you will run into the train.) The implicit method depends on the future boundary conditions so that the flow field variables can be re-adjusted. But if the time step is too large, the link between the future and the old time will fail through a simple two-point link. (3). What I am saying is: with small time step, you can drive a car down the road and then take the train using the explicit method. with large time step, you can hit the train, when driving the car down the road to the train station, using the implicit method. With large time step, using explicit method, you can end up on the road side of the highway rump, at the exit to the train station. with proper time step, using the implicit method, you can not only arrive at the train station safely, but also have extra 15 min to have a cup of coffee. (4). Well, Well, the choice is yours. You don't have to read this messy message at all. (investing in Internet stocks is something like using the implicit method....what you see is not always what you get)

I have reviewed all your messages with great interest. I am doing my PhD on implicit methods for the Navier-Stokes equations under Hypersonic regimes, so maybe I can comment about my experience.

The implict methods are not restricted by the CFL condition or similar, so in theory they can handle infinite time steps. Obviously the implicit methods are only useful looking for steady state solutions.

Problems found --------- Physics At large dt, the boundary condition influence,through its charactersitics, with other boundaries causing problems similar to the explicit one cell problems). Mainly in far from convergent states. This fact would limit(is not a statement) the use of initially large time steps Numerics Depending on the way of perform the approximation of the numerical fluxes in the next time level. I am currently using Taylor's expansion in a fully implicit way. At dt large, there would be a bigger ammount of numerical diffusion. This amount of diffusion can mask natural diffussion due to inherent viscosty. It has been found that the implicit method used works better in the Euler equations than in the N-S ones. In the final system Ax=b, A is less diagonal and takes more time to converge. I used BiCGSTAB to solve the system and I realized that using a large dt, increses the number of iterations of the solver up to 3 orders of magnitude. However while arriving to the steady state, the initial residual b is quite small and the number of iterations also increases due probably to numerical truncation errors.

Apart form all this incovinients, marching gradually for an initial time step with CFL=1 I have been able to arrive to use time steps 10 million times bigger than it explicit counterparts(CFL=10^7), and the program converge, giving the same solution that those obtained with CFL=10 (although is not a good one),. The save in CPU time is huge about several orders of magnitude. I am currently using slope limiters to limit the diffusion and more evolved temporal schemes to prevent advers effects.

Also in my experinece is a key point to treat implicitly the boundaries to guarantee this high time steps.

I do not know if all this stuff will hep you, only to say that stable implicit method can be worked out with huge time steps and the appropiate iterative solver, nevertheless probably it wont converge to the correct solution.

Thanks

I would say the implicit time step restriction is more linked to the non-linearity of the governing equations that to how the BCs are dealt with (implicit versus explicit). I noticed that for BCs that change very little in time (fixed temperature wall condition for example), the maximum allowable time step size is still present and is about of the same size as for BCs that change a lot in time (symmetry condition). I noticed also a very high dependance on the Mach number on the implicit time step restriction. I solved a flat plate at M=0.3 and then at M=25 using a fully implicit method and at M=25, the local time step was limited to a CFL constant of about 1. The reason for this I believe comes from the total energy at M=25 to be almost fully mechanical (more than 99.9%) where the nonlinearity of the terms is great. At M=0.3, the total energy is mostly thermal, hence resulting in equations that are more ``linear'' hence favouring an implicit scheme..

(1). For Mach 25 flow over a flat plate, you have the wall as the boundary condition. So, you said the local time step is limited to CFL = around one. (2). If you remove the wall and change it to the free stream condition, do you still have the same problem(limit)? I think, the limit is coming from the boundary conditions. What do you think? The governing equations are still the same non-linear equations, with or without the wall boundary condition.

I have tested the Isothermal Flat Plate running with Nitrogen at Ma=6.6. I am able to reach CFL up to 10 million, increasing the time step at a ratio of 1.2. And the results agreed well with the theoretical results. Reynolds number was around 10^6/m. I do not see a definitive influence(there is an influence is clear) of Mach number on the convergence of implicit methods, similar results obtained with low Mach numbers. The Mach number did not bound the time step in any case, grid size seem to have bigger influence, not limiting the time step but the ratio of increase.

Hello to all

I read all responses with considerable interest and learnt a few things and also got confirmations of a few things I have only suspected. As for the physical reason sought by R, one heuristic reasoning seems to have been missed out by most contributors although I am sure all of them know it: For a simple linear scalar equation it can be seen from modifed eq.analysis that time step restriction for explicit schemes comes from a need to keep the dissipation coefficient positive, surely a physical as well as mathematical requirement for stability of the solution. The same analysis reveals that the dissipation coefficient is positive and independent of the CFL number for a single step fully implicit scheme. This argument can be extended to the typical heat equation too.

However the argument is indeed heuristic and not a proof of stability. For nonlinear equations, the cell CFL number is solution dependent, but if one has bounds on the solution, a similar argument will go through.

For systems of nonlinear equations, obviously the situation is much more complicated and it is not possible to push the above heuristic analysis through. Strictly speaking, even if the algorithm gurantees bounds on the solution, the same may not exist for the oscillations and any higher order algorithm generates oscillations in the presence of convective terms. If these oscillations take the solution beyond the physical domain of validity of a variable, the algorithm breaksdown even if it is formally stable.

As some contributors have emphasized, conventional stability analysis of the Von Neumann type assumes periodicity of the BCs. Stability with BC has to resort to GKS theory and normal mode analysis. That is why, implicit treatment of BC in an otherwise implicit scheme is important for full stability. Also the DOD argument fails for the simple centered differencing first order explicit scheme for the advection equation because it is not a sufficient condition for stability. One can verify that the numerical viscosity coefficient in this case is unconditionally negative.

Lastly it is difficult to push physical analogies beyond a certain limit to settle questions like stability and convergence. We must however concede that they provide insight and perhaps can even guide hard mathematical analysis.

Ravi

This is definitely interesting. What was your initial CFL number and initial conditions? Did you launch this sim at a CFL of 1 million, or at a CFL of about 1? I can increase the CFL to very high values as well, but only after the solution has developped quite a bit, and is close to the final answer, and which point you might as well stop running anyway ;-)

BTW was the wall pressure implicitly treated? everything else at the wall is constant in time anyway..

``(1). For Mach 25 flow over a flat plate, you have the wall as the boundary condition. So, you said the local time step is limited to CFL = around one. ''

I apologize, I meant the CFL number in the initial stages of convergence is limited to approx. 1 for Mach 25, but is pretty much unlimited for low Mach numbers. After the solution has progressed in time and has almost reached steady-state, then the CFL can be pumped up to high values. But at that point, you might as well just stop the sim: you already have a decent answer.

``(2). If you remove the wall and change it to the free stream condition, do you still have the same problem(limit)? ''

No. But then the initial conditions will be almost exactly equal to the converged results, and (see above) then the CFL can be raised easily.

``I think, the limit is coming from the boundary conditions. What do you think? The governing equations are still the same non-linear equations, with or without the wall boundary condition.''

And the governing equations are also the same ones governing Mach=0.1 flows.. My point is, in certain cases, when the mechanical energy is very high wrt to the thermal energy and when this mechanical energy has a chance to get fully converted to thermal energy somehow (something that wouldn't happen for a Mach 25 flow bounded by freestream conditions), the linearization of the gov. eq. won't point ''as easily'' towards the right root. There's more than one root after all for the NS equations, and it is my belief that such a root is more easily attainable for some flows than others..

Another such flow which limits my CFL condition to 1 in the initial stages of convergence is the one in a quasi one-dimensional nozzle where a strong backpressure exists, creating a stable normal shock in the diverging section. The shock has a strength Ms=3 approximately. For such a flow, unless the solution has advanced in time considerably, I experience a limit to the CFL number, which doesn't exceed one considerably. For the same problem but with a low back pressure, hence resulting in supersonic flow throughout the nozzle, the CFL condition is not as limited in the initial stages of convergence.

But then again.. Maybe there's just some bugs in how I implement my implicit solver. Anyone out there with similar experience to mine? I'd really like to hear from you. Or anyone out there that can solve a backpressure 1D nozzle with a CFL of 1 million in the initial stages of convergence? Let me know.

-bernard

(1). Yes, I understand your problem. (2). The flow is transient in the initial stage. (3). In the final stage, the flow is reaching the steady state and is less sensitive the the time step size or CFL number. In the steady state stage, one can just remove the time derivative term (or using a very large time step to do so), and solve the steady state equations. the trouble is always in the transient state where the flow field and the boundary conditions are still changing rapidly.

At the beginning the gas is at rest at the temperature of the wall and a low density. Then Inflow boundary conditions are applied in one extreme of the domain, amd a shock wave propagates inmediately (the low density of the gas helped).The CFL initial is 1, I have tried with different initial CFL numbers but I had problems with initial numbers above 30. Then I increased the CFL from time step to time step at a ratio of 1.2 (this ratio is very sensible to the gris). The solution converge (density residual below 10^-8) in 100 time steps. The BC was treated with ghost cells, obviously the wall propierties are constant in time but not the ghost cell ones, which depend on the near wall cell. Converged solutions can be obtained also with freestream initial conditions applied, but the results are not as fast (in CPU time) as the rest gas(more "realisitic" initial condition???).

This definitely interests me. Can you be more specific in what your initial conditions are, what kind of turbulence model (or laminar NS equations?) you are using, what kind of implicit solver you use, etc? I'd really, really like to reproduce your results. Also, how big is your mesh? How big is your physical domain?

thanks~

Well my solver is using a Godunov Method using the HLLC Riemann Solver. Extension to second order through MUSCL schemes with different slope limiters. The implicit scheme is based on the work done by:

P. Batten,M. A. Leschziner and U. C. Goldberg "Average-State Jacobians and Implict methods for Compessible Viscous and Turbulent Flows" Journal Of Comp. Phy. 137,38-78(1997)

My BC: LEFT Boundary: INFLOW Nitrogen, at Ma=6.7 and T=57.8K and Pref=166 Pa (Reference Conditions, equal to the experimental conditions )

RIGHT Boundary: OUTFLOW NORTH Boundary: OUTFLOW SOUTH Boundary: Isotermal Wall at T=300K

Initial conditions: rho(t=0)=0.05*rhoREF u(t=0)=UREF v(t=0)=0 T(t=0)=Twall

The Grids tested are 49*48,98*96,196*192( I stopped there, I am interested in compression ramps and was not use to continue), the grids are strechted towards the wall and the leading edge.

The dimensions are L in the flow direction and 0.4 L in the stream direction. Due to the Reynolds number used the equations I am used are laminar, and the method developed is laminar, however the extension to turbulence models is contemplated in the original paper, and also the treatment with gas mixtures.

By the way which are your conditions??, maybe I can run a simulation with similar parameters to contrast results.

Hope to be useful

Hi Salvador,

thank you for your input, I appreciate it. I tried the same case you specified. I am using approximate-factorization as my implicit solver and a Roe approximate riemann solver turned second order accurate through Yee's limiters. The Roe scheme is fully linearized, as well as the diffusion terms.

For your problem, I specified L to be 0.1m, and the convergence history I experienced was similar to what I am used to: a CFL close to 1 until the maximum error reaches the domain exit (which occurs at an iteration count of 100), at which point the CFL number can be increased but very slowly, at a rate of 1.03 approx. per iteration..

I took a look at the paper you referenced: they also mention some limit to the CFL number limit in the initial stages for the Roe linearization, but which is non-existent for the scheme they propose. For this flat plate problem, did you try out the Roe linearization yourself? What kind of behaviour did you get?

Thanks,

-bernard

Sorry I did not try Roe's linearisation. The results obatined in the paper by Batten et al. showed good agreement between Roe and HLLC. I choose HLLC due to my previous experience in explicit solvers and the comment you pointed about the limit on the CFL in Roe's scheme. However the article mentioned that the reason of that phenomena it is not known.(actually I did not risk)

One more comment, just mention that I started with CFL near 1 and I started inmediately to increase the CFL number (different from you), when the shock wave reaches the domain exit, the CFL number is very big already. The ratio of increasing the CFL number is very sensitive to the grid size. With the coarsest grids I can increase the time step at a ratio of 1.2 and for finer grids this ratio decreases to 1.01 or less. It seems reasonable, that when there is a big change in the flow (shock wave propagating) the CFL must be bounded by some relation depending on the boundaries and the domain size. Similar phenomena I have found in 2 and 3D compression ramps. The CFL increase ratio was also very sensitive of the angle of the ramp (seemed reasonable ).

Have you tried in different configurations???

For any flow involving strong shocks propagating through the domain, I experience some limit on the CFL number which does not exceed 1 considerably, which seems to be rather independant on the type of flowfield solved. Only after the solution has stabilized quite a bit can I increase the CFL number. For flows where strong shocks are unlikely to occur (at lower Mach number for example), the limit on the CFL number is far less pronounced, and I can usually start at a CFL of 10 or so and increase it rapidly.

Maybe this behaviour is due entirely to the Roe linearization I am using or maybe I'm just implementing the Roe linearization wrong, which is possible afterall ;-). Could anyone using Roe's linearization confirm (or contradict) my results? I'd definitely like to hear from you.

I am definitely impressed by the stability of the implicit solver you are using and I'll look into it more closely.

Thanks for all your comments,

-bernard