High Mach Number Turbulence can cause NaNs

[BraydenAlbery]

I am running a simulation to model stellar winds. I need an initial turbulence injection to model a realistic ISM.

When injecting turbulence (even with all my physical sources off, just an ambient background with hydro) after a few timesteps I get NaNs for all variables in some cells. This then spreads to the rest of the domain over a few timesteps. At the start, the velocities are on the order of km/s^-1 (realistic for the ISM) and the mach number is between 3-120.
At a certain point. The velocities will explode to 15 orders of magnitude larger and the mach numbers become on the order of 1e10

I have tried various combinations of spatial reconstruction (up to 4c), Riemman solvers (Roe, LLF, HLLC, LHLLC) and time integrators (VL2, RK2, RK3, SSPRK(5,4)) all with the same outcome.

Higher resolutions seem to last longer before the issue. 32^3 tests cause it to happen at around 4e12 s, while 512^3 on a HPC happens around 16.e13 s. Different spatial reconstructions also effect the time it takes. 4c causes it to happen on the first time step.

Has anyone else encountered a similar issue or have a solution? Or any way I could go about solving this?

Thanks for any help!

This is my input file:

\<comment>

problem   = 3D_Hydro_Stellar_wind # name

reference = Brayden # author

configure = python3 configure.py --prob stellarwindv9_thread-check -fft --nghost 4 --nscalars 1 -mpi -hdf5 -h5double --hdf5_path ./hdf5-1.14.3/hdf5 --include ./hdf5-1.14.3/hdf5/include --lib_path ./hdf5-1.14.3/hdf5/lib --fftw_path ./fftw-3.3.10 # problem run command



\<job>

problem_id = stellarwind-production   # problem ID: basename of output filenames



\<output1>

file_type  = hst       # History data dump

id = hist              # file name

dcycle        = 1.0    # time increment between outputs



\<output2>

file_type  = hdf5      # HDF5 data dump

variable   = prim      # variables to be output

id        = primatives # file name

dcycle      =   1.0    # time increment between outputs



\<output3>

file_type = hdf5     # HDF5 data dump

variable = cons      # variable output

id = conservatives   # file name

dcycle = 1.0            # output time step



\<output4>

file_type   = hdf5 # HDF5 data dump

variable    = user_out_var  # Output P/rho, temp, mu, mol atom and ion frac, energy (energy might be redundant)

id          = thermo # file name

dcycle          = 1.0 # time step



#\<output5>

#file_type  = rst        # restart file

#dt         = 1.0e12     # time increment between outputs // 5e9



#\<output6>

file_type = hdf5 #f

#id = slice #f

#variable = prim #f

#x3_slice = 0.0 #f

#dcycle = 1.0 #f

#ghost_zones = 1 #g





\<time>

cfl_number  = 0.1      # The Courant, Friedrichs, & Lewy (CFL) Number

nlim        = -1        # cycle limit

tlim        = 1.57788e14      # time limit

integrator  = rk2     # time integration algorithm

xorder      = 3c          # order of spatial reconstruction

ncycle_out  = 1000         # interval for stdout summary info



\<mesh>

nx1        = 32        # Number of zones in X1-direction

x1min      =  -1.2342710325965468e20       # minimum value of X1

x1max      =  1.2342710325965468e20       # maximum value of X1

ix1_bc     = outflow   # inner-X1 boundary flag

ox1_bc     = outflow   # outer-X1 boundary flag



nx2        = 32        # Number of zones in X2-direction

x2min      = -1.2342710325965468e20      # minimum value of X2

x2max      = 1.2342710325965468e20      # maximum value of X2

ix2_bc     = outflow  # inner-X2 boundary flag

ox2_bc     = outflow   # outer-X2 boundary flag



nx3        = 32        # Number of zones in X3-direction

x3min      =  -1.2342710325965468e20      # minimum value of X3

x3max      =  1.2342710325965468e20      # maximum value of X3

ix3_bc     = outflow   # inner-X3 boundary flag

ox3_bc     = outflow  # outer-X3 boundary flag



refinement  = none    # amr level



\<meshblock>

nx1        = 16        # Number of zones in X1-direction, meshx/meshblockx = maximum number of mpi threads

nx2        = 16    # Number of zones in X2-direction

nx3        = 16       # Number of zones in X3-direction



\<hydro>

gamma = 1.666666666666667 # gamma = C_p/C_v

