Hi, everyone,

I am doing LES of jet diffusion flames. The problem I face now is the reaction terms in the energy equation (finite rate chemistry is adopted in my code), which make the equation a stiff system. Because the low-storage third-order Runge-Kutta method developed by Williamson has been used, it's preferrable for me to stick with higher-order Runge-Kutta method, but perhaps to fully solve this problem, I need to use some semi-implicit Runge-Kutta method, the implicit algorithm to handle the stiff reaction term. Now, I can find some references on semi-implicit RK methods, but cannot find any practical applications. The first thing I am worrying about is the feasibility of such semi-implicit RK method in practical LES, i.e. the additional simulation time needed in comparison with the fully explicit RK methods. For example, my current grid system looks like 40*90*150 (100*100*200 would be expected in future simulations), compact FD schemes are used in spatial discretization. To reach nondimensional time t=100 with 32 processors of a modern pc cluster, I need about 24 hours, which is not bad for me. So how about the situation with the semi-implicit RK method? Do I need 2 days, 1 week or even longer ... , to finish one case? Any experienced statements on this respect is very welcome. And if you have other ideas to handle the reaction term, please also give some advice. Thanks a lot.

Hi,

Pantano (J. Fluid Mech., vol. 514) applied a semi-implicit R-K in DNS of a reacting flow. So, it seems feasible.

Bobby,

Take a look at these:

http://fun3d.larc.nasa.gov/papers/carpenter.pdf

ftp://techreports.larc.nasa.gov/pub/...-2001-2612.pdf

http://www.hsa.lr.tudelft.nl/publica.../2002.2.01.pdf

http://www.math.ntnu.no/preprint/num...004/N1-2004.ps

Also, look into additive Runge-Kutta methods. See

Uri Ascher et al.(1997)

Maripaz Calvo et al. (2001)

Finally look into implicit/explicit multistep methods. See

Uri Ascher et al. (1995)

Thanks a lot, Tom and Runge-Kutta. All the material is quite useful to me.

Here I have aother question, would you please give some comments? The reason why I seek the feasibility of semi-implicit RK is that I want to handle the reaction term in the energy equation, which makes the equation a stiff system. However, other equations, like mass and momentum equations, are not stiff. So could I use the explicit part of the ARK scheme to explicitly advance other nonstiff equations? For example, in CA Kennedy's paper (2003) in Applied Numerical Mathematics, he developed 'ARK3(2)4LSA', which contains two parts, ERK and ESDIRK, to handle the stiff and nonstiff terms in one equation (or one set of equations). In my case, to apply this scheme to the energy equation is straitforward. Then, could I use ARK3(2)4L[2]SA-ERK to advance mass and momentum equations instead of the 3rd-order Williamson's scheme, which is adopted in my code? I think all the coefficents in this ARK2 scheme are designed to possess optimal properties, say, order and stability, in a combined way, with the concern of both explicit and implicit parts. I am not sure if a single ARK3(2)4L[2]SA-ERK scheme can be used (optimally) in a nonstiff equation, like other designed 3rd-order fully explicit RK methods. If not, I feel confused about the implementation of the ARK2 scheme in my case, I mean I have both nonstiff and stiff equations. Of course it may be implemented (with different RK methods for different equations), but with much more evaluation of RHS terms of the governing equations, not to mention the iterative solver for the implicit RK.

Another thing is, do you have any suggestions for the iterative solver? Is Newton iteration okay with my application? Which option is better, taking a free code from some math code library or developing by myself?

Thanks a lot!!!

Bobby,

For COUPLED additive methods, you can use either Runge-Kutta or multistep methods. If you are solving the compressible equations then you have a system of ODEs and both will work fine. The ARK lends itself to varying the timestep much better though. If you are solving the incompressible equations then the multistep methods given by Ascher et al. (1995) make sense because they don't care if your system is ODEs, index-1 DAEs, or index-2 DAEs. All common Runge-Kutta integrators are very sensitive to the index of the equations.

With any partitioned method, there are the individual methods and then there is the interaction between the methods. If the person designing the methods was good then both the individual methods and their coupling was all considered. Obviously, in doing this, the explicit method portion of the ARK might be somewhat compromised relative to what could have been with the coupling constraints. Remember that low-storage methods make the most sense when the memeory of the integrator is a big chunk of what the code uses. Typically, this is only in DNS.

I don't follow why you're confused. Look at some of these IMEX papers for implementation guidance. Take your nonstiff terms and feed them to the explicit method and take the stiff terms and feed them to the implicit method. The jacobian for the method is based only on the stiff terms.

You're on your own for the Newton solver. take a look at this and the references.

