Skip to content
Snippets Groups Projects

turbulence

Last edited by Philipp Grete

Driven turbulence simulations

The new problem generator fmturb uses explicit inverse Fourier transformations (iFTs) on each meshblock in order to reduce communication during the iFT. Thus, it is only efficient if comparatively few modes are used (say < 100).

Quite generally, driven turbulence simulations start from uniform initial conditions (uniform density and pressure, uniform magnetic field in x-direction in case of an MHD setup, and the fluid at rest) and reach a state of stationary, isotropic (or anisotropic depending on the strength of the background magnetic field) turbulence after one to few large eddy turnover times (again depending on the background magnetic field strength). The large eddy turnover time is usually defined as T = L/U with L being the scale of the largest eddies and U the root mean square Mach number in the stationary regime.

The current implementation uses the following forcing spectrum (k/k_peak)^2 * (2 - (k/k_peak)^2). If a different spectrum is required see Making further changes. Here, k_peak is the peak wavenumber of the forcing spectrum. It is related the scales of the largest eddies as L = 1/k_f given that a box size of 1 is currently assumed/hardcoded.

Compilation

Configure with --prob=fmturb, e.g.,

./configure.py --kokkos_arch="Volta70,SKX" --kokkos_devices="Cuda,OpenMP"  --prob=fmturb  --flux=roe -mpi -hdf5 --hdf5_path=$EBROOTHDF5 -b

for K-Athena with Roe Riemann solver, PLM reconstruction, MHD, MPI, and support for HDF5 output.

Problem setup

An example parameter file can be found in inputs/mhd/athinput.fmturb

Following parameters can be changed to control both the initial state

  • rho0 initial mean density
  • p0 initial mean thermal pressure
  • b0 initial mean magnetic field strength in x-direction

as well as the driving field

  • kpeak peak wavenumber of the forcing spectrum. Make sure to update the wavemodes to match kpeak, see below.
  • corr_time autocorrelation time of the acceleration field (in code units). Using delta-in-time forcing, i.e., a very low value, is discouraged, see Grete et al. 2018 ApJL.
  • rseed random seed for the OU process. Only change for new simulation, but keep unchanged for restarting simulations.
  • sol_weight solenoidal weight of the acceleration field. 1.0 is purely solenoidal/rotational and 0.0 is purely dilatational/compressive. Any value between 0.0 and 1.0 is possible. The parameter is related to the resulting rotational power in the 3D acceleration field as 1. - ((1-sol_weight)^2/(1-2*sol_weight+3*sol_weight^2)), see eq (9) in Federrath et al. 2010 A&A.
  • accel_rms root mean square value of the acceleration (controls the "strength" of the forcing field)
  • num_modes number of wavemodes that are specified in the <modes> section of the parameter file. In order to generate a full set of modes run the inputs/mhd/generate_fmturb_modes.py script and replace the corresponding parts of the parameter file with the output of the script. Within the script, the top three variables (k_peak, k_high, and k_low) need to be adjusted in order to generate a complete set (i.e., all) of wavemodes. Alternatively, wavemodes can be chosen/defined manually, e.g., if not all wavemodes are desired or only individual modes should be forced.

Making further changes

Forcing profile

In order to change the forcing profile edit src/fft/few_modes_turbulence.cpp and update the following line (currently line 307)

tmp = std::pow(kmag/kpeak,2.)*(2.-std::pow(kmag/kpeak,2.));

Similarly, inputs/mhd/generate_fmturb_modes.py should be updated, specifically line

if (k_mag/k_peak)**2.*(2.-(k_mag/k_peak)**2.) < 0:

so that the generated wave modes match.

Box sizes

The problem generator/forcing is currently limited to 3D cubic boxes with side length 1. In case different dimensions and/or aspect ratios are required several parts of the code need to be changed, e.g, in calculating the phases of the iFTs. Please get in contact or open an issue for getting directions.

Typical results

The results shown here are obtained from running simulations with the parameters given in the next section.

High level temporal evolution

image

Power spectra

image

Parameters

Verifying that the corr_time, accel_rms, and sol_weigh parameters work as expected:

T_0.25-A_2.00-S_0.00
Forcing correlation time       target: 0.25 actual: 0.21+/-0.035
RMS acceleration               target: 2.00 actual: 2.00+/-0.000
Rel power of sol. modes        target: 0.00 actual: 0.00+/-0.000

