1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
import matplotlib.pyplot as plt
import numpy as np
import yt
"""
Make a turbulent KE power spectrum. Since we are stratified, we use
a rho**(1/3) scaling to the velocity to get something that would
look Kolmogorov (if the turbulence were fully developed).
Ultimately, we aim to compute:
1 ^ ^*
E(k) = integral - V(k) . V(k) dS
2
n ^
where V = rho U is the density-weighted velocity field, and V is the
FFT of V.
(Note: sometimes we normalize by 1/volume to get a spectral
energy density spectrum).
"""
def doit(ds):
# a FFT operates on uniformly gridded data. We'll use the yt
# covering grid for this.
max_level = ds.index.max_level
ref = int(np.prod(ds.ref_factors[0:max_level]))
low = ds.domain_left_edge
dims = ds.domain_dimensions * ref
nx, ny, nz = dims
nindex_rho = 1.0 / 3.0
Kk = np.zeros((nx // 2 + 1, ny // 2 + 1, nz // 2 + 1))
for vel in [("gas", "velocity_x"), ("gas", "velocity_y"), ("gas", "velocity_z")]:
Kk += 0.5 * fft_comp(
ds, ("gas", "density"), vel, nindex_rho, max_level, low, dims
)
# wavenumbers
L = (ds.domain_right_edge - ds.domain_left_edge).d
kx = np.fft.rfftfreq(nx) * nx / L[0]
ky = np.fft.rfftfreq(ny) * ny / L[1]
kz = np.fft.rfftfreq(nz) * nz / L[2]
# physical limits to the wavenumbers
kmin = np.min(1.0 / L)
kmax = np.min(0.5 * dims / L)
kbins = np.arange(kmin, kmax, kmin)
N = len(kbins)
# bin the Fourier KE into radial kbins
kx3d, ky3d, kz3d = np.meshgrid(kx, ky, kz, indexing="ij")
k = np.sqrt(kx3d**2 + ky3d**2 + kz3d**2)
whichbin = np.digitize(k.flat, kbins)
ncount = np.bincount(whichbin)
E_spectrum = np.zeros(len(ncount) - 1)
for n in range(1, len(ncount)):
E_spectrum[n - 1] = np.sum(Kk.flat[whichbin == n])
k = 0.5 * (kbins[0 : N - 1] + kbins[1:N])
E_spectrum = E_spectrum[1:N]
index = np.argmax(E_spectrum)
kmax = k[index]
Emax = E_spectrum[index]
plt.loglog(k, E_spectrum)
plt.loglog(k, Emax * (k / kmax) ** (-5.0 / 3.0), ls=":", color="0.5")
plt.xlabel(r"$k$")
plt.ylabel(r"$E(k)dk$")
plt.savefig("spectrum.png")
def fft_comp(ds, irho, iu, nindex_rho, level, low, delta):
cube = ds.covering_grid(level, left_edge=low, dims=delta, fields=[irho, iu])
rho = cube[irho].d
u = cube[iu].d
nx, ny, nz = rho.shape
# do the FFTs -- note that since our data is real, there will be
# too much information here. fftn puts the positive freq terms in
# the first half of the axes -- that's what we keep. Our
# normalization has an '8' to account for this clipping to one
# octant.
ru = np.fft.fftn(rho**nindex_rho * u)[
0 : nx // 2 + 1, 0 : ny // 2 + 1, 0 : nz // 2 + 1
]
ru = 8.0 * ru / (nx * ny * nz)
return np.abs(ru) ** 2
ds = yt.load("maestro_xrb_lores_23437")
doit(ds)
|
