diff --git a/.gitignore b/.gitignore
index 78cc706149dc24fd17208937694f3ac6deb4caf6..1e37b3381e243aacfcddd94599c6055bfa310ede 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,4 +22,3 @@ wetb/prepost/tests/data/demo_dlc/remote*
/wetb/fatigue_tools/rainflowcounting/compile.py
/docs/api
/htmlcov
-/wetb/wind/tests/test_mann_parameters.py
diff --git a/wetb/wind/turbulence/mann_parameters.py b/wetb/wind/turbulence/mann_parameters.py
index 995afdeff0abd8d146ab356ace9010a83eb04f33..f4e5f29b661e0c1dad756394d936e33d91d90bf9 100644
--- a/wetb/wind/turbulence/mann_parameters.py
+++ b/wetb/wind/turbulence/mann_parameters.py
@@ -11,7 +11,8 @@ from scipy.interpolate import RectBivariateSpline
import numpy as np
from wetb.wind.turbulence.spectra import spectra, logbin_spectra, \
- plot_spectrum
+ plot_spectrum, plot_spectra
+from wetb.utils.test_files import move2test_files
@@ -24,12 +25,7 @@ RBS2 = RectBivariateSpline(xp, yp, sp2)
RBS3 = RectBivariateSpline(xp, yp, sp3)
RBS4 = RectBivariateSpline(xp, yp, sp4)
-def estimate_mann_parameters(sf, u, v, w=None):
- if isinstance(u, (list, tuple)):
- #return fit_mann_model_spectra(*mean_spectra(sf, u, v, w))
- raise NotImplementedError
- else:
- return fit_mann_model_spectra(*spectra(sf, u, v, w))
+
#def mean_spectra(fs, u_ref_lst, u_lst, v_lst=None, w_lst=None):
# if isinstance(fs, (int, float)):
@@ -145,11 +141,13 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
from scipy.optimize import fmin
x = fmin(_local_error, start_vals_for_optimisation, logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size, min_bin_count), disp=False)
- if plt is not False:
+ if plt:
if not hasattr(plt, 'plot'):
import matplotlib.pyplot as plt
- _plot_spectra(k1, uu, vv, ww, uw, plt=plt)
- plot_mann_spectra(*x, plt=plt)
+# plot_spectra(k1, uu, vv, ww, uw, plt=plt)
+# plot_mann_spectra(*x, plt=plt)
+ ae, L, G = x
+ plot_fit(ae, L, G, k1, uu, vv, ww, uw, log10_bin_size=log10_bin_size, plt=plt)
plt.title('ae:%.3f, L:%.1f, G:%.2f' % tuple(x))
plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
@@ -254,26 +252,13 @@ def fit_ae(sf, u, L, G, min_bin_count=None, plt=False):
return ae
-def plot_spectra(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
+def plot_fit(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
# if plt is None:
# import matplotlib.pyplot as plt
- _plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, plt)
+ plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, plt)
plot_mann_spectra(ae, L, G, "-", mean_u, plt)
-def _plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
- if plt is None:
- import matplotlib.pyplot as plt
- bk1, buu, bvv, bww, buw = logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size)
- def plot(xx, label, color, plt):
- plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , '.' + color, label=label)
- plot(buu, 'uu', 'r', plt)
- plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
- plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
- if (bvv) is not None:
- plot(bvv, 'vv', 'g', plt)
- if bww is not None:
- plot(bww, 'ww', 'b', plt)
- plot(buw, 'uw', 'm', plt)
+
def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv', 'ww', 'uw']):
if plt is None:
@@ -312,7 +297,8 @@ if __name__ == "__main__":
nx = 8192
ny, nz = 8, 8
sf = (nx / l)
- fn = os.path.dirname(wind.__file__)+"/tests/test_files/turb/h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat"
+ #fn = os.path.dirname(wind.__file__)+"/tests/test_files/turb/h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat"
+ #fn = os.path.dirname(wind.__file__)+"/tests/test_files/turb/h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat"
u, v, w = [np.fromfile(fn % uvw, np.dtype('<f'), -1).reshape(nx , ny * nz) for uvw in ['u', 'v', 'w']]
ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w), plt=plt)
print (ae, L, G)
diff --git a/wetb/wind/turbulence/spectra.py b/wetb/wind/turbulence/spectra.py
index 40aa0ac3f0d3db47521098ebd2f0825adcae6051..b2376853da9e92791bd0019a7c6d3dfa472b4e05 100644
--- a/wetb/wind/turbulence/spectra.py
+++ b/wetb/wind/turbulence/spectra.py
@@ -35,20 +35,21 @@ def spectrum(x, y=None, k=1):
fftx = fftx[:len(fftx) // 2 ] * 2 # positive half * 2
- if len(fftx.shape) == 2:
- fftx = np.mean(fftx, 1)
+# if len(fftx.shape) == 2:
+# fftx = np.mean(fftx, 1)
return np.real(fftx * len(x) / (2 * k))[1:]
-def spectra(spacial_frq, u, v=None, w=None, detrend=True):
+def spectra(spatial_resolution, u, v=None, w=None, detrend=True):
"""Return the wave number, the uu, vv, ww autospectra and the uw cross spectra
Parameters
----------
- spacial_frq : inf or float
- - For time series: Sample frequency, see Notes\n
- - For turbulence boxes: 1/dx where
- dx is nx/Lx (number of grid points in length direction / box length in meters)
+ spatial_resolution : int, float or array_like
+ Distance between samples in meters
+ - For turbulence boxes: 1/dx = Lx/Nx where dx is distance between points,
+ Nx is number of points and Lx is box length in meters
+ - For time series: Sample frequency / U
u : array_like
The u-wind component\n
- if shape is (r,): One time series with *r* observations\n
@@ -74,36 +75,75 @@ def spectra(spacial_frq, u, v=None, w=None, detrend=True):
- uw: The uw cross spectrum, only if w is provided\n
For two dimensional input, the mean spectra are returned
-
- Notes
- -----
- If the 'mean wind speed' or 'shear compensated reference wind speed' are subtracted from u in advance,
- then the spacial_frq must be divided by the subtracted wind speed
"""
- assert isinstance(spacial_frq, (int, float))
-
- k = 2 * np.pi * spacial_frq
-
- if np.mean(u) > 1:
- k /= np.mean(u)
- u = u - np.mean(u, 0)
- if v is not None:
- v = v - np.mean(v, 0)
- if w is not None:
- w = w - np.mean(w, 0)
+ assert isinstance(spatial_resolution, (int, float))
+
+ k = 2 * np.pi * spatial_resolution
+ if v is not None: assert u.shape==v.shape
+ if w is not None: assert u.shape==w.shape
+
+ if 1 and len(u.shape)==2:
+ assert np.abs(np.mean(u,0)).max()<1
+ if v is not None: assert np.abs(np.mean(v,0)).max()<1
+ if w is not None: assert np.abs(np.mean(w,0)).max()<1
+ if isinstance(k, float):
+ k = np.repeat(k,u.shape[1])
+ else:
+ assert u.shape[1]==k.shape[0]
+ k1_vec = np.array([np.linspace(0, k_ / 2, len(u) / 2)[1:] for k_ in k]).T
+ else:
+ assert np.abs(np.mean(u))<1
+ if v is not None: assert np.abs(np.mean(v))<1
+ if w is not None: assert np.abs(np.mean(w))<1
+ assert isinstance(k,float)
+ k1_vec = np.linspace(0, k / 2, len(u) / 2)[1:]
if detrend:
u, v, w = detrend_wsp(u, v, w)
- k1_vec = np.linspace(0, k / 2, len(u) / 2)[1:]
return [k1_vec] + [spectrum(x1, x2, k=k) for x1, x2 in [(u, u), (v, v), (w, w), (w, u)]]
+def spectra_from_time_series(sample_frq, Uvw_lst):
+ """Return the wave number, the uu, vv, ww autospectra and the uw cross spectra
+
+ Parameters
+ ----------
+ sample_frq : int, float or array_like
+ Sample frequency
+ Uvw_lst : array_like
+ list of U, v and w, [(U1,v1,w1),(U2,v2,w2)...]
+
+ Returns
+ -------
+ k1, uu[, vv[, ww, uw]] : 2-5 x array_like
+ - k1: wave number
+ - uu: The uu auto spectrum\n
+ - vv: The vv auto spectrum, only if v is provided\n
+ - ww: The ww auto spectrum, only if w is provided\n
+ - uw: The uw cross spectrum, only if w is provided\n
+
+ For two dimensional input, the mean spectra are returned
+ """
+ assert isinstance(sample_frq, (int, float))
+ U,v,w = [np.array(Uvw_lst)[:,i,:].T for i in range(3)]
+ k = 2 * np.pi * sample_frq / U.mean(0)
+
+ if v is not None: assert np.abs(np.nanmean(v,0)).max()<1, "Max absolute mean of v is %f"%np.abs(np.nanmean(v,0)).max()
+ if w is not None: assert np.abs(np.nanmean(w,0)).max()<1
+ k1_vec = np.array([np.linspace(0, k_ / 2, U.shape[0] / 2)[1:] for k_ in k]).T
+ u = U-np.nanmean(U, 0)
+ u, v, w = detrend_wsp(u, v, w)
+
+ return [k1_vec] + [spectrum(x1, x2, k=k) for x1, x2 in [(u, u), (v, v), (w, w), (w, u)]]
+
def bin_spectrum(x, y, bin_size, min_bin_count=2):
assert min_bin_count > 0
x = x / bin_size
- low, high = np.floor(np.min(x)), np.ceil(np.max(x))
+ low, high = np.floor(np.nanmin(x)), np.ceil(np.nanmax(x))
bins = int(high - low)
nbr_in_bins = np.histogram(x, bins, range=(low, high))[0]
+ if len(x.shape)==2:
+ min_bin_count*=x.shape[1]
mask = nbr_in_bins >= min_bin_count
return np.histogram(x, bins, range=(low, high), weights=y)[0][mask] / nbr_in_bins[mask], nbr_in_bins
@@ -138,3 +178,22 @@ def detrend_wsp(u, v=None, w=None):
dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend
return dwsp.reshape(wsp.shape)
return [_detrend(wsp) for wsp in [u, v, w]]
+
+
+def plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None, marker_style='.'):
+ if plt is None:
+ import matplotlib.pyplot as plt
+ bk1, buu, bvv, bww, buw = logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size)
+ def plot(xx, label, color, plt):
+ plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , marker_style + color, label=label)
+ plot(buu, 'uu', 'r', plt)
+ plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
+ if mean_u ==1:
+ plt.ylabel(r'Spectral density $k_{1} F(k_{1}) [m^2/s^2]$')
+ else:
+ plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [-]$')
+ if (bvv) is not None:
+ plot(bvv, 'vv', 'g', plt)
+ if bww is not None:
+ plot(bww, 'ww', 'b', plt)
+ plot(buw, 'uw', 'm', plt)
\ No newline at end of file