4
$\begingroup$

I've recently written a code that solves the incompressible/low-Mach number formulation of the Navier-Stokes equation with high-order methods for both time and space. My advisor insisted that I use third order Runge-Kutta method. However, I do not find it to be adapted to our current problem.

In incompressible solvers, a poisson equation for the pressure must be solved, which is very CPU-time consuming and often the bottlenecks of such solvers (for low-Mach number problems, combustion might also be a bottleneck). When using a RK method, the poisson equation has to be solved at every RK-step, which is not efficient. Therefore, I think using a standard $3^{rd}$ order backward difference method might be more appropriate for our case.

To my knowledge, the main advantages of Runge-Kutta methods are that they are relatively simple to implement, self-starting and very stable. In my case, only the last of these of importance.

Are there any other advantages of the RK methods I am not aware of? Will I lose much accuracy/stability by switching to, say a $3^{rd}$ order backward difference method? Also, as mentioned in this post, I've proven than the spatial discretization error might be the leading error in my sims and therefore dictate the global accuracy of my code.

Improve this question
$\endgroup$
3
  • $\begingroup$ A third order RK method vs a third order backward difference method should obviously produce accuracy on similar orders of magnitude wrt time integration. If spatial errors are dominant, then worrying about time integration error isn't likely as big of a deal. Stability is obviously important, but I can't say which scheme is more stable for your problem with what I know. It probably isn't much effort to try running your code using both approaches and just comparing on a few metrics like average runtime for some case, comparing their solutions to see which is more realistic for some case, etc. $\endgroup$
    – spektr
    Feb 16, 2017 at 14:53
  • $\begingroup$ @choward No it is not much effort which is why I am currently doing it ;) Right now I am just testing the code with some manufactured solutions for both formulations. However, the objective of my research is reacting turbulent flows, which are known to be stiff and unstable systems. Even if both methods yield similar results now, I do not know if it will be the case later on. Which is why I wanted someone's output on this, in case there was something I was missing. $\endgroup$
    – solalito
    Feb 16, 2017 at 15:00
  • 2
    $\begingroup$ My suggestion is to go to the work of the guru of Runge-Kutta methods, John Butcher, and see what he has to say on the matter. See scholarpedia.org/article/Runge-Kutta_methods and the references cited there. $\endgroup$
    – Biswajit Banerjee
    Feb 17, 2017 at 0:56

2 Answers 2

Reset to default
2
$\begingroup$

The key issue in this discussion is explicit versus implicit ODE solvers. The traditional Runge-Kutta (RK) methods are explicit and so the time step must be small enough to satisfy the stability requirement. The backward difference methods (BDF) are implicit so the time step can be chosen based, simply, on what is needed for an accurate solution; typically much larger than for an explicit method. I assume that is basically the point of your question about RK versus BDF.

The incompressible Navier-Stokes (NS) equations present a very interesting challenge to both classes of methods. This is discussed in many books on computational fluid dynamics. I found this particular survey, Langtangen, to be an especially nice introduction to the subject.

Specifically, for low-Reynolds number flow, the viscous terms are important which makes implicit methods attractive because these terms dictate a very small time step in explicit solvers.

However, the convective terms can be integrated with an explicit solver and a time step that is much larger. Furthermore, these terms are nonlinear which makes an implicit solver much more unattractive because a set of nonlinear algebraic equations must be solved at every time step. This system includes all the velocity degrees of freedom so can be quite computationally costly.

So, as discussed in, for example, the referenced paper, state of the art algorithms for the incompressible NS equations often use an explicit solver where it works best-- integrating the convective terms-- and an implicit solver where it works best-- integrating the viscous terms.

Normally I would suggest discretizing the equations in space and then using an off-the-shelf ODE solver. But since most standard ODE solvers implement just an explicit or an implicit algorithm, using one of these would not be the optimal approach for this problem.

