Merge pull request #363 from WISDEM/ccblade

Code robustness and cleanup from ccblade repo

wisdem/airfoilprep/airfoilprep.py

@@ -21,11 +21,10 @@
 
 """
 
-from __future__ import print_function
-from math import pi, sin, cos, radians, degrees, tan, ceil, floor
-import numpy as np
 import copy
 
+import numpy as np
+
 # from scipy.interpolate import RectBivariateSpline
 
 
@@ -136,9 +135,9 @@ def correction3D(self, r_over_R, chord_over_r, tsr, alpha_max_corr=30, alpha_lin
         alpha = np.radians(self.alpha)
         cl_2d = self.cl
         cd_2d = self.cd
-        alpha_max_corr = radians(alpha_max_corr)
-        alpha_linear_min = radians(alpha_linear_min)
-        alpha_linear_max = radians(alpha_linear_max)
+        alpha_max_corr = np.radians(alpha_max_corr)
+        alpha_linear_min = np.radians(alpha_linear_min)
+        alpha_linear_max = np.radians(alpha_linear_max)
 
         # parameters in Du-Selig model
         a = 1
@@ -161,7 +160,7 @@ def correction3D(self, r_over_R, chord_over_r, tsr, alpha_max_corr=30, alpha_lin
         )
 
         # not sure where this adjustment comes from (besides AirfoilPrep spreadsheet of course)
-        adj = ((pi / 2 - alpha) / (pi / 2 - alpha_max_corr)) ** 2
+        adj = ((np.pi / 2 - alpha) / (np.pi / 2 - alpha_max_corr)) ** 2
         adj[alpha <= alpha_max_corr] = 1.0
 
         # Du-Selig correction for lift
@@ -227,43 +226,43 @@ def extrapolate(self, cdmax, AR=None, cdmin=0.001, nalpha=15):
         self.cdmax = max(max(self.cd), cdmax)
 
         # extract matching info from ends
-        alpha_high = radians(self.alpha[-1])
+        alpha_high = np.radians(self.alpha[-1])
         cl_high = self.cl[-1]
         cd_high = self.cd[-1]
         cm_high = self.cm[-1]
 
-        alpha_low = radians(self.alpha[0])
+        alpha_low = np.radians(self.alpha[0])
         cl_low = self.cl[0]
         cd_low = self.cd[0]
 
-        if alpha_high > pi / 2:
+        if alpha_high > np.pi / 2:
             raise Exception("alpha[-1] > pi/2")
             return self
-        if alpha_low < -pi / 2:
+        if alpha_low < -np.pi / 2:
             raise Exception("alpha[0] < -pi/2")
             return self
 
         # parameters used in model
-        sa = sin(alpha_high)
-        ca = cos(alpha_high)
+        sa = np.sin(alpha_high)
+        ca = np.cos(alpha_high)
         self.A = (cl_high - self.cdmax * sa * ca) * sa / ca ** 2
         self.B = (cd_high - self.cdmax * sa * sa) / ca
 
         # alpha_high <-> 90
-        alpha1 = np.linspace(alpha_high, pi / 2, nalpha)
+        alpha1 = np.linspace(alpha_high, np.pi / 2, nalpha)
         alpha1 = alpha1[1:]  # remove first element so as not to duplicate when concatenating
         cl1, cd1 = self.__Viterna(alpha1, 1.0)
 
         # 90 <-> 180-alpha_high
-        alpha2 = np.linspace(pi / 2, pi - alpha_high, nalpha)
+        alpha2 = np.linspace(np.pi / 2, np.pi - alpha_high, nalpha)
         alpha2 = alpha2[1:]
-        cl2, cd2 = self.__Viterna(pi - alpha2, -cl_adj)
+        cl2, cd2 = self.__Viterna(np.pi - alpha2, -cl_adj)
 
         # 180-alpha_high <-> 180
-        alpha3 = np.linspace(pi - alpha_high, pi, nalpha)
+        alpha3 = np.linspace(np.pi - alpha_high, np.pi, nalpha)
         alpha3 = alpha3[1:]
-        cl3, cd3 = self.__Viterna(pi - alpha3, 1.0)
-        cl3 = (alpha3 - pi) / alpha_high * cl_high * cl_adj  # override with linear variation
+        cl3, cd3 = self.__Viterna(np.pi - alpha3, 1.0)
+        cl3 = (alpha3 - np.pi) / alpha_high * cl_high * cl_adj  # override with linear variation
 
         if alpha_low <= -alpha_high:
             alpha4 = []
@@ -280,19 +279,19 @@ def extrapolate(self, cdmax, AR=None, cdmin=0.001, nalpha=15):
             alpha5max = -alpha_high
 
         # -90 <-> -alpha_high
-        alpha5 = np.linspace(-pi / 2, alpha5max, nalpha)
+        alpha5 = np.linspace(-np.pi / 2, alpha5max, nalpha)
         alpha5 = alpha5[1:]
         cl5, cd5 = self.__Viterna(-alpha5, -cl_adj)
 
         # -180+alpha_high <-> -90
-        alpha6 = np.linspace(-pi + alpha_high, -pi / 2, nalpha)
+        alpha6 = np.linspace(-np.pi + alpha_high, -np.pi / 2, nalpha)
         alpha6 = alpha6[1:]
-        cl6, cd6 = self.__Viterna(alpha6 + pi, cl_adj)
+        cl6, cd6 = self.__Viterna(alpha6 + np.pi, cl_adj)
 
         # -180 <-> -180 + alpha_high
-        alpha7 = np.linspace(-pi, -pi + alpha_high, nalpha)
-        cl7, cd7 = self.__Viterna(alpha7 + pi, 1.0)
-        cl7 = (alpha7 + pi) / alpha_high * cl_high * cl_adj  # linear variation
+        alpha7 = np.linspace(-np.pi, -np.pi + alpha_high, nalpha)
+        cl7, cd7 = self.__Viterna(alpha7 + np.pi, 1.0)
+        cl7 = (alpha7 + np.pi) / alpha_high * cl_high * cl_adj  # linear variation
 
         alpha = np.concatenate((alpha7, alpha6, alpha5, alpha4, np.radians(self.alpha), alpha1, alpha2, alpha3))
         cl = np.concatenate((cl7, cl6, cl5, cl4, self.cl, cl1, cl2, cl3))
@@ -301,8 +300,8 @@ def extrapolate(self, cdmax, AR=None, cdmin=0.001, nalpha=15):
         cd = np.maximum(cd, cdmin)  # don't allow negative drag coefficients
 
         # Setup alpha and cm to be used in extrapolation
-        cm1_alpha = floor(self.alpha[0] / 10.0) * 10.0
-        cm2_alpha = ceil(self.alpha[-1] / 10.0) * 10.0
+        cm1_alpha = np.floor(self.alpha[0] / 10.0) * 10.0
+        cm2_alpha = np.ceil(self.alpha[-1] / 10.0) * 10.0
         alpha_num = abs(int((-180.0 - cm1_alpha) / 10.0 - 1))
         alpha_cm1 = np.linspace(-180.0, cm1_alpha, alpha_num)
         alpha_cm2 = np.linspace(cm2_alpha, 180.0, int((180.0 - cm2_alpha) / 10.0 + 1))
@@ -355,9 +354,9 @@ def __CMCoeff(self, cl_high, cd_high, cm_high):
             p = -self.cl[0] / (self.cl[1] - self.cl[0])
             cm0 = self.cm[0] + p * (self.cm[1] - self.cm[0])
         self.cm0 = cm0
-        alpha_high = radians(self.alpha[-1])
-        XM = (-cm_high + cm0) / (cl_high * cos(alpha_high) + cd_high * sin(alpha_high))
-        cmCoef = (XM - 0.25) / tan((alpha_high - pi / 2))
+        alpha_high = np.radians(self.alpha[-1])
+        XM = (-cm_high + cm0) / (cl_high * np.cos(alpha_high) + cd_high * np.sin(alpha_high))
+        cmCoef = (XM - 0.25) / np.tan((alpha_high - np.pi / 2))
         return cmCoef
 
     def __getCM(self, i, cmCoef, alpha, cl_ext, cd_ext, alpha_low_deg, alpha_high_deg):
@@ -371,12 +370,15 @@ def __getCM(self, i, cmCoef, alpha, cl_ext, cd_ext, alpha_low_deg, alpha_high_de
                 cm_new = self.cm0
             else:
                 if alpha[i] > 0:
-                    x = cmCoef * tan(radians(alpha[i]) - pi / 2) + 0.25
-                    cm_new = self.cm0 - x * (cl_ext[i] * cos(radians(alpha[i])) + cd_ext[i] * sin(radians(alpha[i])))
+                    x = cmCoef * np.tan(np.radians(alpha[i]) - np.pi / 2) + 0.25
+                    cm_new = self.cm0 - x * (
+                        cl_ext[i] * np.cos(np.radians(alpha[i])) + cd_ext[i] * np.sin(np.radians(alpha[i]))
+                    )
                 else:
-                    x = cmCoef * tan(-radians(alpha[i]) - pi / 2) + 0.25
+                    x = cmCoef * np.tan(-np.radians(alpha[i]) - np.pi / 2) + 0.25
                     cm_new = -(
-                        self.cm0 - x * (-cl_ext[i] * cos(-radians(alpha[i])) + cd_ext[i] * sin(-radians(alpha[i])))
+                        self.cm0
+                        - x * (-cl_ext[i] * np.cos(-np.radians(alpha[i])) + cd_ext[i] * np.sin(-np.radians(alpha[i])))
                     )
         else:
             if alpha[i] == 165:
@@ -421,8 +423,8 @@ def unsteadyparam(self, alpha_linear_min=-5, alpha_linear_max=5):
         cl = self.cl
         cd = self.cd
 
-        alpha_linear_min = radians(alpha_linear_min)
-        alpha_linear_max = radians(alpha_linear_max)
+        alpha_linear_min = np.radians(alpha_linear_min)
+        alpha_linear_max = np.radians(alpha_linear_max)
 
         cn = cl * np.cos(alpha) + cd * np.sin(alpha)
 
@@ -468,7 +470,7 @@ def unsteadyparam(self, alpha_linear_min=-5, alpha_linear_max=5):
 
         # return: control setting, stall angle, alpha for 0 cn, cn slope,
         #         cn at stall+, cn at stall-, alpha for min CD, min(CD)
-        return (0.0, degrees(alphaU), degrees(alpha0), m, cnStallUpper, cnStallLower, alpha[minIdx], cd[minIdx])
+        return (0.0, np.degrees(alphaU), np.degrees(alpha0), m, cnStallUpper, cnStallLower, alpha[minIdx], cd[minIdx])
 
     def plot(self):
         """plot cl/cd/cm polar
