Blow up

[Andrew Cook]

Hi all,

I have been trying to use dedalus to run compressible navier stokes simulations using the equation of state P=C_0^2*rho, P is pressure, C_0 is speed of sound, and rho is density. With the time dependent boundary conditions I am using I am getting blow up, but it seems to fine with simple boundary conditions. I have ran this same exact problem in COMSOL with the same boundary conditions, initial conditions, and equation. In COMSOL there is no blow up and I am getting the expected result. I was wondering if there is something about dedalus that makes the specific scenario I am trying to simulate hard for dedalus to deal with or is there some mistake in how I wrote the problem. I am aware that there are other dedalus scripts that simulate compressible navier stokes, but these one add additional equations for entropy which I would like to avoid if possible.

Attached is my dedalus script, note I added -2*mu*div(E) on both sides of the equation so the solver isn't singular and doesn't change the actual equation.

[Ben Brown]

Andrew,

     Well timed, in that Kyle Augustson and I were just spending a bunch of time today talking about optimal solutions to fully compressible equations.  It's something We've been interested in for a while.

There are aspects of your equations that aren't entirely clear to me.  You have a variable "p" that's used in the continuity equation as well as the momentum equation, but no actual separate pressure evolution equation (or energy equation more generally).  Are these isothermal fully compressible equations, where you assume the temperature is a constant?

Generally, from our experience, I'd strongly recommend solving a log(density) continuity equation.  This has a variety of benefits, including protecting the density rho from crossing through zero (which it cannot physically do).  There are also numerical benefits, and in particular your grad(p)/p nonlinear RHS term would become a grad(ln_rho) term which could be handled linearly on the LHS.

Broadly, you have very little on the LHS to help keep these equations stiffly-stable.  This may not matter if you're operating a very high-Mach flow.  But based on d0 being ~1e-9 in your boundary forcing, I'm guessing you're actually at fairly low Mach.  In that landscape, stabilizing the equations via stiffly stable IMEX solving becomes key.  You can see a variety of examples of how to handle this in work done by Evan Anders and others (e.g., Anders & Brown 2017 PRL, Anders et al 2023 Nature, Powrs et al 2024 PRL, etc. plus early work by various of us in Lecoanet et al 2014 ApJ, Lecoanet et al 2016 MNRAS).

You can also see current "best-in-class" implementations in:

https://github.com/bpbrown/fully_compressible

There, we're focused on convection problems, but the core approaches in the entropy/enthalpy equations (FC_poly_h-s_2.5d.py) seem likely to be some for the best for a variety of reasons, including that these now effectively include log thermodynamics for both density and temperature (enthalpy).  We're working on describing some of those benefits and approaches now.

You're already doing one of the most important things that we've learned (the hard way), which is that your continuity equation also includes the tau polynomials from the velocity equations.  So well done there.  Generally, especially for low Mach number, I'd work to move terms from the RHS to the LHS to provide stiff stability and permit you to timestep on a velocity CFL timescale rather than being strongly constrained by fast (but dynamically unimportant) linear sound waves

--Ben

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/4a23a1af-fb50-4d9b-a551-714c4949e70fn%40googlegroups.com.

[Andrew Cook]

Hi Ben,

Thank you for the quick responses. In my equations I am solving for barotropic flow so constant temperature, that's where I got the constraint P=C^2*rho. and then when I wrote p I was writing pressure.

just for more background on the equations, I am trying to simulate low frequency acoustics directly instead of using the perturbation method which I have done in dedalus. I have done the problem, directly using COMSOL, and for these simulations the only equations I had were compressible navier stokes where the temperature is constant which agreed with the pertubation results I got earlier. Essentially I want to run the simulation I ran in COMSOL in dedalus. Eventually I want to run simulation at higher oscillations, but first I am trying to figure out the problems at the lower frequencies of 1 kHz.

I modified the continuity equation as you suggested and now it does seem to be working in the sense that it now doesn't blow up. I haven't done a throughout analysis, but the velocity seems to be what I would expect. This is great, but I have one last hurdle which is how do I recover the Pressure from Ln_rho? nothing is matching the order I have found in either of my two other simulations. I know for my comsol the ambient pressure was set to 1 atm.

-Andrew