For instructional purposes, Benjamin Seibold, wrote a NS solver based on these principles, Seibold, in, remarkably, a little over 100 lines of MATLAB code! It is well-worth the time studying and experimenting with it.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Saying a library can't do this is just leading someone astray (and is just factually false). It's becoming more standard for "off-the-shelf" ODE solvers like Sundials' ARKODE to have implementations of "IMEX Methods" for splitting explicit and implicit parts. They will allow you to choose what should be implicit, and what should be explicit, and will do so with all the goodies you want. Of course, MATLAB/SciPy can't do that, but that's just a deficiency of MATLAB/SciPy, which is not representative of more comprehensive and performant solver eocsystems. $\endgroup$
    – Chris Rackauckas
    Feb 21, 2017 at 19:47
  • $\begingroup$ You apparently didn't read my post very carefully because the main point is about the characteristics of the different terms in the Navier-Stokes equations. Instead, you focused on a small part of a single sentence. Nevertheless, I modified that sentence to reflect your comment. $\endgroup$
    – Bill Greene
    Feb 21, 2017 at 20:44
  • $\begingroup$ No. The conclusion "using one of these would not be the optimal approach for this problem" still doesn't follow from what you're saying. Yes, the terms of the Navier-Stokes equations have different characteristics. But you can still use an off-the-shelf ODE solver like Sundials in a way that makes use of those characteristics/splittings, the same way as you would in the survey you linked. Though not as straight-forward as "throw it into MATLAB", it's still in the end there's nothing in that survey which isn't doable in something like Sundials (with the banded linear solver). $\endgroup$
    – Chris Rackauckas
    Feb 21, 2017 at 20:59
  • $\begingroup$ Though the staggard grid methods do bring up an interesting thing that to get the full efficiency of a perfectly hand-optimized code, an ODE solver would want to make use of not just a banded linear solver, but also the fact that the Jacobian, if representing the system as a single vector, would have a block structure. Making use of the block structure may make a hand-optimized code faster (if you implemented a full adaptive BDF etc.). $\endgroup$
    – Chris Rackauckas
    Feb 21, 2017 at 21:02
0
$\begingroup$

If you discretize the problem in to ODEs, then you shouldn't need to worry about this because you likely shouldn't be writing the ODE solver anyways. At that point, plug the discretized ODEs into whatever solver is available in your programming language and swap around, try BDF methods or Runge-Kutta methods, and see what's faster. You're not going to beat that by writing a special purpose method, and using an ODE package will have a lot of stuff already built into it (like adaptive timestepping) which will be very helpful for solving the problem.

Improve this answer
$\endgroup$
4
  • $\begingroup$ The problem is that I can test it all I want with any kind of simple 2D manufactured or real solution, it won't help me much. In the end, what matters is the performance/accuracy/stability of the solver for 3D turbulent reacting code. And for that, I can't afford to switch around. Which is why I was looking for some additional insight $\endgroup$
    – solalito
    Feb 20, 2017 at 20:09
  • $\begingroup$ Well there's two things tied up in there. First of all, a 2D problem with the same characteristics will have the same timestepping issues as the 3D problem, so you can use that to find the method which works best. But secondly, for a 3D problem, you don't want to be using hand-written unoptimized code. You want to use a highly optimized high-order adaptive code with high order interpolations to get solve the equations the fastest, which means you might as well switch over to using some package anyways since it'll be faster... $\endgroup$
    – Chris Rackauckas
    Feb 20, 2017 at 21:12
  • $\begingroup$ Writing the ODE solver is the smallest problem in writing a flow solver really. An explicit Runge Kutta loop is just a few lines... Filling some matrices into a format accepted by some library or writing a callback function will take considerably more time. $\endgroup$
    – Vladimir F Героям слава
    Feb 22, 2017 at 14:18
  • $\begingroup$ An explicit Runge-Kutta loop would be really awful though. You'll save a lot of time in the actual calculation (if it takes a long time) if you have proper order adaptivity, PI-controlled stepsizes, etc. It's not hard to find problems where this will be more than hundreds of times faster than a simple RK loop. So yes, you can do this easily if you don't care about speed. But the OP seems to want speed. $\endgroup$
    – Chris Rackauckas
    Feb 22, 2017 at 18:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.