http://www3.interscience.wiley.com/c...63271/ABSTRACT

Hi, Runge-Kutta,

I am grateful for all of your valuable comments and all those important references, which save me loads of time to understand the whole method.

Actually you have exactly answered my doubts. My guess on this is similar, but I need experts' confirmation. Now I have been convinced that the ARK2 method can be used in my case. The only thing left seems the Newton iteration solver ...

As you mentioned, the low-storage property is of great concern in my mind. Although I am doing LES, which is not so stringent on memory usage as in DNS, it's problematic as well in the case of 3D simulations and more governing equations for reactive flows. Yoh JJ and Zhong X (AIAA J. 2004) have developed some low-storage version for their IMEX RK schemes, but unfortunately, their methods have been fully denied by Kennedy CA and Carpenter MH (Applied Numerical Mathematics 2003) for the stiffness leakage. So if you got some comments in this respect, I really appreciate it, because it's a profound further step in practitioners' minds.

Thanks a lot!!!

Bobby,

You must optimize your integrator choice over many attributes. Storage is one that has minimal consequence and you have little control over it. If you stick to an ERK, arguably the single most important matter for any group of reasonably designed methods is the error control/step size selection strategy when the integrator is living at the linear stability boundary. If you are going to use an coupled IMEX method like the ones you mention then your memory required is mostly caught up in the Newton solver. Here, storage reduction via Williamson or van der Houwen strategies are quite nearly pointless. Speaking of pointless, TVD RK methods are another annoying distraction.

If you want to read up on error control strategies for RK methods, look for papers by Soderlind and his students:

http://www.maths.lth.se/na/staff/gustaf/

There are two things to take away from these papers on IMEX Runge-Kutta papers. First is that one needs to consider the proper linear stability function for these methods and zero out the proper terms. It appears that most people have failed this. More subtle, even if you get the linear stability right, the internal stability of the implicit method can bite you. Caveat emptor!! Partitioned methods are complicated.

Could you please elaborate little more why "TVD RK methods are another annoying distraction"? What are their pros and cons. Some time in the past I was planning to use them but later moved to another project. They looked atrractive to me. In future I may still need to use some ode solver for hyperbolic systems, so I would like to make a more knowledgeable choice.

Angen

Angen,

I posted this response several months ago to a similar question:

================================================== ==

Let me refer to these methods, as a group, as SSP/TVD ERK methods. They arose independently by Shu and Osher (1988) and Kraaijevanger (1991) but were looked at in slightly different ways. Shu and Osher were concerned about the property of monotonicity while Kraaijevanger was looking at contractivity. Since then, a cottage industry has sprouted up to make lots of these methods. What has been lacking is a clear demonstration of where these methods offer an advantage over traditional ERK methods. A true test of the methods would be selecting a SSP/TVD ERK method and an ordinary ERK which posess nearly identical leading order truncation errors and linear stability domains and testing both of them on a battery of test problems while simultaneously using an error controller to maintain a fixed local error. Unfortunately, that might jeopardize ones ability to continue publishing papers on these methods. My personal opinion (and experience) is that this property offers little to no true benefit. It is 95% hype. Two of the people who publish this stuff knowingly left Kraaijevanger's work unmentioned and republished his methods. It was the cheapest thing I've ever watched in my career. Kraaijevanger's paper is really amazing. By the way, Kraaijevanger's advisor (Spijker) and another person are addressing the topic in a unifying way now.

http://scholar.google.com/scholar?q=ferracina+spijker&ie=UTF-8&oe=UTF-8&hl=en

It is 2005 now and the fluid dynamics community MUST start controlling their temporal error with error controllers. This stuff has been around for decades but nobody uses it. To agonize over SSP/TVD methods but fail to use an error controller is like building an F1 race car but never putting air in the tires.

Another thing about SSP/TVD ERK methods. One can devise countless measures of contractivity, monotonicity, or stability. Go look at all of the stability measures for Implicit Runge-Kutta (IRK) methods. There are maybe 20. Does this mean that your method needs to posess all 20 of them? The trick is to understand what is needed rather than making a method that satisfies the most demanding stability property. Imagine an ODE with one variable. Integrate that equation from some initial time, Ti, to a final time, Tf. Select an initial condition and plot the solution as time proceeds. Now redo the calculation by slightly perturbing the initial condition. You should have two very similar solution curves. Imagine you have the exact solution. If the distance between the exact solution curves approach each other as time proceeds, the equation is called dissipative. We would like the numerical solution to do the same otherwise any error we commit during the integration will jump us onto an adjacent and diverging solution trajectory. If the distance between the two numerical trajectories decreases as we integrate, the numerical method is called contractive. There are still two issues to be resolved. Is the equation linear or nonlinear and what norm shall we use to decide the distance between the two solution trajectories? Traditional linear stability is for linear problems when the distance is measured with an L2 norm. For implicit methods, algebraic stability is for nonlinear problems when the distance is measured with an L2 norm. SSP/TVD ERK methods are demanding contractivity on nonlinear problems in an L_infinity norm. If one asks for this property in IRKs, the maximum order of accuracy possible is one if stability is requested for the entire complex left half plane. This should tell you that since methods of lesser stability work very well, maybe the criterion is too severe. For stiff problems, L-stability (a linear stability criterion) is more useful than algebraic stability ( a nonlinear stability criterion).

