From 6e46d5780823a96a8fe82a27b711b67d38cfe195 Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Wed, 1 Nov 2017 11:06:06 +0100
Subject: [PATCH] stuff
---
wetb/hawc2/htc_extensions.py | 11 ++-
wetb/hawc2/htc_file.py | 3 +
wetb/hawc2/shear_file.py | 17 +++-
wetb/hawc2/tests/test_shear_file.py | 2 +-
wetb/utils/caching.py | 16 +++-
wetb/utils/tests/test_caching.py | 25 ++++-
wetb/wind/__init__.py | 1 -
wetb/wind/shear.py | 26 ++++--
wetb/wind/stability.py | 1 -
wetb/wind/turbulence/mann_parameters.py | 117 ++++++++++++++++--------
wetb/wind/turbulence/spectra.py | 2 +-
11 files changed, 156 insertions(+), 65 deletions(-)
diff --git a/wetb/hawc2/htc_extensions.py b/wetb/hawc2/htc_extensions.py
index 67b9f7f..72a5900 100644
--- a/wetb/hawc2/htc_extensions.py
+++ b/wetb/hawc2/htc_extensions.py
@@ -15,7 +15,8 @@ from builtins import int
from builtins import str
from future import standard_library
import os
-from wetb.wind.shear import log_shear, power_shear
+
+
standard_library.install_aliases()
@@ -82,10 +83,11 @@ class HTCDefaults(object):
else:
mann.add_line('create_turb_parameters', [L, ae23, Gamma, seed, int(high_frq_compensation)], "L, alfaeps, gamma, seed, highfrq compensation")
if filenames is None:
- fmt = "mann_l%.1f_ae%.2f_g%.1f_h%d_%dx%dx%d_%.3fx%.2fx%.2f_s%04d%c.turb"
+
import numpy as np
dxyz = tuple(np.array(box_dimension) / no_grid_points)
- filenames = ["./turb/" + fmt % ((L, ae23, Gamma, high_frq_compensation) + no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]
+ from wetb.wind.turbulence import mann_turbulence
+ filenames = ["./turb/" + mann_turbulence.name_format % ((L, ae23, Gamma, high_frq_compensation) + no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]
if isinstance(filenames, str):
filenames = ["./turb/%s_s%04d%s.bin" % (filenames, seed, c) for c in ['u', 'v', 'w']]
for filename, c in zip(filenames, ['u', 'v', 'w']):
@@ -153,8 +155,9 @@ class HTCExtensions(object):
z0 = -self.wind.center_pos0[2]
wsp = self.wind.wsp[0]
if shear_type==1: #constant
- return lambda z : parameter
+ return lambda z : wsp
elif shear_type==3:
+ from wetb.wind.shear import power_shear
return power_shear(parameter, z0, wsp)
else:
raise NotImplementedError
\ No newline at end of file
diff --git a/wetb/hawc2/htc_file.py b/wetb/hawc2/htc_file.py
index f783408..655b740 100644
--- a/wetb/hawc2/htc_file.py
+++ b/wetb/hawc2/htc_file.py
@@ -66,6 +66,9 @@ class HTCFile(HTCContents, HTCDefaults, HTCExtensions):
This is a comment
>>> print (hydro_element.wave_breaking[0]) # first value
2
+
+ #Delete element
+ del htc.simulation.logfile #Delete logfile line. Raise keyerror if not exists
"""
filename = None
diff --git a/wetb/hawc2/shear_file.py b/wetb/hawc2/shear_file.py
index 11d6d60..0b2560f 100644
--- a/wetb/hawc2/shear_file.py
+++ b/wetb/hawc2/shear_file.py
@@ -31,7 +31,7 @@ class ShearFile(object):
>>> sf.save('test.dat')
"""
- def __init__(self,v_positions, w_positions, u=None, v=None, w=None, shear=None):
+ def __init__(self,v_positions, w_positions, u=None, v=None, w=None, u0=None, shear=None):
"""
Parameters
----------
@@ -48,6 +48,10 @@ class ShearFile(object):
w : array_like, optional
shear_w component, normalized with U_mean\n
shape must be (#w_positions, #v_positions) or (#w_positions,)
+ u0 : int or float
+ Mean wind speed
+ shear : function
+ Shear function: f(height)->wsp
"""
self.v_positions = v_positions
self.w_positions = w_positions
@@ -66,9 +70,10 @@ class ShearFile(object):
assert uvw[i].shape == shape, (i, uvw[i].shape, shape)
self.u, self.v, self.w = uvw
+ self.u0 = u0
self.shear = shear
- def uvw(self, v,w,shear=None):
+ def uvw(self, v,w,u0=None,shear=None):
"""Calculate u,v,w wind speeds at position(s) (v,w)
Parameters
----------
@@ -76,6 +81,8 @@ class ShearFile(object):
v-coordinate(s)
w : int, float or array_like
w-coordinates(s)
+ u0 : int or float
+ mean wind speed
shear : function or None
if function: f(height)->wsp
if None: self.shear is used, if not None.
@@ -86,10 +93,11 @@ class ShearFile(object):
u,v,w
wind speed(s) or wind speed factor(s) if shear not defined
"""
- shear = shear or self.shear or (lambda z: 1)
+ u0 =u0 or self.u0 or 1
+ shear = shear or self.shear or (lambda z: 0)
from scipy.interpolate import RegularGridInterpolator
wv = np.array([w,v]).T
- return [RegularGridInterpolator((self.w_positions, self.v_positions), uvw)(wv)*shear(w)
+ return [(RegularGridInterpolator((self.w_positions, self.v_positions), uvw)(wv))*u0+shear(w)
for uvw in [self.u, self.v, self.w] if uvw is not None]
@@ -156,6 +164,7 @@ class ShearFile(object):
user_defined_shear_filename = os.path.join(htc_file.modelpath, htc_file.wind.user_defined_shear[0])
shear_file = ShearFile.load(user_defined_shear_filename)
shear_file.shear = htc_file.get_shear()
+ shear_file.u0 = htc_file.wind.wsp[0]
return shear_file
def save(filename, v_coordinates, w_coordinates, u=None, v=None, w=None):
diff --git a/wetb/hawc2/tests/test_shear_file.py b/wetb/hawc2/tests/test_shear_file.py
index 302cae1..f724ca3 100644
--- a/wetb/hawc2/tests/test_shear_file.py
+++ b/wetb/hawc2/tests/test_shear_file.py
@@ -93,7 +93,7 @@ class TestShearFile(unittest.TestCase):
shear_file = ShearFile.load_from_htc(tfp+"htcfiles/test.htc")
np.testing.assert_array_equal(shear_file.w_positions, [30,100,160])
- np.testing.assert_array_almost_equal(shear_file.uvw([0,-55],[65,65])[0],np.array([.85,.9])*8.860807038)
+ np.testing.assert_array_almost_equal(shear_file.uvw([0,-55],[65,65])[0],np.array([.85,.9])*10+8.860807038)
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.test_shearfile']
unittest.main()
diff --git a/wetb/utils/caching.py b/wetb/utils/caching.py
index a4333a2..d8cd22e 100644
--- a/wetb/utils/caching.py
+++ b/wetb/utils/caching.py
@@ -10,9 +10,10 @@ from __future__ import absolute_import
from future import standard_library
import sys
from collections import OrderedDict
+import os
standard_library.install_aliases()
import inspect
-
+import numpy as np
def set_cache_property(obj, name, get_func, set_func=None):
"""Create a cached property
@@ -106,4 +107,15 @@ class cache_method():
if len(self.cache_dict)>self.N:
self.cache_dict.popitem(last=False)
return self.cache_dict[arg_id]
- return wrapped
\ No newline at end of file
+ return wrapped
+
+def cache_binary(f):
+ def wrap(filename,*args,**kwargs):
+ np_filename = os.path.splitext(filename)[0] + ".npy"
+ if os.path.isfile(np_filename) and (not os.path.isfile(filename) or os.path.getmtime(np_filename) > os.path.getmtime(filename)):
+ return np.load(np_filename)
+ else:
+ res = f(filename,*args,**kwargs)
+ np.save(np_filename,res)
+ return res
+ return wrap
\ No newline at end of file
diff --git a/wetb/utils/tests/test_caching.py b/wetb/utils/tests/test_caching.py
index 94530e7..23655a0 100644
--- a/wetb/utils/tests/test_caching.py
+++ b/wetb/utils/tests/test_caching.py
@@ -9,16 +9,18 @@ from __future__ import division
from __future__ import absolute_import
from future import standard_library
from collections import OrderedDict
+import os
standard_library.install_aliases()
import multiprocessing
import time
import unittest
-
+import numpy as np
from wetb.utils.timing import get_time
-from wetb.utils.caching import cache_function, set_cache_property, cache_method
-
+from wetb.utils.caching import cache_function, set_cache_property, cache_method,\
+ cache_binary
+tfp = os.path.dirname(__file__) + "/test_files/"
class Example(object):
def __init__(self, *args, **kwargs):
object.__init__(self, *args, **kwargs)
@@ -54,7 +56,9 @@ class Example(object):
time.sleep(1)
return x*2
-
+@cache_binary
+def open_csv(filename):
+ return np.loadtxt(filename)
def f(x):
return x ** 2
@@ -110,6 +114,19 @@ class TestCacheProperty(unittest.TestCase):
t = time.time()
e.test_cache_property
self.assertAlmostEqual(time.time()-t, 0, places=1)
+
+ def test_cache_binary(self):
+ if os.path.isfile(tfp+"test.npy"):
+ os.remove(tfp+'test.npy')
+ A = open_csv(tfp + "test.csv")
+ self.assertTrue(os.path.isfile(tfp+"test.npy"))
+ np.testing.assert_array_equal(A,np.loadtxt(tfp + "test.csv"))
+ A[0]=-1
+ np.save(tfp+"test.npy",A)
+ B = open_csv(tfp + "test.csv")
+ np.testing.assert_array_equal(A,B)
+ os.remove(tfp+'test.npy')
+
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()
diff --git a/wetb/wind/__init__.py b/wetb/wind/__init__.py
index 36e41c9..e69de29 100644
--- a/wetb/wind/__init__.py
+++ b/wetb/wind/__init__.py
@@ -1 +0,0 @@
-from .turbulence import *
diff --git a/wetb/wind/shear.py b/wetb/wind/shear.py
index dd3281b..de1acbb 100644
--- a/wetb/wind/shear.py
+++ b/wetb/wind/shear.py
@@ -114,7 +114,7 @@ def fit_power_shear_ref(z_u_lst, z_ref, plt=None):
-def log_shear(u_star, z0, z, Obukhov_length=None):
+def log_shear(u_star, z0):
"""logarithmic shear
Parameters
@@ -123,24 +123,28 @@ def log_shear(u_star, z0, z, Obukhov_length=None):
Friction velocity
z0 : int or float
Surface roughness [m]
- z : int, float or array_like
- The heights of interest
Returns
-------
- log_shear : array-like
- Wind speeds at height z
+ log_shear : function f(z,L=None) -> wsp
+ z : int, float or array_like
+ The heights of interest
+ L : int, float or array_like, optional
+ The corresponding Monin-Obukhov length
Example
--------
- >>> log_shear(1, 10, [20, 50, 70])
+ >>> shear = log_shear(1, 10)
+ >>> shear([20, 50, 70])
[ 1.73286795 4.02359478 4.86477537]
"""
K = 0.4 # von Karmans constant
- if Obukhov_length is None:
- return u_star / K * (np.log(np.array(z) / z0))
- else:
- return u_star / K * (np.log(z / z0) - stability_term(z, z0, Obukhov_length))
+ def log_shear(z,Obukhov_length=None):
+ if Obukhov_length is None:
+ return u_star / K * (np.log(np.array(z) / z0))
+ else:
+ return u_star / K * (np.log(z / z0) - stability_term(z, z0, Obukhov_length))
+ return log_shear
def stability_term(z, z0, L):
"""Calculate the stability term for the log shear
@@ -169,6 +173,8 @@ def fit_log_shear(z_u_lst, include_R=False):
- u_z1: Wind speeds or mean wind speed at z1
- z2: another height
- u_z2: Wind speeds or mean wind speeds at z2
+ include_R : boolean, optional
+ If True, the R^2 error is returned
Returns
-------
diff --git a/wetb/wind/stability.py b/wetb/wind/stability.py
index 6826d58..b60a8d3 100644
--- a/wetb/wind/stability.py
+++ b/wetb/wind/stability.py
@@ -32,7 +32,6 @@ def MoninObukhov_length(u,v,w, T):
v = v-np.mean(v)
w = w-np.mean(w)
u_star = (np.mean(u*w)**2+np.mean(v*w)**2)**(1/4)
- t = T-T.mean()
wT = np.mean(w*T)
return -u_star ** 3 * (T.mean() + 273.15) / (K * g * wT)
diff --git a/wetb/wind/turbulence/mann_parameters.py b/wetb/wind/turbulence/mann_parameters.py
index 0dd3c3f..9799b49 100644
--- a/wetb/wind/turbulence/mann_parameters.py
+++ b/wetb/wind/turbulence/mann_parameters.py
@@ -10,9 +10,9 @@ import os
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
+from wetb.wind.turbulence.spectra import spectra, logbin_spectra, plot_spectra,\
+ detrend_wsp
+
@@ -198,13 +198,19 @@ 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 fit_ae2var(variance, L, G):
+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
@@ -226,13 +232,17 @@ def fit_ae2var(variance, L, G):
return ae
-def fit_ae(sf, u, L, G, min_bin_count=None, plt=False):
+
+def fit_ae(spatial_resolution, u, L, G, plt=False):
"""Fit alpha-epsilon to match variance of time series
Parameters
----------
- sf : int or float
- Sample frequency
+ 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
@@ -245,33 +255,37 @@ def fit_ae(sf, u, L, G, min_bin_count=None, plt=False):
ae : float
Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
"""
-# def get_var(k1, ae):
-# uu = get_mann_model_spectra(ae, L, G, k1)[0]
-# return np.trapz(2 * uu, k1)
-# return (np.sum(uu) * 2 * (k1[1] - k1[0]))
- 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[1] / 2)
-
+ #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 = 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)
+ 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, uu = logbin_spectra(*spectra(sf, u), min_bin_count=1)[:2]
+ k1_lb, uu_lb = logbin_spectra(*spectra(sf, u), min_bin_count=1)[:2]
- plt.semilogx(k1, k1 * uu, 'r--', label='uu_logbin')
+ 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()
@@ -312,22 +326,51 @@ if __name__ == "__main__":
ds = gtsdf.Dataset(os.path.dirname(wind.__file__)+"/tests/test_files/WspDataset.hdf5")#'unit_test/test_files/wspdataset.hdf5')
f = 35
u, v = wsp_dir2uv(ds.Vhub_85m, ds.Dir_hub_)
-
+
u_ref = np.mean(u)
u -= u_ref
-
+
sf = f / u_ref
- ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v), plt = plt)
+ ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v), plt = None)
print (ae, L, G)
+ print (u.shape)
+ print (ds.Time[-1])
+ 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 ()
+ 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"""
- l = 16384
- 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"
- 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)
+# """Example of fitting Mann parameters to a "series" of a turbulence box"""
+# l = 16384
+# 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"
+# 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)
+
+
+
+# """Example of fitting Mann parameters to a "series" of a turbulence box"""
+# l = 4778.3936
+# nx = 8192
+# ny, nz = 32,32
+# 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"
+# for s in [1,2,3,4,5,6]:
+# fn = r'C:\mmpe\HAWC2\models\SWT3.6-107\turb/mann_l73.1_ae0.00_g2.0_h1_8192x32x32_0.583x3.44x3.44_s160%d%%s.turb'%s
+# u, v, w = [np.fromfile(fn % uvw, np.dtype('<f'), -1).reshape(nx , ny * nz) for uvw in ['u', 'v', 'w']]
+# u = np.fromfile(fn % 'u', np.dtype('<f'), -1).reshape(nx , ny * nz)
+# #ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w), plt=plt)
+# #print (fit_ae(sf, u, 73.0730383576, 2.01636095317))
+# print (u.std())
+# #print (ae, L, G)
\ No newline at end of file
diff --git a/wetb/wind/turbulence/spectra.py b/wetb/wind/turbulence/spectra.py
index f6aa601..1a226c5 100644
--- a/wetb/wind/turbulence/spectra.py
+++ b/wetb/wind/turbulence/spectra.py
@@ -49,7 +49,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 = Lx/Nx 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
--
GitLab