bladeinfo.py

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)])

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))

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))
 

wetb/hawc2/bladeinfo.py

+'''
+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)
+    
+

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__":

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

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']

wetb/hawc2/tests/test_files/at.dat

+ # 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