Simulation stuck at t=0 when running turbulence with chemistry (gravity disabled)

[YaoguangPei]

Hi all,

I’m running a turbulence simulation with the chemistry module enabled. To avoid extremely high densities that arise from initial fluctuations under self-gravity, I decided to start the run without gravity at the beginning.

However, I seem to have encountered a problem: the simulation gets stuck right after the output step. It only produces results at t = 0, and even after waiting for 5 days, there are no further outputs. The terminal shows the following message:

Setup complete, entering main loop...

cycle=0 time=0.0000000000000000e+00 dt=9.3749999999967193e-03
dt_hyperbolic=9.3749999999967193e-03 ratio=1.000e+00
dt_parabolic=1.7976931348623158e+307 ratio=5.215e-310

Could anyone please help me understand why this might happen and how to avoid it?
Below are my input file and SLURM script for reference.

Thanks in advance for your help!

Best,
Yaoguang Pei
'''

problem = Turbulence with power-law PS
reference =
configure = --prob=turb -fft

problem_id = Turb # problem ID: basename of output filenames file_type = hst # History data dump dt = 0.01 # time increment between outputs file_type = vtk # Binary data dump variable = prim,p,E # variables to be output dt = 2 # time increment between outputs cfl_number = 0.3 # The Courant, Friedrichs, & Lewy (CFL) Number nlim = -1 # cycle limit tlim =40 # time limit integrator = vl2 # time integration algorithm xorder = 2 # order of spatial reconstruction ncycle_out = 1 # interval for stdout summary info nx1 =256 # Number of zones in X1-direction x1min = -2 # minimum value of X1 x1max = 2 # maximum value of X1 ix1_bc = periodic # inner-X1 boundary flag ox1_bc = periodic # outer-X1 boundary flag

nx2 = 256 # Number of zones in X2-direction
x2min = -2 # minimum value of X2
x2max = 2 # maximum value of X2
ix2_bc =periodic # inner-X2 boundary flag
ox2_bc =periodic # outer-X2 boundary flag

nx3 = 256 # Number of zones in X3-direction
x3min = -2 # minimum value of X3
x3max = 2 # maximum value of X3
ix3_bc = periodic # inner-X3 boundary flag
ox3_bc =periodic # outer-X3 boundary flag

refinement = none

nx1 =256 nx2 =32 nx3 =32 gamma = 1.666666666667 # gamma = C_p/C_v sfloor = 0 turb_flag = 3 # 1 for decaying, 2 (impulsive) or 3 (continuous) for driven turbulence dedt = 2.4e2 # Energy injection rate (for driven) or Total energy (for decaying) nlow = 1 # cut-off wavenumber at low-k nhigh = 2 # cut-off wavenumber at high-k expo = 1.667 # power-law exponent tcorr = 0.5 # correlation time for OU process (both impulsive and continuous) dtdrive = 0.1 # time interval between perturbation (impulsive) f_shear = 0.5 # the ratio of the shear component rseed = -1 # if non-negative, seed will be set by hand (slow PS generation) unit_system=ism # chemistry network parameters Zdg = 1 # metallicity, dust and gas metallicity are the same. default = 1. # chemistry solver parameters reltol = 1.0e-3 # relative tolerance, default 1.0e-2 abstol = 1.0e-12 # absolute tolerance, default 1.0e-12 abstol_H2 = 1.0e-10 abstol_E = 1e-4 # for internal energy. Set to 0.01 K user_jac = false # flag for whether use user provided Jacobian. default false maxsteps = 10000 # maximum number of steps in one integration. default 1000 temp_max_heat = 1.0e5 temp_max_rates = 1.0e5 xHe = 0.1 # He per H, default = 0.1 xC = 1.6e-4 # C atom per H, this is Z=1 value, default 1.6e-4 xO = 3.2e-4 # O atom per H, this is Z=1 value, default 3.2e-4 xS = 3.5e-6 # S atom per H, this is Z=1 value, default 3.5e-6 xSi = 1.7e-6 # Si atom per H, this is Z=1 value, default 1.7e-6 maxorder = 2 # maximum order. Default 2. stldet = 0 # stability limit detection. Default 0/false.

<chem_radiation>
G0 = 1
CR = 2e-16 # Cosmic-ray ionization rate per H
n_frequency = 8
output_zone_sec = true # output diagnostic

r_init = 1.0e-6 # initial abundances
"""

"""
#!/bin/bash
#SBATCH -J bla_4pc_256_pre
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=64
#SBATCH --mem=500G
#SBATCH --cpus-per-task=1
#SBATCH -p GPUC
#SBATCH -t 240:00:00

module purge
module load compiler/gcc/12.3.0
module load mpi/openmpi/gnu/4.0.3
module load mathlib/hdf5/gnu/1.12.0
module load mathlib/netcdf/gnu/4.7.4
module load compiler/cmake/3.20.1

export FFTW_DIR=$HOME/fftw-gnu-install
export FFTW_INC=$FFTW_DIR/include
export FFTW_LIB=$FFTW_DIR/lib
export CPATH=$FFTW_INC:$CPATH
export LIBRARY_PATH=$FFTW_LIB:$LIBRARY_PATH
export LD_LIBRARY_PATH=$FFTW_LIB:$LD_LIBRARY_PATH

export SUNDIALS_DIR=$HOME/sundials-gnu-install
export CPATH=$SUNDIALS_DIR/include:$CPATH
export LIBRARY_PATH=$SUNDIALS_DIR/lib64:$LIBRARY_PATH
export LD_LIBRARY_PATH=$SUNDIALS_DIR/lib64:$LD_LIBRARY_PATH

cd ~/athena2

python3 configure.py --prob turb
-fft --fftw_path=$FFTW_DIR
--flux hllc
-mpi
--chemistry gow17
--chem_radiation six_ray
--chem_ode_solver cvode
--cvode_path=$SUNDIALS_DIR
--cflag="-std=c++17"
--cxx="g++"

make clean
make

cd ~/first/bla_4pc_256
export OMPI_MCA_pml=ob1
export OMPI_MCA_btl=tcp,self

srun --mpi=pmix_v3 ~/athena2/bin/athena -i athinput.turb_pre

"""

[tomidakn]

Does it run without chemistry, or with different chem_radiation option?

[YaoguangPei]

When I turn off the chemistry and six_ray chem_radiation, the sim won't stop running. And I have checked the usage rate of CPU when I turn on the chemistry and six_ray chem_radiation, and find CPU is running at full capacity.

[tomidakn]

OK, can you try with chemistry but --chem_radiation=const?

[YaoguangPei]

I have tried the chem_radiation=const, and I find it can keep running.

[YaoguangPei]

Dear Professor,
Could it be that the computational cost of rix_ray is too high, causing my simulation to stop running? If so, should I simply wait for the next output?

[tomidakn]

I rather worry that there is a bug in the implementation of the six-ray module.
@munan , is it possible for you to look into this?

[munan]

Hi Kengo, I am happy to take a look. First, did you pull the newest version of Athena? I updated something about six-ray boundary condition in periodic hydro boundary, and Kyle just merged it in the past week or so. If you still have a problem, please send me your config and input file, and I will try to debug it. Thank you! Munan

…

On Tue, Oct 21, 2025 at 11:27 PM Kengo TOMIDA ***@***.***> wrote: I rather worry that there is a bug in the implementation of the six-ray module. @munan <https://github.com/munan> , is it possible for you to look into this? — Reply to this email directly, view it on GitHub <#697 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAVXJH3OO6IOVLRYXANCPT33Y4IURAVCNFSM6AAAAACJZM4JPSVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTINZUG4YDIMA> . You are receiving this because you were mentioned.Message ID: ***@***.*** com>

[YaoguangPei]

Dear Professor,
It turns out that the issue was on my side — I hadn’t pulled the latest version of Athena. As a result, the periodic boundary condition didn’t work properly. After updating to the newest version, my simulation started running correctly and hasn’t stopped since.
I sincerely appreciate both of your kindness and help.
Best regards,
Yaoguang Pei

[tomidakn]

Great, and thanks @munan for your quick help!