RK Schemes Struggle With Divergence at Low CFL Numbers

[clarkpede]

I've run into some problems using all of the Runge-Kutta algorithms (purely explicit, EDIRK, and SMR91) with a 2D turbulent channel case. Using a convective CFL of around 1, I can advance through time with implicit Euler just fine. That's no surprise, since it's fully implicit. But when I switch to a RK algorithm, using the same timesteps, the code diverges. The residuals rise with each iteration, until the errors are too large for the linear solver.

I've even tried lowering the CFL to 0.01, and I still get divergence after 20 or so iterations. I can understand unstable behavior at high CFL, but CFL=0.01 should be stable.

This is very strange behavior. Any ideas as to what's happening here?

Edit: I should note that I'm using a fixed timestep, not a fixed CFL. I'm estimating the CFL number for a given timestep using the grid spacing and the velocity from the DNS data.

[clarkpede]

Alright. I'm running tests using the 2D mesh in our TestCases repo. According to the grid, CFL = 1 occurs at approximately 0.00363, or 3.63E-3.

In order to check the error, I've been using the following workflow:

1. Get a converged 2D solution (density residuals <= 1E-8)
2. Copy the restart file to form an unsteady restart (e.g. cp restart_flow.dat restart_flow_00000.dat)
3. Do a second run, this time with timestepping for 100 timesteps and a fixed timestep.

I know there's a problem when either the linear solver fails to converge or I get an error because the linear solver exceeds the maximum iterations (when I set LINEAR_SOLVER_MAX_ITER_ERROR= YES). Usually, the residuals start climbing before the linear solver chokes.

With Roe, no slope limiters, and low Mach correction, I can get up to:

• dt = 7E-3 (CFL ~ 2) for Euler Implicit
• CFL ~ 1E-4 for Runge-Kutta LIMEX EDIRK

With Roe, no slope limiters, and no low Mach correction, I can get up to:

• dt = 7E-3 (CFL ~ 2) for Euler Implicit
• CFL ~ 1E-3 for Runge-Kutta LIMEX EDIRK
• CFL ~ 1E-4 for Runge-Kutta Explicit

With JST, k2 = 0.0, k4 = 1e-4, and no low Mach correction, I can get up to:

• dt = 7E-3 (CFL ~ 2) for Euler Implicit
• CFL ~ 1E-3 for Runge-Kutta LIMEX EDIRK
• CFL ~ 1E-4 for Runge-Kutta Explicit

So the results are pretty consistent across spatial numerical schemes. The limit on RK explicit makes sense; it's crashing around the stability limit that SU2 calculates for viscous + convective + sonic stability. But the limit on the RK LIMEX EDIRK isn't much better at all. That is surprising.

I've checked, and the behavior is very similar for SA and SST models.

[clarkpede]

I'm attaching an example set of *.cfg files and an example restart field, in case you're having trouble reproducing the problem.

On a side note: I haven't tested this with other RK LIMEX schemes. SMR91 seems to be broken (see issue #18), and I haven't played around with alternative EDIRK schemes (beyond the default 3-2 stage scheme).

RK_tests.tar.gz

[clarkpede]

A correction: I was running some tests with JST to try to get a hybrid run set up, and I couldn't reproduce my earlier results. Here's my new results:

With JST, k2 = 0.0, k4 = 0.0039, and no low Mach correction, I can get up to:

• dt = 3E-3 (CFL ~ 1E-2) for Euler Implicit
• CFL ~ 1E-2 for Runge-Kutta LIMEX EDIRK
• CFL ~ 1E-3 for Runge-Kutta Explicit

Two important points:

1. With JST and low values of artificial dissipation, the RK schemes and Euler Implicit perform equally poorly. The LIMEX RK schemes aren't that much worse, but we still aren't getting close to the convective stability limit.
2. JST has a lot of trouble solving the linear system and timestepping without divergence without a certain level of artificial dissipation. I can get to much larger timesteps with much higher levels of artificial dissipation (closer to the default values).

[trevilo]

I think the problems described here due to two things. First, I think there are bugs in the EDIRK implementation, which I will fix. Second, and more importantly, when using the MUSCL reconstruction,
we probably shouldn't expect as much of an improvement in the time step as I originally expected due to the fact that the reconstruction is neglected in the Jacobian. This approximation of the Jacobian is problematic because it is wrong at high wavenumber, which is where the most restrictive time step stability constraints come from. The two plots attached below illustrate the issue.

mwn.pdf
imex_explicit_mwn.pdf

The plots show results of modified wavenumber analysis of a standard 2nd order central difference (on uniform mesh with periodic BCs... usual case) compared with a 2nd order central scheme where the gradient is used to reconstruct the face values. This is analogous to SU2 with

SPATIAL_ORDER_FLOW= 1ST_ORDER

versus

SPATIAL_ORDER_FLOW= 2ND_ORDER

but with a central flux (no upwinding or artificial dissipation of any kind). Because the reconstruction is neglected in the Jacobian, the approximate Jacobian corresponds to the standard second order central scheme no matter what. So, when this Jacobian is used to form the implicit part of the linearized IMEX scheme for the spatial scheme with reconstruction, a significant chunk of the eigenvalues corresponding to the high wavenumber modes are left in the explicit part of the scheme, leading to a smaller stable time step than one would like. Specifically, based on this analysis, you'd expect about a factor of 2.4 increase relative to the purely explicit part of the scheme.

I'm working to confirm this hypothesis at the moment. It does seem that the SMR91 scheme is stable at much larger dt than @clarkpede observed above when using ROE and 1ST_ORDER, which is consistent with expectations based on the above reasoning. But, obviously this scheme is also much more dissipative than the 2ND_ORDER version, which muddies the comparison a bit.

Anyway, if it is correct, I don't see a way to avoid this problem in the linearized scheme based on an approximate Jacobian, so we're probably going to have to move to methods that require a nonlinear solve (e.g., BDF or Crank-Nicolson) and use either dual time stepping (as already implemented) or a Newton method (with the approximate Jacobian) as the solver.