Skip to content
Snippets Groups Projects
Commit 22583e0c authored by davcon's avatar davcon
Browse files

Upload New File

parent e5af02de
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
python_scripts/mann_functions.py 0 → 100644
+ 318
− 0
View file @ 22583e0c
'''
Created on 03/06/2014
@author: MMPE
@ adapt from wetb: DAVCON 05/02/2019
'''
import os
from scipy.interpolate import RectBivariateSpline
import numpy as np
from spectra import spectra, logbin_spectra, plot_spectra, detrend_wsp
# Look-Up Tables
sp1, sp2, sp3, sp4 = np.load("C:/Users/davcon/Desktop/0_WakeIndicators/Mann_model/mann_spectra_data.npy")
yp = np.arange(-3, 3.1, 0.1)
xp = np.arange(0, 5.1, 0.1)
RBS1 = RectBivariateSpline(xp, yp, sp1)
RBS2 = RectBivariateSpline(xp, yp, sp2)
RBS3 = RectBivariateSpline(xp, yp, sp3)
RBS4 = RectBivariateSpline(xp, yp, sp4)
def get_mann_model_spectra(ae, L, G, k1):
"""Mann model spectra
Parameters
----------
ae : int or float
Alpha epsilon^(2/3) of Mann model
L : int or float
Length scale of Mann model
G : int or float
Gamma of Mann model
k1 : array_like
Desired wave numbers
Returns
-------
uu : array_like
The u-autospectrum of the wave numbers, k1
vv : array_like
The v-autospectrum of the wave numbers, k1
ww : array_like
The w-autospectrum of the wave numbers, k1
uw : array_like
The u,w cross spectrum of the wave numbers, k1
"""
xq = np.log10(L * k1)
yq = (np.zeros_like(xq) + G)
f = L ** (5 / 3) * ae
uu = f * RBS1.ev(yq, xq)
vv = f * RBS2.ev(yq, xq)
ww = f * RBS3.ev(yq, xq)
uw = f * RBS4.ev(yq, xq)
return uu, vv, ww, uw
def _local_error(x, k1, uu, vv, ww=None, uw=None):
ae, L, G = x
val = 10 ** 99
if ae >= 0 and G >= 0 and G <= 5 and L > 0 and np.log10(k1[0] * L) >= -3 and np.log10(k1[0] * L) <= 3:
tmpuu, tmpvv, tmpww, tmpuw = get_mann_model_spectra(ae, L, G, k1)
val = np.sum((k1 * uu - k1 * tmpuu) ** 2)
if vv is not None:
val += np.sum((k1 * vv - k1 * tmpvv) ** 2)
if ww is not None:
val += np.sum((k1 * ww - k1 * tmpww) ** 2) + np.sum((k1 * uw - k1 * tmpuw) ** 2)
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):
"""Fit a mann model to the spectra
Bins the spectra, into logarithmic sized bins and find the mann model parameters,
that minimize the error between the binned spectra and the Mann model spectra
using an optimization function
Parameters
----------
k1 : array_like
Wave numbers
uu : array_like
The u-autospectrum of the wave numbers, k1
vv : array_like, optional
The v-autospectrum of the wave numbers, k1
ww : array_like, optional
The w-autospectrum of the wave numbers, k1
uw : array_like, optional
The u,w cross spectrum of the wave numbers, k1
log10_bin_size : int or float, optional
Bin size (log 10, based)
start_vals_for_optimization : (ae, L, G), optional
- ae: Alpha epsilon^(2/3) of Mann model\n
- L: Length scale of Mann model\n
- G: Gamma of Mann model
Returns
-------
ae : int or float
Alpha epsilon^(2/3) of Mann model
L : int or float
Length scale of Mann model
G : int or float
Gamma of Mann model
Examples
--------
>>> sf = sample_frq / u_ref
>>> u,v,w = # u,v,w wind components
>>> ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w))
>>> u1,v1 = # u,v wind components
>>> ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v))
"""
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 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)
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]$')
plt.legend()
plt.show()
return x
def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2):
"""Fit a mann model to the spectra
Bins the spectra, into logarithmic sized bins and find the mann model parameters,
that minimize the error between the binned spectra and the Mann model spectra
using an optimization function
Parameters
----------
ae : int or float
Alpha epsilon^(2/3) of Mann model
L : int or float
Length scale of Mann model
G : int or float
Gamma of Mann model
k1 : array_like
Wave numbers
uu : array_like
The u-autospectrum of the wave numbers, k1
vv : array_like, optional
The v-autospectrum of the wave numbers, k1
ww : array_like, optional
The w-autospectrum of the wave numbers, k1
uw : array_like, optional
The u,w cross spectrum of the wave numbers, k1
log10_bin_size : int or float, optional
Bin size (log 10, based)
start_vals_for_optimization : (ae, L, G), optional
- ae: Alpha epsilon^(2/3) of Mann model\n
- L: Length scale of Mann model\n
- G: Gamma of Mann model
Returns
-------
residual : array_like
rms of each spectrum
"""
k1_sp = np.array([sp for sp in logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size) if sp is not None])
bk1, sp_meas = k1_sp[0], k1_sp[1:]
sp_fit = np.array(get_mann_model_spectra(ae, L, G, bk1))[:sp_meas.shape[0]]
return np.sqrt(((bk1 * (sp_meas - sp_fit)) ** 2).mean(1))
def var2ae(variance, spatial_resolution, N, L, G):
"""Fit alpha-epsilon to match variance of time series
Parameters
----------
variance : array-like
variance of u vind component
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,
Nx is number of points and Lx is box length in meters
- For time series: Sample frequency / U
N : int
L : int or float
Length scale of Mann model
G : int or float
Gamma of Mann model
Returns
-------
ae : float
Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
"""
k1 = np.logspace(1,10,1000)/100000000
def get_var(uu):
return np.trapz(2 * uu[:], k1[:])
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
return ae
def fit_ae(spatial_resolution, u, L, G, plt=False):
"""Fit alpha-epsilon to match variance of time series
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,
Nx is number of points and Lx is box length in meters
- For time series: Sample frequency / U
u : array-like
u vind component
L : int or float
Length scale of Mann model
G : int or float
Gamma of Mann model
Returns
-------
ae : float
Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
"""
#if len(u.shape) == 1:
# u = u.reshape(len(u), 1)
# if min_bin_count is None:
# min_bin_count = max(2, 6 - u.shape[0] / 2)
# min_bin_count = 1
def get_var(k1, uu):
l = 0 #128 // np.sqrt(u.shape[1])
return np.trapz(2 * uu[l:], k1[l:])
k1, uu = spectra(spatial_resolution, u)[:2]
v = get_var(k1,uu)
v1 = get_var(k1, get_mann_model_spectra(0.1, L, G, k1)[0])
v2 = get_var(k1, get_mann_model_spectra(0.2, L, G, k1)[0])
ae = (v - v1) / (v2 - v1) * .1 + .1
# print (ae)
#
# k1 = spectra(sf, u)[0]
# v1 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.1, L, G, k1)[0], min_bin_count=min_bin_count)[:2])
# v2 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.2, L, G, k1)[0], min_bin_count=min_bin_count)[:2])
# k1, uu = logbin_spectra(*spectra(sf, u), min_bin_count=2)[:2]
# #variance = np.mean([detrend_wsp(u_)[0].var() for u_ in u.T])
# v = get_var(k1, uu)
# ae = (v - v1) / (v2 - v1) * .1 + .1
# print (ae)
if plt is not 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]
plt.semilogx(k1_lb, k1_lb * uu_lb, 'r--', label='uu_logbin')
muu = get_mann_model_spectra(ae, L, G, k1)[0]
plt.semilogx(k1, k1 * muu, 'g', 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}) [m^2/s^2]$')
plt.show()
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):
# if plt is None:
# import matplotlib.pyplot as 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_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv', 'ww', 'uw']):
if plt is None:
import matplotlib.pyplot as plt
mf = 10 ** (np.linspace(-4, 1, 1000))
muu, mvv, mww, muw = get_mann_model_spectra(ae, L, G, mf)
plt.title("ae: %.3f, L: %.2f, G:%.2f"%(ae,L,G))
if 'uu' in spectra: plt.semilogx(mf, mf * muu * 10 ** 0 / u_ref ** 2, 'r' + style)
if 'vv' in spectra: plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
if 'ww' in spectra: plt.semilogx(mf, mf * mww * 10 ** 0 / u_ref ** 2, 'b' + style)
if 'uw' in spectra: plt.semilogx(mf, mf * muw * 10 ** 0 / u_ref ** 2, 'm' + style)
#
#if __name__ == "__main__":
# import gtsdf
# from geometry import wsp_dir2uv
# #from wetb import wind
# import matplotlib.pyplot as plt
# # """Example of fitting Mann parameters to a "series" of a turbulence box"""
# l = 1800
# nx = 8192; ny, nz = 32,32; lx=0.2197265625;ly,lz=1,1;
# sf = (nx / l)
# s=1;
# fn = r'./turb100%d%%s.bin'%s
# 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= True)
#
# u, v, w = [np.fromfile(fn % uvw, np.dtype('<f'), -1).reshape(nx , ny,nz) for uvw in ['u', 'v', 'w'] ]
# # """ Y and Z vectors """
# Y = np.linspace(-ly/2., ly/2., ny)
# Z = np.linspace(-lz/2., lz/2., nz)
# plt.close('all')
# plt.ion()
#
\ No newline at end of file
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