Skip to content
Snippets Groups Projects
Commit 11f0835e authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files

AtTimeFile test for null-rotor

parent 3188ee93
No related branches found
No related tags found
No related merge requests found
......@@ -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]
......
......@@ -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
......
......@@ -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:
......
......@@ -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):
......
......@@ -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)
#
......@@ -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
......
......@@ -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):
......
......@@ -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()
......@@ -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
......@@ -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()
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment