velocity¶
This is a module designed to process a 2D/3D+1 velocity field with focus on turbulence research.
Basic features: - energy, enstrophy, vorticity fields - energy spectra - n-th order structure functions - two-point velocity autocorrelations - dissipation rate estimation - length scales: integral scales, Taylor micro scales, Kolmogorov scales - Reynolds numbers: turbulent Re, Taylor Re
Features to be added in sight: - probabilistic moments (skewness, kurtosis, …) - triple-point velocity autocorrelations
Philosophy: Prepare a velocity field array “udata”. Then, pass it to a method in this module to obtain any quantities related to turbulence. It should require a single line to obtain the desired quantity from a velocity field unless intermediate steps to obtain the quantity are computationally expensive. The primary example is an autocorrelation function which is used for various quantities like Taylor microscale.
udata = (ux, uy, uz) or (ux, uy) each ui has a shape (height, width, (depth), duration)
If ui-s are individually given, make udata like udata = np.stack((ux, uy))
Overview¶
|
Assumes udata has a shape (d, nrows, ncols, duration) or (d, nrows, ncols) … |
|
Decompose a duidxj tensor into a symmetric and an antisymmetric parts Returns symmetric part (eij) and anti-symmetric part (gij) |
|
Apply the Reynolds decomposition to a velocity field Returns a mean field (time-averaged) and a fluctuating field |
|
Computes divergence of a velocity field |
|
Computes curl of a velocity field using a rate of strain tensor … |
|
Calculate curl of 2D (or 2D+1) field … |
|
Returns energy(vec(x y, z, t) field of udata … |
|
Returns enstrophy(vec(x y, z, t) field of udata … |
|
Returns a time-averaged-energy field … |
|
Returns a time-averaged-enstrophy field … |
|
Return energy averaged over space |
|
Return enstrophy averaged over space |
Turbulence intensity is defined as u/U where u = sqrt((ux**2 + uy**2 + uz**2)/3) # characteristic turbulent velocity U = sqrt((Ux**2 + Uy**2 + Uz**2)) # mean flow speed |
|
|
|
|
Returns nd energy spectrum from velocity data (FFT of a velocity field) |
|
Returns 1D energy spectrum from velocity field data |
|
Returns 1D energy spectrum from velocity field data |
|
Returns dissipation spectrum D(k) = 2 nu k^2 E(k) where E(k) is the 1d energy spectrum |
|
Returns SCALED energy spectrum E(k), its error, and wavenumber. |
|
Returns SCALED 1D energy spectra (E11 and E22) … |
|
Return rescaled dissipation spectra D(k)/(u_eta^3) vs k * eta … |
|
Scales raw energy spectrum by given dissipation rate and viscosity |
Returns a velocity field which satisfies k = sqrt(kx^2 + ky^2) < kmax in the original vel. |
|
|
Returns the dissipation rate computed using the rate-of-strain tensor |
|
Return epsilon computed by isotropic formula involving Taylor microscale |
|
Returns dissipation rate computed by integrated the dissipation specrtrum … |
|
Returns the values of estimated dissipation rate using a long. |
|
Compute spatial autocorrelation function of 2+1 velocity field |
|
Compute spatial autocorrelation function of 2+1 velocity field |
|
Returns two-point velocity autocorrelation tensor, and autocorrelation functions. |
|
Return interpolated functions using the outputs of get_two_point_vel_corr_iso() … |
|
Returns lists of INTERPOLATED autocorrelation function values … |
|
Returns autocorrelation tensor with isotropy assumption … |
|
Structure tensor Dij is essentially the covariance of the two-point velocity difference There is one-to-one correspondence between Dij and Rij. |
|
A method to compute a generalized structure function |
|
Returns the scaled structure functions when raw structure function data are given … |
|
remove nans or infs for a pair of 1D arrays, and returns the compressed arrays with the same length e.g.: arr1 = [0, 0.1, 0.2, np.nan, 0.4, 0.5], arr2 = [-1.2, 19.2. |
|
Returns Taylor microscales as the curvature of the autocorrelation functions at r=0 … |
|
Return Taylor microscales computed by isotropic formulae: lambda_g_iso = (15 nu * u_irms^2 / epsilon) ^ 0.5 |
|
Returns integral scales computed by using autocorrelation functions |
|
Use autocorrelation tensor, Rij to calculate integral length scale Parameters ———- udata Rij: numpy array two-point velocity autocorrelation tensor … |
|
Integral scale defined through energy spectrum. |
|
dissipation rate (epsilon) is known to be proportional to u’^3 / L. |
|
Return integral velocity scale which is identical to u’ (characteristic velocity) See get_characteristic_velocity() |
|
Returns integral scales (related to LARGE EDDIES) |
|
Returns Taylor microscales Parameters ———- udata: numpy array velocity field array r_long: numpy array … |
|
Returns Kolmogorov scales |
|
Returns turbulence reynolds number (Pope 6.59) … |
|
Returns Taylor reynolds number (Pope 6.63) from autocorrelation functions |
|
Reutrns a 2D velocity field with a single Rankine vortex at (x0, y0) |
|
Reutrns a 3D velocity field with a Rankine vortex filament at (x0, y0, z) |
|
Returns udata=(ux, uy, uz) of a slice of isotropic, homogeneous turbulence data (DNS, JHTD) |
Returns values to plot rescaled energy spectrum from Saddoughi (1992) |
|
Returns values to plot energy spectrum from JHTD, computed by Takumi in 2019 Call get_rescaled_energy_spectra_jhtd for scaled energy spectrum. |
|
Returns values to plot rescaled energy spectrum from JHTD, computed by Takumi in 2019 Call get_energy_spectra_jhtd for raw energy spectrum. |
|
Returns the values of rescaled structure function reported in Saddoughi and Veeravalli 1994 paper: r_scaled, dll … |
|
|
General method to get a window with shape (xx.shape[:], duration) or (xx.shape[:]) if duration is None … |
|
Returns the inverse of the signal-loss factor by the window Signal loss factor: 1 (rectangle… |
|
Returns the Hamming window as a function of radius |
|
ONLY WORKS FOR THE 2D data Conducts a cheap bilinear interpolation for missing data. |
|
Returns a mask (N-dim boolean array). |
|
Returns an array whose elements are replaced by fill_value if its mask value is True |
|
Conducts 2d interpolation for udata using scipy.interpolate.griddata … |
|
It is better to always have udata with shape (height, width, depth, duration) (3D) or (height, width, duration) (2D) This method fixes the shape of udata such that its shape is (height, width, depth) or (height, width) |
|
Returns the Hamming window as a function of radius |
|
Returns a equally spaced grid to plot udata |
|
Returns a equally spaced grid to plot FFT of udata |
|
Customizable Kolmogorov Energy spectrum Returns the value(s) of k0 * k^{-5/3} |
|
Universal Kolmogorov Energy spectrum Returns the value(s) of C epsilon^{2/3} k^{-5/3} |
Return Kolmogorov length scale for a given set of dissipation rate and viscosity Parameters ———- epsilon: float, dissipation rate nu: float, viscosity |
|
|
Return 1-component RMS velocity, u’ energy = dim / 2 * u’^2 |
|
Klonecker Delta function: delta_{i, j} |
|
Transformation: Cartesian coord to polar coord |
|
Transformation: cartesian to spherical |
|
natural-sorts elements in a given array alist.sort(key=natural_keys) sorts in human order http://nedbatchelder.com/blog/200712/human_sorting.html (See Toothy’s implementation in the comments) |
Classes and Functions¶
-
ilpm.velocity.get_duidxj_tensor(udata, dx=1.0, dy=1.0, dz=1.0)[source]¶ Assumes udata has a shape (d, nrows, ncols, duration) or (d, nrows, ncols) … one can easily make udata by np.stack((ux, uy))
- Parameters
- udata: numpy array with shape (ux, uy) or (ux, uy, uz)
… assumes ux/uy/uz has a shape (nrows, ncols, duration) or (nrows, ncols, nstacks, duration) … can handle udata without temporal axis
- Returns
- sij: numpy array with shape (nrows, ncols, duration, 2, 2) (dim=2) or (nrows, ncols, nstacks, duration, 3, 3) (dim=3)
- … idea is… sij[spacial coordinates, time, tensor indices]
e.g.- sij(x, y, t) = sij[y, x, t, i, j]
… sij = d ui / dxj
-
ilpm.velocity.decompose_duidxj(sij)[source]¶ Decompose a duidxj tensor into a symmetric and an antisymmetric parts Returns symmetric part (eij) and anti-symmetric part (gij)
- Parameters
- sij, 5d or 6d numpy array (x, y, t, i, j) or (x, y, z, t, i, j)
- Returns
- eij: 5d or 6d numpy array, symmetric part of rate-of-strain tensor.
5d if spatial dimensions are x and y. 6d if spatial dimensions are x, y, and z.
- gij: 5d or 6d numpy array, anti-symmetric part of rate-of-strain tensor.
5d if spatial dimensions are x and y. 6d if spatial dimensions are x, y, and z.
-
ilpm.velocity.reynolds_decomposition(udata)[source]¶ Apply the Reynolds decomposition to a velocity field Returns a mean field (time-averaged) and a fluctuating field
- Parameters
- udata: numpy array
… (ux, uy, uz) or (ux, uy) … ui has a shape (height, width, depth, duration) or (height, width, depth) (3D) … ui may have a shape (height, width, duration) or (height, width) (2D)
- Returns
- u_mean: nd array, mean velocity field
- u_turb: nd array, turbulent velocity field
-
ilpm.velocity.div(udata)[source]¶ Computes divergence of a velocity field
- Parameters
- udata: numpy array
… (ux, uy, uz) or (ux, uy) … ui has a shape (height, width, depth, duration) or (height, width, depth) (3D) … ui may have a shape (height, width, duration) or (height, width) (2D)
- Returns
- div_u: numpy array
… div_u has a shape (height, width, depth, duration) (3D) or (height, width, duration) (2D)
-
ilpm.velocity.curl(udata, dx=1.0, dy=1.0, dz=1.0)[source]¶ Computes curl of a velocity field using a rate of strain tensor … For dim=3, the sign might need to be flipped… not tested … if you already have velocity data as ux = array with shape (m, n) and uy = array with shape (m, n),
udata = np.stack((ugrid1, vgrid1)) omega = vec.curl(udata)
udata: (ux, uy, uz) or (ux, uy)
- Returns
- omega: numpy array
shape: (height, width, duration) (2D) or (height, width, duration) (2D)
-
ilpm.velocity.curl_2d(ux, uy, dx=1.0, dy=1.0)[source]¶ Calculate curl of 2D (or 2D+1) field … A simple method but
- Parameters
- ux: 2D array
x component of a 2D field
- uy: 2D array
y component of a 2D field
- dx: float
data spacing (mm/px)
- dy: float
data spacing (mm/px)
- Returns
- omega: 2D numpy array
vorticity field
-
ilpm.velocity.get_energy(udata)[source]¶ Returns energy(vec(x y, z, t) field of udata … Assumes udata is equally spaced data.
- Parameters
- udata: nd array
- dx: float
data spacing along z (mm/px)
- dy: float
data spacing along z (mm/px)
- dz: float
data spacing along z (mm/px)
- Returns
- energy: nd array
energy
-
ilpm.velocity.get_enstrophy(udata, dx=1.0, dy=1.0, dz=1.0)[source]¶ Returns enstrophy(vec(x y, z, t) field of udata … Assumes udata is equally spaced data.
- Parameters
- udata: nd array
- dx: float
data spacing along z (mm/px)
- dy: float
data spacing along z (mm/px)
- dz: float
data spacing along z (mm/px)
- Returns
- enstrophy: nd array
enstrophy
-
ilpm.velocity.get_time_avg_energy(udata)[source]¶ Returns a time-averaged-energy field … NOT MEAN FLOW ENERGY
- Parameters
- udata: nd array
- Returns
- energy_avg:
time-averaged energy field
-
ilpm.velocity.get_time_avg_enstrophy(udata, dx=1.0, dy=1.0, dz=1.0)[source]¶ Returns a time-averaged-enstrophy field … NOT MEAN FLOW ENSTROPHY
- Parameters
- udata: nd array
- Returns
- enstrophy_avg: nd array
time-averaged enstrophy field
-
ilpm.velocity.get_spatial_avg_energy(udata, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None)[source]¶ Return energy averaged over space
- Parameters
- udata: nd array
- Returns
- energy_vs_t: 1d numpy array
average energy in a field at each time
- energy_vs_t_err
standard deviation of energy in a field at each time
-
ilpm.velocity.get_spatial_avg_enstrophy(udata, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=1.0, dy=1.0, dz=1.0)[source]¶ Return enstrophy averaged over space
- Parameters
- udata: nd array
- Returns
- enstrophy_vs_t: 1d numpy array
average enstrophy in a field at each time
- enstrophy_vs_t_err
standard deviation of enstrophy in a field at each time
-
ilpm.velocity.get_turbulence_intensity_local(udata)[source]¶ Turbulence intensity is defined as u/U where u = sqrt((ux**2 + uy**2 + uz**2)/3) # characteristic turbulent velocity U = sqrt((Ux**2 + Uy**2 + Uz**2)) # mean flow speed
This is ill-defined for turbulence with zero-mean flow !
- Parameters
- udata
- Returns
- ti_local: nd array
turbulent intensity field (scaler field)
-
ilpm.velocity.fft_nd(field, dx=1, dy=1, dz=1)[source]¶ - Parameters
- field: np array, (height, width, depth, duration) or (height, width, duration)
- dx: spacing along x-axis
- dy: spacing along x-axis
- dz: spacing along x-axis
- Returns
- field_fft
- np.asarray([kx, ky, kz])
-
ilpm.velocity.get_energy_spectrum_nd(udata, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=None, dy=None, dz=None, window=None, correct_signal_loss=True)[source]¶ Returns nd energy spectrum from velocity data (FFT of a velocity field)
- Parameters
- udata: nd array
- dx: data spacing in x (units: mm/px)
- dy: data spacing in y (units: mm/px)
- dz: data spacing in z (units: mm/px)
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- nkout: int, default: None
number of bins to compute energy/dissipation spectrum
- window: str
Windowing reduces undesirable effects due to the discreteness of the data. A wideband window such as ‘flattop’ is recommended for turbulent energy spectra.
For the type of applying window function, choose from below: boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann, kaiser (needs beta), gaussian (needs standard deviation), general_gaussian (needs power, width), slepian (needs width), chebwin (needs attenuation), exponential (needs decay scale), tukey (needs taper fraction)
- correct_signal_loss: bool, default: True
If True, it would compensate for the loss of the signals due to windowing. Always recommended to obtain accurate spectral densities.
- Returns
- energy_fft: nd array with shape (height, width, duration) or (height, width, depth, duration)
- ks: nd array with shape (ncomponents, height, width, duration) or (ncomponents, height, width, depth, duration)
-
ilpm.velocity.get_energy_spectrum(udata, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=None, dy=None, dz=None, nkout=None, window=None, correct_signal_loss=True, remove_undersampled_region=True, cc=1.75, notebook=True)[source]¶ Returns 1D energy spectrum from velocity field data
- Parameters
- udata: nd array
- epsilon: nd array or float, default: None
dissipation rate used for scaling energy spectrum If not given, it uses the values estimated using the rate-of-strain tensor
- nu: flaot, viscosity
- x0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- x1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- nkout: int, default: None
number of bins to compute energy/dissipation spectrum
- notebook: bool, default: True
Use tqdm.tqdm_notebook if True. Use tqdm.tqdm otherwise
- window: str
Windowing reduces undesirable effects due to the discreteness of the data. A wideband window such as ‘flattop’ is recommended for turbulent energy spectra.
For the type of applying window function, choose from below: boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann, kaiser (needs beta), gaussian (needs standard deviation), general_gaussian (needs power, width), slepian (needs width), chebwin (needs attenuation), exponential (needs decay scale), tukey (needs taper fraction)
- correct_signal_loss: bool, default: True
If True, it would compensate for the loss of the signals due to windowing. Always recommended to obtain accurate spectral densities.
- remove_undersampled_region: bool, default: True
If True, it will not sample the region with less statistics.
- cc: float, default: 1.75
A numerical factor to compensate for the signal loss due to approximations. … cc=1.75 was obtained from the JHTD data.
- Returns
- ——-
- e_k: numpy array
Energy spectrum with shape (number of data points, duration)
- e_k_err: numpy array
Energy spectrum error with shape (number of data points, duration)
- kk: numpy array
Wavenumber with shape (number of data points, duration)
-
ilpm.velocity.get_1d_energy_spectrum(udata, k='kx', x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=None, dy=None, dz=None, window=None, correct_signal_loss=True)[source]¶ Returns 1D energy spectrum from velocity field data
- Parameters
- udata: nd array
- k: str, default: ‘kx’
string to specify the direction along which the given velocity field is Fourier-transformed
- x0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- x1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- window: str
Windowing reduces undesirable effects due to the discreteness of the data. A wideband window such as ‘flattop’ is recommended for turbulent energy spectra.
For the type of available window function, choose from below: boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann, kaiser (needs beta), gaussian (needs standard deviation), general_gaussian (needs power, width), slepian (needs width), chebwin (needs attenuation), exponential (needs decay scale), tukey (needs taper fraction)
- correct_signal_loss: bool
If True, it would compensate for the loss of the signals due to windowing. Always recommended to obtain accurate spectral densities.
- Returns
- eiis: numpy array
eiis[0] = E11, eiis[1] = E22 … 1D energy spectra with argument k=”k” (kx by default)
- eii_errs: numpy array:
eiis[0] = E11_error, eiis[1] = E22_error
- k: 1d numpy array
Wavenumber with shape (number of data points, ) … Unlike get_energy_spectrum(…), this method NEVER outputs the wavenumber array with shape (number of data points, duration)
-
ilpm.velocity.get_dissipation_spectrum(udata, nu, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=None, dy=None, dz=None, nkout=None, window='flattop', correct_signal_loss=True, notebook=True)[source]¶ Returns dissipation spectrum D(k) = 2 nu k^2 E(k) where E(k) is the 1d energy spectrum
- Parameters
- udata: nd array
- epsilon: nd array or float, default: None
dissipation rate used for scaling energy spectrum If not given, it uses the values estimated using the rate-of-strain tensor
- nu: flaot, viscosity
- x0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- x1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- nkout: int, default: None
number of bins to compute energy/dissipation spectrum
- notebook: bool, default: True
Use tqdm.tqdm_notebook if True. Use tqdm.tqdm otherwise
- window: str
Windowing reduces undesirable effects due to the discreteness of the data. A wideband window such as ‘flattop’ is recommended for turbulent energy spectra.
For the type of available window function, choose from below: boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann, kaiser (needs beta), gaussian (needs standard deviation), general_gaussian (needs power, width), slepian (needs width), chebwin (needs attenuation), exponential (needs decay scale), tukey (needs taper fraction)
- correct_signal_loss: bool
If True, it would compensate for the loss of the signals due to windowing. Always recommended to obtain accurate spectral densities.
- Returns
- D_k: numpy array
Dissipation spectrum with shape (number of data points, duration)
- D_k_err: numpy array
Dissipation spectrum error with shape (number of data points, duration)
- k1d: numpy array
Wavenumber with shape (number of data points, duration)
-
ilpm.velocity.get_rescaled_energy_spectrum(udata, epsilon=100000, nu=1.0034, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=1, dy=1, dz=1, nkout=None, window=None, correct_signal_loss=True, notebook=True)[source]¶ - Returns SCALED energy spectrum E(k), its error, and wavenumber.
E(k) is sometimes called the 3D energy spectrum since it involves 3D FT of a velocity field.
ALWAYS greater than E11(k) and E22(k) which involve 1D FT of the velocity field along a specific direction.
Returns wavenumber with shape (# of data points, duration) instead of (# of data points, ).
… This seems redundant; however, values in the wavenumber array at a given time may fluctuate in principle due to a method how it computes the histogram (spectrum.) It should never fluctuate with the current method but is subject to a change by altering the method of determining the histogram. Therefore, this method outputs the wavenumber array with shape (# of data points, duration).
- Parameters
- udata: nd array
- epsilon: nd array or float, default: None
dissipation rate used for scaling energy spectrum If not given, it uses the values estimated using the rate-of-strain tensor
- nu: flaot, viscosity
- x0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- x1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- nkout: int, default: None
number of bins to compute energy/dissipation spectrum
- notebook: bool, default: True
Use tqdm.tqdm_notebook if True. Use tqdm.tqdm otherwise
- window: str
Windowing reduces undesirable effects due to the discreteness of the data. A wideband window such as ‘flattop’ is recommended for turbulent energy spectra.
For the type of available window function, choose from below: boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann, kaiser (needs beta), gaussian (needs standard deviation), general_gaussian (needs power, width), slepian (needs width), chebwin (needs attenuation), exponential (needs decay scale), tukey (needs taper fraction)
- correct_signal_loss: bool
If True, it would compensate for the loss of the signals due to windowing. Always recommended to obtain accurate spectral densities.
- Returns
- e_k_norm: numpy array
Scaled energy spectrum with shape (number of data points, duration)
- e_k_err_norm: numpy array
Scaled energy spectrum with shape (number of data points, duration)
- k_norm: numpy array
Scaled energy spectrum with shape (number of data points, duration)
-
ilpm.velocity.get_1d_rescaled_energy_spectrum(udata, epsilon=None, nu=1.0034, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=1, dy=1, dz=1, notebook=True, window='flattop', correct_signal_loss=True)[source]¶ Returns SCALED 1D energy spectra (E11 and E22) … Applies the flattop window function by default … Uses dissipation rate estimated by the rate-of-strain tensor unless given]
- Parameters
- udata: nd array
- epsilon: nd array or float, default: None
dissipation rate used for scaling energy spectrum If not given, it uses the values estimated using the rate-of-strain tensor
- nu: flaot, viscosity
- x0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- x1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- notebook: bool, default: True
Use tqdm.tqdm_notebook if True. Use tqdm.tqdm otherwise
- window: str
Windowing reduces undesirable effects due to the discreteness of the data. A wideband window such as ‘flattop’ is recommended for turbulent energy spectra.
For the type of available window function, choose from below: boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann, kaiser (needs beta), gaussian (needs standard deviation), general_gaussian (needs power, width), slepian (needs width), chebwin (needs attenuation), exponential (needs decay scale), tukey (needs taper fraction)
- correct_signal_loss: bool
If True, it would compensate for the loss of the signals due to windowing. Always recommended to obtain accurate spectral densities.
- Returns
- eii_arr_norm: nd array
Scaled spectral densities (Scaled E11 and Scaled E22)
- eii_err_arr_norm: nd array
Scaled errors for (Scaled E11 and Scaled E22)
- k_norm: nd array
Scaled wavenumber
-
ilpm.velocity.get_rescaled_dissipation_spectrum(udata, epsilon=100000, nu=1.0034, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=1, dy=1, dz=1, nkout=None, notebook=True)[source]¶ Return rescaled dissipation spectra D(k)/(u_eta^3) vs k * eta … convention: k = 2pi/ L
- Parameters
- udata: nd array
- nu: flaot, viscosity
- x0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- x1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- nkout: int, default: None
number of bins to compute energy/dissipation spectrum
- notebook: bool, default: True
Use tqdm.tqdm_notebook if True. Use tqdm.tqdm otherwise
- Returns
- D_k_norm: nd array
Scaled dissipation spectrum values with shape (# of data points, duration)
- D_k_err_norm: nd array
Scaled dissipation spectrum error with shape (# of data points, duration)
- k_norm: nd array
Scaled wavenumber with shape (# of data points, duration)
-
ilpm.velocity.scale_energy_spectrum(e_k, kk, epsilon=100000, nu=1.0034, e_k_err=None)[source]¶ Scales raw energy spectrum by given dissipation rate and viscosity
- Parameters
- e_k: numpy array
spectral energy density
- kk: numpy array
wavenumber
- epsilon: numpy array or float
dissipation rate used for scaling. It could be 1D numpy array or float.
- nu numpy array or float
viscosity used for scaling. It could be 1D numpy array or float.
- e_k_err: numpy array
error of dissipation rate used for scaling. It could be 1D numpy array or float.
- Returns
- ——-
- e_k_norm, k_norm if e_k_err is not given
- e_k_norm, e_k_err_norm, k_norm if e_k_err is given
-
ilpm.velocity.get_large_scale_vel_field(kx^2 + ky^2) < kmax in the original vel. field (udata)[source]¶ Returns a velocity field which satisfies k = sqrt(kx^2 + ky^2) < kmax in the original vel. field (udata)
- Parameters
- udata: nd array
- … velocity field data with shape (# of components, physical dimensions (width x height x depth), duration)
- kmax: float
- … value of k below which spectrum is kept. i.e. cutoff k for the low-pass filter
- x0: int
- … index used to specify a region of a vel. field in which spectrum is computed. udata[y0:y1, x0:x1, z0:z1]
- x1: int
- … index used to specify a region of a vel. field in which spectrum is computed. udata[y0:y1, x0:x1, z0:z1]
- y0: int
- … index used to specify a region of a vel. field in which spectrum is computed. udata[y0:y1, x0:x1, z0:z1]
- y1: int
- … index used to specify a region of a vel. field in which spectrum is computed. udata[y0:y1, x0:x1, z0:z1]
- z0: int
- … index used to specify a region of a vel. field in which spectrum is computed. udata[y0:y1, x0:x1, z0:z1]
- z1: int
- … index used to specify a region of a vel. field in which spectrum is computed. udata[y0:y1, x0:x1, z0:z1]
- window: 2d/3d nd array
- … window used to make input data periodic to surpress spectral leakage in FFT
- dx: float
- … x spacing. used to compute k
- dy: float
- … y spacing. used to compute k
- dz: float
- … z spacing. used to compute k
- Returns
- udata_ifft: nd array
- … low-pass filtered velocity field data
- coords: nd array
- … if dim == 2, returns np.asarray([xgrid, ygrid])
- … if dim == 3, returns np.asarray([xgrid, ygrid, zgrid])
-
ilpm.velocity.get_epsilon_using_sij(udata, dx=None, dy=None, dz=None, nu=1.004, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, t0=0, t1=None)[source]¶ Returns the dissipation rate computed using the rate-of-strain tensor
sij: numpy array with shape (nrows, ncols, duration, 2, 2) (dim=2) or (nrows, ncols, nstacks, duration, 3, 3) (dim=3) … idea is… sij[spacial coordinates, time, tensor indices]
e.g.- sij(x, y, t) can be accessed by sij[y, x, t, i, j]
… sij = d ui / dxj Parameters ———- udata: nd array nu: flaot, viscosity x0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- x1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- Returns
- epsilon: numpy array
dissipation rate
-
ilpm.velocity.get_epsilon_iso(udata, lambda_f=None, lambda_g=None, nu=1.004, x=None, y=None, **kwargs)[source]¶ Return epsilon computed by isotropic formula involving Taylor microscale
- Parameters
- udata
- lambda_f: numpy array
long. Taylor microscale
- lambda_g: numpy array
transverse. Taylor microscale
- nu: float
viscosity
- Returns
- epsilon: numpy array
dissipation rate
-
ilpm.velocity.get_epsilon_using_diss_spectrum(udata, nu=1.0034, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, dx=1, dy=1, dz=1, nkout=None, notebook=True)[source]¶ Returns dissipation rate computed by integrated the dissipation specrtrum … must have fully resolved spectrum to yield a reasonable result
- Parameters
- udata: nd array
- nu: flaot, viscosity
- x0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- x1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- y1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t0: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- t1: int
index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1].
- dx: float
spacing in x
- dy: float
spacing in y
- dz: float
spacing in z
- nkout: int, default: None
number of bins to compute energy/dissipation spectrum
- notebook: bool, default: True
Use tqdm.tqdm_notebook if True. Use tqdm.tqdm otherwise
- Returns
- epsilon: numpy array
dissipation rate
-
ilpm.velocity.get_epsilon_using_struc_func(rrs, Dxxs, epsilon_guess=100000, r0=1.0, r1=10.0, p=2, method='Nelder-Mead')[source]¶ Returns the values of estimated dissipation rate using a long. structure function
- Parameters
- rrs: numpy array
separation length for long. transverse function
- Dxxs: numpy array
values of long. transverse function
- epsilon_guess: numpy array
initial guess for dissipation rate
- r0: float
The scaled structure function is expected to have a plateau [r0, r1]
- r1: float
The scaled structure function is expected to have a plateau [r0, r1]
- p: float/int
order of structure function … smaller order is better due to intermittency
- method: str, default: ‘Nelder-Mead’
method used to find the minimum of the test function
- Returns
- epsilons: numpy array
estimated dissipation rate
-
ilpm.velocity.compute_spatial_autocorr(ui, x, y, roll_axis=1, n_bins=None, x0=0, x1=None, y0=0, y1=None, t0=None, t1=None, coarse=1.0, coarse2=0.2, notebook=True)[source]¶ Compute spatial autocorrelation function of 2+1 velocity field Spatial autocorrelation function:
f = <u_j(
ec{x}) u_j( ec{x}+rhat{x_i})> / <u_j( ec{x})^2>
where velocity vector u = (u1, u2). If i = j, f is called longitudinal autocorrelation function. Otherwise, f is called transverse autocorrelation function.
- Example:
u = ux # shape(height, width, duration) x_, y_ = range(u.shape[1]), range(u.shape[0]) x, y = np.meshgrid(x_, y_)
# LONGITUDINAL AUTOCORR FUNC rrs, corrs, corr_errs = compute_spatial_autocorr(u, x, y, roll_axis=1) # roll_axis is i in the description above
# TRANSVERSE AUTOCORR FUNC rrs, corrs, corr_errs = compute_spatial_autocorr(u, x, y, roll_axis=0) # roll_axis is i in the description above
- Parameters
- ui: numpy array, 2 + 1 scalar field. i.e. shape: (height, width, duration)
x: numpy array, 2d grid y: numpy array, 2d grid roll_axis: int, axis index to compute correlation… 0 for y-axis, 1 for x-axis n_bins: int, number of bins to compute statistics x0: int, index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1]. x1: int, index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1]. y0: int, index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1]. y1: int, index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1]. t0: int, index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1]. t1: int, index to specify a portion of data in which autocorrelation funciton is computed. Use data u[y0:y1, x0:x1, t0:t1]. coarse: float (0, 1], Process coarse * possible data points. This is an option to output coarse results. Returns ——- rr: 2d numpy array, (distance, time) corr: 2d numpy array, (autocorrelation values, time) corr_err: 2d numpy array, (std of autocorrelation values, time)
-
ilpm.velocity.compute_spatial_autocorr3d(ui, x, y, z, roll_axis=1, n_bins=None, x0=None, x1=None, y0=None, y1=None, z0=None, z1=None, t0=None, t1=None, coarse=1.0, coarse2=0.2, notebook=True)[source]¶ Compute spatial autocorrelation function of 2+1 velocity field Spatial autocorrelation function:
f = <u_j(
ec{x}) u_j( ec{x}+rhat{x_i})> / <u_j( ec{x})^2>
where velocity vector u = (u1, u2). If i = j, f is called longitudinal autocorrelation function. Otherwise, f is called transverse autocorrelation function.
- Example:
u = ux # shape(height, width, duration) x_, y_ = range(u.shape[1]), range(u.shape[0]) x, y = np.meshgrid(x_, y_)
# LONGITUDINAL AUTOCORR FUNC rrs, corrs, corr_errs = compute_spatial_autocorr(u, x, y, roll_axis=1) # roll_axis is i in the description above
# TRANSVERSE AUTOCORR FUNC rrs, corrs, corr_errs = compute_spatial_autocorr(u, x, y, roll_axis=0) # roll_axis is i in the description above
- Parameters
- ui: numpy array, 3 + 1 scalar field. i.e. shape: (height, width, duration)
x: numpy array, 3d grid y: numpy array, 3d grid z: numpy array, 3d grid roll_axis: int, axis index to compute correlation… 0 for y-axis, 1 for x-axis, 2 for z-axis n_bins: int, number of bins to compute statistics x0: int, index to specify a portion of data in which autocorrelation funciton is computed. Will use data u[y0:y1, x0:x1, z0:z1, t0:t1]. x1: int, index to specify a portion of data in which autocorrelation funciton is computed. y0: int, index to specify a portion of data in which autocorrelation funciton is computed. y1: int, index to specify a portion of data in which autocorrelation funciton is computed. z0: int, index to specify a portion of data in which autocorrelation funciton is computed. z1: int, index to specify a portion of data in which autocorrelation funciton is computed. t0: int, index to specify a portion of data in which autocorrelation funciton is computed. t1: int, index to specify a portion of data in which autocorrelation funciton is computed. coarse: float (0, 1], Process coarse * possible data points. This is an option to output coarse results. coarse2: float (0, 1], Rolls matrix Returns ——- rr: 2d numpy array, (distance, time) corr: 2d numpy array, (autocorrelation values, time) corr_err: 2d numpy array, (std of autocorrelation values, time)
-
ilpm.velocity.get_two_point_vel_corr_iso(udata, x, y, z=None, time=None, n_bins=None, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, t0=None, t1=None, coarse=1.0, coarse2=0.2, notebook=True, return_rij=False, **kwargs)[source]¶ Returns two-point velocity autocorrelation tensor, and autocorrelation functions. Uses the x-component of velocity. (CAUTION required for unisotropic flows)
Pope Eq. 6.44 Parameters ———- udata: 5D or 4D numpy array, 5D if the no. of spatial dimensions is 3. 4D if the no. of spatial dimensions is 2.
… (ux, uy, uz) or (ux, uy) … ui has a shape (height, width, depth, duration) or (height, width, depth) (3D) … ui may have a shape (height, width, duration) or (height, width) (2D)
- Returns
- rij: nd array (r, t, i, j)
… Rij (
- ec{r} , t) = <u_j(
- ec{x}) u_j(
- ec{x}+
- ec{r})> / <u_j(
- ec{x})^2>
- … If system is homogeneous and isotropic,
Rij (
- ec{r} , t) = u_rms^2 [g(r,t) delta_ij + {f(r,t) - g(r,t)} r_i * r_j / r^2]
where f, and g are long. and transverse autocorrelation functions.
-
ilpm.velocity.get_autocorr_functions(r_long, f_long, r_tran, g_tran, time)[source]¶ Return interpolated functions using the outputs of get_two_point_vel_corr_iso() … the returned objects are functions NOT arrays
- Parameters
- r_long: numpy array
output of get_two_point_vel_corr_iso()
- f_long: numpy array
output of get_two_point_vel_corr_iso()
- r_tran: numpy array
output of get_two_point_vel_corr_iso()
- g_tran: numpy array
output of get_two_point_vel_corr_iso()
- time: numpy array
time corresponding to the given autocorrelation functions
- Returns
- f, g: long./trans. autocorrelation functions
-
ilpm.velocity.get_autocorr_functions_int_list(r_long, f_long, r_tran, g_tran)[source]¶ Returns lists of INTERPOLATED autocorrelation function values … outputs of get_two_point_vel_corr_iso() may contain np.nan and np.inf which could be troublesome … this method gets rid of these values, and interpolates the missing values (third-order spline) … the arguments should have a shape (# of data points, duration) … this method conducts intepoltation at given time respectively because 2D interpolation often raises an error due to the discontinuity of the autocorrelation functions along the temporal axis.
- Parameters
- r_long: numpy array
output of get_two_point_vel_corr_iso()
- f_long: numpy array
output of get_two_point_vel_corr_iso()
- r_tran: numpy array
output of get_two_point_vel_corr_iso()
- g_tran: numpy array
output of get_two_point_vel_corr_iso()
- Returns
- fs: list with length = duration = f_long.shape[-1]
list of interpolated longitudinal structure function
- gs: list with length = duration = f_long.shape[-1]
list of interpolated transverse structure function
-
ilpm.velocity.get_autocorrelation_tensor_iso(r_long, f_long, r_tran, g_tran, time)[source]¶ Returns autocorrelation tensor with isotropy assumption … not recommended for practical use due to long convergence time
- Parameters
- r_long: numpy array
output of get_two_point_vel_corr_iso()
- f_long: numpy array
output of get_two_point_vel_corr_iso()
- r_tran: numpy array
output of get_two_point_vel_corr_iso()
- g_tran: numpy array
output of get_two_point_vel_corr_iso()
- time: numpy array
time corresponding to the given autocorrelation functions
- Returns
- rij: autocorrelation tensor
technically a function with arguments: i, j, r, t, udata (tensor indices, separation distance, time index, udata)
-
ilpm.velocity.get_structure_function_long(udata, x, y, z=None, p=2, roll_axis=1, n_bins=None, nu=1.004, u='ux', x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, t0=0, t1=None, coarse=1.0, coarse2=0.2, notebook=True)[source]¶ Structure tensor Dij is essentially the covariance of the two-point velocity difference There is one-to-one correspondence between Dij and Rij. (Pope 6.36) This method returns the LONGITUDINAL STRUCTURE FUNCTION. If p=2, this returns D_LL. … Returns rrs, Dxxs, Dxx_errs, rrs_scaled, Dxxs_scaled, Dxx_errs_scaled Parameters ———- udata x: numpy array
x-coordinate of the spatial grid corresponding to the given udata … it does not have to ve equally spaced
- y: numpy array
y-coordinate of the spatial grid corresponding to the given udata … it does not have to ve equally spaced
- z: numpy array
z-coordinate of the spatial grid corresponding to the given udata … it does not have to ve equally spaced
- p: float or int
order of structure function
- roll_axis: int
“u” and “roll_axis” determines whether this method returns the longitudinal or transverse structure function If you want longitudinal, match “u” and “roll_axis”. e.g.- u=’ux’ and roll_axis=1, u=’uy’ and roll_axis=0 If you want transverse, do not match “u” and “roll_axis”. e.g.- u=’ux’ and roll_axis=0, u=’uy’ and roll_axis=1
- n_bins: int, default=None
number of bins used to take a histogram of velocity difference
- nu: float
viscosity
- u: str, default=’ux’
velocity component used to compute the structure functions. Choices are ‘ux’, ‘uy’, ‘uz’
- x0: int, default: 0
Specified the region of udata to compute the structure function. This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- x1: int, default: None
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- y0: int, default: 0
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- y1 int, default: None
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- z0: int, default: 0
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- z1 int, default: None
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- t0: int, default: 0
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- t1 int, default: None
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- coarse: float, (0, 1], default: 1
the first parameter to save computation time related sampling frequency … The higher “coarse” is, it samples more possible data points. … If “coarse” == 1, it samples all possible data points. … If “coarse” == 0.5, it samples only a half of possible data points. Could overestimate Taylor microscale.
- coarse2: float, (0, 1], default: 0.2
the second parameter to save computation time related to making a histogram … Determines how many sampled data points to be used to make a histogram … If “coarse” == 1, it uses all data points to make a histogram. … If “coarse” == 0.5, it uses only a half of data points to make a histogram
- notebook: bool
… if True, it uses tqdm.tqdm_notebook instead of tqdm.tqdm
- Returns
- rrs: numpy array
two-point separation distance for the structure function
- Dxxs: numpy array
values of the structure function
- Dxx_errs: numpy array
error of the structure function
- rrs_scaled: numpy array
two-point separation distance for the SCALED structure function
- Dxxs_scaled: numpy array
values of the SCALED structure function
- Dxx_errs_scaled: numpy array
error of the SCALED structure function
-
ilpm.velocity.get_structure_function(udata, x, y, z=None, indices=('x', 'x'), roll_axis=1, n_bins=None, nu=1.004, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, t0=0, t1=None, coarse=1.0, coarse2=0.2, notebook=True)[source]¶ A method to compute a generalized structure function
Structure tensor Dij…k is essentially the generalized variance of the two-point velocity difference This method returns the STRUCTURE FUNCTION D_{ij…k}(r x_m) where the tensor indices are “indices” and the subscript of x_m, m, is expressed as roll_axis.
If indices=(1, 1) & roll_axis=1, this returns D_{xx}(r hat{x}). # longitudinal, 2nd order If indices=(0, 0) & roll_axis=0, this returns D_{yy}(r hat{y}). # longitudinal, 2nd order If indices=(0, 0) & roll_axis=0, this returns D_{yy}(r hat{x}). # transverse, 2nd order If indices=(1, 1, 1) & roll_axis=1, this returns D_{xxx}(r hat{x}). # longitudinal, 3nd order If indices=(1, 0, 1) & roll_axis=1, this returns D_{xyx}(r hat{x}). # (1, 0, 1)-component of 3rd order structure function tensor
… Returns rrs, Dijks, Dijk_errs, rrs_scaled, Dijks_scaled, Dijk_errs_scaled
- Parameters
- udata
- x: numpy array
x-coordinate of the spatial grid corresponding to the given udata … it does not have to ve equally spaced
- y: numpy array
y-coordinate of the spatial grid corresponding to the given udata … it does not have to ve equally spaced
- z: numpy array
z-coordinate of the spatial grid corresponding to the given udata … it does not have to ve equally spaced
- indices: array-like, default: (1, 1)
… tensor indices of the structure function tensor: D_{i,j,…, k} … indices=(1,1) corresponds to D_{xx}. … For a 3D spatial field, x, y, z corresponds to 1, 0, 2 respectively. Recall udata.shape = (height, width, depth, duration) … For a 2D spatial field, x, y corresponds to 1, 0 respectively. Recall udata.shape = (height, width, duration)
- roll_axis: int
“u” and “roll_axis” determines whether this method returns the longitudinal or transverse structure function If you want longitudinal, match “u” and “roll_axis”. e.g.- u=’ux’ and roll_axis=1, u=’uy’ and roll_axis=0 If you want transverse, do not match “u” and “roll_axis”. e.g.- u=’ux’ and roll_axis=0, u=’uy’ and roll_axis=1
- n_bins: int, default=None
number of bins used to take a histogram of velocity difference
- nu: float
viscosity
- u: str, default=’ux’
velocity component used to compute the structure functions. Choices are ‘ux’, ‘uy’, ‘uz’
- x0: int, default: 0
Specified the region of udata to compute the structure function. This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- x1: int, default: None
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- y0: int, default: 0
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- y1 int, default: None
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- z0: int, default: 0
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- z1 int, default: None
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- t0: int, default: 0
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- t1 int, default: None
This method uses only udata[y0:y1, x0:x1, z0:z1, t0:t1].
- coarse: float, (0, 1], default: 1
the first parameter to save computation time related sampling frequency … The higher “coarse” is, it samples more possible data points. … If “coarse” == 1, it samples all possible data points. … If “coarse” == 0.5, it samples only a half of possible data points. Could overestimate Taylor microscale.
- coarse2: float, (0, 1], default: 0.2
the second parameter to save computation time related to making a histogram … Determines how many sampled data points to be used to make a histogram … If “coarse” == 1, it uses all data points to make a histogram. … If “coarse” == 0.5, it uses only a half of data points to make a histogram
- notebook: bool
… if True, it uses tqdm.tqdm_notebook instead of tqdm.tqdm
- Returns
- rrs: numpy array
two-point separation distance for the structure function
- Dxxs: numpy array
values of the structure function
- Dxx_errs: numpy array
error of the structure function
- rrs_scaled: numpy array
two-point separation distance for the SCALED structure function
- Dxxs_scaled: numpy array
values of the SCALED structure function
- Dxx_errs_scaled: numpy array
error of the SCALED structure function
-
ilpm.velocity.scale_raw_structure_funciton_long(rrs, Dxxs, Dxx_errs, epsilon, nu=1.004, p=2)[source]¶ Returns the scaled structure functions when raw structure function data are given … This allows users to scale the structure functions with epsilon and nu input by the users.
- Parameters
- rrs: numpy array
separation between two points, output of get_two_point_vel_corr_iso()
- Dxxs: numpy array
longitudinal structure function, output of get_two_point_vel_corr_iso()
- Dxx_errs: numpy array
error of longitudinal structure function, output of get_two_point_vel_corr_iso()
- epsilon: numpy array or float
dissipation rate
- nu: float
viscosity
- p: int/float
order of structure function
- Returns
- rrs_s: numpy array
Scaled r
- Dxxs_s: numpy array
Scaled structure function
- Dxx_errs_s: numpy array
Scaled DLL error
-
ilpm.velocity.remove_nans_for_array_pair(arr1, arr2)[source]¶ remove nans or infs for a pair of 1D arrays, and returns the compressed arrays with the same length e.g.: arr1 = [0, 0.1, 0.2, np.nan, 0.4, 0.5], arr2 = [-1.2, 19.2. 155., np.inf, 0.1]
-> compressed_arr1 = [0, 0.1, 0.2, 0.5], compressed_arr2 = [-1.2, 19.2., 0.1]
- Parameters
- arr1: numpy array
… may include np.nan and np.inf
- arr2: numpy array
… may include np.nan and np.inf
- Returns
- ——-
- compressed_arr1: numpy array
… without any np
- compressed_arr2: numpy array
-
ilpm.velocity.get_taylor_microscales(r_long, f_long, r_tran, g_tran, residual_thd=0.015, deg=2, return_err=False, npt_min=4)[source]¶ Returns Taylor microscales as the curvature of the autocorrelation functions at r=0 … Algorithm: … (1) Polynomial fit (cubic) the long./trans. autocorr functions … (2) Evaluate its second derivate at r=0. … (3) Relate that value to the taylor microscale … Use the first p*100 % of the autocorrelation values. … Designed to use the results of get_two_point_vel_corr_iso().
- Parameters
- r_long: numpy 2d array with shape (no. of elements, duration)
… r for longitudinal autoorrelation function
- f_long: numpy 2d array with shape (no. of elements, duration)
… longitudinal autoorrelation values
- r_tran: numpy 2d array with shape (no. of elements, duration)
… r for longitudinal autoorrelation function
- g_tran: numpy 2d array with shape (no. of elements, duration)
… longitudinal autoorrelation values
- Returns
- lambda_f: numpy 2d array with shape (duration, )
Longitudinal Taylor microscale
- lambda_g: numpy 2d array with shape (duration, )
Transverse Taylor microscale
-
ilpm.velocity.get_taylor_microscales_iso(udata, epsilon, nu=1.004)[source]¶ Return Taylor microscales computed by isotropic formulae: lambda_g_iso = (15 nu * u_irms^2 / epsilon) ^ 0.5
- Parameters
- udata: nd array
- epsilon: float or array with the same length as udata
- nu: float
viscoty
- Returns
- lambda_f_iso: numpy array
longitudinal Taylor microscale
- lambda_g_iso: numpy array
transverse Taylor microscale
-
ilpm.velocity.get_integral_scales(r_long, f_long, r_tran, g_tran, method='trapz')[source]¶ Returns integral scales computed by using autocorrelation functions
- Parameters
- r_long: numpy array
… output of get_two_point_vel_corr_iso(…)
- f_long numpy array
… output of get_two_point_vel_corr_iso(…)
- r_tran numpy array
… output of get_two_point_vel_corr_iso(…)
- g_tran numpy array
… output of get_two_point_vel_corr_iso(…)
- method: string, default: ‘trapz’
… integration method, choose from quadrature or trapezoidal
- Returns
- L11: numpy array
longitudinal integral length scale
- L22: numpy array
transverse integral length scale
-
ilpm.velocity.get_integral_scales_using_rij(udata, Rij, rmax, n=100)[source]¶ Use autocorrelation tensor, Rij to calculate integral length scale Parameters ———- udata Rij: numpy array
two-point velocity autocorrelation tensor … can be obtained by get_two_point_vel_corr_iso(…) but may take a long time
rmax: float n: int
The higher n is, the integrand function becomes more smooth
- Returns
- L11: nd array
… longitudinal integral length scale
- L22: nd array
… transverse integral length scale
-
ilpm.velocity.get_integral_scales_iso_spec(udata, e_k, k)[source]¶ Integral scale defined through energy spectrum. Assumes isotropy and a full 1D energy spectrum. Pope 6.260.
- Parameters
- udata
- e_k: numpy array
output of get_energy_spectrum()
- k: numpy array
output of get_energy_spectrum()
- Returns
- L_iso_spec: 1d array
-
ilpm.velocity.get_integral_scale_large_eddy(udata, epsilon)[source]¶ dissipation rate (epsilon) is known to be proportional to u’^3 / L.
Some just define an integral scale L as u’^3 / epsilon. This method just returns this ratio. It is often interpreted as the characteristic scale of large eddies.
- Parameters
- udata: numpy array
- epsilon: numpy array / float
dissipation rate
- Returns
- L: numpy array
integral length scale (characterisic size of large eddies)
-
ilpm.velocity.get_integral_velocity_scale(udata)[source]¶ Return integral velocity scale which is identical to u’ (characteristic velocity) See get_characteristic_velocity()
- Parameters
- udata: numpy array
- Returns
- u_irms: See get_characteristic_velocity()
-
ilpm.velocity.get_kolmogorov_scale(udata, dx, dy, dz=None, nu=1.004)[source]¶ Returns kolmogorov LENGTH scale … estimates dissipation rate from the rate-of-strain tensor
- Parameters
- udata: numpy array
- dx: float
data spacing in x
- dy: float
data spacing in y
- dz: float
data spacing in z
- nu: float
viscosity
- Returns
- eta: numpy array
kolmogorov length scale
-
ilpm.velocity.get_integral_scales_all(udata, dx, dy, dz=None, nu=1.004)[source]¶ Returns integral scales (related to LARGE EDDIES)
- Parameters
- udata: numpy array
- dx: float
data spacing in x
- dy: float
data spacing in y
- dz: float
data spacing in z
- nu: float
viscosity
- Returns
- L_le: 1d array
- u_L: 1d array
- tau_L: 1d array
-
ilpm.velocity.get_taylor_microscales_all(udata, r_long, f_long, r_tran, g_tran)[source]¶ Returns Taylor microscales Parameters ———- udata: numpy array
velocity field array
- r_long: numpy array
… output of get_two_point_vel_corr_iso(…)
- f_long: numpy array
… output of get_two_point_vel_corr_iso(…)
- r_tran: numpy array
… output of get_two_point_vel_corr_iso(…)
- g_tran: numpy array
… output of get_two_point_vel_corr_iso(…)
- Returns
- lambda_f: 1d array
Taylor microscale (length)
- u_lambda: 1d array
Taylor microscale (velocity)
- tau_lambda: 1d array
Taylor microscale (time)
-
ilpm.velocity.get_kolmogorov_scales_all(udata, dx, dy, dz, nu=1.004)[source]¶ Returns Kolmogorov scales
- Parameters
- udata: numpy array
- dx: float
data spacing in x
- dy: float
data spacing in y
- dz: float
data spacing in z
- nu: float
viscosity
- Returns
- eta: 1d array
- u_eta: 1d array
- tau_eta: 1d array
-
ilpm.velocity.get_turbulence_re(udata, dx, dy, dz=None, nu=1.004)[source]¶ Returns turbulence reynolds number (Pope 6.59) … Integral Reynolds number
- Parameters
- udata: numpy array
- dx: float
data spacing in x
- dy: float
data spacing in y
- dz: float
data spacing in z
- nu: float
viscosity
- Returns
- Re_L: numpy array
Turbulence Reynolds number
-
ilpm.velocity.get_taylor_re(udata, r_long, f_long, r_tran, g_tran, nu=1.004)[source]¶ Returns Taylor reynolds number (Pope 6.63) from autocorrelation functions
- Parameters
- udata: numpy array
velocity field array
- r_long: numpy array
… output of get_two_point_vel_corr_iso(…)
- f_long: numpy array
… output of get_two_point_vel_corr_iso(…)
- r_tran: numpy array
… output of get_two_point_vel_corr_iso(…)
- g_tran: numpy array
… output of get_two_point_vel_corr_iso(…)
- nu: float, default: 1.004 (water, in mm^2/s)
viscosity
- Returns
- Re_lambda: array
Taylor Reynolds number
-
ilpm.velocity.rankine_vortex_2d(xx, yy, x0=0, y0=0, gamma=1.0, a=1.0)[source]¶ Reutrns a 2D velocity field with a single Rankine vortex at (x0, y0)
- Parameters
- xx: numpy array
x-coordinate, 2d grid
- yy: numpy array
y-coordinate, 2d grid
- x0: float
x-coordinate of the position of the rankine vortex
- y0: float
y-coordinate of the position of the rankine vortex
- gamma: float
circulation of the rankine vortex
- a: float
core diameter of the rankine vortex
- Returns
- udata: (ux, uy)
-
ilpm.velocity.rankine_vortex_line_3d(xx, yy, zz, x0=0, y0=0, gamma=1.0, a=1.0, uz0=0)[source]¶ Reutrns a 3D velocity field with a Rankine vortex filament at (x0, y0, z)
- Parameters
- xx: 3d numpy grid
- yy: 3d numpy grid
- zz: 3d numpy grid
- x0: float, location of Rankine vortex filament
- y0: float, location of Rankine vortex filament
- gamma: float, circulation
- a: float, core radius
- uz0: float, constant velocity component in the z-direction
- Returns
- udata: (ux, uy, uz)
-
ilpm.velocity.get_sample_turb_field_3d(return_coord=True)[source]¶ Returns udata=(ux, uy, uz) of a slice of isotropic, homogeneous turbulence data (DNS, JHTD)
- Parameters
- return_coord
- Returns
- udata: numpy array
velocity field
- xx: numpy array
x-coordinate of 3D grid
- yy: numpy array
y-coordinate of 3D grid
- zz: numpy array
z-coordinate of 3D grid
-
ilpm.velocity.get_rescaled_energy_spectrum_saddoughi()[source]¶ Returns values to plot rescaled energy spectrum from Saddoughi (1992)
- Returns
- e: numpy array
spectral energy density: energy stored between [k, k+dk)
- k: numpy array
wavenumber
-
ilpm.velocity.get_energy_spectra_jhtd()[source]¶ Returns values to plot energy spectrum from JHTD, computed by Takumi in 2019 Call get_rescaled_energy_spectra_jhtd for scaled energy spectrum.
- Returns
- datadict: dict
data stored in jhtd_e_specs.h5 is stored: k, ek
-
ilpm.velocity.get_rescaled_energy_spectra_jhtd()[source]¶ Returns values to plot rescaled energy spectrum from JHTD, computed by Takumi in 2019 Call get_energy_spectra_jhtd for raw energy spectrum.
- Returns
- datadict: dict
data stored in jhtd_e_specs.h5 is stored: Scaled k, Scaled ek
-
ilpm.velocity.get_rescaled_structure_function_saddoughi(p=2)[source]¶ Returns the values of rescaled structure function reported in Saddoughi and Veeravalli 1994 paper: r_scaled, dll … this is a curve about a specific Reynolds number! i.e. there is no universal structure function
p: int … order of the structure function
- Returns
- r_scaled: nd array
- dll: nd array
-
ilpm.velocity.get_window_radial(xx, yy, zz=None, wtype='hamming', rmax=None, duration=None, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None, n=500)[source]¶ General method to get a window with shape (xx.shape[:], duration) or (xx.shape[:]) if duration is None … Window types:
boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann, kaiser (needs beta), gaussian (needs standard deviation), general_gaussian (needs power, width), slepian (needs width), chebwin (needs attenuation), exponential (needs decay scale), tukey (needs taper fraction)
- Parameters
- xx: nd array
x-coordinate of the spatial grid of udata
- yy: nd array
y-coordinate of the spatial grid of udata
- zz: nd array
y-coordinate of the spatial grid of udata
- wtype: str
name of window function such as ‘hamming’, ‘flattop’
- rmax: float
window function returns zero for r > rmax
- duration: int, default: None
specifies the temporal dimension of the returning window function.
- x0: int, default: 0
used to specify the region of the returning window function. When coordinates outside the specified region is given to the window function, it returns 0.
- x1: int, default: None
used to specify the region of the returning window function. When coordinates outside the specified region is given to the window function, it returns 0.
- y0: int, default: 0
- y1: int, default: None
used to specify the region of the returning window function. When coordinates outside the specified region is given to the window function, it returns 0.
- z0: int, default: 0
- z1: int, default: None
used to specify the region of the returning window function. When coordinates outside the specified region is given to the window function, it returns 0.
- n: int, default: 500
number of data points when 1D window function is called for the first time. The higher n is, the returning windowing function is more accurate.
- Returns
- window/windows: nd array
- … window: hamming window with the shape as xx
- … window: hamming window with the shape (xx.shape[:], duration)
-
ilpm.velocity.compute_signal_loss_due_to_windowing(xx, yy, window, x0=0, x1=None, y0=0, y1=None)[source]¶ Returns the inverse of the signal-loss factor by the window Signal loss factor: 1 (rectangle… no loss), 0.2-ish (flattop)
- Parameters
- xx: 2d/3d numpy array
x grid
- yy: 2d/3d numpy array
y grid
- window: str
name of the window: flattop, hamming, blackman, etc.
- x0: int
index to specify the region to which the window applies: xx[y0:y1, x0:x1], yy[y0:y1, x0:x1]
- x1: int
index to specify the region to which the window applies: xx[y0:y1, x0:x1], yy[y0:y1, x0:x1]
- y0: int
index to specify the region to which the window applies: xx[y0:y1, x0:x1], yy[y0:y1, x0:x1]
- y1: int
index to specify the region to which the window applies: xx[y0:y1, x0:x1], yy[y0:y1, x0:x1]
- Returns
- gamma: float
inverse of the signal-loss factor
-
ilpm.velocity.get_hamming_window_radial(xx, yy, zz=None, rmax=None, duration=None, x0=0, x1=None, y0=0, y1=None, z0=0, z1=None)[source]¶ - Returns the Hamming window as a function of radius
… r = 0 corresponds to the center of the window.
- Parameters
- r: nd array
- … radius from the center point
- rmax: float
- …
- duration: int
- … if None, this returns the hamming window with shape (height, width), (height, width, depth)
- … Otherwise, this returns (height, width, duration), (height, width, depth, duration)
- … duration = udata.shape[-1] should work for most of purposes.
- Returns
- window/windows: nd array
- … window: hamming window with the shape as xx
- … window: hamming window with the shape (xx.shape[:], duration)
-
ilpm.velocity.clean_udata_cheap(udata, cutoffU=2000, fill_value=nan, verbose=True)[source]¶ ONLY WORKS FOR THE 2D data Conducts a cheap bilinear interpolation for missing data. … literally, computes the average of the values interpolated in the x- and y-directions … griddata performs a better interpolation but this method is much faster. … values near the edges must not be trusted. Parameters ———- udata cutoffU fill_value verbose
-
ilpm.velocity.get_mask_for_unphysical(U, cutoffU=2000.0, fill_value=99999.0, verbose=True)[source]¶ Returns a mask (N-dim boolean array). If elements were below/above a cutoff, np.nan, or np.inf, then they get masked.
- Parameters
- U: array-like
- cutoffU: float
if |value| > cutoff, this method considers those values unphysical.
- fill_value:
- Returns
- mask: nd boolean array
-
ilpm.velocity.fill_unphysical_with_sth(U, mask, fill_value=nan)[source]¶ Returns an array whose elements are replaced by fill_value if its mask value is True
- Parameters
- U: array-like
- mask: multidimensional boolean array
- fill_value: value that replaces masked values
- Returns
- U_filled: numpy array
-
ilpm.velocity.clean_udata(udata, xx, yy, cutoffU=2000, fill_value=nan, verbose=True, method='linear', notebook=True)[source]¶ Conducts 2d interpolation for udata using scipy.interpolate.griddata … applies a mask to ignore the values whose magnitudes are greater than the cutoff … interpolates the missing data (np.nan or np.inf)
- Parameters
- udata
- xx
- yy
- cutoffU
- fill_value
- verbose
- method
- notebook
- Returns
- udata_i: interpolated udata which includes neither np.nan nor np.iinf
-
ilpm.velocity.fix_udata_shape(udata)[source]¶ It is better to always have udata with shape (height, width, depth, duration) (3D) or (height, width, duration) (2D) This method fixes the shape of udata such that its shape is (height, width, depth) or (height, width)
- Parameters
- udata: nd array,
… with shape (height, width, depth) (3D) or (height, width, duration) (2D) … OR shape (height, width, depth, duration) (3D) or (height, width, duration) (2D)
- Returns
- udata: nd array, with shape (height, width, depth, duration) (3D) or (height, width, duration) (2D)
-
ilpm.velocity.get_equally_spaced_grid(udata, spacing=1)[source]¶ Returns a equally spaced grid to plot udata
- Parameters
- udata
- spacing: spacing of the grid
- Returns
- xx, yy, (zz): 2D or 3D numpy arrays
-
ilpm.velocity.get_equally_spaced_kgrid(udata, dx=1)[source]¶ Returns a equally spaced grid to plot FFT of udata
- Parameters
- udata
- dx: spacing of the grid in the real space
- Returns
- kxx, kyy, (kzz): 2D or 3D numpy arrays
-
ilpm.velocity.kolmogorov_53(k, k0=50)[source]¶ Customizable Kolmogorov Energy spectrum Returns the value(s) of k0 * k^{-5/3}
- Parameters
- k: array-like, wavenumber: convention is k= 1/L NOT 2pi/L
- k0: float, coefficient
- Returns
- e_k: power-law spectrum with exponent -5/3 for a given k and k0
-
ilpm.velocity.kolmogorov_53_uni(k, epsilon, c=1.6)[source]¶ Universal Kolmogorov Energy spectrum Returns the value(s) of C epsilon^{2/3} k^{-5/3}
- Parameters
- k: array-like, wavenumber
- epsilon: float, dissipation rate
- c: float, Kolmogorov constant c=1.6 (default)
- … E(k) = c epsilon^(2/3) k^(-5/3)
- … E11(k) = c1 epsilon^(2/3) k^(-5/3)
- … E22(k) = c2 epsilon^(2/3) k^(-5/3)
- … c1:c2:c = 1: 4/3: 55/18
- … If c = 1.6, c1 = 0.491, c2 = 1.125
- … Exp. values: c = 1.5, c1 = 0.5, c2 = ?
- Returns
- e_k: array-like, Kolmogorov energy spectrum for a given range of k
-
ilpm.velocity.compute_kolmogorov_lengthscale_simple(epsilon, nu)[source]¶ Return Kolmogorov length scale for a given set of dissipation rate and viscosity Parameters ———- epsilon: float, dissipation rate nu: float, viscosity
- Returns
- float, Kolmogorov length scale
-
ilpm.velocity.get_characteristic_velocity(udata)[source]¶ Return 1-component RMS velocity, u’ energy = dim / 2 * u’^2
- Parameters
- udata
- Returns
- u_irms
-
ilpm.velocity.klonecker_delta(i, j)[source]¶ Klonecker Delta function: delta_{i, j}
- Parameters
- i: tensor index
- j: tensor index
- Returns
- 1 if i==j, 0 otherwise
-
ilpm.velocity.cart2pol(x, y)[source]¶ Transformation: Cartesian coord to polar coord
- Parameters
- x: numpy array
- y: numpy array
- Returns
- r: numpy array
- phi: numpy array
-
ilpm.velocity.cart2sph(x, y, z)[source]¶ Transformation: cartesian to spherical
- Parameters
- x
- y
- z
- Returns
- r: radius
- theta: elevetaion angle [-pi/2, pi/2]
- phi: azimuthal angle [-pi, pi]
-
ilpm.velocity.natural_sort(arr)[source]¶ natural-sorts elements in a given array alist.sort(key=natural_keys) sorts in human order http://nedbatchelder.com/blog/200712/human_sorting.html (See Toothy’s implementation in the comments)
e.g.- arr = [‘a28’, ‘a01’, ‘a100’, ‘a5’] … WITHOUT natural sorting,
-> [‘a01’, ‘a100’, ‘a28’, ‘a5’]
- … WITH natural sorting,
-> [‘a01’, ‘a5’, ‘a28’, ‘a100’]
- Parameters
- arr: list or numpy array of strings
- Returns
- sorted_array: natural-sorted