@@ -581,6 +583,9 @@ def initFromAerodynFile(cls, aerodynFile, polarType=Polar):
                 if "EOT" in line:
                     break
                 data = [float(s) for s in line.split()]
+                if len(data) < 4:
+                    raise ValueError(f"Error: Expected 4 columns of data but found, {data}")
+
                 alpha.append(data[0])
                 cl.append(data[1])
                 cd.append(data[2])

wisdem/ccblade/ccblade.py

@@ -25,10 +25,9 @@
 import multiprocessing as mp
 
 import numpy as np
+import wisdem.ccblade._bem as _bem
 from scipy.optimize import brentq
 from scipy.interpolate import RectBivariateSpline, bisplev
-
-import wisdem.ccblade._bem as _bem
 from wisdem.airfoilprep import Airfoil
 
 # ------------------
@@ -79,6 +78,9 @@ def __init__(self, alpha, Re, cl, cd, cm=[], x=[], y=[], AFName="DEFAULTAF"):
 
             self.one_Re = True
 
+        if len(alpha) < 2:
+            raise ValueError(f"Need at least 2 angles of attack, but found {alpha}")
+
         kx = min(len(alpha) - 1, 3)
         ky = min(len(Re) - 1, 3)
 

wisdem/test/test_ccblade/test_gradients.py

