diff --git a/wetb/gtsdf/__init__.py b/wetb/gtsdf/__init__.py index 0c84199f52b04b0e90da6e7ac32a120f6e28a62b..4d48fa17a46720365d1de4717a393db9771bf8a9 100644 --- a/wetb/gtsdf/__init__.py +++ b/wetb/gtsdf/__init__.py @@ -67,7 +67,9 @@ class Dataset(object): for i, n in enumerate(self.info['attribute_names']): print (i,n) raise e - + + def __contains__(self, name): + return name in self.info['attribute_names'] __all__ = sorted([m for m in set(dir()) - set(d)]) diff --git a/wetb/hawc2/ae_file.py b/wetb/hawc2/ae_file.py index cee13d9f107fea4ad249ad514e480e70d61d2233..d3628fee19eebc8bc152588b7c572a703e01c179 100644 --- a/wetb/hawc2/ae_file.py +++ b/wetb/hawc2/ae_file.py @@ -45,14 +45,24 @@ class AEFile(object): def _value(self, radius, column, set_nr=1): ae_data = self.ae_sets[set_nr] - return np.interp(radius, ae_data[:, 0], ae_data[:, column]) + if radius: + return np.interp(radius, ae_data[:, 0], ae_data[:, column]) + else: + return ae_data[:,column] - def chord(self, radius, set_nr=1): + def chord(self, radius=None, set_nr=1): return self._value(radius, 1, set_nr) - def thickness(self, radius, set_nr=1): + def thickness(self, radius=None, set_nr=1): return self._value(radius, 2, set_nr) - + + def radius_ae(self, radius=None, set_nr=1): + radii = self.ae_sets[set_nr][:,0] + if radius: + return radii[np.argmin(np.abs(radii-radius))] + else: + return radii + def pc_set_nr(self, radius, set_nr=1): ae_data = self.ae_sets[set_nr] index = np.searchsorted(ae_data[:, 0], radius) @@ -68,6 +78,7 @@ class AEFile(object): if __name__ == "__main__": ae = AEFile(r"tests/test_files/NREL_5MW_ae.txt") - print (ae.thickness(36)) + print (ae.radius_ae(36)) + print (ae.thickness()) print (ae.chord(36)) print (ae.pc_set_nr(36)) diff --git a/wetb/hawc2/at_time_file.py b/wetb/hawc2/at_time_file.py index 10ed8583dd11b1bb280db37a1aa216d43bfcbde1..c84c75e6fdd38bb6646aac366556b261e40744fa 100644 --- a/wetb/hawc2/at_time_file.py +++ b/wetb/hawc2/at_time_file.py @@ -16,22 +16,37 @@ import numpy as np class AtTimeFile(object): """Loads an at time file generated by HAWC2 - Note that the radius in this context is the curved radius - - >>> atfile = AtTimeFile("at_time.dat") # load file - >>> atfile.attribute_names # Attribute names + Functions are generated for each sensor in the file from the template: + xxx(radius=None, curved_length=None) + + E.g. + twist(radius=36) # Twist at radius 36m (Straight line radius, blade_radius must be specified in constructor) + twist(curved_length=36) # Twist 36m from the root along the (curved) center line + twist() # Twist at all aerodynamic calculation points + + + >>> at = AtTimeFile("at.dat", 86.3655) # load file + >>> at.attribute_names # Attribute names ['radius_s', 'twist', 'chord'] - >>> atfile[:3,1]) # first 3 twist rows - [ 0. -0.775186 -2.91652 ] - >>> atfile.twist()[:3]) # first 3 twist rows - [ 0. -0.775186 -2.91652 ] - >>> atfile.twist(curved_length=10) # Twist at curved_length = 10 (interpolated) - -5.34743208242399 - >>> atfile.twist(radius=10) # Twist at curved_length = 10 (interpolated) - -5.34743208242399 + >>> at[:3,1]) # first 3 twist rows + [-14.5 -14.5002 -14.5007] + >>> at.twist()[:3]) # Twist first 3 twist rows + [-14.5 -14.5002 -14.5007] + >>> at.twist(curved_length=36) # Twist at curved_length = 10 (interpolated) + -5.172195702789108 + >>> at.twist(radius=36) # Twist at curved_length = 10 (interpolated) + -5.162084567646019 """ - def __init__(self, filename, blade_radius=None): - self.blade_radius = blade_radius + def __init__(self, filename, bladetip_radius=None): + """ + Parameters + ---------- + filename : string + Filename + bladetip_radius : int, float or None, optional + Radius of blade tip. Used to convert from curved length to radius + """ + self.blade_radius = bladetip_radius with open(filename, encoding='utf-8') as fid: lines = fid.readlines() self.attribute_names = lines[2].lower().replace("#", "").split() @@ -43,7 +58,7 @@ class AtTimeFile(object): if radius is None and curved_length is None: return self.data[:, column] elif radius is not None: - assert self.blade_radius is not None, "blade_radius must be specified in __init__ when requesting value of radius (alternatively you can request for curved_length)" + assert self.blade_radius is not None, "bladetip_radius must be specified in __init__ when requesting value of radius (alternatively you can request for curved_length)" return np.interp(radius/self.blade_radius, self.data[:, 0]/self.data[-1, 0], self.data[:, column]) else: return np.interp(curved_length, self.data[:, 0], self.data[:, column]) @@ -51,7 +66,7 @@ class AtTimeFile(object): for column, att_name in enumerate(self.attribute_names): setattr(self, att_name, func_factory(column)) - def ac_radius(self, radius=None): + def radius_ac(self, radius=None): """Radius (curved distance) of aerodynamic calculation point(s) Parameters @@ -66,7 +81,7 @@ class AtTimeFile(object): Radius of calculation points or radius of calculation point nearest radius """ if radius is None: - return self.radius_s(radius) + return self.radius_s() else: return self.radius_s()[np.argmin(np.abs(self.radius_s() - radius))] @@ -78,8 +93,18 @@ class AtTimeFile(object): if __name__ == "__main__": - at = AtTimeFile(r"tests/test_files/at_time.dat") - print (at.attribute_names) - print (at.twist(36)) - print (at.chord(36)) + at = AtTimeFile(r"tests/test_files/at.dat", 86.3655) # load file + print (at.attribute_names) # Attribute names + print (at[:3,1]) # first 3 twist rows + print (at.twist()[:3]) # first 3 twist rows + print (at.twist(curved_length=36)) # Twist at curved_length = 10 (interpolated) + print (at.twist(radius=36)) # Twist at curved_length = 10 (interpolated) + +# at = AtTimeFile(r"tests/test_files/at.dat", 86.3655) +# print (at.attribute_names) +# print (at.radius_s()) +# print (at.twist(curved_length=36)) +# print (at.chord(curved_length=36)) +# print (at.twist(curved_length=36)) +# print (at.twist(radius=36)) diff --git a/wetb/hawc2/bladeinfo.py b/wetb/hawc2/bladeinfo.py new file mode 100644 index 0000000000000000000000000000000000000000..745685f31f1550761127592954c20890b2f658dc --- /dev/null +++ b/wetb/hawc2/bladeinfo.py @@ -0,0 +1,210 @@ +''' +Created on 3. maj 2017 + +@author: mmpe +''' +import numpy as np +from wetb.hawc2.at_time_file import AtTimeFile +from wetb.hawc2.pc_file import PCFile +from wetb.hawc2.htc_file import HTCFile +from scipy.interpolate.interpolate import interp1d +import os +from wetb.utils.geometry import rad +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 + + 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): + """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_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): + assert htcfile.lower().endswith('.htc') + htcfile = HTCFile(htcfile) + + blade_name = blade_name or htcfile.aero.link[2] + s = htcfile.new_htc_structure + at_time_filename = at_time_filename or 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]) + + 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.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)) + + 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] + + def akima(x, y): + n = len(x) + var = np.zeros((n + 3)) + z = np.zeros((n)) + co = np.zeros((n, 4)) + for i in range(n - 1): + var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) + var[n + 1] = 2 * var[n] - var[n - 1] + var[n + 2] = 2 * var[n + 1] - var[n] + var[1] = 2 * var[2] - var[3] + var[0] = 2 * var[1] - var[2] + + for i in range(n): + wi1 = abs(var[i + 3] - var[i + 2]) + wi = abs(var[i + 1] - var[i]) + if (wi1 + wi) == 0: + z[i] = (var[i + 2] + var[i + 1]) / 2 + else: + z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi) + + for i in range(n - 1): + dx = x[i + 1] - x[i] + a = (z[i + 1] - z[i]) * dx + b = y[i + 1] - y[i] - z[i] * dx + co[i, 0] = y[i] + co[i, 1] = z[i] + co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx + co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2 + co[n - 1, 0] = y[n - 1] + co[n - 1, 1] = z[n - 1] + co[n - 1, 2] = 0 + co[n - 1, 3] = 0 + return co + + def coef2spline(s, co): + x, y = [], [] + for i, c in enumerate(co.tolist()[:-1]): + p = np.poly1d(c[::-1]) + 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 + + 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 + + def c2def_twist(self, radius=None): + if radius is None: + return self.c2def[:, 3] + else: + return np.interp(radius, self.c2def[:, 2], self.c2def[:, 3]) + + + +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) + + diff --git a/wetb/hawc2/pc_file.py b/wetb/hawc2/pc_file.py index f42cc1076518a42f14b38cb307b65c6847a28a39..bd0da0a7403545b893cbea7ac1a2b4dd5f432689 100644 --- a/wetb/hawc2/pc_file.py +++ b/wetb/hawc2/pc_file.py @@ -43,7 +43,7 @@ class PCFile(AEFile): with open (filename) as fid: lines = fid.readlines() nsets = int(lines[0].split()[0]) - self.sets = {} + self.pc_sets = {} lptr = 1 for nset in range(1, nsets + 1): nprofiles = int(lines[lptr].split()[0]) @@ -59,12 +59,12 @@ class PCFile(AEFile): thicknesses.append(thickness) profiles.append(data) lptr += n_rows - self.sets[nset] = (np.array(thicknesses), profiles) + self.pc_sets[nset] = (np.array(thicknesses), profiles) def _Cxxx(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.sets[pc_set_nr] + thicknesses, profiles = self.pc_sets[pc_set_nr] index = np.searchsorted(thicknesses, thickness) if index == 0: index = 1 @@ -112,7 +112,6 @@ class PCFile(AEFile): return self._Cxxx(radius, alpha, 2, ae_set_nr) def CM(self, radius, alpha, ae_set_nr=1): - return self._Cxxx(radius, alpha, 3, ae_set_nr) if __name__ == "__main__": diff --git a/wetb/hawc2/st_file.py b/wetb/hawc2/st_file.py index f03efafb96d5b93b58c330e427abaed48f09e2cb..09ae72407aa278b03c2a5164187bfb82c9a406e4 100644 --- a/wetb/hawc2/st_file.py +++ b/wetb/hawc2/st_file.py @@ -91,7 +91,7 @@ class StFile(object): radius = self.radius(None, mset_nr, set_nr) return np.interp(radius, st_data[:, 0], st_data[:, column]) - def radius(self, radius=None, mset=1, set=1): + def radius_st(self, radius=None, mset=1, set=1): r = self.main_data_sets[mset][set][:, 0] if radius is None: return r diff --git a/wetb/hawc2/tests/test_AtTimeFile.py b/wetb/hawc2/tests/test_AtTimeFile.py index 6f17562bab8741498c04106893e0541d9f50a17f..09573fd33ba82feb1f5dc5d52cc5ba1ec1ea0b21 100644 --- a/wetb/hawc2/tests/test_AtTimeFile.py +++ b/wetb/hawc2/tests/test_AtTimeFile.py @@ -15,7 +15,7 @@ import numpy as np import os - +tfp = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path class TestAtTimeFile(unittest.TestCase): def setUp(self): @@ -24,7 +24,7 @@ class TestAtTimeFile(unittest.TestCase): def test_doc_examples(self): - atfile = AtTimeFile(self.testfilepath + "at_time.dat", blade_radius=20.501) # load file + atfile = AtTimeFile(self.testfilepath + "at_time.dat", bladetip_radius=20.501) # load file self.assertEqual(atfile.attribute_names, ['radius_s', 'twist', 'chord']) np.testing.assert_array_equal(atfile[:3, 1], [ 0., -0.775186, -2.91652 ]) np.testing.assert_array_equal(atfile.twist()[:3], [ 0. , -0.775186 , -2.91652 ]) @@ -48,7 +48,7 @@ class TestAtTimeFile(unittest.TestCase): self.assertEqual(atfile.chord(curved_length=9), 1.3888996578373045) def test_at_time_file_at_radius(self): - atfile = AtTimeFile(self.testfilepath + "at_time.dat", blade_radius=20.501/2) + atfile = AtTimeFile(self.testfilepath + "at_time.dat", bladetip_radius=20.501/2) self.assertEqual(atfile.radius_s(radius=9/2), 9) self.assertEqual(atfile.twist(radius=9/2), -6.635983309665461) self.assertEqual(atfile.chord(radius=9/2), 1.3888996578373045) @@ -56,9 +56,9 @@ class TestAtTimeFile(unittest.TestCase): def test_at_time_file_radius(self): atfile = AtTimeFile(self.testfilepath + "at_time.dat") - self.assertEqual(atfile.ac_radius()[12], 10.2505) - self.assertEqual(atfile.ac_radius(10), 10.2505) - self.assertEqual(atfile.ac_radius(10.5), 10.2505) + self.assertEqual(atfile.radius_ac()[12], 10.2505) + self.assertEqual(atfile.radius_ac(10), 10.2505) + self.assertEqual(atfile.radius_ac(10.5), 10.2505) if __name__ == "__main__": #import sys;sys.argv = ['', 'Test.testName'] diff --git a/wetb/hawc2/tests/test_files/at.dat b/wetb/hawc2/tests/test_files/at.dat new file mode 100644 index 0000000000000000000000000000000000000000..f93af10f4f3c18596f1dd613503a5bc82847fffe --- /dev/null +++ b/wetb/hawc2/tests/test_files/at.dat @@ -0,0 +1,53 @@ + # Version ID : HAWC2MB 12.4 + # Radial output at time: 0.120000000000000 + # Radius_s # twist # Chord + 0.00000E+00 -1.45000E+01 5.38000E+00 + 8.88594E-02 -1.45002E+01 5.38000E+00 + 3.55072E-01 -1.45007E+01 5.38000E+00 + 7.97544E-01 -1.45012E+01 5.38000E+00 + 1.41446E+00 -1.45013E+01 5.38000E+00 + 2.20328E+00 -1.45008E+01 5.38000E+00 + 3.16076E+00 -1.44998E+01 5.38000E+00 + 4.28297E+00 -1.44959E+01 5.38000E+00 + 5.56530E+00 -1.44880E+01 5.38150E+00 + 7.00248E+00 -1.44609E+01 5.40757E+00 + 8.58860E+00 -1.43497E+01 5.47316E+00 + 1.03171E+01 -1.40434E+01 5.57870E+00 + 1.21810E+01 -1.34021E+01 5.71677E+00 + 1.41726E+01 -1.24242E+01 5.87200E+00 + 1.62836E+01 -1.11781E+01 6.02007E+00 + 1.85054E+01 -9.82623E+00 6.13078E+00 + 2.08289E+01 -8.68972E+00 6.18999E+00 + 2.32446E+01 -7.84255E+00 6.19714E+00 + 2.57424E+01 -7.19893E+00 6.15543E+00 + 2.83122E+01 -6.63223E+00 6.06940E+00 + 3.09433E+01 -6.10664E+00 5.94438E+00 + 3.36250E+01 -5.60732E+00 5.78666E+00 + 3.63463E+01 -5.10875E+00 5.60337E+00 + 3.90959E+01 -4.60713E+00 5.40154E+00 + 4.18626E+01 -4.07749E+00 5.18393E+00 + 4.46350E+01 -3.51118E+00 4.95580E+00 + 4.74017E+01 -2.93766E+00 4.72176E+00 + 5.01513E+01 -2.36236E+00 4.48569E+00 + 5.28726E+01 -1.79469E+00 4.25127E+00 + 5.55543E+01 -1.25288E+00 4.02160E+00 + 5.81854E+01 -7.42980E-01 3.79894E+00 + 6.07552E+01 -2.63044E-01 3.58510E+00 + 6.32530E+01 1.76184E-01 3.38191E+00 + 6.56687E+01 5.69023E-01 3.19012E+00 + 6.79922E+01 9.24087E-01 3.01062E+00 + 7.02140E+01 1.25492E+00 2.84386E+00 + 7.23251E+01 1.55945E+00 2.68993E+00 + 7.43167E+01 1.84168E+00 2.54872E+00 + 7.61805E+01 2.10535E+00 2.42040E+00 + 7.79091E+01 2.34909E+00 2.29766E+00 + 7.94952E+01 2.57177E+00 2.17211E+00 + 8.09324E+01 2.77156E+00 2.04065E+00 + 8.22147E+01 2.94258E+00 1.90297E+00 + 8.33369E+01 3.08259E+00 1.75760E+00 + 8.42945E+01 3.19421E+00 1.59874E+00 + 8.50832E+01 3.28140E+00 1.41325E+00 + 8.56999E+01 3.34684E+00 1.20310E+00 + 8.61425E+01 3.39238E+00 1.00803E+00 + 8.64090E+01 3.41914E+00 8.77215E-01 + 8.64979E+01 3.42796E+00 8.33540E-01