Turbulence simulation runtime issue

[YaoguangPei]

When running turbulence simulations, I have encountered persistent issues with the time step (dt) inevitably dropping to around 1e-9, and in some cases even reaching 1e-10, across multiple test runs. I have attempted several approaches to address this:

1）Lowering the resolution, from 1024 to 512, and then to 128. However, this showed little to no significant improvement.
2）Applying pencil decomposition within each MeshBlock. This yielded some improvement, but dt still dropped to 1e-7 by simulation time t = 1.2, causing the run to stall.
3）Enabling super-time-stepping (STS). Unfortunately, this worsened the situation, and dt dropped to 1e-10 even earlier in the run.
4）Reducing the turbulence speed, which I control via the dedt parameter. Even after cutting the driving strength to one-fourth of the original, the problem persisted.

Therefore, I would like to ask: Are there any more effective strategies to improve the simulation runtime?
For reference, I have included below my simulation configuration and SLURM job script.

athinput.turb:

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 = 0.1 # time increment between outputs file_type = rst dt = 0.1 cfl_number = 0.3 # The Courant, Friedrichs, & Lewy (CFL) Number nlim =-1 # cycle limit tlim = 10 # time limit integrator = vl2 # time integration algorithm xorder = 2 # order of spatial reconstruction ncycle_out = 1 # interval for stdout summary info sts_integrator=rkl2 #If not enabled, the following three lines will not appear. dt_diagnostics = -1 sts_max_dt_ratio=100 nx1 = 128 # 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 = 128 # 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 = 128 # 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 =128 nx2 =16 nx3 =16 x1min = -2 x1max = 2 x2min = -2 x2max = 2 x3min = -2 x3max = 2 level = 1

iso_sound_speed = 0.19 # equivalent to sqrt(gamma*p/d) for p=0.1, d=1
dfloor=28.5

turb_flag = 3 # 1 for decaying, 2 (impulsive) or 3 (continuous) for driven turbulence dedt = 914.28# Energy injection rate (for driven) or Total energy (for decaying) nlow = 0 # cut-off wavenumber at low-k nhigh =128 # cut-off wavenumber at high-k expo = 5/3 # power-law exponent tcorr = 0.1 # 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 four_pi_G=1 nu_iso = 1e-3

Slurm:
#!/bin/bash
#SBATCH -J iso2
#SBATCH --nodes=2
#SBATCH --ntasks-per-node=32
#SBATCH --mem=200G
#SBATCH --cpus-per-task=1
#SBATCH -p GPUC
#SBATCH -t 240:00:00

cd ~/athena
module purge
module load compiler/intel/2021.3.0
module load mpi/intelmpi/2021.3.0
module load mathlib/fftw/intelmpi/3.3.9_double
module load mathlib/hdf5/intel/1.12.0

export FFTW_INC=/public/software/mathlib/fftw/intelmpi/3.3.9_double/include
export FFTW_LIB=/public/software/mathlib/fftw/intelmpi/3.3.9_double/lib
export HDF5_INC=/public/software/mathlib/hdf5/intel/1.12.0/include
export HDF5_LIB=/public/software/mathlib/hdf5/intel/1.12.0/lib

export CPATH=$FFTW_INC:$HDF5_INC:$CPATH
export LIBRARY_PATH=$FFTW_LIB:$HDF5_LIB:$LIBRARY_PATH
export LD_LIBRARY_PATH=$FFTW_LIB:$HDF5_LIB:$LD_LIBRARY_PATH

python3 configure.py --prob turb -fft --fftw_path=$FFTW_DIR --flux roe --eos isothermal --grav fft -mpi -sts##I have enabled super-time-stepping (STS) in this case

sed -i -E
-e "s|^(CXX\s*:=\s*).|\1mpiicpc|"
-e "s|^(CXXFLAGS\s:=\s*).|\1-O3 -std=c++11 -I${FFTW_INC} -D__INTEL_COMPILER -fopenmp|"
-e "s|^(LDFLAGS\s:=\s*).|\1-L${FFTW_LIB} -fopenmp| "
-e "s|^(LDLIBS\s:=\s*).*|\1-lfftw3 -lfftw3_mpi |"
Makefile
make clean
make
cd ~/first/iso2
mpirun -np 64 ~/athena/bin/athena -i athinput.turb

[tianninglyu]

Hi Yaoguang,

In turbulent self-gravity simulations, the standard practice is to first evolve the system for several crossing times with turbulence only (without self-gravity), so that the turbulence can fully develop and reach a statistical steady state. In your setup, it seems that self-gravity is turned on from the very beginning, which tends to produce overly dense structures prematurely, leading to very small timesteps. This approach can also easily violate the Truelove criterion (Truelove et al. 1997, i.e., requiring at least 4 cells per Jeans length), which may cause unphysical fragmentation. You can refer to Brucy et al. 2024 for the general procedure of turbulence simulation.

In addition, your dedt value looks too large. By default, the turbulence problem generator in Athena++ initializes with a uniform density of 1. From my numerical experiments, for a box with density = 1, box size = 1, and isothermal sound speed = 0.1, driving turbulence to Mach 10 requires dedt ≈ 1.5. I don't know the background of your research, but generally, such simulations inject energy at the largest scale, for example, k=[1,3], and your k=[0,128] does make the dedt larger. But I still don't think dedt>900 is a reasonable value. I suggest you to estimate a reasonable dedt value using Equation(7) in Mac Low (1999).

Could you share your hst file? That would greatly help in diagnosing the issue. When sharing, please enclose the file content with triple backticks (```), which will make the format much clearer.

Best regards,
Tianning

[YaoguangPei]

Dear Tianning,
Thank you very much for your response; it has been extremely helpful to me. I have identified the source of the issue: the drop in dedt was caused by the parameter four_pi_G, which I mistakenly set about 1000 times larger than it should be. As a result, in the later stages of the simulation, gravity contributed an excessive amount of energy to the kinetic energy, making the CFL condition particularly restrictive.
In addition, I have not yet tried the strategy of disabling self-gravity in the initial phase, but I plan to test this approach. Thank you very much for suggesting it.
Your advice regarding the dedt setting has also been very helpful. In my current simulation with density = 40, box_size = 4, and isothermal_speed = 0.19, I set dedt = 30. I also adjusted k to the range [1, 3].
Thank you again for your help
Best regards,
Yaoguang Pei