Skip to content

Commit

Permalink
Fix bug in wave drag derivatives (#432)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
A-CGray authored Jul 16, 2024
1 parent 0448bff commit 97be635
Show file tree
Hide file tree
Showing 8 changed files with 76 additions and 37 deletions.
4 changes: 3 additions & 1 deletion .github/workflows/oas.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand Down
5 changes: 4 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
# Ignore OpenMDAO reports folders
**/reports/

*.pyc
*.py~
*.db
Expand All @@ -19,6 +22,6 @@ openaerostruct/docs/_srcdocs/
*.msg
src/adjoint/tempReverse/
src/adjoint/tempForward/
openaerostruct.egg-info/
openaerostruct.egg-info/
.vscode
*.plt
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 |
Expand Down
2 changes: 1 addition & 1 deletion openaerostruct/__init__.py
Original file line number Diff line number Diff line change
@@ -1 +1 @@
__version__ = "2.7.2"
__version__ = "2.7.3"
44 changes: 29 additions & 15 deletions openaerostruct/aerodynamics/wave_drag.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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"]
Expand All @@ -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) / (
Expand All @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion openaerostruct/docs/installation.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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)
Expand All @@ -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__":
Expand Down
36 changes: 27 additions & 9 deletions openaerostruct/tests/test_scaneagle.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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__":
Expand Down

0 comments on commit 97be635

Please sign in to comment.