Fix bug in wave drag derivatives (#432)

* Fix bug in wave drag derivatives

* black

* Ignore OpenMDAO reports folders

* tests

* tweak scaneagle test

* bump patch version

* Bump oldest supported versions of things

* Forgot to bump `OPENMDAO_VERSION_OLDEST`

* Need openmdao 3.22 to work with numpy 1.24

* Try running tests on our docker images

* Revert "Try running tests on our docker images"

This reverts commit 9ef3db2.

* Put it all back how it was

* Use older pip for "oldest" test

* Revert "Bump oldest supported versions of things"

This reverts commit 5bf3e8e.

* Add `PIP_VERSION_OLDEST` variable

* Change optimisation failure test to work with older OpenMDAO

* Modify scaneagle tolerance

* Remove commented out old code

Author: A-CGray

.github/workflows/oas.yml

@@ -20,12 +20,13 @@ jobs:
   test:
     runs-on: ubuntu-latest
     strategy:
-      fail-fast: false  # continue other jobs even if one of the jobs in matrix fails 
+      fail-fast: false  # continue other jobs even if one of the jobs in matrix fails
       matrix:
         dep-versions: ["oldest", "latest"]
         # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Set versions to test here ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         PYTHON_VERSION_OLDEST: ['3.8']
         PYTHON_VERSION_LATEST: ['3.11']
+        PIP_VERSION_OLDEST: ['23.3.2'] # OpenMDAO 3.15 doesn't install with pip>=24
         WHEEL_VERSION_OLDEST: ['0.38.4']  # latest wheel cannot build the oldest OpenMDAO (3.17 or older)
         NUMPY_VERSION_OLDEST: ['1.20']  # latest is most recent on PyPI
         SCIPY_VERSION_OLDEST: ['1.6.0']  # latest is most recent on PyPI
@@ -94,6 +95,7 @@ jobs:
     - name: Install OAS and its dependencies (oldest versions)
       if: ${{ matrix.dep-versions == 'oldest' }}
       run: |
+        python -m pip install pip==${{ matrix.PIP_VERSION_OLDEST }}
         python -m pip install wheel==${{ matrix.WHEEL_VERSION_OLDEST }}
         pip install numpy==${{ matrix.NUMPY_VERSION_OLDEST }} scipy==${{ matrix.SCIPY_VERSION_OLDEST }} openmdao==${{ matrix.OPENMDAO_VERSION_OLDEST }} mphys==${{ matrix.MPHYS_VERSION_OLDEST }}
         pip install -e .[test]

.gitignore

@@ -1,3 +1,6 @@
+# Ignore OpenMDAO reports folders
+**/reports/
+
 *.pyc
 *.py~
 *.db
@@ -19,6 +22,6 @@ openaerostruct/docs/_srcdocs/
 *.msg
 src/adjoint/tempReverse/
 src/adjoint/tempForward/
-openaerostruct.egg-info/ 
+openaerostruct.egg-info/
 .vscode
 *.plt

README.md

@@ -94,7 +94,7 @@ Version Information
 The oldest and latest versions of the dependencies that we test regularly are the following (other versions may work, but no guarantees):
 
 | Dependency         | oldest | latest |
-|--------------------|--------| ------ |
+| ------------------ | ------ | ------ |
 | Python             | 3.8    | 3.11   |
 | NumPy              | 1.20   | latest |
 | SciPy              | 1.6.0  | latest |

openaerostruct/__init__.py

@@ -1 +1 @@
-__version__ = "2.7.2"
+__version__ = "2.7.3"

openaerostruct/aerodynamics/wave_drag.py

@@ -67,14 +67,15 @@ def compute(self, inputs, outputs):
             chords = inputs["chords"]
             CL = inputs["CL"]
 
-            mean_chords = (chords[:-1] + chords[1:]) / 2.0
-            panel_areas = mean_chords * widths
-            avg_cos_sweep = np.sum(cos_sweep * panel_areas) / np.sum(panel_areas)  # weighted average of 1/4 chord sweep
-            avg_t_over_c = np.sum(t_over_c * panel_areas) / np.sum(panel_areas)  # weighted average of streamwise t/c
+            panel_mid_chords = (chords[:-1] + chords[1:]) / 2.0
+            panel_areas = panel_mid_chords * widths
+            sum_panel_areas = np.sum(panel_areas)
+            avg_cos_sweep = np.sum(cos_sweep * panel_areas) / sum_panel_areas  # weighted average of 1/4 chord sweep
+            avg_t_over_c = np.sum(t_over_c * panel_areas) / sum_panel_areas  # weighted average of streamwise t/c
             MDD = self.ka / avg_cos_sweep - avg_t_over_c / avg_cos_sweep**2 - CL / (10 * avg_cos_sweep**3)
             Mcrit = MDD - (0.1 / 80.0) ** (1.0 / 3.0)
 
-            if M > Mcrit:
+            if np.real(M) > np.real(Mcrit):
                 outputs["CDw"] = 20 * (M - Mcrit) ** 4
             else:
                 outputs["CDw"] = 0.0
@@ -86,6 +87,14 @@ def compute(self, inputs, outputs):
 
     def compute_partials(self, inputs, partials):
         """Jacobian for wave drag."""
+        # Explicitly zero out the partials to begin, we can't assume the input partials arrays contain zeros already
+        partials["CDw", "CL"][:] = 0.0
+        partials["CDw", "lengths_spanwise"][:] = 0.0
+        partials["CDw", "widths"][:] = 0.0
+        partials["CDw", "Mach_number"][:] = 0.0
+        partials["CDw", "chords"][:] = 0.0
+        partials["CDw", "t_over_c"][:] = 0.0
+
         if self.with_wave:
             ny = self.surface["mesh"].shape[1]
             t_over_c = inputs["t_over_c"]
@@ -96,16 +105,16 @@ def compute_partials(self, inputs, partials):
             chords = inputs["chords"]
             CL = inputs["CL"]
 
-            chords = (chords[:-1] + chords[1:]) / 2.0
-            panel_areas = chords * widths
+            panel_mid_chords = (chords[:-1] + chords[1:]) / 2.0
+            panel_areas = panel_mid_chords * widths
             sum_panel_areas = np.sum(panel_areas)
             avg_cos_sweep = np.sum(cos_sweep * panel_areas) / sum_panel_areas
             avg_t_over_c = np.sum(t_over_c * panel_areas) / sum_panel_areas
 
-            MDD = 0.95 / avg_cos_sweep - avg_t_over_c / avg_cos_sweep**2 - CL / (10 * avg_cos_sweep**3)
+            MDD = self.ka / avg_cos_sweep - avg_t_over_c / avg_cos_sweep**2 - CL / (10 * avg_cos_sweep**3)
             Mcrit = MDD - (0.1 / 80.0) ** (1.0 / 3.0)
 
-            if M > Mcrit:
+            if np.real(M) > np.real(Mcrit):
                 dCDwdMDD = -80 * (M - Mcrit) ** 3
                 dMDDdCL = -1.0 / (10 * avg_cos_sweep**3)
                 dMDDdavg = (-10 * self.ka * avg_cos_sweep**2 + 20 * avg_t_over_c * avg_cos_sweep + 3 * CL) / (
@@ -114,14 +123,19 @@ def compute_partials(self, inputs, partials):
                 dMDDdtoc = -1.0 / (avg_cos_sweep**2)
                 dtocavgdtoc = panel_areas / sum_panel_areas
 
-                ccos = np.sum(widths * chords)
-                ccos2w = np.sum(chords * widths**2 / lengths_spanwise)
+                ccos = np.sum(widths * panel_mid_chords)
+                ccos2w = np.sum(panel_mid_chords * widths**2 / lengths_spanwise)
 
-                davgdcos = 2 * chords * widths / lengths_spanwise / ccos - chords * ccos2w / ccos**2
-                dtocdcos = chords * t_over_c / ccos - chords * np.sum(chords * widths * t_over_c) / ccos**2
-                davgdw = -1 * chords * widths**2 / lengths_spanwise**2 / ccos
+                davgdcos = (
+                    2 * panel_mid_chords * widths / lengths_spanwise / ccos - panel_mid_chords * ccos2w / ccos**2
+                )
+                dtocdcos = (
+                    panel_mid_chords * t_over_c / ccos
+                    - panel_mid_chords * np.sum(panel_mid_chords * widths * t_over_c) / ccos**2
+                )
+                davgdw = -1 * panel_mid_chords * widths**2 / lengths_spanwise**2 / ccos
                 davgdc = widths**2 / lengths_spanwise / ccos - widths * ccos2w / ccos**2
-                dtocdc = t_over_c * widths / ccos - widths * np.sum(chords * widths * t_over_c) / ccos**2
+                dtocdc = t_over_c * widths / ccos - widths * np.sum(panel_mid_chords * widths * t_over_c) / ccos**2
 
                 dcdchords = np.zeros((ny - 1, ny))
                 i, j = np.indices(dcdchords.shape)

openaerostruct/docs/installation.rst

@@ -70,7 +70,7 @@ To run the tests on your machine, use the [test] option. This will install the p
 Then run the tests from the OpenAeroStruct root directory by calling:
 
 .. code-block:: bash
-  
+
     testflo -v .
 
 To install the dependencies to build the documentation locally, run:

openaerostruct/tests/test_aerostruct_wingbox_wave_fuel_vol_constraint_opt.py

@@ -289,7 +289,7 @@ def test(self):
             # Airfoil properties for viscous drag calculation
             "k_lam": 0.05,  # percentage of chord with laminar
             # flow, used for viscous drag
-            "t_over_c_cp": np.array([0.08, 0.08, 0.08, 0.10, 0.10, 0.08]),
+            "t_over_c_cp": 0.12 * np.ones(6),
             "original_wingbox_airfoil_t_over_c": 0.12,
             "c_max_t": 0.38,  # chordwise location of maximum thickness
             "with_viscous": True,
@@ -317,10 +317,11 @@ def test(self):
 
         # Add problem information as an independent variables component
         indep_var_comp = om.IndepVarComp()
-        indep_var_comp.add_output("v", val=0.85 * 295.07, units="m/s")
+        Mach = 0.77
+        indep_var_comp.add_output("v", val=Mach * 295.07, units="m/s")
         indep_var_comp.add_output("alpha", val=0.0, units="deg")
-        indep_var_comp.add_output("Mach_number", val=0.85)
-        indep_var_comp.add_output("re", val=0.348 * 295.07 * 0.85 * 1.0 / (1.43 * 1e-5), units="1/m")
+        indep_var_comp.add_output("Mach_number", val=Mach)
+        indep_var_comp.add_output("re", val=0.348 * 295.07 * Mach * 1.0 / (1.43 * 1e-5), units="1/m")
         indep_var_comp.add_output("rho", val=0.348, units="kg/m**3")
         indep_var_comp.add_output("CT", val=0.53 / 3600, units="1/s")
         indep_var_comp.add_output("R", val=14.307e6, units="m")
@@ -412,7 +413,7 @@ def test(self):
         prob.model.add_design_var("wing.twist_cp", lower=-15.0, upper=15.0, scaler=0.1)
         prob.model.add_design_var("wing.spar_thickness_cp", lower=0.003, upper=0.1, scaler=1e2)
         prob.model.add_design_var("wing.skin_thickness_cp", lower=0.003, upper=0.1, scaler=1e2)
-        prob.model.add_design_var("wing.geometry.t_over_c_cp", lower=0.07, upper=0.2, scaler=10.0)
+        prob.model.add_design_var("wing.geometry.t_over_c_cp", lower=0.06, upper=0.2, scaler=10.0)
 
         prob.model.add_constraint("AS_point_0.CL", equals=0.5)
         prob.model.add_constraint("AS_point_0.wing_perf.failure", upper=0.0)
@@ -427,10 +428,11 @@ def test(self):
         # Set up the problem
         prob.setup()
 
-        prob.run_driver()
+        optFailed = prob.run_driver()
 
-        assert_near_equal(prob["AS_point_0.fuelburn"][0], 78292.730861421, 1e-5)
-        assert_near_equal(prob["wing.structural_mass"][0] / 1.25, 16168.330307591294, 1e-5)
+        self.assertFalse(optFailed)
+        assert_near_equal(prob["AS_point_0.fuelburn"][0], 85348.88283214, 1e-5)
+        assert_near_equal(prob["wing.structural_mass"][0], 13029.71120634, 1e-5)
 
 
 if __name__ == "__main__":

openaerostruct/tests/test_scaneagle.py

@@ -8,10 +8,11 @@
 from openaerostruct.geometry.utils import generate_mesh
 from openaerostruct.integration.aerostruct_groups import AerostructGeometry, AerostructPoint
 from openaerostruct.utils.constants import grav_constant
+from openaerostruct.utils.testing import assert_check_totals
 
 
 class Test(unittest.TestCase):
-    def test(self):
+    def setUp(self):
         # Total number of nodes to use in the spanwise (num_y) and
         # chordwise (num_x) directions. Vary these to change the level of fidelity.
         num_y = 21
@@ -86,7 +87,7 @@ def test(self):
             "c_max_t": 0.303,  # chordwise location of maximum (NACA0015)
             # thickness
             "with_viscous": True,
-            "with_wave": False,  # if true, compute wave drag
+            "with_wave": True,  # if true, compute wave drag
             # Material properties taken from http://www.performance-composites.com/carbonfibre/mechanicalproperties_2.asp
             "E": 85.0e9,
             "G": 25.0e9,
@@ -197,16 +198,33 @@ def test(self):
         # We're trying to minimize fuel burn
         prob.model.add_objective("AS_point_0.fuelburn", scaler=0.1)
 
-        # Set up the problem
-        prob.setup()
+        self.prob = prob
 
+    def test_opt(self):
+        # Set up the problem
+        self.prob.setup()
         # Actually run the optimization problem
-        prob.run_driver()
+        self.prob.run_driver()
+
+        assert_near_equal(self.prob["AS_point_0.fuelburn"][0], 4.6365011384888275, 1e-5)
+        assert_near_equal(self.prob["wing.twist_cp"], np.array([2.25819837, 10.39881572, 5.0]), 1e-5)
+        assert_near_equal(self.prob["wing.sweep"][0], 18.964409030629632, 1e-5)
+        assert_near_equal(self.prob["alpha"][0], 2.0366563718492547, 1e-5)
 
-        assert_near_equal(prob["AS_point_0.fuelburn"][0], 4.6365011384888275, 1e-5)
-        assert_near_equal(prob["wing.twist_cp"], np.array([2.25819837, 10.39881572, 5.0]), 1e-5)
-        assert_near_equal(prob["wing.sweep"][0], 18.964409030629632, 1e-5)
-        assert_near_equal(prob["alpha"][0], 2.0366563718492547, 1e-5)
+    def test_totals(self):
+        # Set up the problem
+        self.prob.setup()
+        self.prob.run_model()
+        totals = self.prob.check_totals(
+            method="fd",
+            form="central",
+            step=5e-4,
+            step_calc="rel",
+            compact_print=True,
+            abs_err_tol=1e-4,
+            rel_err_tol=1e-4,
+        )
+        assert_check_totals(totals, atol=1e-4, rtol=1e-4)
 
 
 if __name__ == "__main__":