[yt-svn] commit/yt: 4 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Sep 2 16:11:37 PDT 2014
4 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/6a9ad8480134/
Changeset: 6a9ad8480134
Branch: yt
User: mzingale
Date: 2014-09-02 20:18:50
Summary: add a recipe showing how to do a power spectrum
Affected #: 2 files
diff -r 2c3d2e84074ac8260feb8f056b0720c3a7442e82 -r 6a9ad84801349749ce49e43f6465501e3baf5b2b doc/source/cookbook/calculating_information.rst
--- a/doc/source/cookbook/calculating_information.rst
+++ b/doc/source/cookbook/calculating_information.rst
@@ -90,3 +90,14 @@
See :ref:`filtering-particles` for more information.
.. yt_cookbook:: particle_filter_sfr.py
+
+Making a Turbulent Kinetic Energy Power Spectrum
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This recipe shows how to use `yt` to read data and put it on a uniform
+grid to interface with the NumPy FFT routines and create a turbulent
+kinetic energy power spectrum. (Note: the dataset use here is sort
+of low resolution, so the turbulence is not very well-developed. The
+spike at high wavenumbers is due to non-periodicity in the z-direction).
+
+.. yt_cookbook:: power_spectrum_example.py
diff -r 2c3d2e84074ac8260feb8f056b0720c3a7442e82 -r 6a9ad84801349749ce49e43f6465501e3baf5b2b doc/source/cookbook/power_spectrum_example.py
--- /dev/null
+++ b/doc/source/cookbook/power_spectrum_example.py
@@ -0,0 +1,117 @@
+import numpy as np
+import matplotlib.pyplot as plt
+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.
+
+(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.product(ds.ref_factors[0:max_level]))
+
+ low = ds.domain_left_edge
+ dims = ds.domain_dimensions*ref
+
+ nx, ny, nz = dims
+
+ nindex_rho = 1./3.
+
+ 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.max(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./3.), 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
+
+
+if __name__ == "__main__":
+
+ ds = yt.load("maestro_xrb_lores_23437")
+ doit(ds)
https://bitbucket.org/yt_analysis/yt/commits/c98fe077b293/
Changeset: c98fe077b293
Branch: yt
User: mzingale
Date: 2014-09-02 20:54:44
Summary: fix typo
Affected #: 1 file
diff -r 6a9ad84801349749ce49e43f6465501e3baf5b2b -r c98fe077b29375a79634e7d11feb6e5455246ec7 doc/source/cookbook/calculating_information.rst
--- a/doc/source/cookbook/calculating_information.rst
+++ b/doc/source/cookbook/calculating_information.rst
@@ -96,8 +96,8 @@
This recipe shows how to use `yt` to read data and put it on a uniform
grid to interface with the NumPy FFT routines and create a turbulent
-kinetic energy power spectrum. (Note: the dataset use here is sort
-of low resolution, so the turbulence is not very well-developed. The
-spike at high wavenumbers is due to non-periodicity in the z-direction).
+kinetic energy power spectrum. (Note: the dataset used here is of low
+resolution, so the turbulence is not very well-developed. The spike
+at high wavenumbers is due to non-periodicity in the z-direction).
.. yt_cookbook:: power_spectrum_example.py
https://bitbucket.org/yt_analysis/yt/commits/f6f3c6119843/
Changeset: f6f3c6119843
Branch: yt
User: mzingale
Date: 2014-09-02 20:59:31
Summary: add a note explaining \hat{V}
Affected #: 1 file
diff -r c98fe077b29375a79634e7d11feb6e5455246ec7 -r f6f3c6119843f16d7a6b3be639de1209f6f464bc doc/source/cookbook/power_spectrum_example.py
--- a/doc/source/cookbook/power_spectrum_example.py
+++ b/doc/source/cookbook/power_spectrum_example.py
@@ -13,8 +13,9 @@
E(k) = integral - V(k) . V(k) dS
2
- n
-where V = rho U is the density-weighted velocity field.
+ 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).
https://bitbucket.org/yt_analysis/yt/commits/35f68c6d1808/
Changeset: 35f68c6d1808
Branch: yt
User: ngoldbaum
Date: 2014-09-03 01:11:30
Summary: Merged in mzingale/yt-new (pull request #1181)
add a recipe showing how to do a power spectrum
Affected #: 2 files
diff -r 99232fd7d61dac91f4d2e5bbc368559be178883e -r 35f68c6d1808a977a5bfbb996b106bc3f074d990 doc/source/cookbook/calculating_information.rst
--- a/doc/source/cookbook/calculating_information.rst
+++ b/doc/source/cookbook/calculating_information.rst
@@ -90,3 +90,14 @@
See :ref:`filtering-particles` for more information.
.. yt_cookbook:: particle_filter_sfr.py
+
+Making a Turbulent Kinetic Energy Power Spectrum
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This recipe shows how to use `yt` to read data and put it on a uniform
+grid to interface with the NumPy FFT routines and create a turbulent
+kinetic energy power spectrum. (Note: the dataset used here is of low
+resolution, so the turbulence is not very well-developed. The spike
+at high wavenumbers is due to non-periodicity in the z-direction).
+
+.. yt_cookbook:: power_spectrum_example.py
diff -r 99232fd7d61dac91f4d2e5bbc368559be178883e -r 35f68c6d1808a977a5bfbb996b106bc3f074d990 doc/source/cookbook/power_spectrum_example.py
--- /dev/null
+++ b/doc/source/cookbook/power_spectrum_example.py
@@ -0,0 +1,118 @@
+import numpy as np
+import matplotlib.pyplot as plt
+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.product(ds.ref_factors[0:max_level]))
+
+ low = ds.domain_left_edge
+ dims = ds.domain_dimensions*ref
+
+ nx, ny, nz = dims
+
+ nindex_rho = 1./3.
+
+ 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.max(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./3.), 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
+
+
+if __name__ == "__main__":
+
+ ds = yt.load("maestro_xrb_lores_23437")
+ doit(ds)
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list