@@ -12,12 +12,14 @@
 limitations under the License.
 """
 
+import os
 import unittest
-from os import path
 
 import numpy as np
 from wisdem.ccblade.ccblade import CCBlade, CCAirfoil
 
+basepath = os.path.join(os.path.dirname(os.path.realpath(__file__)), "../../../examples/_airfoil_files/")
+
 
 class TestGradients(unittest.TestCase):
     def setUp(self):
@@ -96,18 +98,17 @@ def setUp(self):
         self.mu = 1.81206e-5
 
         afinit = CCAirfoil.initFromAerodynFile  # just for shorthand
-        basepath = path.join(path.dirname(path.realpath(__file__)), "../../../examples/_airfoil_files/")
 
         # load all airfoils
         airfoil_types = [0] * 8
-        airfoil_types[0] = afinit(basepath + "Cylinder1.dat")
-        airfoil_types[1] = afinit(basepath + "Cylinder2.dat")
-        airfoil_types[2] = afinit(basepath + "DU40_A17.dat")
-        airfoil_types[3] = afinit(basepath + "DU35_A17.dat")
-        airfoil_types[4] = afinit(basepath + "DU30_A17.dat")
-        airfoil_types[5] = afinit(basepath + "DU25_A17.dat")
-        airfoil_types[6] = afinit(basepath + "DU21_A17.dat")
-        airfoil_types[7] = afinit(basepath + "NACA64_A17.dat")
+        airfoil_types[0] = afinit(basepath + os.sep + "Cylinder1.dat")
+        airfoil_types[1] = afinit(basepath + os.sep + "Cylinder2.dat")
+        airfoil_types[2] = afinit(basepath + os.sep + "DU40_A17.dat")
+        airfoil_types[3] = afinit(basepath + os.sep + "DU35_A17.dat")
+        airfoil_types[4] = afinit(basepath + os.sep + "DU30_A17.dat")
+        airfoil_types[5] = afinit(basepath + os.sep + "DU25_A17.dat")
+        airfoil_types[6] = afinit(basepath + os.sep + "DU21_A17.dat")
+        airfoil_types[7] = afinit(basepath + os.sep + "NACA64_A17.dat")
 
         # place at appropriate radial stations
         af_idx = [0, 0, 1, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7]
@@ -3684,18 +3685,17 @@ def setUp(self):
         self.mu = 1.81206e-5
 
         afinit = CCAirfoil.initFromAerodynFile  # just for shorthand
-        basepath = path.join(path.dirname(path.realpath(__file__)), "../../../examples/_airfoil_files/")
 
         # load all airfoils
         airfoil_types = [0] * 8
-        airfoil_types[0] = afinit(basepath + "Cylinder1.dat")
-        airfoil_types[1] = afinit(basepath + "Cylinder2.dat")
-        airfoil_types[2] = afinit(basepath + "DU40_A17.dat")
-        airfoil_types[3] = afinit(basepath + "DU35_A17.dat")
-        airfoil_types[4] = afinit(basepath + "DU30_A17.dat")
-        airfoil_types[5] = afinit(basepath + "DU25_A17.dat")
-        airfoil_types[6] = afinit(basepath + "DU21_A17.dat")
-        airfoil_types[7] = afinit(basepath + "NACA64_A17.dat")
+        airfoil_types[0] = afinit(basepath + os.sep + "Cylinder1.dat")
+        airfoil_types[1] = afinit(basepath + os.sep + "Cylinder2.dat")
+        airfoil_types[2] = afinit(basepath + os.sep + "DU40_A17.dat")
+        airfoil_types[3] = afinit(basepath + os.sep + "DU35_A17.dat")
+        airfoil_types[4] = afinit(basepath + os.sep + "DU30_A17.dat")
+        airfoil_types[5] = afinit(basepath + os.sep + "DU25_A17.dat")
+        airfoil_types[6] = afinit(basepath + os.sep + "DU21_A17.dat")
+        airfoil_types[7] = afinit(basepath + os.sep + "NACA64_A17.dat")
 
         # place at appropriate radial stations
         af_idx = [0, 0, 1, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7]
@@ -4646,18 +4646,17 @@ def setUp(self):
         self.mu = 1.81206e-5
 
         afinit = CCAirfoil.initFromAerodynFile  # just for shorthand
-        basepath = path.join(path.dirname(path.realpath(__file__)), "../../../examples/_airfoil_files/")
 
         # load all airfoils
         airfoil_types = [0] * 8
-        airfoil_types[0] = afinit(basepath + "Cylinder1.dat")
-        airfoil_types[1] = afinit(basepath + "Cylinder2.dat")
-        airfoil_types[2] = afinit(basepath + "DU40_A17.dat")
-        airfoil_types[3] = afinit(basepath + "DU35_A17.dat")
-        airfoil_types[4] = afinit(basepath + "DU30_A17.dat")
-        airfoil_types[5] = afinit(basepath + "DU25_A17.dat")
-        airfoil_types[6] = afinit(basepath + "DU21_A17.dat")
-        airfoil_types[7] = afinit(basepath + "NACA64_A17.dat")
+        airfoil_types[0] = afinit(basepath + os.sep + "Cylinder1.dat")
+        airfoil_types[1] = afinit(basepath + os.sep + "Cylinder2.dat")
+        airfoil_types[2] = afinit(basepath + os.sep + "DU40_A17.dat")
+        airfoil_types[3] = afinit(basepath + os.sep + "DU35_A17.dat")
+        airfoil_types[4] = afinit(basepath + os.sep + "DU30_A17.dat")
+        airfoil_types[5] = afinit(basepath + os.sep + "DU25_A17.dat")
+        airfoil_types[6] = afinit(basepath + os.sep + "DU21_A17.dat")
+        airfoil_types[7] = afinit(basepath + os.sep + "NACA64_A17.dat")
 
         # place at appropriate radial stations
         af_idx = [0, 0, 1, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7]
@@ -5113,18 +5112,17 @@ def setUp(self):
         self.mu = 1.81206e-5
 
         afinit = CCAirfoil.initFromAerodynFile  # just for shorthand
-        basepath = path.join(path.dirname(path.realpath(__file__)), "../../../examples/_airfoil_files/")
 
         # load all airfoils
         airfoil_types = [0] * 8
-        airfoil_types[0] = afinit(basepath + "Cylinder1.dat")
-        airfoil_types[1] = afinit(basepath + "Cylinder2.dat")
-        airfoil_types[2] = afinit(basepath + "DU40_A17.dat")
-        airfoil_types[3] = afinit(basepath + "DU35_A17.dat")
-        airfoil_types[4] = afinit(basepath + "DU30_A17.dat")
-        airfoil_types[5] = afinit(basepath + "DU25_A17.dat")
-        airfoil_types[6] = afinit(basepath + "DU21_A17.dat")
-        airfoil_types[7] = afinit(basepath + "NACA64_A17.dat")
+        airfoil_types[0] = afinit(basepath + os.sep + "Cylinder1.dat")
+        airfoil_types[1] = afinit(basepath + os.sep + "Cylinder2.dat")
+        airfoil_types[2] = afinit(basepath + os.sep + "DU40_A17.dat")
+        airfoil_types[3] = afinit(basepath + os.sep + "DU35_A17.dat")
+        airfoil_types[4] = afinit(basepath + os.sep + "DU30_A17.dat")
+        airfoil_types[5] = afinit(basepath + os.sep + "DU25_A17.dat")
+        airfoil_types[6] = afinit(basepath + os.sep + "DU21_A17.dat")
+        airfoil_types[7] = afinit(basepath + os.sep + "NACA64_A17.dat")
 
         # place at appropriate radial stations
         af_idx = [0, 0, 1, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7]

wisdem/test/test_ccblade/test_om_gradients.py

@@ -110,8 +110,9 @@ def test_ccblade_loads(self):
                 if "airfoil" not in input_name and "rho" not in input_name and "mu" not in input_name:
                     new_check[comp_name][(output_name, input_name)] = check[comp_name][(output_name, input_name)]
 
-        assert_check_partials(new_check, rtol=5e-5, atol=10.)
+        assert_check_partials(new_check, rtol=5e-5, atol=10.0)
 
+    @unittest.skip("Not useful and now OpenMDAO complains")
     def test_ccblade_twist(self):
         """
         Right now this just compares fd to fd so it is not a meaningful test.
@@ -194,15 +195,15 @@ def test_ccblade_twist(self):
 
         prob.run_model()
 
-        check = prob.check_partials(out_stream=None, compact_print=True)
+        check = prob.check_partials(out_stream=None, compact_print=True, step=1e-7)
 
         # Manually filter some entries out of the assert_check_partials call.
         # Will want to add this functionality to OpenMDAO itself at some point.
         new_check = {}
         for comp_name in check:
             new_check[comp_name] = {}
             for (output_name, input_name) in check[comp_name]:
-                if "airfoil" not in input_name and "rho" not in input_name and "mu" not in input_name:
+                if not input_name in ["airfoil", "rho", "mu"]:
                     new_check[comp_name][(output_name, input_name)] = check[comp_name][(output_name, input_name)]
 
         assert_check_partials(new_check)  # , rtol=5e-5, atol=1e-4)
@@ -290,7 +291,7 @@ def test_ccblade_standalone(self):
                 if "airfoil" not in input_name and "rho" not in input_name and "mu" not in input_name:
                     new_check[comp_name][(output_name, input_name)] = check[comp_name][(output_name, input_name)]
 
-        assert_check_partials(new_check, rtol=5e-4, atol=50.)
+        assert_check_partials(new_check, rtol=5e-4, atol=50.0)
 
 
 def suite():