From 43fc8cd93510073aea5ddd65666cbf7aefc18af1 Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Fri, 21 Jul 2017 08:49:49 +0200
Subject: [PATCH] updates in wetb.wind.turbulence
---
.gitignore | 1 +
wetb/wind/turbulence/mann_parameters.py | 40 ++++++---
wetb/wind/turbulence/spectra.py | 109 ++++++------------------
3 files changed, 53 insertions(+), 97 deletions(-)
diff --git a/.gitignore b/.gitignore
index 1e37b338..78cc7061 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,3 +22,4 @@ 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 f4e5f29b..995afdef 100644
--- a/wetb/wind/turbulence/mann_parameters.py
+++ b/wetb/wind/turbulence/mann_parameters.py
@@ -11,8 +11,7 @@ from scipy.interpolate import RectBivariateSpline
import numpy as np
from wetb.wind.turbulence.spectra import spectra, logbin_spectra, \
- plot_spectrum, plot_spectra
-from wetb.utils.test_files import move2test_files
+ plot_spectrum
@@ -25,7 +24,12 @@ 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)):
@@ -141,13 +145,11 @@ 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:
+ if plt is not False:
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)
- ae, L, G = x
- plot_fit(ae, L, G, k1, uu, vv, ww, uw, log10_bin_size=log10_bin_size, plt=plt)
+ _plot_spectra(k1, uu, vv, ww, uw, plt=plt)
+ plot_mann_spectra(*x, 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]$')
@@ -252,13 +254,26 @@ def fit_ae(sf, u, L, G, min_bin_count=None, plt=False):
return ae
-def plot_fit(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
+def plot_spectra(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:
@@ -297,8 +312,7 @@ 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 b2376853..40aa0ac3 100644
--- a/wetb/wind/turbulence/spectra.py
+++ b/wetb/wind/turbulence/spectra.py
@@ -35,21 +35,20 @@ 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(spatial_resolution, u, v=None, w=None, detrend=True):
+def spectra(spacial_frq, u, v=None, w=None, detrend=True):
"""Return the wave number, the uu, vv, ww autospectra and the uw cross spectra
Parameters
----------
- 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
+ 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)
u : array_like
The u-wind component\n
- if shape is (r,): One time series with *r* observations\n
@@ -75,75 +74,36 @@ def spectra(spatial_resolution, 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(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:]
+ 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)
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.nanmin(x)), np.ceil(np.nanmax(x))
+ low, high = np.floor(np.min(x)), np.ceil(np.max(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
@@ -178,22 +138,3 @@ 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
--
GitLab