dfloor     = 1.0e-30    # required otherwise Athena++ uses a default value

pfloor     = 1.0e-30    # required otherwise Athena++ uses a default value

sfloor     = 0.0        # required otherwise Athena++ uses a default value



\<turbulence>

dedt =  2.567118990202951e+51 # total energy injection for turbulence

expo = 6 # k expoentent

nlow = 2 # low n = k*l_box/2pi

nhigh = 64 # high n = k*l_box/2pi

rseed = 219 # random seed

f_shear = 1 # solenoidal turbulence

tcorr = 0 # not used for turb flag 1

dtdrive = 0 # not used for turb flag 1



\<problem>

turb_flag = 1 # turb at start 

n_inj = 5 # Radius of injection cells

rho_0 = 2.525309218890173e-22 # Initial background density

v_wind = 2.0e8 # Stellar wind speed

rho_wind = 6.298768747466559e20 # Stellar mass loss rate cgs

xpos = 0 # star x

ypos = 0 # star y

zpos = 0 # star z

[c-white]

The first thing that comes to mind is units. It looks like you're trying to run the code in CGS. While this might work, it can also lead to numerical issues. Codes as large as this are rarely purely scale-free -- at some point a number may be considered "effectively 0" or "effectively infinity" by the algorithm, based on its size relative to 1. Even if not, floating-point numbers have finite range, and some of your values already exceed single-precision.

Given the perfectly scale-free nature of the equations of hydrodynamics, it's always possible (and I usually suggest) to scale your values by some typical normalizations for your problem. For this, you might want to end up with a system where rho_0, v_wind, and x1max are all 1; then tlim would probably also end up closer to 1 as well.

Other than that, you could try turning on different pieces of physics one at a time. It looks like you have turbulent driving, and also some sort of custom mass injection. Can the code run fine with one of them turned off?

[c-white]

Also (thinking out loud here) -- turb_flag=1 is for punching the simulation really hard right at the start. I don't know why the turbulence community decided this was a reasonable thing to do, but here we are. If lowering dedt has a positive effect on robustness, it's possible that using turb_flag=3 and spreading the injection over an interval will allow you to get the same energy injected without the crashing.

And just to mention, when changing the algorithm choices, the simpler, more diffusive, less accurate, lower-order direction will almost always increase robustness. That is, xorder=2 (or even 1), integrator=vl2 (or even rk1), flux=hlle (or even llf) is a good place to seek robustness. If that doesn't work, 4c/rk4/roe will almost certainly fail.

[BraydenAlbery]

Thanks for getting back to me! Yes, turning down the turbulence did help with robustness. I will look into using turb flag 2 or 3 to initialise my background medium. Is there a simple way to turn off turbulence driving after a certain time/energy injection? Or will it be a case of doing something like using an rst file after rebuilding the binary with turbulence commented out in the problem generator?

[c-white]

I don't think anyone ever wrote a parameter for stopping turbulent driving after a certain time, so the easiest way would be a restart. Instead of recompiling, it may be enough to change the input parameter turb_flag to 0.

[BraydenAlbery]

I've been doing tests with turb flag 3, and it appears to be behaving itself! One last question, is there a way to force the output of an rst file? Otherwise, I lose some progress when I have to stop the simulation when the turbulence is finished injecting.

[c-white]

Usually it's best to use the input parameters to control when outputs are made. But if there's no better way, the code checks for SIGTERM, SIGINT, and SIGALRM signals, any one of which should signal it to finish its current timestep, make a restart file, and exit. How to trigger these signals and whether the code ever sees them depends on your operating system.

[ltlancas]

This sounds like a problem I'm familiar with! What I have done in a similar simulation is to drive turbulence with an instantaneous energy injection turb_flag=1 and allow that to evolve over some time set by the <time>/tlim parameter. To be clear on what Chris already implicitly said, to force a restart you need to add an output block with file_type = rst (looks like your <output5> would do) with dt=tlim, though it should output a restart file at the end of the simulation anyway as long as you have an output block with the rst type.

[BraydenAlbery]

My makeshift solution is to inject the turbulence with turb_flag = 3. I have a function in Mesh::UserWorkInLoop that tracks the total kinetic energy, it is only called when a boolean TurbInj = false. Once this reaches my dedt total that I was injecting with turb_flag =1, I kill the simulations and restart with turb_flag = 0 and TurbInj = true. It appears to be working so far!