T_0.25-A_2.00-S_0.50
Forcing correlation time       target: 0.25 actual: 0.29+/-0.031
RMS acceleration               target: 2.00 actual: 2.00+/-0.000
Rel power of sol. modes        target: 0.67 actual: 0.68+/-0.042

T_0.25-A_2.00-S_1.00
Forcing correlation time       target: 0.25 actual: 0.30+/-0.053
RMS acceleration               target: 2.00 actual: 2.00+/-0.000
Rel power of sol. modes        target: 1.00 actual: 1.00+/-0.000

T_0.50-A_1.00-S_1.00
Forcing correlation time       target: 0.50 actual: 0.54+/-0.102
RMS acceleration               target: 1.00 actual: 1.00+/-0.000
Rel power of sol. modes        target: 1.00 actual: 1.00+/-0.000

T_1.00-A_0.50-S_1.00
Forcing correlation time       target: 1.00 actual: 1.06+/-0.140
RMS acceleration               target: 0.50 actual: 0.50+/-0.000
Rel power of sol. modes        target: 1.00 actual: 1.00+/-0.000

T_0.2-A_1.25-S_0.30
Forcing correlation time       target: 0.20 actual: 0.22+/-0.045
RMS acceleration               target: 1.25 actual: 1.25+/-0.000
Rel power of sol. modes        target: 0.27 actual: 0.28+/-0.038

T_0.2-A_1.25-S_0.30-NPROC_8
Forcing correlation time       target: 0.20 actual: 0.22+/-0.045
RMS acceleration               target: 1.25 actual: 1.25+/-0.000
Rel power of sol. modes        target: 0.27 actual: 0.28+/-0.038

T_0.2-A_1.25-S_0.30-RST
Forcing correlation time       target: 0.20 actual: 0.22+/-0.044
RMS acceleration               target: 1.25 actual: 1.25+/-0.000
Rel power of sol. modes        target: 0.27 actual: 0.28+/-0.038

T_0.2-A_1.25-S_0.30-MHD
Forcing correlation time       target: 0.20 actual: 0.22+/-0.044
RMS acceleration               target: 1.25 actual: 1.25+/-0.000
Rel power of sol. modes        target: 0.27 actual: 0.28+/-0.038

T_0.2-A_1.25-S_0.30-GCC
Forcing correlation time       target: 0.20 actual: 0.22+/-0.044
RMS acceleration               target: 1.25 actual: 1.25+/-0.000
Rel power of sol. modes        target: 0.27 actual: 0.28+/-0.038

Restarts and MPI

image

GPU and CPU

image

Consistency of acceleration field

As each meshblock does a full iFT of all modes the following slices from a run with 8 meshblocks illustrate that there is no discontinuities at the meshblock boundary.

Plot shows x-, y-, and z-acceleration (in rows top to bottom) slices in the x-, y-, and z-direction (in columns from left to right).

image

Scripts used to create/analyze test data

Running the simulations

#!/bin/bash

KATHENADIR=/mnt/home/gretephi/src/kathena
N=64 # numerical resolution


# first test for Ms approx 0.5 solenoidal
ACCRMS=0.50 # with k=2 -> T = 1
TCORR=1.00
SOLW=1.00
TLIM=10.0 # should be 10 turnover times
DTHST=0.05
DTHDF=0.1

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW
mkdir $DIR
cd $DIR
srun -N 1 -n 1 $KATHENADIR/bin/athena.cuda.hydro -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$N  meshblock/nx2=$N meshblock/nx3=$N problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..


# first test for Ms approx 1.0 solenoidal
ACCRMS=1.00 # with k=2 -> T = 0.5
TCORR=0.50
SOLW=1.00
TLIM=5.0 # should be 10 turnover times
DTHST=0.025
DTHDF=0.05

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW
mkdir $DIR
cd $DIR
srun -N 1 -n 1 $KATHENADIR/bin/athena.cuda.hydro -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$N  meshblock/nx2=$N meshblock/nx3=$N problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..

