Shock tube example + shock capturing

[tt-aiken]

I would like to set up a new problem in examples/fluids/ based on the Sod shock tube, to be used as a benchmark for an implementation of the YZB shock capturing term introduced by Tezduyar and Senga in 2005. There are more modern DC operators, but I chose this one because it's well proven and has a pretty simple formulation. A shock tube is good because there is an analytical solution, so the effectiveness of the operator can be quantitatively evaluated.

[jedbrown]

Cool! To clarify, would this be a quasi-1D experiment (perhaps implemented as a "pencil" domain with slip or periodic boundary conditions)? With a shock capturing scheme, we'd also be able to solve the richer benchmark problems such as shock-bubble interaction. (I'm mentioning this, but not trying to create scope creep.)

This would be a good chance to better decouple the physics (compressible Navier-Stokes with gravity) from the "density current" set of initial and boundary conditions. Cc: @LeilaGhaffari

[tt-aiken]

Yes, my plan is to run the case as 1D using a long line of elements with slip BCs on the sides. I think that a BC to capture the end wall (slip + no penetration) would be more interesting, to see how the reflection is resolved, but I'm open to suggestions on what would be most useful. Of course there could be a command line option to use wall or periodic, but I only want to choose one for now.

I was planning to implement this using only the Euler equations, since that's what the analytical solution is for. Would this be of too limited use? A full Navier-Stokes implementation would provide a more difficult test of the YZB term, particularly looking at its impact where the shock and boundary layer intersect (3D), but I don't know how that would be verified.

[jedbrown]

So what currently lives in eulervortex.h, except for the initial condition, is generic Euler and I think you can use it, dropping in your shock capturing scheme with a parameter that makes it possible to turn off. Alternatively, you could start from densitycurrent.h (which has SUPG for N-S). Do you want explicit or implicit time stepping?

[tt-aiken]

I plan to do both implicit and explicit. How about I start with Euler, do implicit and explicit, and then move to full Navier-Stokes. This way the 1D Euler shock tube can be used as verification, and the 3D Navier-Stokes shock tube can help with exploring the impact of tuning parameters for the artificial viscosity on overall flowfield behavior (ex. boundary layer growth, shock attenuation). I think this is well scoped given the scaffolding that's already in place with the RHS and IFunctions in eulervortex and densitycurrent.

[jedbrown]

Sounds good to me.

[knepley]

What kind of implicit solver?

[LeilaGhaffari]

This would be a good chance to better decouple the physics (compressible Navier-Stokes with gravity) from the "density current" set of initial and boundary conditions. Cc: @LeilaGhaffari

I am happy to work on the decoupling if that helps with Tim's project. Just to clarify, when you say "compressible Navier-Stokes with gravity", do you mean excluding gravitational potential from total energy (E)?

[LeilaGhaffari]

What kind of implicit solver?

We use PETSc's TSSolve() in our fluid examples.

[jedbrown]

I prefer leaving it inside total energy because that makes it easier to check that energy is conserved. The point is we retain the ability to have g nonzero, though some of Tim's experiments may select g=0.

@knepley We've focused on time-accurate solvers so we use RK for explicit time integration and either alpha or RosW for implicit integration, as with SUPG.

[knepley]

@jedbrown Just to help me understand, why do implicit if you are resolving everything?

[tt-aiken]

Can this decoupling of the Navier-Stokes (and Euler?) solvers be accomplished simply using a new header file? I've copied the eulervortex RHS and IFunctions into the shock tube header file for the time being, which doesn't feel like a long term solution.

[jedbrown]

@knepley Depending on the problem, viscous stiffness, acoustic stiffness, and time accuracy with SUPG (can't "lump" the stabilization term, which includes a gradient of the test function).

@tt-aiken @LeilaGhaffari How does this proposal sound:

1. Rename problems/eulervortex.c to problems/euler.c and qfunctions/eulervortex.c to qfunctions/euler.c.
2. Use PetscOptionsFList to select -euler_problem vortex or -euler_problem shocktube. These would call a function EulerSetUp_Vortex or EulerSetUp_ShockTube, which sets the qfunctions for initial and boundary conditions. So you'll have Exact_Euler_Vortex and Exact_Euler_ShockTube. The latter might only know the solution at time 0, which is okay.
3. The IFunction and RHSFunction for Euler can be the same (for vortex and shocktube, or any other problem using the Euler equations), but we need to add stabilization. My inclination would be to add a scalar parameter for which 0 means no stabilization, just clearly delineate it in the code so the reader can see what is stabilization versus what is the physics of the problem.

@tt-aiken Will your shock capturing formulation be combined with SUPG or will it be independent (each could in principle be turned on/off separately)?

[tt-aiken]

@jedbrown This sounds reasonable based on the current organization I've seen. How would this apply to the Navier Stokes routines contained within density_current, given that navierstokes.c already exists?

As for the YZB operator and SU/SUPG, my plan was to have them be decoupled. This gives more flexibility to mix different stabilization schemes if they are added down the line.

[LeilaGhaffari]

@jedbrown , this sounds good to me too. This would help us avoid repeating the Euler solver for each of those problems. As for stabilization, our current SU/SUPG is not functional in density_current nor in euler_vortex. I was planning to debug that with MMS using Enzyme after merging #834 . Can Tim put this off for now?

[jedbrown]

Heh, I think we should rename navierstokes.c to main.c and create an executable called fluids or compressible, but that's just me.

But for now, let's just split up the Euler part so you can focus on your project. We can apply lessons learned to viscous models later.

If YZB can be used without SU/SUPG, then I'd just proceed with it. If this stabilization is needed, I'd add pure SU to Euler -- it's a simpler operator and can still be used with explicit RK time stepping, though it'll blur the contact wave in particular more than SUPG.