Commit 74d51824 authored by Mads M. Pedersen's avatar Mads M. Pedersen

Improve var2ae by including dt and T

parent b6709093
Pipeline #20206 passed with stages
in 2 minutes and 23 seconds
Subproject commit cf64c48ffb671d35a457d0539e078bfe3eb07f0b
Subproject commit 167b16d8f1a26161ebc5a258c62389777d133489
......@@ -87,7 +87,8 @@ def _local_error(x, k1, uu, vv, ww=None, uw=None):
return val
def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2, min_bin_count=2, start_vals_for_optimisation=(0.01, 50, 3.3), plt=False):
def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
min_bin_count=2, start_vals_for_optimisation=(0.01, 50, 3.3), plt=False):
"""Fit a mann model to the spectra
Bins the spectra, into logarithmic sized bins and find the mann model parameters,
......@@ -140,7 +141,7 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
# 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_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]$')
......@@ -192,7 +193,7 @@ def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2):
return np.sqrt(((bk1 * (sp_meas - sp_fit)) ** 2).mean(1))
def var2ae(variance, L, G):
def var2ae(variance, L, G, U, T=600, sample_frq=10, plt=False):
"""Fit alpha-epsilon to match variance of time series
Parameters
......@@ -203,14 +204,22 @@ def var2ae(variance, L, G):
Length scale of Mann model
G : int or float
Gamma of Mann model
U : int or float
Mean wind speed
T: int or float
Length [s] of signal, from which the variance is calculated
sample_frq: int or float
Sample frequency [Hz] of signal from which the variance is calculated
Returns
-------
ae : float
Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
Alpha epsilon^(2/3) of Mann model that makes the energy of the model in the
frequency range [1/length, sample_frq] equal to the variance of u
"""
k1 = np.logspace(1, 10, 1000) / 100000000
k_low, k_high = 2 * np.pi / (U * np.array([T, 1 / sample_frq]))
k1 = 10 ** (np.linspace(np.log10(k_low), np.log10(k_high), 1000))
def get_var(uu):
return np.trapz(2 * uu[:], k1[:])
......@@ -218,6 +227,14 @@ def var2ae(variance, L, G):
v1 = get_var(get_mann_model_spectra(0.1, L, G, k1)[0])
v2 = get_var(get_mann_model_spectra(0.2, L, G, k1)[0])
ae = (variance - v1) / (v2 - v1) * .1 + .1
if plt is not False:
if not hasattr(plt, 'plot'):
import matplotlib.pyplot as plt
muu = get_mann_model_spectra(ae, L, G, k1)[0]
plt.semilogx(k1, k1 * muu, label='ae:%.3f, L:%.1f, G:%.2f' % (ae, L, G))
plt.legend()
plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
return ae
......@@ -227,8 +244,8 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
Parameters
----------
spatial_resolution : int, float or array_like
Distance between samples in meters
- For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
Number of points pr meterDistance between samples in meters
- For turbulence boxes: 1/dx = Nx/Lx 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
......@@ -271,7 +288,7 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
if not hasattr(plt, 'plot'):
import matplotlib.pyplot as plt
plt.semilogx(k1, k1 * uu, 'b-', label='uu')
k1_lb, uu_lb = logbin_spectra(*spectra(sf, u), min_bin_count=1)[:2]
k1_lb, uu_lb = logbin_spectra(*spectra(spatial_resolution, u), min_bin_count=1)[:2]
plt.semilogx(k1_lb, k1_lb * uu_lb, 'r--', label='uu_logbin')
muu = get_mann_model_spectra(ae, L, G, k1)[0]
......@@ -329,11 +346,11 @@ if __name__ == "__main__":
plt.plot(u)
plt.plot(detrend_wsp(u)[0])
plt.show()
print(fit_ae(sf, detrend_wsp(u)[0], L, G, plt))
print(var2ae(detrend_wsp(u)[0].var(), L, G,))
print(fit_ae(sf, detrend_wsp(u)[0], L, G, plt))
print(var2ae(detrend_wsp(u)[0].var(), L, G,))
print()
print(fit_ae(sf, u[:21000], L, G))
print(var2ae(u[:21000].var(), L, G,))
print(fit_ae(sf, u[:21000], L, G))
print(var2ae(u[:21000].var(), L, G,))
# """Example of fitting Mann parameters to a "series" of a turbulence box"""
......
......@@ -53,7 +53,7 @@ def spectra(spatial_resolution, u, v=None, w=None, detrend=True):
----------
spatial_resolution : int, float or array_like
Distance between samples in meters
- For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
- For turbulence boxes: 1/dx = Nx/Lx 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
......@@ -91,11 +91,11 @@ def spectra(spatial_resolution, u, v=None, w=None, detrend=True):
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
# assert np.abs(np.mean(u, 0)).max() < 2
# 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:
......@@ -146,15 +146,13 @@ def spectra_from_time_series(sample_frq, Uvw_lst):
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
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
......@@ -194,7 +192,8 @@ def detrend_wsp(u, v=None, w=None):
t = np.arange(dwsp.shape[0])
A = np.vstack([t, np.ones(len(t))]).T
for i in range(dwsp.shape[1]):
trend, offset = np.linalg.lstsq(A, dwsp[:, i], rcond=None)[0]
m = ~np.isnan(dwsp[:, i])
trend, offset = np.linalg.lstsq(A[:m.sum()], dwsp[:, i][m], rcond=None)[0]
dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend
return dwsp.reshape(wsp.shape)
return [_detrend(wsp) for wsp in [u, v, w]]
......
'''
Created on 05/06/2012
@author: Mads
'''
import os
import sys
import unittest
import matplotlib.pyplot as plt
import numpy as np
from wetb.utils.test_files import get_test_file, move2test_files
from wetb.wind.turbulence.mann_turbulence import load_uvw, fit_mann_parameters,\
fit_mann_parameters_from_time_series
from wetb.utils.timing import print_time
from wetb.wind.turbulence.spectra import spectra_from_time_series, spectra,\
plot_spectra, logbin_spectra
import warnings
from wetb.wind.turbulence import mann_turbulence
from wetb.wind.turbulence.mann_parameters import var2ae, fit_ae
from numpy import spacing
tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
class TestMannTurbulence(unittest.TestCase):
def test_fit_mann_parameters_turbulence_box(self):
# for uvw in 'uvw':
# move2test_files(tfp + 'h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat'%uvw)
fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
u, v, w = load_uvw(fn_lst, (8192, 8, 8))
dx = 16384 / 8192
fx = 1 / dx # spatial resolution
plt = None
ae, L, G = fit_mann_parameters(fx, u, v, w, plt=plt)
self.assertAlmostEqual(ae, .15, delta=0.01)
self.assertAlmostEqual(L, 10, delta=0.3)
self.assertAlmostEqual(G, 3.3, delta=0.06)
def test_spectra_from_timeseries(self):
fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
u, v, w = load_uvw(fn_lst, (8192, 8, 8))
dx = 16384 / 8192
fx = 1 / dx # spatial resolution
k1, uu, vv, ww, uw = logbin_spectra(*spectra(fx, u, v, w))
U = u + 4
sample_frq = 2
k12, uu2, vv2, ww2, uw2 = logbin_spectra(
*spectra_from_time_series(sample_frq, [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]))
np.testing.assert_allclose(uu, uu2, 0.02)
U = u + 8
sample_frq = 2
k13, uu3, vv3, ww3, uw3 = logbin_spectra(
*spectra_from_time_series(sample_frq, [(U_[::2], v_[::2], w_[::2]) for U_, v_, w_ in zip(U.T, v.T, w.T)]))
np.testing.assert_allclose(uu[:-3], uu3[:-2], 0.1)
# One set of time series with U=4
U = u + 4
Uvw_lst = [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]
# Another set of time series with U=8 i.e. only every second point to have
# same sample_frq. (nan added to have same length)
U = u + 4
Uvw_lst.extend([(np.r_[U_[::2], U_[::2] + np.nan], np.r_[v_[::2], v_[::2] + np.nan],
np.r_[w_[::2], w_[::2] + np.nan]) for U_, v_, w_ in zip(U.T, v.T, w.T)])
sample_frq = 2
with warnings.catch_warnings():
warnings.simplefilter("ignore")
k14, uu4, vv4, ww4, uw4 = logbin_spectra(*spectra_from_time_series(sample_frq, Uvw_lst))
np.testing.assert_allclose(uu[:-3], uu3[:-2], rtol=0.1)
if 0:
import matplotlib.pyplot as plt
plt.semilogx(k1, uu * k1)
plt.semilogx(k12, uu2 * k12)
plt.semilogx(k13, uu3 * k13)
plt.semilogx(k14, uu4 * k14)
plt.show()
def test_fit_mann_parameters_from_timeseries(self):
fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
u, v, w = load_uvw(fn_lst, (8192, 8, 8))
dx = 16384 / 8192
fx = 1 / dx # spatial resolution
ae, L, G = fit_mann_parameters(fx, u, v, w)
self.assertAlmostEqual(ae, .15, delta=0.01)
self.assertAlmostEqual(L, 10, delta=0.3)
self.assertAlmostEqual(G, 3.3, delta=0.06)
#import matplotlib.pyplot as plt
plt = None
U = u + 4
sample_frq = 2
ae, L, G = fit_mann_parameters_from_time_series(
sample_frq, [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)], plt=plt)
self.assertAlmostEqual(ae, .15, delta=0.01)
self.assertAlmostEqual(L, 10, delta=0.3)
self.assertAlmostEqual(G, 3.3, delta=0.06)
# One set of time series with U=4
U = u + 4
Uvw_lst = [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]
# Another set of time series with U=8 i.e. only every second point to have
# same sample_frq. (nan added to have same length)
U = u + 4
Uvw_lst.extend([(np.r_[U_[::2], U_[::2] + np.nan], np.r_[v_[::2], v_[::2] + np.nan],
np.r_[w_[::2], w_[::2] + np.nan]) for U_, v_, w_ in zip(U.T, v.T, w.T)])
sample_frq = 2
ae, L, G = fit_mann_parameters_from_time_series(sample_frq, Uvw_lst, plt=plt)
self.assertAlmostEqual(ae, .15, delta=0.01)
self.assertAlmostEqual(L, 10, delta=0.3)
self.assertAlmostEqual(G, 3.3, delta=0.06)
def test_var2ae_U(self):
u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
dx = 2
for U in [1, 10, 100]:
# should be independent of U
dt = dx / U
T = 16384 / U
self.assertAlmostEqual(var2ae(variance=u.var(), L=10, G=3.3, U=U, T=T, sample_frq=1 / dt), .15, delta=.021)
def test_var2ae_T(self):
u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
dx = 2
U = 10
dt = dx / U
for i in np.arange(6):
# reshape to more and shorter series. Variance should decrease while ae should be roughly constant
n = 2**i
u_ = u.T.reshape((u.T.shape * np.array([n, 1 / n])).astype(int)).T
var = u_.var(0).mean()
ae = var2ae(variance=var, L=10, G=3.3, U=U, T=dx * u_.shape[0] / U, sample_frq=1 / dt)
self.assertAlmostEqual(ae, .15, delta=.025)
def test_var2ae_dt(self):
u = mann_turbulence.load(get_test_file("h2a16384_8_8_65536_32_32_0.15_40_4.0u.dat"), (16384, 8, 8))
dx = 4
U = 10
T = u.shape[0] * dx / U
for i in np.arange(9):
# average every neighbouring samples to decrease dt.
# Variance should decrease while ae should be roughly constant
n = 2**i
u_ = u.reshape(u.shape[0] // n, n, u.shape[1]).mean(1)
var = u_.var(0).mean()
ae = var2ae(variance=var, L=40, G=4, U=U, T=T, sample_frq=1 / (n * dx / U), plt=False)
#print(u_.shape, var, ae)
self.assertAlmostEqual(ae, .15, delta=.04)
def test_fit_ae2var(self):
u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
self.assertAlmostEqual(fit_ae(spatial_resolution=2, u=u, L=10, G=3.3), .15, delta=.02)
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment