Dual time stepping for fluid dynamics

Question

I'm attempting to implement the Weiss and Smith preconditioner in an existing finite volume code and I am struggling with the idea of dual time stepping. My inner time steps are predictor-corrector, second order explicit with local time step acceleration. The time step is based on the stable CFL number based on the convective speed since it is using the preconditioned values.

My outer time step is a second order backward difference implicit method. This is all straight from the Weiss and Smith 1995 AIAA paper. The outer time step is freely chosen to resolve whatever features of interest are desired. The inner time steps over pseudo-time are performed until convergence, in theory when the pseudo-time is marched to infinity. When this happens, the solution should be at the physical time step used for the outer scheme.

I have two issues with this that are related. Right now in my code, the flow evolves to the time determined by the number of iterations times the inner time step rather than simply becoming my outer time step. For example, if I run a vortex convection case convecting at 1 m/s with an outer time step of 1 s and an inner time step of, say, 0.001 s, I would expect the vortex to move 1 meter each outer timestep. However, it actually convects 0.001*<number of inner steps> meters.

So I feel like I'm missing something with the whole concept. How does my inner step converge to a physical time of 1s per outer time step rather than a physical time of the number of inner steps times the inner time step?

Answer 1

I believe you're switching around the definitions of "inner" and "outer" stepping (this nomenclature gets worse in segregated schemes where you have at least 3 separate iterative schemes going on.)

The core of what you're missing is that "pseudo-time" isn't a time at all; it's just a convenient way to achieve an iterative solution. Take a 1st order backward difference in time equation (spatial discretization doesn't matter, except that we're evaluating it at time n+1):

$\frac{\partial u}{\partial t} + f(u)\approx \frac{u^{n+1}-u^n}{\Delta t} + f(u)^{n+1} \equiv R(u) = 0$

Now that expression is an implicit equation and you generally need some kind of iterative method to solve it for time $n+1$. Here you can use any (non-)linear algebraic solver to find $u$ such that $R(u)=0$. An alternative is to cast that problem as an unsteady one $\frac{du}{d\tau}=R(u)$ and step forward until the residual $R$ is 0 (I think of this as equivalent to using Jacobi iteration for solving R(u)=0). You may have seen this as a first method for solving laplace's elliptic equation $\nabla^2 u=0$ by posing it as a parabolic equation $\frac{du}{dt} = \nabla^2 u$ and stepping till convergence.

We can modify that expression to this:

$R^{\star}(u) = \frac{u-u^n}{\Delta t} + f(u)$

and write:

$ \frac{u^{\star}-u}{\Delta \tau} + R^{\star}(u) = 0$

Where $u^{\star}$ is the value at some pseudo time and $\Delta \tau$ is the pseudo-timestep. This is equal to the above whenever $u^{\star}=u$ which implies $u=u^{n+1}$. That we can step forward explicitly in pseudo-time. The subscripts/superscripts for dual time-stepping always end up getting a little cloudy - see "Multigrid Unsteady Navier-Stokes Calculations with Aeroeleastic Applications" for example of 2nd order.

In implementation, dual time-stepping is really easy once you wrap your head around how it works. It's algorithmically equivalent to your original explicit time-stepping scheme except now you have an additional source term representing the implicit unsteady term, and for each time step you iterate until the residual is near-zero. Another way of looking at this is to realize that each physical unsteady timestep is equivalent to solving a single "steady" problem.