diff --git a/wetb/dlc/high_level.py b/wetb/dlc/high_level.py
index 0edfa43b109cd4eba748e6ba881f2a663129a525..a7ad916e0e992545568c845a50d86efe7bd2851e 100644
--- a/wetb/dlc/high_level.py
+++ b/wetb/dlc/high_level.py
@@ -236,16 +236,20 @@ class DLCHighLevel(object):
pass
elif not hasattr(self, "res_folder") or self.res_folder == "":
files = glob.glob(os.path.join(self.res_path, "*"+ext)) + glob.glob(os.path.join(self.res_path, "*/*"+ext))
+ if len(files)==0:
+ raise Exception('No *%s files found in:\n%s or\n%s'%(ext, self.res_path, os.path.join(self.res_path, "*/")))
else:
files = []
-
for dlc_id in fatigue_dlcs:
dlc_id = str(dlc_id)
if "%" in self.res_folder:
folder = self.res_folder % dlc_id
else:
folder = self.res_folder
- files.extend(glob.glob(os.path.join(self.res_path , folder, "*"+ext)))
+ dlc_files = (glob.glob(os.path.join(self.res_path , folder, "*"+ext)))
+ if len(dlc_files)==0:
+ raise Exception('DLC%s included in fatigue analysis, but no *%s files found in:\n%s'%(dlc_id, ext, os.path.join(self.res_path , folder)))
+ files.extend(dlc_files)
keys = list(zip(*self.dist_value_keys))[1]
fmt = self.format_tag_value
tags = [[fmt(tag.replace(key, "")) for tag, key in zip(os.path.basename(f).split("_"), keys)] for f in files]
diff --git a/wetb/dlc/tests/test_files/res/dlc31_iec61400-1ed3/dlc31_wsp25_wdir000_s0000.int b/wetb/dlc/tests/test_files/res/dlc31_iec61400-1ed3/dlc31_wsp25_wdir000_s0000.int
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/wetb/fatigue_tools/fatigue.py b/wetb/fatigue_tools/fatigue.py
index 0aa87132fd3e5474d2cfda7b60c38fd8378ce0dd..85625df0f793e1c6254627bd8ef293c38b7f0f75 100644
--- a/wetb/fatigue_tools/fatigue.py
+++ b/wetb/fatigue_tools/fatigue.py
@@ -19,6 +19,7 @@ from __future__ import print_function
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
+import warnings
standard_library.install_aliases()
import numpy as np
from wetb.fatigue_tools.rainflowcounting import rainflowcount
@@ -156,11 +157,12 @@ def cycle_matrix(signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_win
if isinstance(ampl_bins, int):
ampl_bins = np.linspace(0, 1, num=ampl_bins + 1) * ampls[weights>0].max()
cycles, ampl_edges, mean_edges = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights)
-
- ampl_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * ampls)[0]
- ampl_bin_mean = np.nanmean(ampl_bin_sum / np.where(cycles,cycles,np.nan),1)
- mean_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * means)[0]
- mean_bin_mean = np.nanmean(mean_bin_sum / np.where(cycles, cycles, np.nan), 1)
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore")
+ ampl_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * ampls)[0]
+ ampl_bin_mean = np.nanmean(ampl_bin_sum / np.where(cycles,cycles,np.nan),1)
+ mean_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * means)[0]
+ mean_bin_mean = np.nanmean(mean_bin_sum / np.where(cycles, cycles, np.nan), 1)
cycles = cycles / 2 # to get full cycles
return cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges
diff --git a/wetb/fatigue_tools/rainflowcounting/rainflowcount.py b/wetb/fatigue_tools/rainflowcounting/rainflowcount.py
index bac7c262bbe9582c28961aaabea600b9fd6d9d6b..146cf9c07356b72fdc179ce6b0fb01997c247aa6 100644
--- a/wetb/fatigue_tools/rainflowcounting/rainflowcount.py
+++ b/wetb/fatigue_tools/rainflowcounting/rainflowcount.py
@@ -64,7 +64,8 @@ def rainflow_windap(signal, levels=255., thresshold=(255 / 50)):
check_signal(signal)
#type <double> is required by <find_extreme> and <rainflow>
signal = signal.astype(np.double)
-
+ if np.all(np.isnan(signal)):
+ return None
offset = np.nanmin(signal)
signal -= offset
if np.nanmax(signal) > 0:
diff --git a/wetb/gtsdf/__init__.py b/wetb/gtsdf/__init__.py
index 4d48fa17a46720365d1de4717a393db9771bf8a9..9023647877f345cb26b6e33194d5964478e61973 100644
--- a/wetb/gtsdf/__init__.py
+++ b/wetb/gtsdf/__init__.py
@@ -64,8 +64,8 @@ class Dataset(object):
try:
return self(name)
except:
- for i, n in enumerate(self.info['attribute_names']):
- print (i,n)
+# for i, n in enumerate(self.info['attribute_names']):
+# print (i,n)
raise e
def __contains__(self, name):
diff --git a/wetb/hawc2/bladeinfo.py b/wetb/hawc2/bladeinfo.py
index 12c1d4aa556aa99c6facaabe0a2609d1711a98ce..e89bd9ab10c5d3619f38b442870bf0e115250cd8 100644
--- a/wetb/hawc2/bladeinfo.py
+++ b/wetb/hawc2/bladeinfo.py
@@ -14,53 +14,54 @@ from wetb.hawc2.st_file import StFile
-class BladeInfo(object):
- def twist(self, radius=None):
- """Aerodynamic twist [deg]. Negative when leading edge is twisted towards wind(opposite to normal definition)\n
+# class BladeInfo(object):
+# def twist(self, radius=None):
+# """Aerodynamic twist [deg]. Negative when leading edge is twisted towards wind(opposite to normal definition)\n
+#
+# Parameters
+# ---------
+# radius : float or array_like, optional
+# radius [m] of interest
+# If None (default) the twist of all points are returned
+#
+# Returns
+# -------
+# twist [deg] : float
+# """
+# raise NotImplementedError()
+#
+# def chord(self, radius=None):
+# """Aerodynamic chord
+#
+# Parameters
+# ---------
+# radius : float or array_like, optional
+# radius [m] of interest
+# If None (default) the twist of all points are returned
+#
+# Returns
+# -------
+# chord [m] : float
+# """
+# raise NotImplementedError()
+#
+#
+#
+# def c2nd(self, radius_nd):
+# """Return center line position
+#
+# Parameters
+# ---------
+# radius_nd : float or array_like, optional
+# non dimensional radius
+#
+# Returns
+# -------
+# x,y,z : float
+# """
- Parameters
- ---------
- radius : float or array_like, optional
- radius [m] of interest
- If None (default) the twist of all points are returned
-
- Returns
- -------
- twist [deg] : float
- """
- raise NotImplementedError()
-
- def chord(self, radius=None):
- """Aerodynamic chord
-
- Parameters
- ---------
- radius : float or array_like, optional
- radius [m] of interest
- If None (default) the twist of all points are returned
-
- Returns
- -------
- chord [m] : float
- """
- raise NotImplementedError()
-
-
- def c2nd(self, radius_nd):
- """Return center line position
-
- Parameters
- ---------
- radius_nd : float or array_like, optional
- non dimensional radius
-
- Returns
- -------
- x,y,z : float
- """
-
-class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
+class H2aeroBladeInfo(PCFile, AtTimeFile):
"""Provide HAWC2 info about a blade
From AE file:
@@ -85,45 +86,64 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
"""
- def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
+ def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, blade_name=None):
if isinstance(htcfile, str):
assert htcfile.lower().endswith('.htc')
htcfile = HTCFile(htcfile)
-
+ self.htcfile = htcfile
blade_name = blade_name or htcfile.aero.link[2]
- s = htcfile.new_htc_structure
+ #s = htcfile.new_htc_structure
at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
- mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
- mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
- st_filename = st_filename or os.path.join(htcfile.modelpath, mainbody_blade.timoschenko_input.filename[0])
+ #mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
+ #self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
PCFile.__init__(self, pc_filename, ae_filename)
blade_radius = self.ae_sets[1][-1,0]
- if os.path.isfile(st_filename):
- StFile.__init__(self, st_filename)
if os.path.isfile(at_time_filename):
AtTimeFile.__init__(self, at_time_filename, blade_radius)
-
-
- self.c2def = np.array([v.values[1:5] for v in mainbody_blade.c2_def if v.name_ == "sec"])
+ self.curved_length = self.radius_curved_ac()[-1]
+ else:
+ raise NotImplementedError
+ #z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
+ #self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
- #self.c2nd = lambda x : interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')(np.max([np.min([x, np.ones_like(x)], 0), np.zeros_like(x) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
- x, y, z, twist = self.hawc2_splines()
- def interpolate(r):
- r = (np.max([np.min([r, np.ones_like(r)], 0), np.zeros_like(r) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
- return np.array([np.interp(r, np.array(z) / z[-1], xyz) for xyz in [x,y,z, twist]]).T
-
- self.c2nd = interpolate #lambda r : interp1d(np.array(z) / z[-1], np.array([x, y, z]).T, axis=0, kind=1)(np.max([np.min([r, np.ones_like(r)], 0), np.zeros_like(r) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
+
+ self.hawc2_splines_data = self.hawc2_splines()
+
+
+ @property
+ def c2def(self):
+ if not hasattr(self, "_c2def"):
+ self._c2def = np.array([v.values[1:5] for v in self.htcfile.blade_c2_def if v.name_ == "sec"])
+ return self._c2def
+
+
+ def c2nd(self, l_nd, curved_length=True):
+ curve_l_nd, x, y, z, twist = self.hawc2_splines_data
+ if curved_length:
+ l_nd = (np.max([np.min([l_nd, np.ones_like(l_nd)], 0), np.zeros_like(l_nd) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
+ return np.array([np.interp(l_nd, curve_l_nd, xyz) for xyz in [x,y,z, twist]]).T
+ else:
+ l_nd = (np.max([np.min([l_nd, np.ones_like(l_nd)], 0), np.zeros_like(l_nd)], 0))
+ return np.array([np.interp(l_nd, z/z[-1], xyz) for xyz in [x,y,z, twist]]).T
+
+ def c2(self, l, curved_length=True):
+ if curved_length:
+ L = self.curved_length
+ else:
+ L = self.c2def[-1,2]
+ return self.c2nd(l/L, curved_length)
+
def hawc2_splines(self):
- curve_z = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
- curve_z_nd = curve_z / curve_z[-1]
+ curve_l = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
+ curve_l_nd = curve_l / curve_l[-1]
def akima(x, y):
n = len(x)
@@ -166,10 +186,11 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
z = np.linspace(0, s[i + 1] - s[i ], 10, endpoint=i >= co.shape[0] - 2)
x.extend(s[i] + z)
y.extend(p(z))
- return y
+ return x,y
- x, y, z, twist = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(4)]
- return x, y, z, twist
+ x, y, z, twist = [coef2spline(curve_l_nd, akima(curve_l_nd, self.c2def[:, i]))[1] for i in range(4)]
+ curve_l_nd = coef2spline(curve_l_nd, akima(curve_l_nd, self.c2def[:, 0]))[0]
+ return curve_l_nd, x, y, z, twist
def xyztwist(self, l=None, curved_length=False):
"""Return splined x,y,z and twist
@@ -188,44 +209,120 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
x,y,z,twist
"""
if l is None:
- return self.c2def[:, 3]
+ return self.c2def[:, :3]
else:
- r_nd = np.linspace(0,1,100)
+ r_nd = np.linspace(0,1,1000)
if curved_length:
- curved_length = np.cumsum(np.sqrt((np.diff(self.c2nd(np.linspace(0,1,100)),1,0)[:,:3]**2).sum(1)))
- assert np.all(l>=curved_length[0]) and np.all(l<=curved_length[-1])
- return self.c2nd(r_nd[np.argmin(np.abs(curved_length-l))+1])
+ #curved_length = np.cumsum(np.sqrt((np.diff(self.c2nd(np.linspace(0,1,100)),1,0)[:,:3]**2).sum(1)))
+ return self.c2nd(l/self.radius_curved_ac()[-1])
+# z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
+# assert np.all(l>=curved_length[0]) and np.all(l<=curved_length[-1])
+# return self.c2nd(r_nd[np.argmin(np.abs(curved_length-l))+1])
else:
assert np.all(l>=self.c2def[0,2]) and np.all(l<=self.c2def[-1,2])
return self.c2nd(l/self.c2def[-1, 2])
-
-
-class H2aeroBladeInfo(H2BladeInfo):
-
- def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
- """
- at_time_filename: file name of at time file containing twist and chord data
- """
- PCFile.__init__(self, pc_filename, ae_filename)
- blade_radius = self.ae_sets[1][-1,0]
- AtTimeFile.__init__(self, at_time_filename, blade_radius)
-
- assert('twist' in self.attribute_names)
- htcfile = HTCFile(htc_filename)
-
-
- self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
- #self._c2nd = interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')
-
- ac_radii = self.ac_radius()
- c2_axis_length = np.r_[0, np.cumsum(np.sqrt((np.diff(self.c2def[:, :3], 1, 0) ** 2).sum(1)))]
- self._c2nd = interp1d(c2_axis_length / c2_axis_length[-1], self.c2def[:, :3], axis=0, kind='cubic')
- #self._c2nd = interp1d(self.c2def[:,2]/self.c2def[-1,2], self.c2def[:, :3], axis=0, kind='cubic')
-
- def c2nd(self, r_nd):
- r_nd_min = np.zeros_like(r_nd) + self.c2def[0, 2] / self.c2def[-1, 2]
- r_nd_max = np.ones_like(r_nd)
- r_nd = np.max([np.min([r_nd, r_nd_max], 0), r_nd_min], 0)
- return self._c2nd(r_nd)
+
+
+class H2BladeInfo(H2aeroBladeInfo):
+ """Provide HAWC2 info about a blade
+
+ From AE file:
+ - chord(radius=None, set_nr=1):
+ - thickness(radius=None, set_nr=1)
+ - radius_ae(radius=None, set_nr=1)
+
+ From PC file
+ - CL(radius, alpha, ae_set_nr=1)
+ - CD(radius, alpha, ae_set_nr=1)
+ - CM(radius, alpha, ae_set_nr=1)
+ From at_time_filename
+ - attribute_names
+ - xxx(radius=None, curved_length=None) # xxx for each attribute name
+ - radius_curved_ac(radius=None) # Curved length of nearest/all aerodynamic calculation points
+
+ From ST file
+ - radius_st(radius=None, mset=1, set=1)
+ - xxx(radius=None, mset=1, set=1) # xxx for each of r, m, x_cg,y_cg, ri_x, ri_y, xs, ys, E, G, Ix, Iy, K, kx, ky, A, pitch, xe, ye
+ with template
+
+ """
+ def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
+ if isinstance(htcfile, str):
+ htcfile = HTCFile(htcfile)
+ s = htcfile.new_htc_structure
+# at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
+# pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
+# ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
+#
+ mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
+ if blade_name is None:
+ blade_name = htcfile.aero.link[2]
+ self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
+ H2aeroBladeInfo.__init__(self, htcfile, ae_filename=ae_filename, pc_filename=pc_filename, at_time_filename=at_time_filename, blade_name=blade_name)
+
+# def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
+#
+# if isinstance(htcfile, str):
+# assert htcfile.lower().endswith('.htc')
+# htcfile = HTCFile(htcfile)
+#
+# blade_name = blade_name or htcfile.aero.link[2]
+ st_filename = st_filename or os.path.join(htcfile.modelpath, self.mainbody_blade.timoschenko_input.filename[0])
+#
+# if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
+# PCFile.__init__(self, pc_filename, ae_filename)
+# blade_radius = self.ae_sets[1][-1,0]
+#
+ if os.path.isfile(st_filename):
+ StFile.__init__(self, st_filename)
+# if os.path.isfile(at_time_filename):
+# AtTimeFile.__init__(self, at_time_filename, blade_radius)
+# self.curved_length = self.radius_curved_ac()[-1]
+# else:
+# raise NotImplementedError
+# #z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
+# #self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
+#
+# self.c2def = np.array([v.values[1:5] for v in mainbody_blade.c2_def if v.name_ == "sec"])
+#
+# self.hawc2_splines_data = self.hawc2_splines()
+ @property
+ def c2def(self):
+ if not hasattr(self, "_c2def"):
+ self._c2def = np.array([v.values[1:5] for v in self.mainbody_blade.c2_def if v.name_ == "sec"])
+ return self._c2def
+
+#
+# class H2aeroBladeInfo(H2BladeInfo):
+#
+# def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
+# """
+# at_time_filename: file name of at time file containing twist and chord data
+# """
+# PCFile.__init__(self, pc_filename, ae_filename)
+# blade_radius = self.ae_sets[1][-1,0]
+# AtTimeFile.__init__(self, at_time_filename, blade_radius)
+#
+# assert('twist' in self.attribute_names)
+# htcfile = HTCFile(htc_filename)
+#
+#
+# self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
+# #self._c2nd = interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')
+#
+# ac_radii = self.radius_curved_ac()
+# self.hawc2_splines_data = self.hawc2_splines()
+# self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
+#
+# # c2_axis_length = np.r_[0, np.cumsum(np.sqrt((np.diff(self.c2def[:, :3], 1, 0) ** 2).sum(1)))]
+# # self._c2nd = interp1d(c2_axis_length / c2_axis_length[-1], self.c2def[:, :3], axis=0, kind='cubic')
+# #self._c2nd = interp1d(self.c2def[:,2]/self.c2def[-1,2], self.c2def[:, :3], axis=0, kind='cubic')
+# #
+# # def c2nd(self, r_nd):
+# # r_nd_min = np.zeros_like(r_nd) + self.c2def[0, 2] / self.c2def[-1, 2]
+# # r_nd_max = np.ones_like(r_nd)
+# # r_nd = np.max([np.min([r_nd, r_nd_max], 0), r_nd_min], 0)
+# # return self._c2nd(r_nd)
+#
diff --git a/wetb/hawc2/htc_contents.py b/wetb/hawc2/htc_contents.py
index bd40ea4e2b004e22c953d6b472d774a7e53ee865..2df59c6014deb1c32c413f7b1863f21297463a71 100644
--- a/wetb/hawc2/htc_contents.py
+++ b/wetb/hawc2/htc_contents.py
@@ -14,6 +14,7 @@ from builtins import zip
from builtins import int
from builtins import str
from future import standard_library
+import os
standard_library.install_aliases()
from collections import OrderedDict
import collections
@@ -119,8 +120,8 @@ class HTCContents(object):
ending = "__%d" % (1 + float("0%s" % ending.replace("__", "")))
self[contents.name_ + ending] = contents
- def add_section(self, name, allow_duplicate_section=False):
- if name in self and allow_duplicate_section is False:
+ def add_section(self, name, allow_duplicate=False):
+ if name in self and allow_duplicate is False:
return self[name]
section = HTCSection(name)
self._add_contents(section)
@@ -363,8 +364,9 @@ class HTCDefaults(object):
def add_mann_turbulence(self, L=29.4, ae23=1, Gamma=3.9, seed=1001, high_frq_compensation=True,
filenames=None,
- no_grid_points=(4096, 32, 32), box_dimension=(6000, 100, 100),
- std_scaling=(1, .8, .5)):
+ no_grid_points=(16384, 32, 32), box_dimension=(6000, 100, 100),
+ dont_scale=False,
+ std_scaling=None):
wind = self.add_section('wind')
wind.turb_format = 1
mann = wind.add_section('mann')
@@ -383,15 +385,38 @@ class HTCDefaults(object):
setattr(mann, 'filename_%s' % c, filename)
for c, n, dim in zip(['u', 'v', 'w'], no_grid_points, box_dimension):
setattr(mann, 'box_dim_%s' % c, "%d %.4f" % (n, dim / (n - 1)))
- if std_scaling is None:
+ if dont_scale:
mann.dont_scale = 1
else:
try:
del mann.dont_scale
except KeyError:
pass
+ if std_scaling is not None:
mann.std_scaling = "%f %f %f" % std_scaling
+ else:
+ try:
+ del mann.std_scaling
+ except KeyError:
+ pass
+
+
+
+ def add_turb_export(self, filename="export_%s.turb", samplefrq = None):
+ exp = self.wind.add_section('turb_export', allow_duplicate=True)
+ for uvw in 'uvw':
+ exp.add_line('filename_%s'%uvw, [filename%uvw])
+ sf = samplefrq or max(1,int( self.wind.mann.box_dim_u[1]/(self.wind.wsp[0] * self.deltat())))
+ exp.samplefrq = sf
+ if "time" in self.output:
+ exp.time_start = self.output.time[0]
+ else:
+ exp.time_start = 0
+ exp.nsteps = (self.simulation.time_stop[0]-exp.time_start[0]) / self.deltat()
+ for vw in 'vw':
+ exp.add_line('box_dim_%s'%vw, self.wind.mann['box_dim_%s'%vw].values)
+
def import_dtu_we_controller_input(self, filename):
dtu_we_controller = [dll for dll in self.dll if dll.name[0] == 'dtu_we_controller'][0]
diff --git a/wetb/hawc2/htc_file.py b/wetb/hawc2/htc_file.py
index 2ae9e14518d7c0f06903f92904e35d594bdcdcb3..087fb242f3e9aefa4eaa9a9a24573daebe11891d 100644
--- a/wetb/hawc2/htc_file.py
+++ b/wetb/hawc2/htc_file.py
@@ -197,7 +197,7 @@ class HTCFile(HTCContents, HTCDefaults):
with open(filename, 'w', encoding='cp1252') as fid:
fid.write(str(self))
- def set_name(self, name, htc_folder="htc", log_folder="log", res_folder="res"):
+ def set_name(self, name, htc_folder="htc", log_folder="log", res_folder="res", animation_folder='animation', visualization_folder="visualization"):
#if os.path.isabs(folder) is False and os.path.relpath(folder).startswith("htc" + os.path.sep):
self.contents #load if not loaded
fmt_folder = lambda folder : "./" + os.path.relpath(folder).replace("\\", "/")
@@ -207,6 +207,10 @@ class HTCFile(HTCContents, HTCDefaults):
self.simulation.logfile = os.path.join(fmt_folder(log_folder), "%s.log" % name).replace("\\", "/")
elif 'test_structure' in self and 'logfile' in self.test_structure: # hawc2aero
self.test_structure.logfile = os.path.join(fmt_folder(log_folder), "%s.log" % name).replace("\\", "/")
+ if 'simulation' in self and 'animation' in self.simulation:
+ self.simulation.animation = os.path.join(fmt_folder(animation_folder), "%s.dat" % name).replace("\\", "/")
+ if 'simulation' in self and 'visualization' in self.simulation:
+ self.simulation.visualization = os.path.join(fmt_folder(visualization_folder), "%s.hdf5" % name).replace("\\", "/")
self.output.filename = os.path.join(fmt_folder(res_folder), "%s" % name).replace("\\", "/")
def set_time(self, start=None, stop=None, step=None):
@@ -354,6 +358,20 @@ class HTCFile(HTCContents, HTCDefaults):
if errorcode or 'Elapsed time' not in log:
raise Exception (str(stdout) + str(stderr))
+
+ def deltat(self):
+ return self.simulation.newmark.deltat[0]
+
+ def get_body(self, name):
+ lst = [b for b in self.new_htc_structure if b.name_=="main_body" and b.name[0]==name]
+ if len(lst)==1:
+ return lst[0]
+ else:
+ if len(lst)==0:
+ raise ValueError("Body '%s' not found"%name)
+ else:
+ raise NotImplementedError()
+
class H2aeroHTCFile(HTCFile):
def __init__(self, filename=None, modelpath=None):
diff --git a/wetb/hawc2/pc_file.py b/wetb/hawc2/pc_file.py
index bd0da0a7403545b893cbea7ac1a2b4dd5f432689..301a5b81ab93840138ae5a6326e3f198ceb0c093 100644
--- a/wetb/hawc2/pc_file.py
+++ b/wetb/hawc2/pc_file.py
@@ -74,6 +74,26 @@ class PCFile(AEFile):
Cx1 = np.interp(alpha, Cx1[:, 0], Cx1[:, column])
th0, th1 = thicknesses[index - 1:index + 1]
return Cx0 + (Cx1 - Cx0) * (thickness - th0) / (th1 - th0)
+
+ def _CxxxH2(self, radius, alpha, column, ae_set_nr=1):
+ thickness = self.thickness(radius, ae_set_nr)
+ pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+ thicknesses, profiles = self.pc_sets[pc_set_nr]
+ index = np.searchsorted(thicknesses, thickness)
+ if index == 0:
+ index = 1
+
+ Cx0, Cx1 = profiles[index - 1:index + 1]
+
+ Cx0 = np.interp(np.arange(360), Cx0[:,0]+180, Cx0[:,column])
+ Cx1 = np.interp(np.arange(360), Cx1[:,0]+180, Cx1[:,column])
+ #Cx0 = np.interp(alpha, Cx0[:, 0], Cx0[:, column])
+ #Cx1 = np.interp(alpha, Cx1[:, 0], Cx1[:, column])
+ th0, th1 = thicknesses[index - 1:index + 1]
+ cx = Cx0 + (Cx1 - Cx0) * (thickness - th0) / (th1 - th0)
+ return np.interp(alpha+180, np.arange(360), cx)
+
+
def CL(self, radius, alpha, ae_set_nr=1):
"""Lift coefficient
@@ -93,6 +113,10 @@ class PCFile(AEFile):
"""
return self._Cxxx(radius, alpha, 1, ae_set_nr)
+
+ def CL_H2(self, radius, alpha, ae_set_nr=1):
+ return self._CxxxH2(radius, alpha, 1, ae_set_nr)
+
def CD(self, radius, alpha, ae_set_nr=1):
"""Drag coefficient
diff --git a/wetb/hawc2/st_file.py b/wetb/hawc2/st_file.py
index 7859144fab6a442ff110c681297e46fd4a32690a..d999671cf3ef2bc7763d03c202078d1fee84b595 100644
--- a/wetb/hawc2/st_file.py
+++ b/wetb/hawc2/st_file.py
@@ -89,7 +89,7 @@ class StFile(object):
def _value(self, radius, column, mset_nr=1, set_nr=1):
st_data = self.main_data_sets[mset_nr][set_nr]
if radius is None:
- radius = self.radius(None, mset_nr, set_nr)
+ radius = self.radius_st(None, mset_nr, set_nr)
return np.interp(radius, st_data[:, 0], st_data[:, column])
def radius_st(self, radius=None, mset=1, set=1):
diff --git a/wetb/hawc2/tests/test_AtTimeFile.py b/wetb/hawc2/tests/test_AtTimeFile.py
index d6ae0cf24f018976aa39e2726f88157303e2af64..909151b2bb01bf27342b844c01ff51a7104206ac 100644
--- a/wetb/hawc2/tests/test_AtTimeFile.py
+++ b/wetb/hawc2/tests/test_AtTimeFile.py
@@ -60,6 +60,11 @@ class TestAtTimeFile(unittest.TestCase):
self.assertEqual(atfile.radius_curved_ac(10), 10.2505)
self.assertEqual(atfile.radius_curved_ac(10.5), 10.2505)
+ def test_rotor_name_null(self):
+ atfile = AtTimeFile(self.testfilepath + "at_time/test_rotor_name_null.dat", bladetip_radius=20.501) # load file
+ print (atfile.radius_s())
+
+
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()
diff --git a/wetb/hawc2/tests/test_files/at_time/test_rotor_name_null.dat b/wetb/hawc2/tests/test_files/at_time/test_rotor_name_null.dat
new file mode 100644
index 0000000000000000000000000000000000000000..778834e44e09c58137cbe825905a028add1013c5
Binary files /dev/null and b/wetb/hawc2/tests/test_files/at_time/test_rotor_name_null.dat differ
diff --git a/wetb/hawc2/tests/test_htc_file.py b/wetb/hawc2/tests/test_htc_file.py
index 8d2515f26998fb4fe134a7935f8d278c76a73d82..0326807eb6ea8865fa59285a3810ecb471075353 100644
--- a/wetb/hawc2/tests/test_htc_file.py
+++ b/wetb/hawc2/tests/test_htc_file.py
@@ -148,17 +148,35 @@ class TestHtcFile(unittest.TestCase):
htcfile.add_mann_turbulence(30.1, 1.1, 3.3, 102, False)
s = """begin mann;
create_turb_parameters\t30.1 1.1 3.3 102 0;\tL, alfaeps, gamma, seed, highfrq compensation
- filename_u\t./turb/mann_l30.1_ae1.10_g3.3_h0_4096x32x32_1.465x3.12x3.12_s0102u.turb;
- filename_v\t./turb/mann_l30.1_ae1.10_g3.3_h0_4096x32x32_1.465x3.12x3.12_s0102v.turb;
- filename_w\t./turb/mann_l30.1_ae1.10_g3.3_h0_4096x32x32_1.465x3.12x3.12_s0102w.turb;
- box_dim_u\t4096 1.4652;
+ filename_u\t./turb/mann_l30.1_ae1.10_g3.3_h0_16384x32x32_0.366x3.12x3.12_s0102u.turb;
+ filename_v\t./turb/mann_l30.1_ae1.10_g3.3_h0_16384x32x32_0.366x3.12x3.12_s0102v.turb;
+ filename_w\t./turb/mann_l30.1_ae1.10_g3.3_h0_16384x32x32_0.366x3.12x3.12_s0102w.turb;
+ box_dim_u\t16384 0.3662;
box_dim_v\t32 3.2258;
- box_dim_w\t32 3.2258;
- std_scaling\t1 0.8 0.5;"""
+ box_dim_w\t32 3.2258;"""
for a, b in zip(s.split("\n"), str(htcfile.wind.mann).split("\n")):
self.assertEqual(a.strip(), b.strip())
self.assertEqual(htcfile.wind.turb_format[0], 1)
self.assertEqual(htcfile.wind.turb_format.comments, "")
+
+ def test_add_turb_export(self):
+ htc = HTCFile()
+ htc.add_mann_turbulence(30.1, 1.1, 3.3, 102, False)
+ htc.set_time(100,700,0.01)
+ htc.add_turb_export()
+ s = """begin turb_export;
+ filename_u\texport_u.turb;
+ filename_v\texport_v.turb;
+ filename_w\texport_w.turb;
+ samplefrq\t3;
+ time_start\t100;
+ nsteps\t60000.0;
+ box_dim_v\t32 3.2258;
+ box_dim_w\t32 3.2258;
+end turb_export;"""
+ for a, b in zip(s.split("\n"), str(htc.wind.turb_export).split("\n")):
+ self.assertEqual(a.strip(), b.strip())
+
def test_sensors(self):
diff --git a/wetb/signal/filters/_differentiation.py b/wetb/signal/filters/_differentiation.py
index 94406ca98696259225cb45d60256edaadb872a99..241fe62bab9945d6647ef08742ecf364058fa0e2 100644
--- a/wetb/signal/filters/_differentiation.py
+++ b/wetb/signal/filters/_differentiation.py
@@ -40,8 +40,10 @@ def differentiation(x, type='center', sample_frq=None, cutoff_frq=None):
x = np.array(x)
if type=="left":
dy = np.r_[np.nan, x[1:]-x[:-1]]
- if type=="center":
+ elif type=="center":
dy = np.r_[x[1]-x[0], (x[2:]-x[:-2])/2, x[-1]-x[-2]]
- if type=="right":
+ elif type=="right":
dy = np.r_[x[1:]-x[:-1], np.nan]
+ else:
+ raise NotImplementedError()
return dy
\ No newline at end of file
diff --git a/wetb/utils/process_exec.py b/wetb/utils/process_exec.py
index eeaf03497d7ad6d738ba1d96e49990c14d9a8c9b..589024b8b5880acd2997dce5f061e7047c44c103 100644
--- a/wetb/utils/process_exec.py
+++ b/wetb/utils/process_exec.py
@@ -32,7 +32,7 @@ def pexec(args, cwd=None):
args[i] = str(args[i]).replace('/', os.path.sep).replace('\\', os.path.sep).replace('"', '')
cmd = "%s" % '{} /c "{}"'.format (os.environ.get("COMSPEC", "cmd.exe"), subprocess.list2cmdline(args))
- if os.path.isfile(cwd):
+ if cwd and os.path.isfile(cwd):
cwd = os.path.dirname(cwd)
proc = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, cwd=cwd)
stdout, stderr = proc.communicate()