# first test for Ms approx 1.0 solenoidal
ACCRMS=2.00 # with k=2 -> T = 0.25
TCORR=0.25
SOLW=1.00
TLIM=2.5 # should be 10 turnover times
DTHST=0.0125
DTHDF=0.025

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW
mkdir $DIR
cd $DIR
srun -N 1 -n 1 $KATHENADIR/bin/athena.cuda.hydro -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$N  meshblock/nx2=$N meshblock/nx3=$N problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..



# first test for Ms approx 2 mixed
ACCRMS=2.00 # with k=2 -> T = 0.25
TCORR=0.25
SOLW=0.50
TLIM=2.50 # should be 10 turnover times
DTHST=0.0125
DTHDF=0.025

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW
mkdir $DIR
cd $DIR
srun -N 1 -n 1 $KATHENADIR/bin/athena.cuda.hydro -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$N  meshblock/nx2=$N meshblock/nx3=$N problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..

# first test for Ms approx 2 compressive
ACCRMS=2.00 # with k=2 -> T = 0.25
TCORR=0.25
SOLW=0.00
TLIM=2.50 # should be 10 turnover times
DTHST=0.0125
DTHDF=0.025

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW
mkdir $DIR
cd $DIR
srun -N 1 -n 1 $KATHENADIR/bin/athena.cuda.hydro -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$N  meshblock/nx2=$N meshblock/nx3=$N problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..



# first test for some parameters in between 
ACCRMS=1.25 # with k=2 -> T = 0.4
TCORR=0.2 # -> 0.5T
SOLW=0.30
TLIM=4.0 # should be 10 turnover times
DTHST=0.02
DTHDF=0.04

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW
mkdir $DIR
cd $DIR
srun -N 1 -n 1 $KATHENADIR/bin/athena.cuda.hydro -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$N  meshblock/nx2=$N meshblock/nx3=$N problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..

# same run with more  proc 
ACCRMS=1.25 # with k=2 -> T = 0.4
TCORR=0.2 # -> 0.5T
SOLW=0.30
TLIM=4.0 # should be 10 turnover times
DTHST=0.02
DTHDF=0.04

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW-NPROC_8
mkdir $DIR
cd $DIR
srun -N 1 -n 8 $KATHENADIR/bin/athena.cuda.hydro -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$((N/2))  meshblock/nx2=$((N/2)) meshblock/nx3=$((N/2)) problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..

# parameters in between restarted
ACCRMS=1.25 # with k=2 -> T = 0.4
TCORR=0.2 # -> 0.5T
SOLW=0.30
TLIM=4.0 # should be 10 turnover times
DTHST=0.02
DTHDF=0.04

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW-RST
mkdir $DIR
cd $DIR
cp ../T_$TCORR-A_$ACCRMS-S_$SOLW/Turb.00001.rst .
srun -N 1 -n 1 $KATHENADIR/bin/athena.cuda.hydro -r Turb.00001.rst  |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..

# parameters in between MHD
ACCRMS=1.25 # with k=2 -> T = 0.4
TCORR=0.2 # -> 0.5T
SOLW=0.30
TLIM=4.0 # should be 10 turnover times
DTHST=0.02
DTHDF=0.04

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW-MHD
mkdir $DIR
cd $DIR
srun -N 1 -n 1 $KATHENADIR/bin/athena.cuda.mhd -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$N  meshblock/nx2=$N meshblock/nx3=$N problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..


# parameters in between MHD
ACCRMS=1.25 # with k=2 -> T = 0.4
TCORR=0.2 # -> 0.5T
SOLW=0.30
TLIM=4.0 # should be 10 turnover times
DTHST=0.02
DTHDF=0.04

DIR=T_$TCORR-A_$ACCRMS-S_$SOLW-GCC
mkdir $DIR
cd $DIR
srun -N 1 -n 1 $KATHENADIR/bin/athena.gcc.mhd -i $KATHENADIR/inputs/mhd/athinput.fmturb output1/dt=$DTHST output2/dt=$DTHDF output4/dt=$DTHDF output3/dt=$DTHDF time/tlim=$TLIM mesh/nx1=$N mesh/nx2=$N mesh/nx3=$N meshblock/nx1=$N  meshblock/nx2=$N meshblock/nx3=$N problem/accel_rms=$ACCRMS problem/corr_time=$TCORR problem/sol_weight=$SOLW |tee ath.out
for X in `seq -w 0001 0100`; do srun -n 2 python ~/src/energy-transfer-analysis/run_analysis.py --res $N --data_path 0$X.athdf  --data_type AthenaPP --type flow --eos adiabatic --gamma 1.0001 --outfile flow-$X.hdf5 -forced; done
cd ..