My advise is to look for an ERK method with low leading order error that accepts an error controller well. That is more important than whether the method is SSP/TVD. By the way, here is Kraaijevanger's method WITH an embedded method. I'm not saying this is the best but it is partly what you want and partly what you need!!

a(2,1) = 0.39175222657188905833d0

a(3,1) = 0.21766909626116921036d0

a(3,2) = 0.36841059305037202075d0

a(4,1) = 0.08269208665781075441d0

a(4,2) = 0.13995850219189573938d0

a(4,3) = 0.25189177427169263984d0

a(5,1) = 0.06796628363711496324d0

a(5,2) = 0.11503469850463199467d0

a(5,3) = 0.20703489859738471851d0

a(5,4) = 0.54497475022851992204d0

b(1) = 0.14681187608478644956d0

b(2) = 0.24848290944497614757d0

b(3) = 0.10425883033198029567d0

b(4) = 0.27443890090134945681d0

b(5) = 0.22600748323690765039d0

bh(1) = 3298630537163.d0/18795049059328.d0

bh(2) = 1771072646821.d0/14103719399933.d0

bh(3) = 25277699738943.d0/87361585864100.d0

bh(4) = 1305111367099.d0/5907422403405.d0

bh(5) = 3614287973383.d0/19159025146192.d0

================================================== =============

Let me add a few other things. First, the region of nonlinear contractivity is a subset of the traditional linear stability domain. This means that if you are anywhere near the maximum stepsize of these TVD/Contractive ERKs then you are likely outside of the TVD/Contractive region of the method. Also the analysis applies to real eigenvalues - viscous terms. Next, here's an attempt by practitioners to assess the value of TVD methods.

http://www3.interscience.wiley.com/c...73119/ABSTRACT

Thanks for the detailed answer and references. If I understand it well the main points you are making are

1. Authors of SSP/TVD (i.e. Shu and Osher, 1988) did not mention Kraaijevanger (1991) work. Should they? (1988 vs 1991)

2. There are many variants of these methods and they offer very little over traditional ERK with time step control.

3. There is not time step control in Shu and Osher SSP/TVD.

4. Nobody is using time integration error controlers in CFD community.

I will comment only on points three and four. I did not read original Shu and Osher (1988) paper, but I did read later reviews. I noticed there is no error control strategy defined there. At first, me too I considered it to be a drawback, especially since I did use time error controlers combined with traditional ERK mainly with implicit methods or high order explicit pseudo-spectral methods applied to parabolic systems. However, I am wonedering if we really need time error controlers for explicit methods applied to hyperbolic systems. Afterall, here the main time step limitation is Cu=1 (or something very close to it like Cu=1.6 or 0.8) and not error limited time steps. If someone needs a superaccuracy the explicit methods can be combined with spatial grid adaptivity, where time step control is very nicely combined with error control, grid control, and Cu limitation. Moreover, it can be done at less cost because grids can be refined only locally.

Angen

Angen,

I wasn't referring to Shu and Osher but let's forget about naming names. Yes, I believe they offer little over traditional methods WHEN the local error is controlled properly. Nonlinear contractivity is an excessive requirement on these methods. It is my impression that few people use error control in CFD.

Error can come from many places; sometimes nicely predictable and other times not. The stability boundary is not necessarily always going to be where the uniform grid, periodic case suggests. Why integrate blindly right on top of the linear stability boundary? Everytime you accidentally cross it, the global error grows geometrically. It has been standard pratice in the ODE community to use error control for many decades. Why does CFD spend 99% of its time making better spatial discretizations and then use integration techniques from a long time ago? My impression is that most explicit integrations in CFD are done with a blindfold on and the gas pedal on the floor.

Maybe I'm being a bit harsh ...

In any event, there is still useful work to be done to make error control strategies for ERKs work better. They are not completely ready for CFD prime time. Too bad there is no funding for any of this stuff ...