diff --git a/geom_files/aircraft.py b/geom_files/aircraft.py new file mode 100644 index 0000000..fbf0647 --- /dev/null +++ b/geom_files/aircraft.py @@ -0,0 +1,61 @@ +"""Auto-generated geometry file +Generated by make_aircraft.py - Aircraft configuration with wing and tail +""" +import numpy as np + +input_dict = { + "title": "AIRCRAFT 1", + "mach": np.float64(0.1), + "iysym": np.int32(0), + "izsym": np.int32(0), + "zsym": np.float64(0.0), + "Sref": np.float64(1.5825), + "Cref": np.float64(0.35), + "Bref": np.float64(4.0), + "XYZref": np.array([0.084, 0.0, 0.0], dtype=np.float64), + "CDp": np.float64(0.0116), + "surfaces": { + "Wing": { + "yduplicate": np.float64(0.0), + "scale": np.array([1.0, 1.0, 1.0], dtype=np.float64), + "translate": np.array([0.0, 0.0, 0.0], dtype=np.float64), + "mesh": np.array([[[0.0, -0.0, 0.0], [0.0, 0.15000000000000002, 0.0], [0.0, 0.30000000000000004, 0.0], [0.0, 0.45, 0.0], [0.0, 0.6, 0.0], [0.0, 0.75, 0.0], [0.0, 0.86, 0.0], [0.0, 0.97, 0.0], [0.0, 1.08, 0.0], [0.0, 1.19, 0.0], [0.0, 1.3, 0.0], [0.0, 1.3900000000000001, 0.0], [0.0, 1.48, 0.0], [0.0, 1.57, 0.0], [0.0, 1.66, 0.0], [0.0, 1.75, 0.0], [0.0, 1.8, 0.020000000000000004], [0.0, 1.85, 0.04000000000000001], [0.0, 1.9, 0.060000000000000005], [0.0, 1.95, 0.08], [0.0, 2.0, 0.1]], [[0.0642857142857143, -0.0, 0.0], [0.0642857142857143, 0.15000000000000002, 0.0], [0.0642857142857143, 0.30000000000000004, 0.0], [0.0642857142857143, 0.45, 0.0], [0.0642857142857143, 0.6, 0.0], [0.0642857142857143, 0.75, 0.0], [0.06285600000000001, 0.86, 0.0], [0.061426285714285725, 0.97, 0.0], [0.059996571428571434, 1.08, 0.0], [0.05856685714285715, 1.19, 0.0], [0.057137142857142866, 1.3, 0.0], [0.05428028571428572, 1.3900000000000001, 0.0], [0.05142342857142858, 1.48, 0.0], [0.04856657142857143, 1.57, 0.0], [0.04570971428571429, 1.66, 0.0], [0.04285285714285715, 1.75, 0.0], [0.039995428628571424, 1.8, 0.020000000000000004], [0.03713800011428571, 1.85, 0.04000000000000001], [0.034280571600000004, 1.9, 0.060000000000000005], [0.031423143085714283, 1.95, 0.08], [0.02856571457142857, 2.0, 0.1]], [[0.1285714285714286, -0.0, 0.0], [0.1285714285714286, 0.15000000000000002, 0.0], [0.1285714285714286, 0.30000000000000004, 0.0], [0.1285714285714286, 0.45, 0.0], [0.1285714285714286, 0.6, 0.0], [0.1285714285714286, 0.75, 0.0], [0.12571200000000002, 0.86, 0.0], [0.12285257142857145, 0.97, 0.0], [0.11999314285714287, 1.08, 0.0], [0.1171337142857143, 1.19, 0.0], [0.11427428571428573, 1.3, 0.0], [0.10856057142857144, 1.3900000000000001, 0.0], [0.10284685714285716, 1.48, 0.0], [0.09713314285714286, 1.57, 0.0], [0.09141942857142858, 1.66, 0.0], [0.0857057142857143, 1.75, 0.0], [0.07999085725714285, 1.8, 0.020000000000000004], [0.07427600022857142, 1.85, 0.04000000000000001], [0.06856114320000001, 1.9, 0.060000000000000005], [0.06284628617142857, 1.95, 0.08], [0.05713142914285714, 2.0, 0.1]], [[0.1928571428571429, -0.0, 0.0], [0.1928571428571429, 0.15000000000000002, 0.0], [0.1928571428571429, 0.30000000000000004, 0.0], [0.1928571428571429, 0.45, 0.0], [0.1928571428571429, 0.6, 0.0], [0.1928571428571429, 0.75, 0.0], [0.18856800000000004, 0.86, 0.0], [0.1842788571428572, 0.97, 0.0], [0.1799897142857143, 1.08, 0.0], [0.17570057142857146, 1.19, 0.0], [0.1714114285714286, 1.3, 0.0], [0.16284085714285718, 1.3900000000000001, 0.0], [0.15427028571428575, 1.48, 0.0], [0.1456997142857143, 1.57, 0.0], [0.1371291428571429, 1.66, 0.0], [0.12855857142857144, 1.75, 0.0], [0.11998628588571429, 1.8, 0.020000000000000004], [0.11141400034285713, 1.85, 0.04000000000000001], [0.10284171480000003, 1.9, 0.060000000000000005], [0.09426942925714288, 1.95, 0.08], [0.08569714371428572, 2.0, 0.1]], [[0.2571428571428572, -0.0, 0.0], [0.2571428571428572, 0.15000000000000002, 0.0], [0.2571428571428572, 0.30000000000000004, 0.0], [0.2571428571428572, 0.45, 0.0], [0.2571428571428572, 0.6, 0.0], [0.2571428571428572, 0.75, 0.0], [0.25142400000000004, 0.86, 0.0], [0.2457051428571429, 0.97, 0.0], [0.23998628571428574, 1.08, 0.0], [0.2342674285714286, 1.19, 0.0], [0.22854857142857146, 1.3, 0.0], [0.21712114285714287, 1.3900000000000001, 0.0], [0.2056937142857143, 1.48, 0.0], [0.19426628571428572, 1.57, 0.0], [0.18283885714285716, 1.66, 0.0], [0.1714114285714286, 1.75, 0.0], [0.1599817145142857, 1.8, 0.020000000000000004], [0.14855200045714284, 1.85, 0.04000000000000001], [0.13712228640000002, 1.9, 0.060000000000000005], [0.12569257234285713, 1.95, 0.08], [0.11426285828571428, 2.0, 0.1]], [[0.32142857142857145, -0.0, 0.0], [0.32142857142857145, 0.15000000000000002, 0.0], [0.32142857142857145, 0.30000000000000004, 0.0], [0.32142857142857145, 0.45, 0.0], [0.32142857142857145, 0.6, 0.0], [0.32142857142857145, 0.75, 0.0], [0.31428, 0.86, 0.0], [0.3071314285714286, 0.97, 0.0], [0.29998285714285716, 1.08, 0.0], [0.29283428571428577, 1.19, 0.0], [0.2856857142857143, 1.3, 0.0], [0.27140142857142857, 1.3900000000000001, 0.0], [0.2571171428571429, 1.48, 0.0], [0.24283285714285716, 1.57, 0.0], [0.22854857142857146, 1.66, 0.0], [0.21426428571428574, 1.75, 0.0], [0.19997714314285714, 1.8, 0.020000000000000004], [0.18569000057142854, 1.85, 0.04000000000000001], [0.17140285800000002, 1.9, 0.060000000000000005], [0.15711571542857145, 1.95, 0.08], [0.14282857285714284, 2.0, 0.1]], [[0.3857142857142858, -0.0, 0.0], [0.3857142857142858, 0.15000000000000002, 0.0], [0.3857142857142858, 0.30000000000000004, 0.0], [0.3857142857142858, 0.45, 0.0], [0.3857142857142858, 0.6, 0.0], [0.3857142857142858, 0.75, 0.0], [0.3771360000000001, 0.86, 0.0], [0.3685577142857144, 0.97, 0.0], [0.3599794285714286, 1.08, 0.0], [0.3514011428571429, 1.19, 0.0], [0.3428228571428572, 1.3, 0.0], [0.32568171428571435, 1.3900000000000001, 0.0], [0.3085405714285715, 1.48, 0.0], [0.2913994285714286, 1.57, 0.0], [0.2742582857142858, 1.66, 0.0], [0.2571171428571429, 1.75, 0.0], [0.23997257177142858, 1.8, 0.020000000000000004], [0.22282800068571426, 1.85, 0.04000000000000001], [0.20568342960000005, 1.9, 0.060000000000000005], [0.18853885851428576, 1.95, 0.08], [0.17139428742857143, 2.0, 0.1]], [[0.45, -0.0, 0.0], [0.45, 0.15000000000000002, 0.0], [0.45, 0.30000000000000004, 0.0], [0.45, 0.45, 0.0], [0.45, 0.6, 0.0], [0.45, 0.75, 0.0], [0.439992, 0.86, 0.0], [0.42998400000000003, 0.97, 0.0], [0.419976, 1.08, 0.0], [0.40996800000000005, 1.19, 0.0], [0.39996000000000004, 1.3, 0.0], [0.379962, 1.3900000000000001, 0.0], [0.35996400000000006, 1.48, 0.0], [0.33996600000000005, 1.57, 0.0], [0.31996800000000003, 1.66, 0.0], [0.29997, 1.75, 0.0], [0.27996800039999997, 1.8, 0.020000000000000004], [0.2599660008, 1.85, 0.04000000000000001], [0.23996400120000005, 1.9, 0.060000000000000005], [0.2199620016, 1.95, 0.08], [0.199960002, 2.0, 0.1]]], dtype=np.float64), + "afiles": "../geom_files/airfoils/ag40d.dat", + }, + "Horizontal Tail": { + "yduplicate": np.float64(0.0), + "claf": 1.1078, + "scale": np.array([1.0, 1.0, 1.0], dtype=np.float64), + "translate": np.array([0.0, 0.0, 0.0], dtype=np.float64), + "mesh": np.array([[[1.60247523857, -0.0, 0.0], [1.60247523857, 0.052000000000000005, 0.0], [1.60247523857, 0.10400000000000001, 0.0], [1.60247523857, 0.15600000000000003, 0.0], [1.60247523857, 0.20800000000000002, 0.0], [1.60247523857, 0.26, 0.0], [1.60247523857, 0.31200000000000006, 0.0], [1.60247523857, 0.36400000000000005, 0.0], [1.60247523857, 0.41600000000000004, 0.0], [1.60247523857, 0.468, 0.0], [1.60247523857, 0.52, 0.0]], [[1.6104285546720376, -0.0, 0.0], [1.6104285546720376, 0.052000000000000005, 0.0], [1.6104285546720376, 0.10400000000000001, 0.0], [1.6104285546720376, 0.15600000000000003, 0.0], [1.6104285546720376, 0.20800000000000002, 0.0], [1.6104285546720376, 0.26, 0.0], [1.6104285546720376, 0.31200000000000006, 0.0], [1.6104285546720376, 0.36400000000000005, 0.0], [1.6104285546720376, 0.41600000000000004, 0.0], [1.6104285546720376, 0.468, 0.0], [1.6104285546720376, 0.52, 0.0]], [[1.633509976984071, -0.0, 0.0], [1.633509976984071, 0.052000000000000005, 0.0], [1.633509976984071, 0.10400000000000001, 0.0], [1.633509976984071, 0.15600000000000003, 0.0], [1.633509976984071, 0.20800000000000002, 0.0], [1.633509976984071, 0.26, 0.0], [1.633509976984071, 0.31200000000000006, 0.0], [1.633509976984071, 0.36400000000000005, 0.0], [1.633509976984071, 0.41600000000000004, 0.0], [1.633509976984071, 0.468, 0.0], [1.633509976984071, 0.52, 0.0]], [[1.6694601350724732, -0.0, 0.0], [1.6694601350724732, 0.052000000000000005, 0.0], [1.6694601350724732, 0.10400000000000001, 0.0], [1.6694601350724732, 0.15600000000000003, 0.0], [1.6694601350724732, 0.20800000000000002, 0.0], [1.6694601350724732, 0.26, 0.0], [1.6694601350724732, 0.31200000000000006, 0.0], [1.6694601350724732, 0.36400000000000005, 0.0], [1.6694601350724732, 0.41600000000000004, 0.0], [1.6694601350724732, 0.468, 0.0], [1.6694601350724732, 0.52, 0.0]], [[1.7147599769840711, -0.0, 0.0], [1.7147599769840711, 0.052000000000000005, 0.0], [1.7147599769840711, 0.10400000000000001, 0.0], [1.7147599769840711, 0.15600000000000003, 0.0], [1.7147599769840711, 0.20800000000000002, 0.0], [1.7147599769840711, 0.26, 0.0], [1.7147599769840711, 0.31200000000000006, 0.0], [1.7147599769840711, 0.36400000000000005, 0.0], [1.7147599769840711, 0.41600000000000004, 0.0], [1.7147599769840711, 0.468, 0.0], [1.7147599769840711, 0.52, 0.0]], [[1.76497523857, -0.0, 0.0], [1.76497523857, 0.052000000000000005, 0.0], [1.76497523857, 0.10400000000000001, 0.0], [1.76497523857, 0.15600000000000003, 0.0], [1.76497523857, 0.20800000000000002, 0.0], [1.76497523857, 0.26, 0.0], [1.76497523857, 0.31200000000000006, 0.0], [1.76497523857, 0.36400000000000005, 0.0], [1.76497523857, 0.41600000000000004, 0.0], [1.76497523857, 0.468, 0.0], [1.76497523857, 0.52, 0.0]], [[1.815190500155929, -0.0, 0.0], [1.815190500155929, 0.052000000000000005, 0.0], [1.815190500155929, 0.10400000000000001, 0.0], [1.815190500155929, 0.15600000000000003, 0.0], [1.815190500155929, 0.20800000000000002, 0.0], [1.815190500155929, 0.26, 0.0], [1.815190500155929, 0.31200000000000006, 0.0], [1.815190500155929, 0.36400000000000005, 0.0], [1.815190500155929, 0.41600000000000004, 0.0], [1.815190500155929, 0.468, 0.0], [1.815190500155929, 0.52, 0.0]], [[1.860490342067527, -0.0, 0.0], [1.860490342067527, 0.052000000000000005, 0.0], [1.860490342067527, 0.10400000000000001, 0.0], [1.860490342067527, 0.15600000000000003, 0.0], [1.860490342067527, 0.20800000000000002, 0.0], [1.860490342067527, 0.26, 0.0], [1.860490342067527, 0.31200000000000006, 0.0], [1.860490342067527, 0.36400000000000005, 0.0], [1.860490342067527, 0.41600000000000004, 0.0], [1.860490342067527, 0.468, 0.0], [1.860490342067527, 0.52, 0.0]], [[1.896440500155929, -0.0, 0.0], [1.896440500155929, 0.052000000000000005, 0.0], [1.896440500155929, 0.10400000000000001, 0.0], [1.896440500155929, 0.15600000000000003, 0.0], [1.896440500155929, 0.20800000000000002, 0.0], [1.896440500155929, 0.26, 0.0], [1.896440500155929, 0.31200000000000006, 0.0], [1.896440500155929, 0.36400000000000005, 0.0], [1.896440500155929, 0.41600000000000004, 0.0], [1.896440500155929, 0.468, 0.0], [1.896440500155929, 0.52, 0.0]], [[1.9195219224679625, -0.0, 0.0], [1.9195219224679625, 0.052000000000000005, 0.0], [1.9195219224679625, 0.10400000000000001, 0.0], [1.9195219224679625, 0.15600000000000003, 0.0], [1.9195219224679625, 0.20800000000000002, 0.0], [1.9195219224679625, 0.26, 0.0], [1.9195219224679625, 0.31200000000000006, 0.0], [1.9195219224679625, 0.36400000000000005, 0.0], [1.9195219224679625, 0.41600000000000004, 0.0], [1.9195219224679625, 0.468, 0.0], [1.9195219224679625, 0.52, 0.0]], [[1.92747523857, -0.0, 0.0], [1.92747523857, 0.052000000000000005, 0.0], [1.92747523857, 0.10400000000000001, 0.0], [1.92747523857, 0.15600000000000003, 0.0], [1.92747523857, 0.20800000000000002, 0.0], [1.92747523857, 0.26, 0.0], [1.92747523857, 0.31200000000000006, 0.0], [1.92747523857, 0.36400000000000005, 0.0], [1.92747523857, 0.41600000000000004, 0.0], [1.92747523857, 0.468, 0.0], [1.92747523857, 0.52, 0.0]]], dtype=np.float64), + "naca": "0012", + "control_assignments": { + "Elevator": { + "assignment": np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=np.int64), + "xhinged": 0.0, + "vhinged": np.array([0, 1, 0], dtype=np.int64), + "gaind": -1.0, + "refld": 1.0, + }, + }, + }, + "Vertical Tail": { + "component": np.int32(2), + "claf": 1.1078, + "scale": np.array([1.0, 1.0, 1.0], dtype=np.float64), + "translate": np.array([0.0, 0.0, 0.0], dtype=np.float64), + "mesh": np.array([[[1.60247523857, 0.0, 0.296], [1.60247523857, 0.0, 0.2368], [1.60247523857, 0.0, 0.1776], [1.60247523857, 0.0, 0.1184], [1.60247523857, 0.0, 0.0592], [1.60247523857, 0.0, -0.0]], [[1.6104285546720376, 0.0, 0.296], [1.6104285546720376, 0.0, 0.2368], [1.6104285546720376, 0.0, 0.1776], [1.6104285546720376, 0.0, 0.1184], [1.6104285546720376, 0.0, 0.0592], [1.6104285546720376, 0.0, -0.0]], [[1.633509976984071, 0.0, 0.296], [1.633509976984071, 0.0, 0.2368], [1.633509976984071, 0.0, 0.1776], [1.633509976984071, 0.0, 0.1184], [1.633509976984071, 0.0, 0.0592], [1.633509976984071, 0.0, -0.0]], [[1.6694601350724732, 0.0, 0.296], [1.6694601350724732, 0.0, 0.2368], [1.6694601350724732, 0.0, 0.1776], [1.6694601350724732, 0.0, 0.1184], [1.6694601350724732, 0.0, 0.0592], [1.6694601350724732, 0.0, -0.0]], [[1.7147599769840711, 0.0, 0.296], [1.7147599769840711, 0.0, 0.2368], [1.7147599769840711, 0.0, 0.1776], [1.7147599769840711, 0.0, 0.1184], [1.7147599769840711, 0.0, 0.0592], [1.7147599769840711, 0.0, -0.0]], [[1.76497523857, 0.0, 0.296], [1.76497523857, 0.0, 0.2368], [1.76497523857, 0.0, 0.1776], [1.76497523857, 0.0, 0.1184], [1.76497523857, 0.0, 0.0592], [1.76497523857, 0.0, -0.0]], [[1.815190500155929, 0.0, 0.296], [1.815190500155929, 0.0, 0.2368], [1.815190500155929, 0.0, 0.1776], [1.815190500155929, 0.0, 0.1184], [1.815190500155929, 0.0, 0.0592], [1.815190500155929, 0.0, -0.0]], [[1.860490342067527, 0.0, 0.296], [1.860490342067527, 0.0, 0.2368], [1.860490342067527, 0.0, 0.1776], [1.860490342067527, 0.0, 0.1184], [1.860490342067527, 0.0, 0.0592], [1.860490342067527, 0.0, -0.0]], [[1.896440500155929, 0.0, 0.296], [1.896440500155929, 0.0, 0.2368], [1.896440500155929, 0.0, 0.1776], [1.896440500155929, 0.0, 0.1184], [1.896440500155929, 0.0, 0.0592], [1.896440500155929, 0.0, -0.0]], [[1.9195219224679625, 0.0, 0.296], [1.9195219224679625, 0.0, 0.2368], [1.9195219224679625, 0.0, 0.1776], [1.9195219224679625, 0.0, 0.1184], [1.9195219224679625, 0.0, 0.0592], [1.9195219224679625, 0.0, -0.0]], [[1.92747523857, 0.0, 0.296], [1.92747523857, 0.0, 0.2368], [1.92747523857, 0.0, 0.1776], [1.92747523857, 0.0, 0.1184], [1.92747523857, 0.0, 0.0592], [1.92747523857, 0.0, -0.0]]], dtype=np.float64), + "naca": "0012", + "control_assignments": { + "Rudder": { + "assignment": np.array([0, 1, 2, 3, 4, 5], dtype=np.int64), + "xhinged": 0.5, + "vhinged": np.array([0, 0, 1], dtype=np.int64), + "gaind": 1.0, + "refld": -1.0, + }, + }, + }, + }, + "dname": ['Elevator', 'Rudder'], +} diff --git a/geom_files/make_ffd_wing.py b/geom_files/make_ffd_wing.py new file mode 100644 index 0000000..c347e97 --- /dev/null +++ b/geom_files/make_ffd_wing.py @@ -0,0 +1,24 @@ +# ============================================================================= +# Standard Python Modules +# ============================================================================= +import os + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np + +# ============================================================================= +# Local modules +# ============================================================================= +from optvl.utils.ffd_utils import write_FFD_file +from wing_mesh import mesh + +filename = "wing_test_ffd" + +def create_ffd(): + write_FFD_file(mesh,12,6, filename) + + +if __name__ == "__main__": + create_ffd() \ No newline at end of file diff --git a/geom_files/make_pkls/make_aircraft_L1.py b/geom_files/make_pkls/make_aircraft_L1.py new file mode 100644 index 0000000..3c0e145 --- /dev/null +++ b/geom_files/make_pkls/make_aircraft_L1.py @@ -0,0 +1,130 @@ +# ============================================================================= +# Extension modules +# ============================================================================= +from optvl import OVLSolver + +# ============================================================================= +# Standard Python Modules +# ============================================================================= +import os + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np +import pickle +from openaerostruct.meshing.mesh_generator import generate_mesh +from openaerostruct.geometry.utils import taper + + +mesh_dict = { + "num_y": 5, + "num_x": 3, + "wing_type": "rect", + "symmetry": True, + "span": 2.0, + "root_chord": 0.4, + "span_cos_spacing": 0.0, + "chord_cos_spacing": 1.0, + } +mesh_wing = generate_mesh(mesh_dict) + +taper(mesh_wing,0.35,True,ref_axis_pos=0.0) + + +mesh_dict = { + "num_y": 3, + "num_x": 11, + "wing_type": "rect", + "symmetry": True, + "span": 1.0, + "root_chord": 0.325, + "span_cos_spacing": 0.0, + "chord_cos_spacing": 1.0, + } +mesh_tail = generate_mesh(mesh_dict) + +mesh_wing[:,:,1] = -mesh_wing[:,:,1] +mesh_wing = mesh_wing[:,::-1,:] + +mesh_tail[:,:,1] = -mesh_tail[:,:,1] +mesh_tail = mesh_tail[:,::-1,:] + +mesh_tail[:,:,0] += 1.5 + + +surf = { + "Wing": { + # General + # "component": np.int32(0), # logical surface component index (for grouping interacting surfaces, see AVL manual) + # "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh_wing), # (nx,ny,3) numpy array containing mesh coordinates + + }, + "Horizontal Tail": { + # General + # "component": np.int32(1), # logical surface component index (for grouping interacting surfaces, see AVL manual) + # "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh_tail), # (nx,ny,3) numpy array containing mesh coordinates + # Control Surface Specification + "control_assignments": { + "Elevator" : {"assignment":np.arange(0,mesh_tail.shape[1]), + "xhinged": 0.0, # x/c location of hinge + "vhinged": np.array([0,1,0]), # vector giving hinge axis about which surface rotates + "gaind": -1.0, # control surface gain + "refld": 1.0 # control surface reflection, sign of deflection for duplicated surface + } + }, + + } +} + + +input_dict = { + "title": "MACH MDAO AVL", + "mach": np.float64(0.0), + "iysym": np.int32(0), + "izsym": np.int32(0), + "zsym": np.float64(0.0), + "Sref": np.float64(1.13), + "Cref": np.float64(0.4), + "Bref": np.float64(1.0), + "XYZref": np.array([0.0, 0.0, 0.0],dtype=np.float64), + "CDp": np.float64(0.0), + "surfaces": surf, + # "bodies": fuselage, + # Global Control and DV info + "dname": ["Elevator"], # Name of control input for each corresonding index +} + +# For verification +# base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder +# geom_dir = os.path.join(base_dir, '..', 'geom_files') +# rect_file = os.path.join(geom_dir, 'aircraft_L1.avl') + + +# solver = OVLSolver(input_dict=input_dict,debug=True) +# solver = OVLSolver(geo_file=rect_file,debug=True) + +# solver.set_variable("alpha", 25.0) +# solver.set_variable("beta", 5.0) +# solver.execute_run() + +with open("aircraft_L1.pkl", 'wb') as f: + pickle.dump(input_dict, f) diff --git a/geom_files/make_pkls/make_rect_with_body.py b/geom_files/make_pkls/make_rect_with_body.py new file mode 100644 index 0000000..11bc622 --- /dev/null +++ b/geom_files/make_pkls/make_rect_with_body.py @@ -0,0 +1,103 @@ +# ============================================================================= +# Extension modules +# ============================================================================= +from optvl import OVLSolver + +# ============================================================================= +# Standard Python Modules +# ============================================================================= +import os + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np +import pickle + + +mesh = np.zeros((2,2,3)) +mesh[1,:,0] = 1.0 +mesh[:,1,1] = 1.0 +mesh[:,:,2] = 1.0 + + +surf = { + "Wing": { + # General + # "component": np.int32(1), # logical surface component index (for grouping interacting surfaces, see AVL manual) + # "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh), # (nx,ny,3) numpy array containing mesh coordinates + # Control Surface Specification + "control_assignments": { + "Elevator" : {"assignment":np.arange(0,mesh.shape[1]), + "xhinged": 0.5, # x/c location of hinge + "vhinged": np.array([0,1,0]), # vector giving hinge axis about which surface rotates + "gaind": -1.0, # control surface gain + "refld": 1.0 # control surface reflection, sign of deflection for duplicated surface + } + }, + + } +} + +fuselage = {"Fuse pod": { + # General + # 'yduplicate': np.float64(0), # body is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0] + ), # scaling factors applied to all x,y,z coordinates (chords areal so scaled by Xscale) + "translate": np.array([0.0, 0.0, 0.0]), # offset added on to all X,Y,Z values in this surface + # Discretization + "nvb": np.int32(4), # number of source-line nodes + "bspace": np.float64(2.0), # lengthwise node spacing parameter + "bfile": "../geom_files/fuseSimple.dat", # body oml file name +} +} + +input_dict = { + "title": "MACH MDAO AVL", + "mach": np.float64(0.0), + "iysym": np.int32(0), + "izsym": np.int32(0), + "zsym": np.float64(0.0), + "Sref": np.float64(1.123), + "Cref": np.float64(0.25), + "Bref": np.float64(6.01), + "XYZref": np.array([0.0, 0, 0],dtype=np.float64), + "CDp": np.float64(0.0), + "surfaces": surf, + "bodies": fuselage, + # Global Control and DV info + "dname": ["Elevator"], # Name of control input for each corresonding index +} + +# For verification +# base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder +# geom_dir = os.path.join(base_dir, '..', 'geom_files') +# rect_file = os.path.join(geom_dir, 'rect_with_body.avl') + + +# solver = OVLSolver(input_dict=input_dict,debug=True) +# solver = OVLSolver(geo_file=rect_file,debug=True) + +# solver.set_variable("alpha", 25.0) +# solver.set_variable("beta", 5.0) +# solver.execute_run() + +# solver.plot_geom() + +with open("rect_with_body.pkl", 'wb') as f: + pickle.dump(input_dict, f) + + + + diff --git a/geom_files/make_py/make_aircraft.py b/geom_files/make_py/make_aircraft.py new file mode 100644 index 0000000..0ab7eff --- /dev/null +++ b/geom_files/make_py/make_aircraft.py @@ -0,0 +1,212 @@ +# ============================================================================= +# Standard Python Modules +# ============================================================================= +import os +import copy + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np +from openaerostruct.meshing.mesh_generator import generate_mesh +from openaerostruct.meshing.section_mesh_generator import generate_mesh as generate_section_mesh +from openaerostruct.geometry.utils import shear_z + +# ============================================================================= +# Local modules +# ============================================================================= +from write_dict_to_py import write_dict_to_py +from optvl import OVLSolver + + +mesh_dict_wing = { + "num_sections": 4, + "sec_name": ["sec0", "sec1", "sec2", "sec3"], + "symmetry": True, + "ref_axis_pos": 0.0, + "root_section": 0, + # Geometry Parameters + "taper": np.array([0.6666, 0.75, 0.8888, 1.0]), # Wing taper for each section + "span": np.array([0.25, 0.45, 0.55, 0.75]), # Wing span for each section + "sweep": np.array([0.0, 0.0, 0.0, 0,0]), # Wing sweep for each section + "root_chord": 0.45, # Wing root chord for each section + # Mesh Parameters + "nx": 8, # Number of chordwise panels. Same for all sections + "ny": np.array([6, 6, 6, 6]), # Number of spanwise panels for each section +} + + +_, sections_wing = generate_section_mesh(mesh_dict_wing) + + +mesh_dict_h_tail = { + "num_y": 21, + "num_x": 11, + "wing_type": "rect", + "symmetry": True, + "span": 0.52*2, + "root_chord": 0.325, + "span_cos_spacing": 0.0, + "chord_cos_spacing": 1.0, + } +mesh_h_tail = generate_mesh(mesh_dict_h_tail) + +mesh_dict_v_tail = { + "num_y": 11, + "num_x": 11, + "wing_type": "rect", + "symmetry": True, + "span": 0.296*2, + "root_chord": 0.325, + "span_cos_spacing": 0.0, + "chord_cos_spacing": 1.0, + } +mesh_v_tail = generate_mesh(mesh_dict_v_tail) + + +shear_z(sections_wing[0],np.linspace(0.1,0.0,6)) + +for i in range(4): + if i == 0: + mesh_wing = sections_wing[i][:,:-1,:] + elif i == 3: + mesh_wing = np.hstack([mesh_wing,sections_wing[i]]) + else: + mesh_wing = np.hstack([mesh_wing,sections_wing[i][:,:-1,:]]) +# mesh_wing = np.hstack(sections_wing) + +mesh_wing[:,:,1] = -mesh_wing[:,:,1] +mesh_wing = mesh_wing[:,::-1,:] + +mesh_h_tail[:,:,1] = -mesh_h_tail[:,:,1] +mesh_h_tail = mesh_h_tail[:,::-1,:] + +mesh_h_tail[:,:,0] += 1.60247523857 + +mesh_v_tail[:,:,2] = -mesh_v_tail[:,:,1] +mesh_v_tail[:,:,1] = 0.0 + +mesh_v_tail[:,:,0] += 1.60247523857 + + +surf = { + "Wing": { + # General + # "component": np.int32(0), # logical surface component index (for grouping interacting surfaces, see AVL manual) + "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh_wing), # (nx,ny,3) numpy array containing mesh coordinates + # Airfoils + # "afiles": np.array(["../geom_files/airfoils/A_1.dat","../geom_files/airfoils/A_2.dat","../geom_files/airfoils/A_3.dat","../geom_files/airfoils/A_4.dat","../geom_files/airfoils/A_5.dat"]), + # "afiles": np.concatenate([np.repeat("../airfoils/A_1.dat",5),np.repeat("../airfoils/A_2.dat",5),np.repeat("../airfoils/A_3.dat",5),np.repeat("../airfoils/A_4.dat",5),["../airfoils/A_5.dat"]]), + "afiles": "../geom_files/airfoils/ag40d.dat", + + }, + "Horizontal Tail": { + # General + # "component": np.int32(1), # logical surface component index (for grouping interacting surfaces, see AVL manual) + "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + "claf": 1.1078, + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh_h_tail), # (nx,ny,3) numpy array containing mesh coordinates + # Airfoils + "naca": "0012", + # Control Surface Specification + "control_assignments": { + "Elevator" : {"assignment":np.arange(0,mesh_h_tail.shape[1]), + "xhinged": 0.0, # x/c location of hinge + "vhinged": np.array([0,1,0]), # vector giving hinge axis about which surface rotates + "gaind": -1.0, # control surface gain + "refld": 1.0 # control surface reflection, sign of deflection for duplicated surface + } + }, + + }, + "Vertical Tail": { + # General + "component": np.int32(2), # logical surface component index (for grouping interacting surfaces, see AVL manual) + # "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + "claf": 1.1078, + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh_v_tail), # (nx,ny,3) numpy array containing mesh coordinates + # Airfoils + "naca": "0012", + # Control Surface Specification + "control_assignments": { + "Rudder" : {"assignment":np.arange(0,mesh_v_tail.shape[1]), + "xhinged": 0.5, # x/c location of hinge + "vhinged": np.array([0,0,1]), # vector giving hinge axis about which surface rotates + "gaind": 1.0, # control surface gain + "refld": -1.0 # control surface reflection, sign of deflection for duplicated surface + } + }, + + } +} + + +input_dict = { + "title": "AIRCRAFT 1", + "mach": np.float64(0.1), + "iysym": np.int32(0), + "izsym": np.int32(0), + "zsym": np.float64(0.0), + "Sref": np.float64(1.5825), + "Cref": np.float64(0.35), + "Bref": np.float64(4.0), + "XYZref": np.array([0.084, 0.0, 0.0],dtype=np.float64), + "CDp": np.float64(0.011600), + "surfaces": surf, + # "bodies": fuselage, + # Global Control and DV info + "dname": ["Elevator", "Rudder"], # Name of control input for each corresonding index +} + +# For verification +# base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder +# geom_dir = os.path.join(base_dir, '..', 'geom_files') +# rect_file = os.path.join(geom_dir, 'aircraft.avl') + + +# solver = OVLSolver(input_dict=input_dict,debug=True) +# solver = OVLSolver(geo_file=rect_file,debug=True) + +# solver.set_variable("alpha", 25.0) +# solver.set_variable("beta", 5.0) +# solver.execute_run() + +# solver.plot_geom() +# solver.plot_cp() + +# Write as a plain text Python file instead of pickle +output_file = "../aircraft.py" + +write_dict_to_py( + output_file, + input_dict, + var_name='input_dict', + description='Generated by make_aircraft.py - Aircraft configuration with wing and tail' +) + +print(f"Generated {output_file}") diff --git a/geom_files/rect_out.avl b/geom_files/rect_out.avl new file mode 100644 index 0000000..420353e --- /dev/null +++ b/geom_files/rect_out.avl @@ -0,0 +1,41 @@ +# generated using OptVL v2.4.0.dev0 +#=============================================================================== +#------------------------------------ Header ----------------------------------- +#=============================================================================== +MACH MDAO AVL +#Mach +0.12341234 +#IYsym IZsym Zsym +0 0 0.0 +#Sref Cref Bref +1.0 2.0 3.0 +#Xref Yref Zref +4.0 5.0 6.0 +#CD0 +0.0 +#=============================================================================== +#------------------------------------- Wing ------------------------------------ +#=============================================================================== +SURFACE +Wing +#Nchordwise Cspace [Nspanwise Sspace] +1 1.0 1 -2.0 +SCALE +1.0 1.0 1.0 +TRANSLATE +0.0 0.0 0.0 +ANGLE +0.0 +#--------------------------------------- +SECTION +#Xle Yle Zle | Chord Ainc Nspan Sspace + 0.000000 0.000000 0.000000 1.000000 0.000000 + CONTROL +#surface gain xhinge hvec SgnDup + Elevator -1.0 0.5 0.000000 1.000000 0.000000 1.0 +SECTION +#Xle Yle Zle | Chord Ainc Nspan Sspace + 0.000000 1.000000 0.000000 1.000000 0.000000 + CONTROL +#surface gain xhinge hvec SgnDup + Elevator -1.0 0.5 0.000000 1.000000 0.000000 1.0 diff --git a/geom_files/wing_test_ffd.xyz b/geom_files/wing_test_ffd.xyz new file mode 100644 index 0000000..59b4d89 --- /dev/null +++ b/geom_files/wing_test_ffd.xyz @@ -0,0 +1,50 @@ +1 +12 2 6 +-0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 + -0.05 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 + 0.54545455 0.63636364 0.72727273 0.81818182 0.90909091 1.05 +-5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 + -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 -5.05 + -4. -4. -4. -4. -4. -4. -4. -4. -4. -4. -4. -4. + -4. -4. -4. -4. -4. -4. -4. -4. -4. -4. -4. -4. + -3. -3. -3. -3. -3. -3. -3. -3. -3. -3. -3. -3. + -3. -3. -3. -3. -3. -3. -3. -3. -3. -3. -3. -3. + -2. -2. -2. -2. -2. -2. -2. -2. -2. -2. -2. -2. + -2. -2. -2. -2. -2. -2. -2. -2. -2. -2. -2. -2. + -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. + -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. + 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 + 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 +-0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 + 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 + -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 + 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 + -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 + 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 + -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 + 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 + -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 + 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 + -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 + 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 diff --git a/optvl/om_wrapper.py b/optvl/om_wrapper.py index 9a43395..fdcd9ac 100644 --- a/optvl/om_wrapper.py +++ b/optvl/om_wrapper.py @@ -18,13 +18,15 @@ class OVLGroup(om.Group): output_dir: the output directory for the files generated input_param_vals: flag to turn on the flght parameters (Mach, Velocity, etc.) as inputs input_ref_val: flag to turn on the geometric reference values (Sref, Cref, Bref) as inputs + input_mesh_dim: flag that sets the mesh input shape dimention count output_stability_derivs: flag to turn on the output of stability derivatives output_body_axis_derivs: flag to turn on the output of body axis derivatives output_con_surf_derivs: flag to turn on the output of control surface deflections """ def initialize(self): - self.options.declare("geom_file", types=str) + self.options.declare("geom_file", types=str, default=None) + self.options.declare("input_dict", types=dict, default=None) self.options.declare("mass_file", default=None) self.options.declare("write_grid", types=bool, default=False) self.options.declare("write_grid_sol_time", types=bool, default=False) @@ -33,6 +35,7 @@ def initialize(self): self.options.declare("input_param_vals", types=bool, default=False) self.options.declare("input_ref_vals", types=bool, default=False) self.options.declare("input_airfoil_geom", types=bool, default=False) + self.options.declare("input_mesh_dim", types=int, default=1) self.options.declare("output_stability_derivs", types=bool, default=False) self.options.declare("output_body_axis_derivs", types=bool, default=False) @@ -40,17 +43,19 @@ def initialize(self): def setup(self): geom_file = self.options["geom_file"] + input_dict = self.options["input_dict"] mass_file = self.options["mass_file"] input_param_vals = self.options["input_param_vals"] input_ref_vals = self.options["input_ref_vals"] input_airfoil_geom = self.options["input_airfoil_geom"] + input_mesh_dim = self.options["input_mesh_dim"] output_stability_derivs = self.options["output_stability_derivs"] output_body_axis_derivs = self.options["output_body_axis_derivs"] output_con_surf_derivs = self.options["output_con_surf_derivs"] - self.ovl = OVLSolver(geo_file=geom_file, mass_file=mass_file, debug=False) + self.ovl = OVLSolver(geo_file=geom_file, input_dict=input_dict, mass_file=mass_file, debug=False) self.add_subsystem( "solver", @@ -59,6 +64,7 @@ def setup(self): input_param_vals=input_param_vals, input_ref_vals=input_ref_vals, input_airfoil_geom=input_airfoil_geom, + input_mesh_dim=input_mesh_dim, ), promotes=["*"], ) @@ -72,6 +78,7 @@ def setup(self): output_stability_derivs=output_stability_derivs, output_body_axis_derivs=output_body_axis_derivs, output_con_surf_derivs=output_con_surf_derivs, + input_mesh_dim=input_mesh_dim, ), promotes=["*"], ) @@ -91,6 +98,7 @@ def setup(self): AIRFOIL_GEOM_VARS = ["xasec", "casec", "tasec", "xuasec", "xlasec", "zuasec", "zlasec"] +AVL_GEOM_VARS = ["xles","yles","zles","chords"] # helper functions used by the AVL components @@ -107,25 +115,49 @@ def add_ovl_geom_vars(self, ovl, add_as="inputs", include_airfoil_geom=False): surf_data = ovl.get_surface_params() for surf in surf_data: + idx_surf = ovl.get_surface_index(surf) for key in surf_data[surf]: if key in AIRFOIL_GEOM_VARS: if not include_airfoil_geom: continue + if key in AVL_GEOM_VARS: + if ovl.avl.SURF_MESH_L.LSURFMSH[idx_surf]: + continue geom_key = f"{surf}:{key}" if add_as == "inputs": self.add_input(geom_key, val=surf_data[surf][key], tags="geom") elif add_as == "outputs": self.add_output(geom_key, val=surf_data[surf][key], tags="geom") -def add_ovl_mesh_out_as_output(self, ovl): - surf_data = ovl.get_surface_params() +def add_ovl_mesh_vars(self, ovl, add_as="inputs", input_mesh_dim=1): + # add the geometric parameters as inputs + surfs = ovl.unique_surface_names - meshes,_ = ovl.get_cp_data() + for surf in surfs: + idx_surf = ovl.get_surface_index(surf) + if ovl.avl.SURF_MESH_L.LSURFMSH[idx_surf]: + mesh_key = f"{surf}:mesh" + mesh = ovl.get_mesh(idx_surf) - for surf in surf_data: - idx_surf = ovl.surface_names.index(surf) - out_name = f"{surf}:mesh" - self.add_output(out_name, val=meshes[idx_surf], tags="geom_mesh") + if input_mesh_dim == 1: + mesh = mesh.transpose((1,0,2)).flatten() + elif input_mesh_dim == 2: + mesh = mesh.transpose((1,0,2)).reshape((mesh.shape[0]*mesh.shape[1],3)) + + if add_as == "inputs": + self.add_input(mesh_key, val=mesh, tags="geom_mesh") + elif add_as == "outputs": + self.add_output(mesh_key, val=mesh, tags="geom_mesh") + +# def add_ovl_mesh_out_as_output(self, ovl): +# surf_data = ovl.get_surface_params() + +# meshes,_ = ovl.get_cp_data() + +# for surf in surf_data: +# idx_surf = ovl.surface_names.index(surf) +# out_name = f"{surf}:mesh" +# self.add_output(out_name, val=meshes[idx_surf], tags="geom_mesh") def add_ovl_conditions_as_inputs(sys, ovl): # TODO: add all the condition constraints @@ -151,11 +183,14 @@ def add_ovl_refs_as_inputs(sys, ovl): def om_input_to_surf_dict(sys, inputs): - geom_inputs = sys.list_inputs(tags="geom", val=False, out_stream=None) + geom_inputs = sys.list_inputs(tags=["geom"], val=False, out_stream=None) + geom_mesh_inputs = sys.list_inputs(tags=["geom_mesh"], val=False, out_stream=None) # convert to a list witout tuples geom_inputs = [x[0] for x in geom_inputs] + geom_mesh_inputs = [x[0] for x in geom_mesh_inputs] surf_data = {} + surf_mesh_data = {} for input_var in inputs: if input_var in geom_inputs: # split the input name into surface name and parameter name @@ -169,7 +204,15 @@ def om_input_to_surf_dict(sys, inputs): surf_data[surf][param] = inputs[input_var][0] else: surf_data[surf][param] = inputs[input_var] - return surf_data + elif input_var in geom_mesh_inputs: + # split the input name into surface name and parameter name + surf, param = input_var.split(":") + + # update the corresponding parameter in the surface data + if surf not in surf_mesh_data: + surf_mesh_data[surf] = {} + surf_mesh_data[surf][param] = inputs[input_var] + return surf_data, surf_mesh_data def om_surf_dict_to_input(surf_dict): @@ -216,6 +259,7 @@ def initialize(self): self.options.declare("input_param_vals", types=bool, default=False) self.options.declare("input_ref_vals", types=bool, default=False) self.options.declare("input_airfoil_geom", types=bool, default=False) + self.options.declare("input_mesh_dim", types=int, default=1) def setup(self): self.ovl = self.options["ovl"] @@ -226,6 +270,7 @@ def setup(self): self.num_states = self.ovl.get_mesh_size() self.num_cs = self.ovl.get_num_control_surfs() self.num_vel = self.ovl.NUMAX + self.input_mesh_dim = self.options["input_mesh_dim"] self.add_output("gamma", val=np.zeros(self.num_states)) self.add_output("gamma_d", val=np.zeros((self.num_cs, self.num_states))) @@ -241,16 +286,28 @@ def setup(self): self.control_names = add_ovl_controls_as_inputs(self, self.ovl) add_ovl_geom_vars(self, self.ovl, add_as="inputs", include_airfoil_geom=input_airfoil_geom) + add_ovl_mesh_vars(self,self.ovl,add_as="inputs",input_mesh_dim=self.input_mesh_dim) self.res_slice = (slice(0, self.num_states),) self.res_d_slice = (slice(0, self.num_cs), slice(0, self.num_states)) self.res_u_slice = (slice(0, self.num_vel), slice(0, self.num_states)) def apply_nonlinear(self, inputs, outputs, residuals): + # Set case parameters om_set_avl_inputs(self, inputs) - surf_data = om_input_to_surf_dict(self, inputs) + surf_data, surf_mesh_data = om_input_to_surf_dict(self, inputs) self.ovl.set_surface_params(surf_data) + # set the meshes + for surf in surf_mesh_data: + if "mesh" in surf_mesh_data[surf].keys(): + idx_surf = self.ovl.get_surface_index(surf) + # Handle all possible mesh shapes + if self.input_mesh_dim < 3: + nx = self.ovl.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + surf_mesh_data[surf]["mesh"] = surf_mesh_data[surf]["mesh"].reshape((ny,nx,3)).transpose((1,0,2)) + self.ovl.set_mesh(idx_surf,surf_mesh_data[surf]["mesh"]) gam_arr = outputs["gamma"] gam_d_arr = outputs["gamma_d"] @@ -276,10 +333,21 @@ def apply_nonlinear(self, inputs, outputs, residuals): def solve_nonlinear(self, inputs, outputs): start_time = time.time() + # set case params om_set_avl_inputs(self, inputs) # update the surface parameters - surf_data = om_input_to_surf_dict(self, inputs) + surf_data, surf_mesh_data = om_input_to_surf_dict(self, inputs) + # set the meshes + for surf in surf_mesh_data: + if "mesh" in surf_mesh_data[surf].keys(): + idx_surf = self.ovl.get_surface_index(surf) + # Handle all possible mesh shapes + if self.input_mesh_dim < 3: + nx = self.ovl.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + surf_mesh_data[surf]["mesh"] = surf_mesh_data[surf]["mesh"].reshape((ny,nx,3)).transpose((1,0,2)) + self.ovl.set_mesh(idx_surf,surf_mesh_data[surf]["mesh"]) self.ovl.set_surface_params(surf_data) # def_dict = self.ovl.get_control_deflections() @@ -313,7 +381,19 @@ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode): if con_key in d_inputs: con_seeds[con_key] = d_inputs[con_key] - geom_seeds = om_input_to_surf_dict(self, d_inputs) + geom_seeds, mesh_seeds = om_input_to_surf_dict(self, d_inputs) + + # Reshape mesh seeds + if self.input_mesh_dim != 2: + for surface in mesh_seeds: + if "mesh" in mesh_seeds[surface].keys(): + idx_surf = self.ovl.get_surface_index(surf_name=surface) + if self.input_mesh_dim == 1: + mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seeds[surface]["mesh"].flatten()) + elif self.input_mesh_dim == 3: + nx = self.ovl.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seeds[surface]["mesh"].reshape((ny,nx,3)).transpose((1,0,2))) param_seeds = {} for param in self.ovl.param_idx_dict: @@ -325,8 +405,8 @@ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode): if ref in d_inputs: ref_seeds[ref] = d_inputs[ref] - _, res_seeds, _, _, res_d_seeds, res_u_seeds = self.ovl._execute_jac_vec_prod_fwd( - con_seeds=con_seeds, geom_seeds=geom_seeds, param_seeds=param_seeds, ref_seeds=ref_seeds + _, res_seeds, _, _, _, res_d_seeds, res_u_seeds = self.ovl._execute_jac_vec_prod_fwd( + con_seeds=con_seeds, geom_seeds=geom_seeds, mesh_seeds=mesh_seeds, param_seeds=param_seeds, ref_seeds=ref_seeds ) d_residuals["gamma"] += res_seeds @@ -340,12 +420,23 @@ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode): res_d_seeds = d_residuals["gamma_d"] res_u_seeds = d_residuals["gamma_u"] - con_seeds, geom_seeds, mesh_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = ( + con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = ( self.ovl._execute_jac_vec_prod_rev( res_seeds=res_seeds, res_d_seeds=res_d_seeds, res_u_seeds=res_u_seeds ) ) + if self.input_mesh_dim != 2: + for surface in mesh_seeds: + if "mesh" in mesh_seeds[surface].keys(): + idx_surf = self.ovl.get_surface_index(surf_name=surface) + if self.input_mesh_dim == 1: + mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seeds[surface]["mesh"].flatten()) + elif self.input_mesh_dim == 3: + nx = self.ovl.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seeds[surface]["mesh"].reshape((ny,nx,3)).transpose((1,0,2))) + if "gamma" in d_outputs: d_outputs["gamma"] += gamma_seeds @@ -356,10 +447,13 @@ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode): d_outputs["gamma_u"] += gamma_u_seeds d_input_geom = om_surf_dict_to_input(geom_seeds) + d_input_mesh = om_surf_dict_to_input(mesh_seeds) for d_input in d_inputs: if d_input in d_input_geom: d_inputs[d_input] += d_input_geom[d_input] + elif d_input in d_input_mesh: + d_inputs[d_input] += d_input_mesh[d_input] elif d_input in ["alpha", "beta"]: d_inputs[d_input] += con_seeds[d_input] elif d_input in self.control_names: @@ -402,12 +496,14 @@ def initialize(self): self.options.declare("input_param_vals", types=bool, default=False) self.options.declare("input_ref_vals", types=bool, default=False) self.options.declare("input_airfoil_geom", types=bool, default=False) + self.options.declare("input_mesh_dim", types=int, default=1) def setup(self): self.ovl = self.options["ovl"] self.num_states = self.ovl.get_mesh_size() self.num_cs = self.ovl.get_num_control_surfs() self.num_vel = self.ovl.NUMAX + self.input_mesh_dim = self.options["input_mesh_dim"] input_param_vals = self.options["input_param_vals"] input_ref_vals = self.options["input_ref_vals"] input_airfoil_geom = self.options["input_airfoil_geom"] @@ -426,6 +522,7 @@ def setup(self): self.control_names = add_ovl_controls_as_inputs(self, self.ovl) add_ovl_geom_vars(self, self.ovl, add_as="inputs", include_airfoil_geom=input_airfoil_geom) + add_ovl_mesh_vars(self, self.ovl, add_as="inputs",input_mesh_dim=self.input_mesh_dim) # add the outputs for func_key in self.ovl.case_var_to_fort_var: @@ -455,10 +552,21 @@ def compute(self, inputs, outputs): # TODO: set_constraint does not correctly do derives yet start_time = time.time() + # set the case variables om_set_avl_inputs(self, inputs) # update the surface parameters - surf_data = om_input_to_surf_dict(self, inputs) + surf_data, surf_mesh_data = om_input_to_surf_dict(self, inputs) + # set the meshes + for surf in surf_mesh_data: + if "mesh" in surf_mesh_data[surf].keys(): + idx_surf = self.ovl.get_surface_index(surf) + # Handle all possible mesh shapes + if self.input_mesh_dim < 3: + nx = self.ovl.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + surf_mesh_data[surf]["mesh"] = surf_mesh_data[surf]["mesh"].reshape((ny,nx,3)).transpose((1,0,2)) + self.ovl.set_mesh(idx_surf,surf_mesh_data[surf]["mesh"]) self.ovl.set_surface_params(surf_data) gam_arr = inputs["gamma"] @@ -538,11 +646,24 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): if ref in d_inputs: ref_seeds[ref] = d_inputs[ref] - geom_seeds = self.om_input_to_surf_dict(self, d_inputs) + geom_seeds, mesh_seeds = self.om_input_to_surf_dict(self, d_inputs) + + # Reshape mesh seeds + if self.input_mesh_dim != 2: + for surface in mesh_seeds: + if "mesh" in mesh_seeds[surface].keys(): + idx_surf = self.ovl.get_surface_index(surf_name=surface) + if self.input_mesh_dim == 1: + mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seeds[surface]["mesh"].flatten()) + elif self.input_mesh_dim == 3: + nx = self.ovl.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seeds[surface]["mesh"].reshape((ny,nx,3)).transpose((1,0,2))) func_seeds, _, csd_seeds, stab_derivs_seeds, body_axis_seeds, _, _ = self.ovl._execute_jac_vec_prod_fwd( con_seeds=con_seeds, geom_seeds=geom_seeds, + mesh_seeds=mesh_seeds, gamma_seeds=gamma_seeds, gamma_d_seeds=gamma_d_seeds, gamma_u_seeds=gamma_u_seeds, @@ -610,7 +731,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): # print(var_name, body_axis_seeds[func_key]) print(f" running rev mode derivs for {func_key}") - con_seeds, geom_seeds, mesh_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = ( + con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = ( self.ovl._execute_jac_vec_prod_rev( func_seeds=func_seeds, consurf_derivs_seeds=csd_seeds, @@ -619,6 +740,17 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): ) ) + if self.input_mesh_dim != 2: + for surface in mesh_seeds: + if "mesh" in mesh_seeds[surface].keys(): + idx_surf = self.ovl.get_surface_index(surf_name=surface) + if self.input_mesh_dim == 1: + mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seeds[surface]["mesh"].flatten()) + elif self.input_mesh_dim == 3: + nx = self.ovl.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seeds[surface]["mesh"].reshape((ny,nx,3)).transpose((1,0,2))) + if "gamma" in d_inputs: d_inputs["gamma"] += gamma_seeds @@ -629,10 +761,13 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): d_inputs["gamma_u"] += gamma_u_seeds d_input_geom = om_surf_dict_to_input(geom_seeds) + d_input_mesh = om_surf_dict_to_input(mesh_seeds) for d_input in d_inputs: if d_input in d_input_geom: d_inputs[d_input] += d_input_geom[d_input] + elif d_input in d_input_mesh: + d_inputs[d_input] += d_input_mesh[d_input] elif d_input in ["alpha", "beta"] or d_input in self.control_names: d_inputs[d_input] += con_seeds[d_input] elif d_input in param_seeds: @@ -720,20 +855,23 @@ class OVLMeshReader(om.ExplicitComponent): """ def initialize(self): - self.options.declare("geom_file", types=str) + self.options.declare("geom_file", types=str, default=None) + self.options.declare("input_dict", types=dict, default=None) self.options.declare("mass_file", default=None) - self.options.declare("mesh_output",default=False) + self.options.declare("input_mesh_dim", types=int, default=1) def setup(self): geom_file = self.options["geom_file"] + input_dict = self.options["input_dict"] mass_file = self.options["mass_file"] - mesh_output = self.options["mesh_output"] + input_mesh_dim = self.options["input_mesh_dim"] - avl = OVLSolver(geo_file=geom_file, mass_file=mass_file, debug=False) + avl = OVLSolver(geo_file=geom_file, input_dict=input_dict, mass_file=mass_file, debug=False) add_ovl_geom_vars(self, avl, add_as="outputs", include_airfoil_geom=True) + add_ovl_mesh_vars(self, avl, add_as="outputs",input_mesh_dim=input_mesh_dim) - if mesh_output: - add_ovl_mesh_out_as_output(self,avl) + # if mesh_output: + # add_ovl_mesh_out_as_output(self,avl) class Differencer(om.ExplicitComponent): def setup(self): diff --git a/optvl/optvl_class.py b/optvl/optvl_class.py index 8ab2a53..508fdf4 100644 --- a/optvl/optvl_class.py +++ b/optvl/optvl_class.py @@ -403,10 +403,11 @@ def __init__( deriv_key = self._get_deriv_key(var, func) self.case_body_derivs_to_fort_var[deriv_key] = ["CASE_R", f"{func_to_prefix[func]}TOT_U_BA", idx_var] - # In the case where we used a file then we have to initialize these before _init_map_data so ad seeds work correctly + # In the case where we used a file then we have to initialize these before _init_map_data so ad seeds work correclty if not input_dict: self.mesh_idx_first = np.zeros(self.get_num_surfaces(),dtype=np.int32) self.y_offsets = np.zeros(self.get_num_surfaces(),dtype=np.float64) + self.point_sets = {} # the case parameters are stored in a 1d array, # these indices correspond to the position of each parameter in that arra @@ -415,6 +416,9 @@ def __init__( # set the default solver tolerance self.set_avl_fort_arr("CASE_R", "EXEC_TOL", 2e-5) + # init the DVGeo variable + self.DVGeo = None + if timing: print(f"AVL init took {time.time() - start_time} seconds") @@ -788,6 +792,9 @@ def check_type(key, avl_vars, given_val, cast_type=True): # a duplicated mesh as its normally only passed as a dummy argument into the SDUPL subroutine and not stored in the fortran layer. self.y_offsets = np.zeros(self.get_num_surfaces(),dtype=np.float64) + # Class dictionary to store the point set names if a DVGeo object is used + self.point_sets = {} + # Load surfaces if num_surfs > 0: surf_names = list(input_dict["surfaces"].keys()) @@ -1058,8 +1065,10 @@ def check_type(key, avl_vars, given_val, cast_type=True): surf_dict["flatten mesh"] = True self.set_mesh(idx_surf, surf_dict["mesh"],flatten=surf_dict["flatten mesh"],update_nvs=True,update_nvc=True) # set_mesh handles the Fortran indexing and ordering self.avl.makesurf_mesh(idx_surf + 1) #+1 for Fortran indexing + self.point_sets[idx_surf] = "optvl_%s_coords" % surf_name # store the pointset name for DVGeo later else: self.avl.makesurf(idx_surf + 1) # +1 to convert to 1 based indexing + self.point_sets[idx_surf] = None # No pointset available if no mesh if "yduplicate" in surf_dict.keys(): self.avl.sdupl(idx_surf + 1, surf_dict["yduplicate"], "YDUP") @@ -3149,6 +3158,82 @@ def _write_data(key_list: List[str], newline: bool = True): fid.write("#surface gain\n") fid.write(f" {design_var_names[idx_des_var - 1]} ") fid.write(f" {data['gaing'][idx_sec][idx_local_des_var]}\n") + + # region --- pyGeo API + def set_DVGeo(self, DVGeo, pointSetKwargs=None, customPointSetFamilies=None): + """ + Set the DVGeometry object that will manipulate 'geometry' in + this object. Note that does not **strictly** need a + DVGeometry object, but if optimization with geometric + changes is desired, then it is required. + + Parameters + ---------- + DVGeo : A DVGeometry object. + Object responsible for manipulating the geometry. + + pointSetKwargs : dict + Keyword arguments to be passed to the DVGeo addPointSet call. + Useful for DVGeometryMulti, specifying FFD projection tolerances, etc. + These arguments are used for all point sets added by this solver. + + customPointSetFamilies : dict of dicts + This argument is used to split up the surface points added to the DVGeo by the solver into potentially + multiple subsets. The keys of the dictionary will be used to determine what families should be + added to the dvgeo object as separate point sets. The values of each key is another dictionary, which can be empty. + If desired, the inner dictionaries can contain custom kwargs for the addPointSet call for each surface family, + specified by the keys of the top level dictionary. + The surface families need to be all part of the designSurfaceFamily. + Useful for DVGeometryMulti, specifying FFD projection tolerances, etc. + If this is provided together with pointSetKwargs, the regular pointSetKwargs + will be appended to each component's dictionary. If the same argument + is also provided in pointSetKwargs, the value specified in customPointSetFamilies + will be used. + + """ + + self.DVGeo = DVGeo + + # save the common kwargs dict. default is empty + if pointSetKwargs is None: + self.pointSetKwargs = {} + else: + self.pointSetKwargs = pointSetKwargs + + # save if we have customPointSetFamilies. this default is not mutable so we can just set it as is. + self.customPointSetFamilies = customPointSetFamilies + + def update_DVGeo(self): + """If DVGeo is present this function will embed all the meshes into it if needed and then perform an update + of the pointset and set the meshes back into AVL. + """ + # Check if we have an DVGeo object to deal with: + if self.DVGeo is not None: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a mesh, skip it + + # Embed the points if they haven't been already + if point_set_name not in self.DVGeo.points: + mesh = self.get_mesh(idx_surf) + coords0 = mesh.transpose((1,0,2)).reshape((mesh.shape[0]*mesh.shape[1],3)) + self.DVGeo.addPointSet(coords0, point_set_name, **self.pointSetKwargs) + # print(f"Embedeed point set {point_set_name}!") + + # Check if our point-set is up to date and update the mesh accordingly + if not self.DVGeo.pointSetUpToDate(point_set_name): + coords = self.DVGeo.update(point_set_name) + mesh_old = self.get_mesh(idx_surf) + mesh_new = copy.deepcopy(coords.reshape((mesh_old.shape[1],mesh_old.shape[0],3)).transpose((1,0,2))) + self.set_mesh(idx_surf,mesh_new) + + # Now update all of our surfaces in AVL + self.avl.update_surfaces() # region --- Utility functions def get_num_surfaces(self) -> int: @@ -3828,6 +3913,7 @@ def _execute_jac_vec_prod_fwd( con_seeds: Optional[Dict[str, float]] = None, geom_seeds: Optional[Dict[str, Dict[str, any]]] = None, mesh_seeds: Optional[Dict[str, Dict[str, np.ndarray]]] = None, + dvgeo_seeds: Optional[Dict[str, Dict[str, np.ndarray]]] = None, param_seeds: Optional[Dict[str, float]] = None, ref_seeds: Optional[Dict[str, float]] = None, gamma_seeds: Optional[np.ndarray] = None, @@ -3842,6 +3928,7 @@ def _execute_jac_vec_prod_fwd( con_seeds: Case constraint AD seeds geom_seeds: Geometric AD seeds in the same format as the geometric data mesh_seeds: Mesh geometry AD seeds in the same format as the mesh data + dvgeo_seeds: DVGeo DV seeds for a given surface and then design variable param_seeds: Case parameter AD seeds ref_seeds: Reference condition AD seeds gamma_seeds: Circulation AD seeds @@ -3875,6 +3962,9 @@ def _execute_jac_vec_prod_fwd( if mesh_seeds is None: mesh_seeds = {} + if self.DVGeo is None: + dvgeo_seeds = {} + if gamma_seeds is None: gamma_seeds = np.zeros(mesh_size) @@ -3899,13 +3989,38 @@ def _execute_jac_vec_prod_fwd( # self.clear_ad_seeds() self.set_variable_ad_seeds(con_seeds) self.set_geom_ad_seeds(geom_seeds) - self.set_mesh_ad_seeds(mesh_seeds) + # self.set_mesh_ad_seeds(mesh_seeds) self.set_gamma_ad_seeds(gamma_seeds) self.set_gamma_d_ad_seeds(gamma_d_seeds) self.set_gamma_u_ad_seeds(gamma_u_seeds) self.set_parameter_ad_seeds(param_seeds) self.set_reference_ad_seeds(ref_seeds) + # Since DVGeo seeds operate entirely within the python layer we set them here + if self.DVGeo is not None and dvgeo_seeds is not None: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a pointset, skip it + + # If no mesh seed was provided for the surface we will need to start it at zero + if surface not in mesh_seeds.keys(): + nx = self.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + mesh_seeds[surface] = {} + mesh_seeds[surface]["mesh"] = np.zeros((nx*ny,3)) + + # Loop over the design variables and accumulate the sensitivity product into the mesh_seeds + mesh_seeds[surface]["mesh"] += self.DVGeo.totalSensitivityProd(dvgeo_seeds[surface], point_set_name).reshape( + mesh_seeds[surface]["mesh"].shape + ) + + self.set_mesh_ad_seeds(mesh_seeds) + self.avl.update_surfaces_d() self.avl.get_res_d() self.avl.velsum_d() @@ -3942,6 +4057,37 @@ def _execute_jac_vec_prod_fwd( self.set_parameter_ad_seeds(param_seeds, mode="FD", scale=step) self.set_reference_ad_seeds(ref_seeds, mode="FD", scale=step) + + # Since DVGeo operates entirely within the python layer we have have to do this + if self.DVGeo is not None and dvgeo_seeds is not None: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a pointset, skip it + + # Apply the FD step to the DV seeds + for dv in dvgeo_seeds[surface].keys(): + # current design var values + currentDV = self.DVGeo.getValues()[dv] + + # Set the updated DVs + self.DVGeo.setDesignVars({dv: currentDV + dvgeo_seeds[surface][dv]*step}) + + # get mesh size + nx = self.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + + # Compute the perturbed mesh values + coords = self.DVGeo.update(point_set_name) + mesh_pertub = copy.deepcopy(coords.reshape((ny,nx,3)).transpose((1,0,2))) + + self.set_mesh(idx_surf,mesh_pertub) + + # propogate the seeds through without resolving self.avl.update_surfaces() self.avl.get_res() @@ -3966,6 +4112,36 @@ def _execute_jac_vec_prod_fwd( self.set_parameter_ad_seeds(param_seeds, mode="FD", scale=-1 * step) self.set_reference_ad_seeds(ref_seeds, mode="FD", scale=-1 * step) + # Set the mesh seeds back + if self.DVGeo is not None and dvgeo_seeds is not None: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a pointset, skip it + + # Restore the orignal dvgeo seeds + for dv in dvgeo_seeds[surface].keys(): + # current design var values + currentDV = self.DVGeo.getValues()[dv] + + # Set the updated DVs + self.DVGeo.setDesignVars({dv: currentDV - dvgeo_seeds[surface][dv]*step}) + + # get mesh size + nx = self.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + + # Compute the original mesh values + coords = self.DVGeo.update(point_set_name) + mesh_orig = copy.deepcopy(coords.reshape((ny,nx,3)).transpose((1,0,2))) + + self.set_mesh(idx_surf,mesh_orig) + + self.avl.update_surfaces() self.avl.get_res() self.avl.velsum() @@ -4130,6 +4306,30 @@ def _execute_jac_vec_prod_rev( print(f" Time to extract seeds: {time.time() - time_last}") time_last = time.time() + # Create dv_geo seeds a empty dict of surface keys + dvgeo_seeds = {} + for surf_key in self.unique_surface_names: + dvgeo_seeds[surf_key] = {} + + # If a DVGeo is present then propagate the mesh seeds all the way back to the DVs + if self.DVGeo is not None and self.DVGeo.getNDV() > 0: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a mesh, skip it + dvgeo_seeds[surface] = {} + + # Get the sensitivities + dvgeo_seeds[surface].update( + self.DVGeo.totalSensitivity(mesh_seeds[surface]["mesh"], point_set_name) + ) + + + self.set_function_ad_seeds(func_seeds, scale=0.0) self.set_residual_ad_seeds(res_seeds, scale=0.0) self.set_residual_d_ad_seeds(res_d_seeds, scale=0.0) @@ -4144,7 +4344,7 @@ def _execute_jac_vec_prod_rev( if print_timings: print(f" Total Time: {time.time() - time_start}") - return con_seeds, geom_seeds, mesh_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds + return con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds def execute_run_sensitivities( self, @@ -4178,7 +4378,7 @@ def execute_run_sensitivities( # TODO: remove seeds if it doesn't effect accuracy # self.clear_ad_seeds() time_last = time.time() - _, _, _, pfpU, _, _, _, _ = self._execute_jac_vec_prod_rev(func_seeds={func: 1.0}) + _, _, _, _, pfpU, _, _, _, _ = self._execute_jac_vec_prod_rev(func_seeds={func: 1.0}) if print_timings: print(f"Time to get RHS: {time.time() - time_last}") time_last = time.time() @@ -4195,7 +4395,7 @@ def execute_run_sensitivities( # get the resulting adjoint vector (dfunc/dRes) from fortran dfdR = self.get_residual_ad_seeds() # self.clear_ad_seeds() - con_seeds, geom_seeds, mesh_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( + con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( func_seeds={func: 1.0}, res_seeds=dfdR ) if print_timings: @@ -4205,7 +4405,7 @@ def execute_run_sensitivities( sens[func].update(con_seeds) # I don't know if it's worth combining geom_seeds and mesh_seeds into one just to make this one part less nasty for key in geom_seeds: - sens[func][key] = geom_seeds[key] | mesh_seeds[key] + sens[func][key] = geom_seeds[key] | mesh_seeds[key] | dvgeo_seeds[key] # sens[func].update(geom_seeds) # sens[func].update(mesh_seeds) sens[func].update(param_seeds) @@ -4222,7 +4422,7 @@ def execute_run_sensitivities( # get the RHS of the adjoint equation (pFpU) # TODO: remove seeds if it doesn't effect accuracy - _, _, _, pfpU, pf_pU_d, _, _, _ = self._execute_jac_vec_prod_rev(consurf_derivs_seeds={func_key: 1.0}) + _, _, _, _, pfpU, pf_pU_d, _, _, _ = self._execute_jac_vec_prod_rev(consurf_derivs_seeds={func_key: 1.0}) if print_timings: print(f"Time to get RHS: {time.time() - time_last}") time_last = time.time() @@ -4242,7 +4442,7 @@ def execute_run_sensitivities( dfdR = self.get_residual_ad_seeds() dfdR_d = self.get_residual_d_ad_seeds() # self.clear_ad_seeds() - con_seeds, geom_seeds, mesh_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( + con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( consurf_derivs_seeds={func_key: 1.0}, res_seeds=dfdR, res_d_seeds=dfdR_d ) if print_timings: @@ -4251,7 +4451,7 @@ def execute_run_sensitivities( sens[func_key].update(con_seeds) for key in geom_seeds: - sens[func_key][key] = geom_seeds[key] | mesh_seeds[key] + sens[func_key][key] = geom_seeds[key] | mesh_seeds[key] | dvgeo_seeds[key] # sens[func_key].update(geom_seeds) # sens[func_key].update(mesh_seeds) sens[func_key].update(param_seeds) @@ -4268,7 +4468,7 @@ def execute_run_sensitivities( # get the RHS of the adjoint equation (pFpU) # TODO: remove seeds if it doesn't effect accuracy - _, _, _, pfpU, _, pf_pU_u, _, _ = self._execute_jac_vec_prod_rev(stab_derivs_seeds={func_key: 1.0}) + _, _, _, _, pfpU, _, pf_pU_u, _, _ = self._execute_jac_vec_prod_rev(stab_derivs_seeds={func_key: 1.0}) if print_timings: print(f"Time to get RHS: {time.time() - time_last}") time_last = time.time() @@ -4288,7 +4488,7 @@ def execute_run_sensitivities( dfdR = self.get_residual_ad_seeds() dfdR_u = self.get_residual_u_ad_seeds() # self.clear_ad_seeds() - con_seeds, geom_seeds, mesh_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( + con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( stab_derivs_seeds={func_key: 1.0}, res_seeds=dfdR, res_u_seeds=dfdR_u ) @@ -4297,8 +4497,10 @@ def execute_run_sensitivities( time_last = time.time() sens[func_key].update(con_seeds) - for surf_key in geom_seeds: - sens[func_key][surf_key] = geom_seeds[surf_key] | mesh_seeds[surf_key] + for key in geom_seeds: + sens[func_key][key] = geom_seeds[key] | mesh_seeds[key] | dvgeo_seeds[key] + # sens[func_key].update(geom_seeds) + # sens[func_key].update(mesh_seeds) sens[func_key].update(param_seeds) sens[func_key].update(ref_seeds) # sd_deriv_seeds[func_key] = 0.0 @@ -4314,7 +4516,7 @@ def execute_run_sensitivities( # get the RHS of the adjoint equation (pFpU) # TODO: remove seeds if it doesn't effect accuracy - _, _, _, pfpU, _, pf_pU_u, _, _ = self._execute_jac_vec_prod_rev(body_axis_derivs_seeds={func_key: 1.0}) + _, _, _, _, pfpU, _, pf_pU_u, _, _ = self._execute_jac_vec_prod_rev(body_axis_derivs_seeds={func_key: 1.0}) if print_timings: print(f"Time to get RHS: {time.time() - time_last}") time_last = time.time() @@ -4334,7 +4536,7 @@ def execute_run_sensitivities( dfdR = self.get_residual_ad_seeds() dfdR_u = self.get_residual_u_ad_seeds() # self.clear_ad_seeds() - con_seeds, geom_seeds, mesh_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( + con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( body_axis_derivs_seeds={func_key: 1.0}, res_seeds=dfdR, res_u_seeds=dfdR_u ) @@ -4343,8 +4545,10 @@ def execute_run_sensitivities( time_last = time.time() sens[func_key].update(con_seeds) - for surf_key in geom_seeds: - sens[func_key][surf_key] = geom_seeds[surf_key] | mesh_seeds[surf_key] + for key in geom_seeds: + sens[func_key][key] = geom_seeds[key] | mesh_seeds[key] | dvgeo_seeds[key] + # sens[func_key].update(geom_seeds) + # sens[func_key].update(mesh_seeds) sens[func_key].update(param_seeds) sens[func_key].update(ref_seeds) diff --git a/optvl/utils/ffd_utils.py b/optvl/utils/ffd_utils.py new file mode 100644 index 0000000..e3cdc7e --- /dev/null +++ b/optvl/utils/ffd_utils.py @@ -0,0 +1,76 @@ +# ============================================================================= +# Standard Python Modules +# ============================================================================= +import os + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np + +# ============================================================================= +# Local modules +# ============================================================================= + + +def write_FFD_file(mesh, mx, my, filename="ffd", cushion=0.05): + """Utility functions that generates a box FFD around a wing mesh. + Based on the utility function in OpenAeroStruct. + + Args: + mesh: mesh array for the surface to be embedded in the FFD np.ndarray (nx,ny,3) + mx: number of control points in the x-direction + my: number of control points in the y-direction + filename: name of the ffd file to write without extension.Defaults to "ffd.dat". + cushion: Amount of cushion to apply from the LE/TE and tips. Defaults to 0.05. + + Returns: + filename of the written FFD file as a str + """ + nx, ny = mesh.shape[:2] + + half_ffd = np.zeros((mx, my, 3)) + + LE = mesh[0, :, :] + TE = mesh[-1, :, :] + + half_ffd[0, :, 0] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), LE[:, 0]) + half_ffd[0, :, 1] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), LE[:, 1]) + half_ffd[0, :, 2] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), LE[:, 2]) + + half_ffd[-1, :, 0] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), TE[:, 0]) + half_ffd[-1, :, 1] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), TE[:, 1]) + half_ffd[-1, :, 2] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), TE[:, 2]) + + for i in range(my): + half_ffd[:, i, 0] = np.linspace(half_ffd[0, i, 0], half_ffd[-1, i, 0], mx) + half_ffd[:, i, 1] = np.linspace(half_ffd[0, i, 1], half_ffd[-1, i, 1], mx) + half_ffd[:, i, 2] = np.linspace(half_ffd[0, i, 2], half_ffd[-1, i, 2], mx) + + half_ffd[0, :, 0] -= cushion + half_ffd[-1, :, 0] += cushion + half_ffd[:, 0, 1] -= cushion + half_ffd[:, -1, 1] += cushion + + bottom_ffd = half_ffd.copy() + bottom_ffd[:, :, 2] -= cushion + + top_ffd = half_ffd.copy() + top_ffd[:, :, 2] += cushion + + ffd = np.vstack((bottom_ffd, top_ffd)) + + filename = filename + ".xyz" + + with open(filename, "w") as f: + f.write("1\n") + f.write("{} {} {}\n".format(mx, 2, my)) + x = np.array_str(ffd[:, :, 0].flatten(order="F"))[1:-1] + "\n" + y = np.array_str(ffd[:, :, 1].flatten(order="F"))[1:-1] + "\n" + z = np.array_str(ffd[:, :, 2].flatten(order="F"))[1:-1] + "\n" + + f.write(x) + f.write(y) + f.write(z) + + return filename \ No newline at end of file diff --git a/tests/test_body_axis_derivs_partial_derivs.py b/tests/test_body_axis_derivs_partial_derivs.py index 7359402..a471e97 100644 --- a/tests/test_body_axis_derivs_partial_derivs.py +++ b/tests/test_body_axis_derivs_partial_derivs.py @@ -196,7 +196,7 @@ def test_rev_gamma_u(self): # for var_key in bd_d_fwd[deriv_func]: bd_d_rev = {deriv_func: 1.0} - gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(body_axis_derivs_seeds=bd_d_rev)[5] + gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(body_axis_derivs_seeds=bd_d_rev)[6] rev_sum = np.sum(gamma_u_seeds_rev * gamma_u_seeds_fwd) @@ -247,7 +247,7 @@ def test_rev_ref(self): for deriv_func, var_dict in self.ovl_solver.case_body_derivs_to_fort_var.items(): body_axis_deriv_seeds_rev[deriv_func] = np.random.rand(1)[0] - ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(body_axis_derivs_seeds=body_axis_deriv_seeds_rev)[7] + ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(body_axis_derivs_seeds=body_axis_deriv_seeds_rev)[8] self.ovl_solver.clear_ad_seeds_fast() diff --git a/tests/test_consurf_partial_derivs.py b/tests/test_consurf_partial_derivs.py index f1a0939..bb82ec2 100644 --- a/tests/test_consurf_partial_derivs.py +++ b/tests/test_consurf_partial_derivs.py @@ -168,7 +168,7 @@ def test_rev_gamma_d(self): res_d_seeds_fwd = self.ovl_solver._execute_jac_vec_prod_fwd(gamma_d_seeds=gamma_d_seeds_fwd)[5] self.ovl_solver.clear_ad_seeds_fast() - gamma_d_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_d_seeds=res_d_seeds_rev)[4] + gamma_d_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_d_seeds=res_d_seeds_rev)[5] gamma_sum = np.sum(gamma_d_seeds_rev * gamma_d_seeds_fwd) res_sum = np.sum(res_d_seeds_rev * res_d_seeds_fwd) @@ -359,8 +359,7 @@ def test_rev_gamma_d(self): for deriv_func in cs_d_fwd: cs_d_rev = {deriv_func: 1.0} - gamma_d_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(consurf_derivs_seeds=cs_d_rev)[4] - + gamma_d_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(consurf_derivs_seeds=cs_d_rev)[5] rev_sum = np.sum(gamma_d_seeds_rev * gamma_d_seeds_fwd) fwd_sum = np.sum(cs_d_fwd[deriv_func]) diff --git a/tests/test_om_wrapper.py b/tests/test_om_wrapper.py index 2a589e1..e20f493 100644 --- a/tests/test_om_wrapper.py +++ b/tests/test_om_wrapper.py @@ -3,10 +3,12 @@ # ============================================================================= from optvl import OVLSolver, OVLGroup + # ============================================================================= # Standard Python Modules # ============================================================================= import os +import sys # ============================================================================= # External Python modules @@ -28,6 +30,11 @@ geom_file = os.path.join(geom_dir, "aircraft.avl") mass_file = os.path.join(geom_dir, "aircraft.mass") +# Add geom_files to path for importing +sys.path.insert(0, geom_dir) + +# Import input_dict from the .py files instead of loading from .pkl +from aircraft import input_dict class TestOMWrapper(unittest.TestCase): def setUp(self): @@ -220,6 +227,162 @@ def test_OM_total_derivs(self): err_msg=f"deriv of {key[0]} wrt {key[1]} does not agree with FD to rtol={rtol}", ) +class TestOMWrapperMesh(unittest.TestCase): + def setUp(self): + self.ovl_solver = OVLSolver(input_dict=input_dict, mass_file=mass_file) + + model = om.Group() + model.add_subsystem( + "ovlsolver", + OVLGroup( + input_dict=input_dict, + mass_file=mass_file, + output_stability_derivs=True, + output_body_axis_derivs=True, + input_param_vals=True, + input_ref_vals=True, + input_airfoil_geom=True, + input_mesh_dim=3, + ), + ) + + self.prob = om.Problem(model) + + def test_aero_coef(self): + self.ovl_solver.execute_run() + run_data = self.ovl_solver.get_total_forces() + + prob = self.prob + prob.setup(mode="rev") + prob.run_model() + + for func in run_data: + om_val = prob.get_val(f"ovlsolver.{func}") + assert om_val == run_data[func] + + def test_surface_mesh_setting(self): + prob = self.prob + prob.setup(mode="rev") + + np.random.seed(111) + + for surf in self.ovl_solver.unique_surface_names: + idx_surf = self.ovl_solver.get_surface_index(surf) + arr = self.ovl_solver.get_mesh(idx_surf) + arr += np.random.rand(*arr.shape) * 0.1 + # print(f'setting {surf_key}:{geom_key} to {arr}') + # #set surface data + self.ovl_solver.set_mesh(idx_surf, arr) + self.ovl_solver.avl.update_surfaces() + self.ovl_solver.execute_run() + + # set om surface data + prob.set_val(f"ovlsolver.{surf}:mesh", arr) + prob.run_model() + + run_data = self.ovl_solver.get_total_forces() + for func in run_data: + om_val = prob.get_val(f"ovlsolver.{func}") + assert om_val == run_data[func] + + stab_derivs = self.ovl_solver.get_stab_derivs() + for func in stab_derivs: + om_val = prob.get_val(f"ovlsolver.{func}") + assert om_val == stab_derivs[func] + + body_axis_derivs = self.ovl_solver.get_body_axis_derivs() + for func in body_axis_derivs: + om_val = prob.get_val(f"ovlsolver.{func}") + assert om_val == body_axis_derivs[func] + + def test_CL_solve(self): + prob = self.prob + cl_star = 0.9 + prob.model.add_design_var("ovlsolver.alpha", lower=-10, upper=10) + prob.model.add_constraint("ovlsolver.CL", equals=cl_star) + prob.model.add_objective("ovlsolver.CD", scaler=1e3) + prob.setup(mode="rev") + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["debug_print"] = ["desvars", "ln_cons", "nl_cons", "objs"] + prob.driver.options["tol"] = 1e-6 + prob.driver.options["disp"] = True + + prob.setup(mode="rev") + prob.run_driver() + om.n2(prob, show_browser=False, outfile="vlm_opt.html") + + om_val = prob.get_val("ovlsolver.alpha") + self.ovl_solver.set_trim_condition("CL", cl_star) + self.ovl_solver.execute_run() + alpha = self.ovl_solver.get_parameter("alpha") + + np.testing.assert_allclose( + om_val, + alpha, + rtol=1e-5, + err_msg="solved alpha", + ) + + def test_CM_solve(self): + prob = self.prob + prob.model.add_design_var("ovlsolver.alpha", lower=-10, upper=10) + prob.model.add_constraint("ovlsolver.Cm", equals=0.0, scaler=1e3) + prob.model.add_objective("ovlsolver.CD", scaler=1e3) + prob.setup(mode="rev") + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["debug_print"] = ["desvars", "ln_cons", "nl_cons", "objs"] + prob.driver.options["tol"] = 1e-6 + prob.driver.options["disp"] = True + + prob.setup(mode="rev") + prob.run_driver() + om.n2(prob, show_browser=False, outfile="vlm_opt.html") + + om_val = prob.get_val(f"ovlsolver.alpha") + + self.ovl_solver.set_constraint("alpha", "Cm", 0.00) + self.ovl_solver.execute_run() + alpha = self.ovl_solver.get_parameter("alpha") + + np.testing.assert_allclose( + om_val, + alpha, + rtol=1e-5, + err_msg="solved alpha", + ) + + def test_OM_total_derivs(self): + prob = self.prob + cl_star = 0.9 + dcl_dalpha_star = -0.05 + # prob.model.add_design_var("ovlsolver.Wing:mesh") # This is a really costly test + # prob.model.add_design_var("ovlsolver.Wing:aincs") # dCL/dalpha on the first element fails + prob.model.add_design_var("ovlsolver.Elevator", lower=-10, upper=10) + prob.model.add_design_var("ovlsolver.alpha", lower=-10, upper=10) + prob.model.add_design_var("ovlsolver.Sref") + prob.model.add_design_var("ovlsolver.Mach") + prob.model.add_design_var("ovlsolver.X cg") + prob.model.add_constraint("ovlsolver.CL", equals=cl_star) + prob.model.add_constraint("ovlsolver.dCL/dalpha", equals=-dcl_dalpha_star) + prob.model.add_constraint("ovlsolver.dCl/dp", equals=-dcl_dalpha_star) + prob.model.add_objective("ovlsolver.CD", scaler=1e3) + prob.model.add_objective("ovlsolver.Cm", scaler=1e3) + prob.setup(mode="rev") + prob.run_model() + om.n2(prob, show_browser=False, outfile="vlm_opt.html") + deriv_err = prob.check_totals(step=1e-7, form="central") + rtol = 1e-2 + for key, data in deriv_err.items(): + np.testing.assert_allclose( + data["J_fd"], + data["J_rev"], + rtol=rtol, + atol=1e-9, + err_msg=f"deriv of {key[0]} wrt {key[1]} does not agree with FD to rtol={rtol}", + ) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_partial_derivs.py b/tests/test_partial_derivs.py index d1527ca..4869b84 100644 --- a/tests/test_partial_derivs.py +++ b/tests/test_partial_derivs.py @@ -16,6 +16,7 @@ import numpy as np + base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder geom_dir = os.path.join(base_dir, '..', 'geom_files') @@ -215,7 +216,7 @@ def test_rev_gamma(self): self.ovl_solver.clear_ad_seeds_fast() for func_key in self.ovl_solver.case_var_to_fort_var: - gamma_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[3] + gamma_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[4] rev_sum = np.sum(gamma_seeds_rev * gamma_seeds_fwd) fwd_sum = np.sum(func_seeds_fwd[func_key]) @@ -262,7 +263,7 @@ def test_rev_param(self): self.ovl_solver.clear_ad_seeds_fast() for func_key in self.ovl_solver.case_var_to_fort_var: - param_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[6] + param_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[7] # print(f"{func_key} wrt {param_key}", "fwd ", func_seeds_fwd[func_key], "rev", param_seeds_rev[param_key]) tol = 1e-14 @@ -316,7 +317,7 @@ def test_rev_ref(self): self.ovl_solver.clear_ad_seeds_fast() for func_key in self.ovl_solver.case_var_to_fort_var: - ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[7] + ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[8] # print(f"{func_key} wrt {ref_key}", "fwd ", func_seeds_fwd[func_key], "rev", ref_seeds_rev[ref_key]) tol = 1e-14 @@ -465,7 +466,7 @@ def test_fwd_param(self): def test_rev_param(self): num_res = self.ovl_solver.get_mesh_size() res_seeds_rev = np.random.rand(num_res) - param_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_seeds=res_seeds_rev)[6] + param_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_seeds=res_seeds_rev)[7] self.ovl_solver.clear_ad_seeds_fast() @@ -498,7 +499,7 @@ def test_fwd_ref(self): def test_rev_ref(self): num_res = self.ovl_solver.get_mesh_size() res_seeds_rev = np.random.rand(num_res) - ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_seeds=res_seeds_rev)[7] + ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_seeds=res_seeds_rev)[8] self.ovl_solver.clear_ad_seeds_fast() diff --git a/tests/test_pygeo.py b/tests/test_pygeo.py new file mode 100644 index 0000000..1e802d6 --- /dev/null +++ b/tests/test_pygeo.py @@ -0,0 +1,466 @@ +# ============================================================================= +# Standard modules +# ============================================================================= +import os +import copy + +# ============================================================================= +# External Python modules +# ============================================================================= +import unittest +import numpy as np +import sys +import psutil + +# ============================================================================= +# Extension modules +# ============================================================================= +from optvl import OVLSolver +from optvl.utils.ffd_utils import write_FFD_file + +try: + from pygeo import DVGeometry + HAS_PYGEO = True +except ImportError: + HAS_PYGEO = False + +# Add geom_files to path for importing +base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder +geom_dir = os.path.join(base_dir, "..", "geom_files") +sys.path.insert(0, geom_dir) + +from wing_mesh import mesh # (nx, ny, 3) numpy array + +# --------------------------------------------------------------------------- +# Path setup — reuse the same mesh used in test_mesh_input.py +# --------------------------------------------------------------------------- +base_dir = os.path.dirname(os.path.abspath(__file__)) +geom_dir = os.path.join(base_dir, '..', 'geom_files') +sys.path.insert(0, geom_dir) + + +surf = { + "Wing": { + # General + "component": np.int32(1), # logical surface component index (for grouping interacting surfaces, see AVL manual) + "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + "claf": 1.0, # CL alpha (dCL/da) scaling factor per section (provide a single entry and OptVL applies to all strips, otherwise provide a vector corresponding to each strip) + + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + "angle": np.float64(0.0), # offset added on to the Ainc values for all the defining sections in this surface + "aincs": np.ones(mesh.shape[1]), # incidence angle vector (provide a single entry and OptVL applies to all strips, otherwise provide a vector corresponding to each strip) + + # Geometry: Mesh + "mesh": np.float64(mesh), # (nx,ny,3) numpy array containing mesh coordinates + "flatten mesh": True, # True by default so can be turned off or just excluded (not recommended) + + # Control Surface Specification + "control_assignments": { + "flap" : {"assignment":np.arange(0,mesh.shape[1]), + "xhinged": 0.8, # x/c location of hinge + "vhinged": np.zeros(3), # vector giving hinge axis about which surface rotates + "gaind": 1.0, # control surface gain + "refld": 1.0 # control surface reflection, sign of deflection for duplicated surface + } + }, + + # Design Variables (AVL) Specification + "design_var_assignments": { + "des" : {"assignment":np.arange(0,mesh.shape[1]), + "gaing":1.0} + }, + } +} + +geom = { + "title": "Aircraft", + "mach": np.float64(0.0), + "iysym": np.int32(0), + "izsym": np.int32(0), + "zsym": np.float64(0.0), + "Sref": np.float64(10.0), + "Cref": np.float64(1.0), + "Bref": np.float64(10.0), + "XYZref": np.array([0.25, 0, 0],dtype=np.float64), + "CDp": np.float64(0.0), + "surfaces": surf, + # Global Control and DV info + "dname": ["flap"], # Name of control input for each corresonding index + "gname": ["des"], # Name of design var for each corresonding index +} + +@unittest.skipUnless(HAS_PYGEO, "pygeo is not installed — skipping FFD tests") +class TestFFDIntegration(unittest.TestCase): + """Tests for the pygeo FFD integration in OVLSolver + + Each test creates its own OVLSolver instance so that solver state is + isolated between test cases. + """ + + # ------------------------------------------------------------------ + # Shared FFD file (written once for the whole test class) + # ------------------------------------------------------------------ + ffd_path = os.path.join(geom_dir, "wing_test_ffd.xyz") + + def setUp(self): + self._make_solver() + self._make_dvgeo() + + def tearDown(self): + # Get the memory usage of the current process using psutil + process = psutil.Process() + mb_memory = process.memory_info().rss / (1024 * 1024) # Convert bytes to MB + print(f"{self.id():80} Memory usage: {mb_memory:.2f} MB") + + def _make_solver(self): + """Sets a freshly initialised OVLSolver with the mesh geometry.""" + self.ovl_solver = OVLSolver(input_dict=copy.deepcopy(geom)) + + def _make_dvgeo(self): + """Sets DVGeometry object built from the class-level FFD file to OptVL""" + self.DVGeo = DVGeometry(self.ffd_path) + self.ovl_solver.set_DVGeo(self.DVGeo) + # self.ovl_solver.update_DVGeo() + + def finite_dif(self, dvgeo_seeds, step=1e-7): + + # Loop over all surfaces + for surface in self.ovl_solver.unique_surface_names: + # Get the pointset name + idx_surf = self.ovl_solver.get_surface_index(surf_name=surface) + if idx_surf in self.ovl_solver.point_sets.keys(): + point_set_name = self.ovl_solver.point_sets[idx_surf] + else: + continue # This surface doesn't have a pointset, skip it + + # Apply the FD step to the DV seeds + for dv in dvgeo_seeds[surface].keys(): + # current design var values + currentDV = self.ovl_solver.DVGeo.getValues()[dv] + + # Set the updated DVs + self.ovl_solver.DVGeo.setDesignVars({dv: currentDV + dvgeo_seeds[surface][dv]*step}) + + # get mesh size + nx = self.ovl_solver.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl_solver.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + + # Compute the perturbed mesh values + coords = self.ovl_solver.DVGeo.update(point_set_name) + mesh_pertub = copy.deepcopy(coords.reshape((ny,nx,3)).transpose((1,0,2))) + + self.ovl_solver.set_mesh(idx_surf,mesh_pertub) + + self.ovl_solver.avl.update_surfaces() + self.ovl_solver.avl.get_res() + self.ovl_solver.avl.exec_rhs() + self.ovl_solver.avl.get_res() + self.ovl_solver.avl.velsum() + self.ovl_solver.avl.aero() + # self.ovl_solver.execute_run() + coef_data_peturb = self.ovl_solver.get_total_forces() + consurf_derivs_peturb = self.ovl_solver.get_control_stab_derivs() + stab_deriv_derivs_peturb = self.ovl_solver.get_stab_derivs() + body_axis_deriv_petrub = self.ovl_solver.get_body_axis_derivs() + body_forces_peturb = self.ovl_solver.get_body_forces() + + for surface in self.ovl_solver.unique_surface_names: + # Get the pointset name + idx_surf = self.ovl_solver.get_surface_index(surf_name=surface) + if idx_surf in self.ovl_solver.point_sets.keys(): + point_set_name = self.ovl_solver.point_sets[idx_surf] + else: + continue # This surface doesn't have a pointset, skip it + + # Apply the FD step to the DV seeds + for dv in dvgeo_seeds[surface].keys(): + # current design var values + currentDV = self.ovl_solver.DVGeo.getValues()[dv] + + # Set the updated DVs + self.ovl_solver.DVGeo.setDesignVars({dv: currentDV - dvgeo_seeds[surface][dv]*step}) + + # get mesh size + nx = self.ovl_solver.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.ovl_solver.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + + # Compute the perturbed mesh values + coords = self.ovl_solver.DVGeo.update(point_set_name) + mesh_orig = copy.deepcopy(coords.reshape((ny,nx,3)).transpose((1,0,2))) + + self.ovl_solver.set_mesh(idx_surf,mesh_orig) + + self.ovl_solver.avl.update_surfaces() + self.ovl_solver.avl.get_res() + self.ovl_solver.avl.exec_rhs() + self.ovl_solver.avl.get_res() + self.ovl_solver.avl.velsum() + self.ovl_solver.avl.aero() + # self.ovl_solver.execute_run() + + coef_data = self.ovl_solver.get_total_forces() + consurf_derivs = self.ovl_solver.get_control_stab_derivs() + stab_deriv_derivs = self.ovl_solver.get_stab_derivs() + body_axis_deriv = self.ovl_solver.get_body_axis_derivs() + body_forces = self.ovl_solver.get_body_forces() + + body_func_seeds = {} + for body in body_forces: + body_func_seeds[body] = {} + for key in body_forces[body]: + body_func_seeds[body][key] = (body_forces_peturb[body][key] - body_forces[body][key]) / step + + + func_seeds = {} + for func_key in coef_data: + func_seeds[func_key] = (coef_data_peturb[func_key] - coef_data[func_key]) / step + + consurf_derivs_seeds = {} + for func_key in consurf_derivs: + consurf_derivs_seeds[func_key] = (consurf_derivs_peturb[func_key] - consurf_derivs[func_key]) / step + + stab_derivs_seeds = {} + for func_key in stab_deriv_derivs: + stab_derivs_seeds[func_key] = (stab_deriv_derivs_peturb[func_key] - stab_deriv_derivs[func_key]) / step + + body_axis_derivs_seeds = {} + for deriv_func in body_axis_deriv: + body_axis_derivs_seeds[deriv_func] = ( + body_axis_deriv_petrub[deriv_func] - body_axis_deriv[deriv_func] + ) / step + + return func_seeds, consurf_derivs_seeds, stab_derivs_seeds, body_axis_derivs_seeds + + def test_mesh_deformation_propagates(self): + self._make_solver() + self._make_dvgeo() + + idx_surf = self.ovl_solver.get_surface_index("Wing") + mesh_orig = copy.deepcopy(self.ovl_solver.get_mesh(idx_surf)) # (nx, ny, 3) + DVGeo = self.ovl_solver.DVGeo + + DVGeo.addLocalDV("local_z", lower=-1.0, upper=1.0, axis="z", scale=1.0) + dz = 0.1 + + # Perturb all local z DVs by dz + current_dvs = DVGeo.getValues() + for key in current_dvs: + current_dvs[key] = current_dvs[key] + dz + DVGeo.setDesignVars(current_dvs) + + self.ovl_solver.update_DVGeo() + mesh_perturb = self.ovl_solver.get_mesh(idx_surf) # (nx, ny, 3) + + # The z-coordinates should have shifted by approximately dz + dz_actual = mesh_perturb[:, :, 2] - mesh_orig[:, :, 2] + np.testing.assert_allclose( + dz_actual, dz, + atol=1e-3, + err_msg="Z-coordinates of the mesh should shift by the applied FFD dz offset", + ) + + # Y-coordinates should be unchanged (we only moved z) + dy_actual = mesh_perturb[:, :, 1] - mesh_orig[:, :, 1] + np.testing.assert_allclose( + dy_actual, 0.0, + atol=1e-10, + err_msg="Y-coordinates should not change when only a Z FFD deformation is applied", + ) + + + def test_totals(self): + + self._make_solver() + self._make_dvgeo() + + # Add some DVs + nrefaxpts = self.ovl_solver.DVGeo.addRefAxis("c4", xFraction=0.25, alignIndex="k") + + def twist(val, geo): + for i in range(nrefaxpts): + geo.rot_y["c4"].coef[i] = val[i] + + def sweep(val, geo): + # the extractCoef method gets the unperturbed ref axis control points + C = geo.extractCoef("c4") + C_orig = C.copy() + # we will sweep the wing about the first point in the ref axis + sweep_ref_pt = C_orig[0, :] + + theta = -val[0] * np.pi / 180 + # rot_mtx = np.array([[np.cos(theta), 0.0, -np.sin(theta)], [0.0, 1.0, 0.0], [np.sin(theta), 0.0, np.cos(theta)]]) + rot_mtx = np.array([[np.cos(theta), -np.sin(theta), 0.0], [np.sin(theta), np.cos(theta), 0.0], [0.0, 0.0, 1.0]]) + + # modify the control points of the ref axis + # by applying a rotation about the first point in the x-z plane + for i in range(nrefaxpts): + # get the vector from each ref axis point to the wing root + vec = C[i, :] - sweep_ref_pt + # need to now rotate this by the sweep angle and add back the wing root loc + C[i, :] = sweep_ref_pt + rot_mtx @ vec + # use the restoreCoef method to put the control points back in the right place + geo.restoreCoef(C, "c4") + + self.ovl_solver.DVGeo.addGlobalDV("twist", func=twist, value=np.zeros(nrefaxpts), lower=-10, upper=10, scale=0.05) + self.ovl_solver.DVGeo.addGlobalDV("sweep", func=sweep, value=0.0, lower=0, upper=45, scale=0.05) + self.ovl_solver.DVGeo.addLocalDV("shape", lower=-0.25, upper=0.25, axis="z", scale=1.0) + + # Execute the solver + self.ovl_solver.update_DVGeo() + self.ovl_solver.set_variable("alpha", 5.0) + self.ovl_solver.set_variable("beta", 0.0) + self.ovl_solver.set_parameter("Mach", 0.0) + self.ovl_solver.execute_run() + + # compare the analytical gradients with finite difference for dvgeo + surf_key = list(self.ovl_solver.surf_mesh_to_fort_var.keys())[0] + dvgeo_vars = self.ovl_solver.DVGeo.getVarNames() + cs_names = self.ovl_solver.get_control_names() + + consurf_vars = [] + for func_key in self.ovl_solver.case_derivs_to_fort_var: + consurf_vars.append(self.ovl_solver._get_deriv_key(cs_names[0], func_key)) + + func_vars = self.ovl_solver.case_var_to_fort_var + stab_derivs = self.ovl_solver.case_stab_derivs_to_fort_var + body_axis_derivs = self.ovl_solver.case_body_derivs_to_fort_var + + sens = self.ovl_solver.execute_run_sensitivities( + func_vars, + consurf_derivs=consurf_vars, + stab_derivs=stab_derivs, + body_axis_derivs=body_axis_derivs, + print_timings=False, + ) + + # for con_key in self.ovl_solver.con_var_to_fort_var: + sens_FD = {} + for surf_key in self.ovl_solver.surf_geom_to_fort_var: + sens_FD[surf_key] = {} + for dvgeo_var_key in dvgeo_vars: + # arr = self.ovl_solver.get_mesh(self.ovl_solver.get_surface_index(surf_key)).reshape(-1,3) + arr = self.ovl_solver.DVGeo.getValues()[dvgeo_var_key] + np.random.seed(arr.size) + rand_arr = np.random.rand(*arr.shape) + rand_arr /= np.linalg.norm(rand_arr) + + func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif({surf_key: {dvgeo_var_key: rand_arr}}, step=1.0e-7) + + for func_key in func_vars: + dvgeo_var_dot = np.sum(sens[func_key][surf_key][dvgeo_var_key] * rand_arr) + func_dot = func_seeds[func_key] + + rel_err = np.abs(dvgeo_var_dot - func_dot) / np.abs(func_dot + 1e-20) + + # print( + # f"{func_key:5} wrt {surf_key}:{dvgeo_var_key:10} | AD:{dvgeo_var_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" + # ) + tol = 1e-7 + if np.abs(dvgeo_var_dot) < tol or np.abs(func_dot) < tol: + # If either value is basically zero, use an absolute tolerance + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + atol=1e-4, + err_msg=f"{func_key:5} wrt {surf_key}:{dvgeo_var_key:10}", + ) + else: + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + rtol=5e-3, + err_msg=f"{func_key:5} wrt {surf_key}:{dvgeo_var_key:10}", + ) + + for func_key in consurf_vars: + # for cs_key in consurf_vars[func_key]: + dvgeo_var_dot = np.sum(sens[func_key][surf_key][dvgeo_var_key] * rand_arr) + func_dot = consurf_deriv_seeds[func_key] + + # rel_err = np.abs(dvgeo_var_dot - func_dot) / np.abs(func_dot + 1e-20) + # print( + # f"{func_key} wrt {surf_key}:{mesh_key:10} | AD:{dvgeo_var_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" + # ) + + tol = 1e-8 + if np.abs(dvgeo_var_dot) < tol or np.abs(func_dot) < tol: + # If either value is basically zero, use an absolute tolerance + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + atol=1e-4, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + else: + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + rtol=6e-3, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + + for func_key in stab_derivs_seeds: + if func_key == "spiral parameter" or func_key == "lateral parameter": + continue # These derivatives blow up in different way for both adjoint and FD + dvgeo_var_dot = np.sum(sens[func_key][surf_key][dvgeo_var_key] * rand_arr) + func_dot = stab_derivs_seeds[func_key] + + rel_err = np.abs(dvgeo_var_dot - func_dot) / np.abs(func_dot + 1e-20) + + # print( + # f"{func_key} wrt {surf_key}:{dvgeo_var_key:10} | AD:{dvgeo_var_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" + # ) + + tol = 5e-7 + if np.abs(dvgeo_var_dot) < tol or np.abs(func_dot) < tol: + # If either value is basically zero, use an absolute tolerance + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + atol=5e-9, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + else: + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + rtol=6e-3, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + + for func_key in body_axis_derivs_seeds: + dvgeo_var_dot = np.sum(sens[func_key][surf_key][dvgeo_var_key] * rand_arr) + func_dot = body_axis_derivs_seeds[func_key] + + rel_err = np.abs(dvgeo_var_dot - func_dot) / np.abs(func_dot + 1e-20) + + # print( + # f"{func_key} wrt {surf_key}:{dvgeo_var_key:10} | AD:{dvgeo_var_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" + # ) + + tol = 1e-6 + if np.abs(dvgeo_var_dot) < tol or np.abs(func_dot) < tol: + # If either value is basically zero, use an absolute tolerance + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + atol=5e-8, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + else: + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + rtol=6e-3, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/tests/test_stab_derivs_partial_derivs.py b/tests/test_stab_derivs_partial_derivs.py index e404a84..4491f1c 100644 --- a/tests/test_stab_derivs_partial_derivs.py +++ b/tests/test_stab_derivs_partial_derivs.py @@ -165,7 +165,7 @@ def test_rev_gamma_u(self): res_u_seeds_fwd = self.ovl_solver._execute_jac_vec_prod_fwd(gamma_u_seeds=gamma_u_seeds_fwd)[6] self.ovl_solver.clear_ad_seeds_fast() - gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_u_seeds=res_u_seeds_rev)[5] + gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_u_seeds=res_u_seeds_rev)[6] gamma_sum = np.sum(gamma_u_seeds_rev * gamma_u_seeds_fwd) res_sum = np.sum(res_u_seeds_rev * res_u_seeds_fwd) @@ -378,7 +378,7 @@ def test_rev_gamma_u(self): # for var_key in sd_d_fwd[deriv_func]: sd_d_rev = {deriv_func: 1.0} - gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(stab_derivs_seeds=sd_d_rev)[5] + gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(stab_derivs_seeds=sd_d_rev)[6] rev_sum = np.sum(gamma_u_seeds_rev * gamma_u_seeds_fwd) @@ -425,7 +425,7 @@ def test_rev_ref(self): for deriv_func, var_dict in self.ovl_solver.case_stab_derivs_to_fort_var.items(): stab_deriv_seeds_rev[deriv_func] = np.random.rand(1)[0] - ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(stab_derivs_seeds=stab_deriv_seeds_rev)[7] + ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(stab_derivs_seeds=stab_deriv_seeds_rev)[8] self.ovl_solver.clear_ad_seeds_fast() diff --git a/tests/test_total_derivs.py b/tests/test_total_derivs.py index 5d562ad..e7b87c3 100644 --- a/tests/test_total_derivs.py +++ b/tests/test_total_derivs.py @@ -211,6 +211,7 @@ def test_geom(self): surf_key = list(self.ovl_solver.surf_geom_to_fort_var.keys())[0] geom_vars = self.ovl_solver.surf_geom_to_fort_var[surf_key] + # geom_vars += self.ovl_solver.surf_mesh_to_fort_var[surf_key] cs_names = self.ovl_solver.get_control_names() consurf_vars = [] diff --git a/tests/wing_mesh.npy b/tests/wing_mesh.npy deleted file mode 100644 index df1bec4..0000000 Binary files a/tests/wing_mesh.npy and /dev/null differ