Analyzing the data

Note, that the follow scripts are just snippets that cover some basics. Please get in contact for more complete/feature rich scripts.

Calculating the integral time of the acceleration field:

def get_integral_time(Id):
    ds_arr = yt.load(RootDir +  Id + '/Turb.acc.000*.athdf')
    
    acc_arr = []
    times = []
    for ds in ds_arr:    
        cg = ds.covering_grid(0,ds.domain_left_edge,ds.domain_dimensions)
        acc_arr.append(np.array([
            cg[('athena_pp', 'acceleration_x')],
            cg[('athena_pp', 'acceleration_y')],
            cg[('athena_pp', 'acceleration_z')],
        ]))
        times.append(ds.current_time)
    
    acc_arr = np.array(acc_arr)
    num_snaps = acc_arr.shape[0]
    corr = np.zeros((4,num_snaps,num_snaps))
    for i in range(num_snaps):
        for j in range(num_snaps):
            if j < i:
                continue
            corr[0,i,j-i] = scipy.stats.pearsonr(acc_arr[i,:,:,:,:].reshape(-1),acc_arr[j,:,:,:,:].reshape(-1))[0]
            corr[1,i,j-i] = scipy.stats.pearsonr(acc_arr[i,0,:,:,:].reshape(-1),acc_arr[j,0,:,:,:].reshape(-1))[0]
            corr[2,i,j-i] = scipy.stats.pearsonr(acc_arr[i,1,:,:,:].reshape(-1),acc_arr[j,1,:,:,:].reshape(-1))[0]
            corr[3,i,j-i] = scipy.stats.pearsonr(acc_arr[i,2,:,:,:].reshape(-1),acc_arr[j,2,:,:,:].reshape(-1))[0]
            
    return corr

Check power in acceleration field, solenoidal weight, and correlation time:

for Id in Ids:
    id_split = Id.split('-')
    
    t_corr = float(id_split[0].split('_')[1])
    t_corr_actuals = []
    this_data = autocorr[Id][0] # use the all components of the acc vector
    num_points = this_data.shape[1]
    for i in range(num_points):
        this_slice = this_data[i,:]
        # find the first zero crossing - we only integrate till that point as it's noise afterwards anyway
        # this also ensures that the t_corr from later snapshots is not included as too few snapshots would follow
        # to actually integrate for a full t_corr
        idx_0 = np.argwhere(np.array(this_slice) < 0)
        if len(idx_0) == 0:
            continue
        t_corr_actuals.append(np.trapz(this_slice[:idx_0[0][0]],dx=sim_dict[Id]['code_time_between_dumps']))

    t_corr_actual = (np.mean(t_corr_actuals),np.std(t_corr_actuals))
    
    a_rms = float(id_split[1].split('_')[1])
    a_rms_actual = get_mean('a' + '/moments/' + 'rms')
    
    zeta = float(id_split[2].split('_')[1])
    sol_weight_actual = get_mean_squared_ratio('a_s_mag' + '/moments/' + 'rms','a' + '/moments/' + 'rms')
    sol_weight = 1. - ((1-zeta)**2/(1-2*zeta+3*zeta**2)) # see eq 9 of Federrath et al 2010 http://dx.doi.org/10.1051/0004-6361/200912437
    
    msg = (
        f"\n{Id}\n"
        f"{'Forcing correlation time':30} target: {t_corr:.2f} actual: {t_corr_actual[0]:.2f}+/-{t_corr_actual[1]:.3f}\n"
        f"{'RMS acceleration':30} target: {a_rms:.2f} actual: {a_rms_actual[0]:.2f}+/-{a_rms_actual[1]:.3f}\n"
        f"{'Rel power of sol. modes':30} target: {sol_weight:.2f}"
        f" actual: {sol_weight_actual[0]:.2f}+/-{sol_weight_actual[1]:.3f}"
    )
    print(msg)

Comments

    Please register or sign in to add a comment.