steady Navier–Stokes by natural parameter continuation

[gdmcbain]

This is extends the creeping flow over the backward-facing step from ex24 to finite Reynolds number using pacopy.natural parameter continuation in the Reynolds number.

This example is twice the length of the previous longest (ex04 and ex24) and still requires substantial rewriting and tweaking, but it's starting to take form and starting to work so I thought I'd share it early.

It also motivates the proposal to add ‘milestones’ to sigma-py/pacopy#5 (so that the solution can be obtained at specified values of Reynolds number).

Natural parameter continuation is preferred to arc-length continuation (pacopy.euler_newton) as used in ex23 and discussed in #119 because there's no attempt here to track bifurcation, so the solution is monotonic in the parameter. Physically, the bifurcation involves three-dimensional effects; see, e.g.

• Barkley, D., Gomes, M. G. M. & Henderson, R. D. (2002). Three-dimensional instability in flow over a backward-facing step. Journal of Fluid Mechanics, 473, 167–190. doi: 10.1017/s002211200200232x

[kinnala]

An interesting way to structure the example. I guess it's related to pacopy and its interface? I'll read it more carefully later.

[gdmcbain]

An interesting way to structure the example. I guess it's related to pacopy and its interface?

Yes, for pacopy

The only thing the user must provide is a class with some simple methods, e.g., a function evaluation f(u, lmbda), a Jacobian a solver jac_solver(u, lmbda, rhs) etc.

I figured that since a class was needed anyway, I might as well pack the other data and methods in it too. It can be judged at the end whether this is clearer or more cumbersome.

(While working on this, it also occurred to me that a similar structure might be useful for abstracting time-stepping, where the residual function f would also depend on the rate of change of u and the jacobians with respect to u and its rate of change would be the stiffness and mass matrices, respectively, but that's for later.)

[gdmcbain]

Once this is done, quite a lot of code could and should be shared with the creeping version, ex24.

[gdmcbain]

Ah, I had the indices around the wrong way on the vector-gradient! Is the order of indices documented somewhere? (I'm not sure that it shows up in the elastic examples because the tensors there are symmetric.) Anyway, after 1aaa7f0, I'm getting the right answers now, like Barkley, Gomes, & Henderson (2002). Something like $u_j v_{i,j}$ would be written u[j] * dv[i][j]which makes sense, but I had been thinking more like $u_j \partial_j v_i$ which led me astray.

[gdmcbain]

It doesn't take too long to run now (about a minute and a half, depending on fineness of mesh and parameters for Newton iteration and continuation, etc., the hardware), but here are some results.

The messy Reynolds numbers are an artefact of pacopy.natural, hopefully to be fixed by sigma-py/pacopy#6.

[kinnala]

Impressive work. I believe this starts to be ready for merge?

[gdmcbain]

Thank you. (I confess that I am quite excited about this one.)

Before merging, I would like to sort out the steps in the Reynolds number. They really should be at round values, as, e.g., by Gartling (1990, §3):

The backward step flow was formulated and solved as a steady flow problem. In order to obtain a solution for the specified Re of 800, a series of steady state solutions were obtained at intermediate Re-values of 200, 400 and 600. Zeroth-order continuation was used to advance from one Re solution to another. The indicated sequence of steps in Re does not represent an optimized path to the final required state but was selected for convenience and with a view towards comparing solutions at other standard conditions.

The adoption of pacopy.natural might have been premature; its use here will probably involve lots of tweaks and that might be unreasonably disruptive upstream.

Locally, there's also:

• refactoring with the creeping version in ex24
  • (I haven't decided what the best way to do that is, or about what was perhaps charitably described as an ‘interesting way to structure the example’; the big class makes this example seem more different to those preceding than it really is.)
• organizing the sequence of viscous incompressible flow examples in the documentation (ex20, ex18, ex 24, this), as mentioned in Lagrange multipliers etc for Dirichlet conditions #150

but they could be deferred till after merging.

[gdmcbain]

O. K., sigma-py/pacopy#6 has been accepted so the natural parameter continuation hits a list of specified Reynolds numbers (150, 450, 750) and plots the stream-lines there. I think this is ready to merge.

[gdmcbain]

In the course of developing this example, I learnt numpy.einsum and was impressed. It occurred to me that it might be possible to use it to write a higher-level interface to the bilinear_form and linear_form decorators in skfem.assembly.asm, something like FEniCS's UFL or nutils.expression.