diff --git a/geom_files/aircraft.py b/geom_files/aircraft.py new file mode 100644 index 00000000..fbf0647f --- /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 00000000..c347e97d --- /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 00000000..3c0e1453 --- /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 00000000..11bc6225 --- /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 00000000..0ab7effb --- /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 00000000..420353ef --- /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 00000000..59b4d895 --- /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 9a433951..0c9c8023 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,24 +35,29 @@ 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_case_strip_vars", types=bool, default=False) self.options.declare("output_stability_derivs", types=bool, default=False) self.options.declare("output_body_axis_derivs", types=bool, default=False) self.options.declare("output_con_surf_derivs", types=bool, default=False) 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_case_strip_vars = self.options["output_case_strip_vars"] 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 +66,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=["*"], ) @@ -69,9 +77,11 @@ def setup(self): input_param_vals=input_param_vals, input_ref_vals=input_ref_vals, input_airfoil_geom=input_airfoil_geom, + output_case_strip_vars=output_case_strip_vars, 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 +101,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 +118,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 +186,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 +207,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 +262,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 +273,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 +289,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 +336,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 +384,22 @@ 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: + 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((nx*ny,3))) + 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"].transpose((1,0,2)).reshape((nx*ny,3))) + # 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 +411,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 +426,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 +453,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: @@ -396,18 +496,21 @@ class OVLFuncsComp(om.ExplicitComponent): def initialize(self): self.options.declare("ovl", types=OVLSolver, recordable=False) + self.options.declare("output_case_strip_vars", types=bool, default=False) self.options.declare("output_stability_derivs", types=bool, default=False) self.options.declare("output_body_axis_derivs", types=bool, default=False) self.options.declare("output_con_surf_derivs", types=bool, default=False) 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,11 +529,22 @@ 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: self.add_output(func_key) + if self.options["output_case_strip_vars"]: + for func_key in self.ovl.case_strip_var_to_fort_var: + if func_key in ["CF strip", "Cm strip"]: #skip it for now + continue + for surface in self.ovl.unique_surface_names: + idx_surf = self.ovl.get_surface_index(surf_name=surface) + num_strips = self.ovl.get_surface_num_strips(idx_surf) + surf_func_name = f"{surface}:{func_key}" + self.add_output(surf_func_name,shape=(num_strips,)) + if self.options["output_con_surf_derivs"]: for func_key in self.ovl.case_derivs_to_fort_var: self.add_output(func_key) @@ -455,10 +569,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"] @@ -484,6 +609,15 @@ def compute(self, inputs, outputs): outputs[func_key] = run_data[func_key] # print(f" CD {run_data['CD']} CL {run_data['CL']}") + if self.options["output_case_strip_vars"]: + strip_data = self.ovl.get_strip_forces() + for func_key in self.ovl.case_strip_var_to_fort_var: + if func_key in ["CF strip", "Cm strip"]: #skip it for now + continue + for surface in self.ovl.unique_surface_names: + surf_func_name = f"{surface}:{func_key}" + outputs[surf_func_name] = strip_data[surface][func_key] + if self.options["output_con_surf_derivs"]: consurf_derivs_seeds = self.ovl.get_control_stab_derivs() for func_key in consurf_derivs_seeds: @@ -538,11 +672,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) - - func_seeds, _, csd_seeds, stab_derivs_seeds, body_axis_seeds, _, _ = self.ovl._execute_jac_vec_prod_fwd( + 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,strip_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, @@ -553,6 +700,11 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): for func_key in func_seeds: d_outputs[func_key] += func_seeds[func_key] + for surf in strip_func_seeds: + for strip_func_key in strip_func_seeds[surf]: + out_name = f"{surf}:{strip_func_key}" + d_outputs[out_name] += strip_func_seeds[surf][strip_func_key] + for func_key in csd_seeds: for con_name in csd_seeds[func_key]: var_name = f"d{func_key}_d{con_name}" @@ -579,6 +731,17 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): if np.abs(func_seeds[func_key]) > 0.0: print(f" running rev mode derivs for {func_key}") + strip_func_seeds = {} + for surf in self.ovl.unique_surface_names: + for strip_func_key in self.ovl.case_strip_var_to_fort_var: + out_name = f"{surf}:{strip_func_key}" + if out_name in d_outputs: + if surf not in strip_func_seeds.keys(): + strip_func_seeds[surf] = {} + strip_func_seeds[surf][strip_func_key] = d_outputs[out_name] + if np.any(np.abs(strip_func_seeds[surf][strip_func_key]) > 0.0): + print(f" running rev mode derivs for {surf}:{strip_func_key}") + csd_seeds = {} con_names = self.ovl.get_control_names() for func_key in self.ovl.case_derivs_to_fort_var: @@ -610,15 +773,27 @@ 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, + strip_func_seeds=strip_func_seeds, consurf_derivs_seeds=csd_seeds, stab_derivs_seeds=stab_derivs_seeds, body_axis_derivs_seeds=body_axis_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_inputs: d_inputs["gamma"] += gamma_seeds @@ -629,10 +804,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 +898,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 8ab2a536..7413e901 100644 --- a/optvl/optvl_class.py +++ b/optvl/optvl_class.py @@ -194,6 +194,49 @@ class OVLSolver(object): #TODO: add CF_LSRF(3,NFMAX), CM_LSRF(3,NFMAX) } + case_strip_var_to_fort_var = { + # geometric quantities + "X LE": ["STRP_R", "RLE", (slice(None), 0)], # control point leading edge coordinates + "Y LE": ["STRP_R", "RLE", (slice(None), 1)], # control point leading edge coordinates + "Z LE": ["STRP_R", "RLE", (slice(None), 2)], # control point leading edge coordinates + "chord": ["STRP_R", "CHORD"], + "width": ["STRP_R", "WSTRIP"], + "twist": ["STRP_R", "AINC"], + + # strip contributions to total lift and drag from strip integration + "CL": ["STRP_R", "CLSTRP"], + "CD": ["STRP_R", "CDSTRP"], + "CDv" : ["STRP_R","CDV_LSTRP"], # strip viscous drag in stability axes + "downwash" : ["STRP_R","DWWAKE"], + + + # strip contributions to non-dimensionalized forces + "CX": ["STRP_R", "CFSTRP", (slice(None), 0)], + "CY": ["STRP_R", "CFSTRP", (slice(None), 1)], + "CZ": ["STRP_R", "CFSTRP", (slice(None), 2)], + + # strip contributions to total moments (body frame) + "Cl": ["STRP_R", "CMSTRP", (slice(None), 0)], # previously CR + "Cm": ["STRP_R", "CMSTRP", (slice(None), 1)], # previously CM + "Cn": ["STRP_R", "CMSTRP", (slice(None), 2)], # previously CN + + + # forces non-dimentionalized by strip quantities + "CL strip" : ["STRP_R", "CL_LSTRP"], + "CD strip" : ["STRP_R", "CD_LSTRP"], + "CF strip" : ["STRP_R", "CF_LSTRP"], # forces in 3 directions + "Cm strip" : ["STRP_R", "CM_LSTRP"], # moments in 3 directions + + # additional forces and moments + "CL perp" : ["STRP_R", "CLT_LSTRP"], # strip CL referenced to Vperp, + "Cm c/4" : ["STRP_R","CMC4_LSTRP"], # strip pitching moment about c/4 and + "Cm LE" : ["STRP_R","CMLE_LSTRP"], # strip pitching moment about LE vector + "spanloading" : ["STRP_R","CNC"], # strip spanloading + + # TODO: add + # & CF_LSTRP(3,NSMAX), CM_LSTRP(3,NSMAX), ! strip forces in body axes referenced to strip area and 1/4 chord + } + body_geom_to_fort_var = { "scale": ["BODY_GEOM_R", "XYZSCAL_B"], "translate": ["BODY_GEOM_R", "XYZTRAN_B"], @@ -403,10 +446,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 +459,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 +835,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 +1108,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") @@ -2166,6 +2218,12 @@ def _get_surface_strip_indices(self, idx_surf: int): idx_srp_end = np.sum(num_strips[: idx_surf + 1]) return idx_srp_beg, idx_srp_end + + def get_surface_num_strips(self, idx_surf: int): + # num_strips = np.trim_zeros(self.get_avl_fort_arr("SURF_I", "NJ")) + num_strips = np.trim_zeros(self.get_avl_fort_arr("SURF_I", "NJ",slicer=idx_surf)) + + return num_strips # region --- modal analysis api def execute_eigen_mode_calc(self): @@ -3149,6 +3207,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: @@ -3503,7 +3637,7 @@ def get_mesh_ad_seeds(self) -> Dict[str, Dict[str, float]]: return mesh_seeds - def set_mesh_ad_seeds(self, mesh_seeds: Dict[str, float], mode: str = "AD", scale=1.0) -> None: + def set_mesh_ad_seeds(self, mesh_seeds: Dict[str, Dict[str, float]], mode: str = "AD", scale=1.0) -> None: for surf_key in mesh_seeds: for mesh_key in mesh_seeds[surf_key]: blk, var, slicer = self.surf_mesh_to_fort_var[surf_key][mesh_key] @@ -3638,6 +3772,49 @@ def set_function_ad_seeds(self, func_seeds: Dict[str, float], scale=1.0): val = func_seeds[_var] * scale self.set_avl_fort_arr(blk, var, val) + def get_strip_function_ad_seeds(self): + strip_func_seeds = {} + for surf_key in self.surface_names: + strip_func_seeds[surf_key] = {} + idx_surf = self.get_surface_index(surf_key) + idx_srp_beg, idx_srp_end = self._get_surface_strip_indices(idx_surf) + for _var in self.case_strip_var_to_fort_var: + vardata = self.case_strip_var_to_fort_var[_var] + if len(vardata) == 2: + blk = vardata[0] + var = vardata[1] + slicer = None + elif len(vardata) == 3: + blk = vardata[0] + var = vardata[1] + slicer = vardata[2] + blk += self.ad_suffix + var += self.ad_suffix + val = self.get_avl_fort_arr(blk, var, slicer) + strip_func_seeds[surf_key][_var] = copy.deepcopy(val[idx_srp_beg:idx_srp_end]) + + return strip_func_seeds + + def set_strip_function_ad_seeds(self, strip_func_seeds: Dict[str, Dict[str, float]], scale=1.0): + for surf_key in strip_func_seeds: + idx_surf = self.get_surface_index(surf_key) + idx_srp_beg, idx_srp_end = self._get_surface_strip_indices(idx_surf) + for _var in strip_func_seeds[surf_key]: + vardata = self.case_strip_var_to_fort_var[_var] + if len(vardata) == 2: + blk = vardata[0] + var = vardata[1] + slicer = slice(idx_srp_beg,idx_srp_end) + elif len(vardata) == 3: + blk = vardata[0] + var = vardata[1] + slicer = vardata[2] + slicer = (slice(idx_srp_beg,idx_srp_end),slicer[1]) + blk += self.ad_suffix + var += self.ad_suffix + val = strip_func_seeds[surf_key][_var] * scale + self.set_avl_fort_arr(blk, var, val,slicer) + def get_consurf_derivs_ad_seeds(self): cs_deriv_seeds = {} consurf_names = self.get_control_names() @@ -3828,6 +4005,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 +4020,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 +4054,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 +4081,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() @@ -3913,6 +4120,7 @@ def _execute_jac_vec_prod_fwd( # extract derivatives seeds and set the output dict of functions func_seeds = self.get_function_ad_seeds() + strip_func_seeds = self.get_strip_function_ad_seeds() res_seeds = self.get_residual_ad_seeds() consurf_derivs_seeds = self.get_consurf_derivs_ad_seeds() stab_derivs_seeds = self.get_stab_derivs_ad_seeds() @@ -3942,6 +4150,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() @@ -3949,6 +4188,7 @@ def _execute_jac_vec_prod_fwd( self.avl.aero() coef_data_peturb = self.get_total_forces() + strip_forces_peturb = self.get_strip_forces() consurf_derivs_petrub = self.get_control_stab_derivs() stab_deriv_petrub = self.get_stab_derivs() body_axis_deriv_petrub = self.get_body_axis_derivs() @@ -3966,12 +4206,43 @@ 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() self.avl.aero() coef_data = self.get_total_forces() + strip_forces = self.get_strip_forces() consurf_derivs = self.get_control_stab_derivs() stab_deriv = self.get_stab_derivs() body_axis_deriv = self.get_body_axis_derivs() @@ -3984,6 +4255,14 @@ def _execute_jac_vec_prod_fwd( for func_key in coef_data: func_seeds[func_key] = (coef_data_peturb[func_key] - coef_data[func_key]) / step + strip_func_seeds = {} + for surface in self.surface_names: + strip_func_seeds[surface] = {} + for strip_func_key in strip_forces[surface]: + if strip_func_key not in self.case_strip_var_to_fort_var: + continue + strip_func_seeds[surface][strip_func_key] = (strip_forces_peturb[surface][strip_func_key] - strip_forces[surface][strip_func_key]) / step + consurf_derivs_seeds = {} for deriv_func in consurf_derivs: consurf_derivs_seeds[deriv_func] = ( @@ -4007,6 +4286,7 @@ def _execute_jac_vec_prod_fwd( # TODO-clean: the way these arrays are returned is a bit of a mess return ( func_seeds, + strip_func_seeds, res_seeds, consurf_derivs_seeds, stab_derivs_seeds, @@ -4018,6 +4298,7 @@ def _execute_jac_vec_prod_fwd( def _execute_jac_vec_prod_rev( self, func_seeds: Optional[Dict[str, float]] = None, + strip_func_seeds: Optional[Dict[str, Dict[str, any]]] = None, res_seeds: Optional[np.ndarray] = None, consurf_derivs_seeds: Optional[Dict[str, float]] = None, stab_derivs_seeds: Optional[Dict[str, float]] = None, @@ -4066,6 +4347,9 @@ def _execute_jac_vec_prod_rev( if func_seeds is None: func_seeds = {} + if strip_func_seeds is None: + strip_func_seeds = {} + if res_seeds is None: res_seeds = np.zeros(mesh_size) @@ -4088,6 +4372,7 @@ def _execute_jac_vec_prod_rev( # self.clear_ad_seeds() time_last = time.time() self.set_function_ad_seeds(func_seeds) + self.set_strip_function_ad_seeds(strip_func_seeds) self.set_residual_ad_seeds(res_seeds) self.set_residual_d_ad_seeds(res_d_seeds) self.set_residual_u_ad_seeds(res_u_seeds) @@ -4130,7 +4415,32 @@ 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_strip_function_ad_seeds(strip_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) self.set_residual_u_ad_seeds(res_u_seeds, scale=0.0) @@ -4144,11 +4454,12 @@ 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, funcs: List[str], + strip_funcs: Optional[List[str]] = None, stab_derivs: Optional[List[str]] = None, body_axis_derivs: Optional[List[str]] = None, consurf_derivs: Optional[List[str]] = None, @@ -4178,7 +4489,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 +4506,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,12 +4516,56 @@ 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) sens[func].update(ref_seeds) + if strip_funcs is not None: + for surf in strip_funcs: + sens[surf] = {} + for func in strip_funcs[surf]: + + sens[surf][func] = {} + # get the RHS of the adjoint equation (pFpU) + # TODO: remove seeds if it doesn't effect accuracy + # self.clear_ad_seeds() + time_last = time.time() + _, _, _, _, pfpU, _, _, _, _ = self._execute_jac_vec_prod_rev(strip_func_seeds={surf:{func: 1.0}}) + if print_timings: + print(f"Time to get RHS: {time.time() - time_last}") + time_last = time.time() + + # self.clear_ad_seeds() + # u solver adjoint equation with RHS + self.set_gamma_ad_seeds(-1 * pfpU) + solve_gamma_u_adj = False + solve_gamma_d_adj = False + self.avl.solve_adjoint(solve_gamma_u_adj, solve_gamma_d_adj) + if print_timings: + print(f"Time to solve adjoint: {time.time() - time_last}") + time_last = time.time() + # 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, dvgeo_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( + strip_func_seeds={surf:{func: 1.0}}, res_seeds=dfdR + ) + if print_timings: + print(f"Time to combine derivs: {time.time() - time_last}") + time_last = time.time() + + sens[surf][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[surf][func][key] = geom_seeds[key] | mesh_seeds[key] | dvgeo_seeds[key] + # sens[func].update(geom_seeds) + # sens[func].update(mesh_seeds) + sens[surf][func].update(param_seeds) + sens[surf][func].update(ref_seeds) + + if consurf_derivs is not None: if print_timings: print("Running consurf derivs") @@ -4222,7 +4577,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 +4597,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 +4606,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 +4623,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 +4643,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 +4652,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 +4671,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 +4691,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 +4700,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 00000000..e3cdc7e0 --- /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/src/ad_src/Makefile_tapenade b/src/ad_src/Makefile_tapenade index 908cfead..cfa39794 100644 --- a/src/ad_src/Makefile_tapenade +++ b/src/ad_src/Makefile_tapenade @@ -55,6 +55,10 @@ fullRoutines = "\ CRSAX_D, CMSAX_D, CNSAX_D,\ CXTOT_U_BA, CYTOT_U_BA, CZTOT_U_BA,\ CRTOT_U_BA, CMTOT_U_BA, CNTOT_U_BA\ + CLSTRP, CDSTRP, CDV_LSTRP, DWWAKE\ + CFSTRP, CMSTRP, CL_LSTRP, CD_LSTRP\ + CF_LSTRP, CF_LSTRP, CLT_LSTRP, CMC4_LSTRP\ + CMLE_LSTRP, CNC\ ) \ " diff --git a/src/ad_src/forward_ad_src/aero_d.f b/src/ad_src/forward_ad_src/aero_d.f index 718de70a..72bb28f6 100644 --- a/src/ad_src/forward_ad_src/aero_d.f +++ b/src/ad_src/forward_ad_src/aero_d.f @@ -14,12 +14,14 @@ C cytot_rx crtot_rx cmtot_rx cntot_rx cdtot_ry cltot_ry C cytot_ry crtot_ry cmtot_ry cntot_ry cdtot_rz cltot_rz C cytot_rz crtot_rz cmtot_rz cntot_rz xnp sm bb -C rr +C rr cdstrp clstrp cfstrp cmstrp cf_lstrp cd_lstrp +C cl_lstrp cdv_lstrp clt_lstrp cmc4_lstrp cmle_lstrp +C cnc dwwake C with respect to varying inputs: alfa vinf vinf_a vinf_b wrot C sref cref bref xyzref mach cdref rle chord rle1 -C chord1 rle2 chord2 wstrip ensy ensz xsref ysref -C zsref rv1 rv2 rv rc gam gam_u gam_d src src_u -C vv vv_u vv_d wv wv_u wv_d +C chord1 rle2 chord2 wstrip ess ensy ensz xsref +C ysref zsref rv1 rv2 rv rc gam gam_u gam_d src +C src_u vv vv_u vv_d wv wv_u wv_d C RW status of diff variables: alfa:in vinf:in vinf_a:in vinf_b:in C wrot:in sref:in cref:in bref:in xyzref:in mach:in C cdref:in clff:out cyff:out cdff:out spanef:out @@ -40,9 +42,13 @@ C cltot_rz:out cytot_rz:out crtot_rz:out cmtot_rz:out C cntot_rz:out xnp:out sm:out bb:out rr:out rle:in C chord:in rle1:in chord1:in rle2:in chord2:in wstrip:in -C ensy:in ensz:in xsref:in ysref:in zsref:in rv1:in -C rv2:in rv:in rc:in gam:in gam_u:in gam_d:in src:in -C src_u:in vv:in vv_u:in vv_d:in wv:in wv_u:in wv_d:in +C ess:in ensy:in ensz:in xsref:in ysref:in zsref:in +C cdstrp:out clstrp:out cfstrp:out cmstrp:out cf_lstrp:out +C cd_lstrp:out cl_lstrp:out cdv_lstrp:out clt_lstrp:out +C cmc4_lstrp:out cmle_lstrp:out cnc:out dwwake:out +C rv1:in rv2:in rv:in rc:in gam:in gam_u:in gam_d:in +C src:in src_u:in vv:in vv_u:in vv_d:in wv:in wv_u:in +C wv_d:in C*********************************************************************** C Module: aero.f C @@ -333,10 +339,13 @@ SUBROUTINE AERO_D() C variations of useful results: cdtot cytot cltot cdtot_a cltot_a C cdtot_u cytot_u cltot_u cdtot_d cytot_d cltot_d C cftot cftot_u cftot_d cmtot cmtot_u cmtot_d cdvtot +C cdstrp clstrp cfstrp cmstrp cf_lstrp cd_lstrp +C cl_lstrp cdv_lstrp clt_lstrp cmc4_lstrp cmle_lstrp +C cnc C with respect to varying inputs: alfa vinf wrot sref cref bref C xyzref rle chord rle1 chord1 rle2 chord2 wstrip -C ensy ensz xsref ysref zsref rv1 rv2 rv gam gam_u -C gam_d vv vv_u vv_d wv wv_u wv_d +C ess ensy ensz xsref ysref zsref rv1 rv2 rv gam +C gam_u gam_d vv vv_u vv_d wv wv_u wv_d C AERO C C @@ -354,6 +363,7 @@ SUBROUTINE SFFORC_D() REAL vrot(3), vrot_u(3), wrot_u(3) REAL vrot_diff(3), vrot_u_diff(3), wrot_u_diff(3) REAL vperp(3) + REAL vperp_diff(3) REAL g(3), r(3), rh(3), mh(3) REAL g_diff(3), r_diff(3) REAL f(3), f_u(3, 6), f_d(3, ndmax), f_g(3, ngmax) @@ -469,12 +479,19 @@ SUBROUTINE SFFORC_D() REAL vsq REAL vsqi REAL vspan + REAL vspan_diff REAL vpsq + REAL vpsq_diff REAL vpsqi + REAL vpsqi_diff REAL delx + REAL delx_diff REAL dely + REAL dely_diff REAL delz + REAL delz_diff REAL dmag + REAL dmag_diff INTEGER is INTEGER jj REAL dcm @@ -487,8 +504,12 @@ SUBROUTINE SFFORC_D() REAL arg1_diff REAL temp INTEGER ii1 - REAL temp0 - REAL(kind=8) temp1 + REAL(kind=8) temp0 + REAL temp1 + REAL temp2 + REAL temp3 + REAL temp4 + REAL temp5 INTEGER ii2 INTEGER ii3 DATA icrs /2, 3, 1/ @@ -582,9 +603,32 @@ SUBROUTINE SFFORC_D() ENDDO ENDDO ENDDO + DO ii1=1,NSTRIP + DO ii2=1,3 + cf_lstrp_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO + DO ii1=1,NSTRIP + cd_lstrp_diff(ii1) = 0.D0 + ENDDO + DO ii1=1,NSTRIP + cl_lstrp_diff(ii1) = 0.D0 + ENDDO DO ii1=1,NSTRIP cdv_lstrp_diff(ii1) = 0.D0 ENDDO + DO ii1=1,NSTRIP + clt_lstrp_diff(ii1) = 0.D0 + ENDDO + DO ii1=1,NSTRIP + cmc4_lstrp_diff(ii1) = 0.D0 + ENDDO + DO ii1=1,NSTRIP + cmle_lstrp_diff(ii1) = 0.D0 + ENDDO + DO ii1=1,NSTRIP + cnc_diff(ii1) = 0.D0 + ENDDO DO ii1=1,numax cmz_u_diff(ii1) = 0.D0 ENDDO @@ -668,6 +712,9 @@ SUBROUTINE SFFORC_D() DO ii1=1,3 rc4_diff(ii1) = 0.D0 ENDDO + DO ii1=1,3 + vperp_diff(ii1) = 0.D0 + ENDDO DO ii1=1,numax veffmag_u_diff(ii1) = 0.D0 ENDDO @@ -829,6 +876,7 @@ SUBROUTINE SFFORC_D() cmx = 0. cmy = 0. cmz = 0. + cnc_diff(j) = 0.D0 cnc(j) = 0. C DO n=1,numax @@ -1129,7 +1177,10 @@ SUBROUTINE SFFORC_D() cmz = cmz + temp C C-------- accumulate strip spanloading = c*CN - cnc(j) = cnc(j) + cr*(ensy(j)*dcfy+ensz(j)*dcfz) + temp0 = ensy(j)*dcfy + ensz(j)*dcfz + cnc_diff(j) = cnc_diff(j) + temp0*cr_diff + cr*(dcfy*ensy_diff + + (j)+ensy(j)*dcfy_diff+dcfz*ensz_diff(j)+ensz(j)*dcfz_diff) + cnc(j) = cnc(j) + cr*temp0 C C-------- freestream and rotation derivatives DO n=1,numax @@ -1464,7 +1515,11 @@ SUBROUTINE SFFORC_D() Ccc write(23,*) I,ddcmx C C---------- accumulate strip spanloading = c*CN - cnc(j) = cnc(j) + cr*(ensy(j)*dcfy+ensz(j)*dcfz) + temp0 = ensy(j)*dcfy + ensz(j)*dcfz + cnc_diff(j) = cnc_diff(j) + temp0*cr_diff + cr*(dcfy* + + ensy_diff(j)+ensy(j)*dcfy_diff+dcfz*ensz_diff(j)+ensz(j) + + *dcfz_diff) + cnc(j) = cnc(j) + cr*temp0 C C---------- freestream and rotation derivatives DO n=1,numax @@ -1749,29 +1804,29 @@ SUBROUTINE SFFORC_D() C DO n=1,numax temp = veff_u(1, n)*veffmag + veff(1)*veffmag_u(n) - temp0 = veff(1)*clv_u(n) + temp1 = veff(1)*clv_u(n) dcvfx_u_diff = cdv*(veffmag*veff_u_diff(1, n)+veff_u(1, n)* + veffmag_diff+veffmag_u(n)*veff_diff(1)+veff(1)* - + veffmag_u_diff(n)) + temp*cdv_diff + temp0*(cdv_clv* + + veffmag_u_diff(n)) + temp*cdv_diff + temp1*(cdv_clv* + veffmag_diff+veffmag*cdv_clv_diff) + veffmag*cdv_clv*( + clv_u(n)*veff_diff(1)+veff(1)*clv_u_diff(n)) - dcvfx_u = temp*cdv + veffmag*cdv_clv*temp0 - temp0 = veff_u(2, n)*veffmag + veff(2)*veffmag_u(n) + dcvfx_u = temp*cdv + veffmag*cdv_clv*temp1 + temp1 = veff_u(2, n)*veffmag + veff(2)*veffmag_u(n) temp = veff(2)*clv_u(n) dcvfy_u_diff = cdv*(veffmag*veff_u_diff(2, n)+veff_u(2, n)* + veffmag_diff+veffmag_u(n)*veff_diff(2)+veff(2)* - + veffmag_u_diff(n)) + temp0*cdv_diff + temp*(cdv_clv* + + veffmag_u_diff(n)) + temp1*cdv_diff + temp*(cdv_clv* + veffmag_diff+veffmag*cdv_clv_diff) + veffmag*cdv_clv*( + clv_u(n)*veff_diff(2)+veff(2)*clv_u_diff(n)) - dcvfy_u = temp0*cdv + veffmag*cdv_clv*temp - temp0 = veff_u(3, n)*veffmag + veff(3)*veffmag_u(n) + dcvfy_u = temp1*cdv + veffmag*cdv_clv*temp + temp1 = veff_u(3, n)*veffmag + veff(3)*veffmag_u(n) temp = veff(3)*clv_u(n) dcvfz_u_diff = cdv*(veffmag*veff_u_diff(3, n)+veff_u(3, n)* + veffmag_diff+veffmag_u(n)*veff_diff(3)+veff(3)* - + veffmag_u_diff(n)) + temp0*cdv_diff + temp*(cdv_clv* + + veffmag_u_diff(n)) + temp1*cdv_diff + temp*(cdv_clv* + veffmag_diff+veffmag*cdv_clv_diff) + veffmag*cdv_clv*( + clv_u(n)*veff_diff(3)+veff(3)*clv_u_diff(n)) - dcvfz_u = temp0*cdv + veffmag*cdv_clv*temp + dcvfz_u = temp1*cdv + veffmag*cdv_clv*temp C cfx_u_diff(n) = cfx_u_diff(n) + dcvfx_u_diff cfx_u(n) = cfx_u(n) + dcvfx_u @@ -1789,21 +1844,21 @@ SUBROUTINE SFFORC_D() ENDDO C DO n=1,ncontrol - temp0 = veff(1)*clv_d(n) - dcvfx_d_diff = temp0*(cdv_clv*veffmag_diff+veffmag* + temp1 = veff(1)*clv_d(n) + dcvfx_d_diff = temp1*(cdv_clv*veffmag_diff+veffmag* + cdv_clv_diff) + veffmag*cdv_clv*(clv_d(n)*veff_diff(1)+ + veff(1)*clv_d_diff(n)) - dcvfx_d = veffmag*cdv_clv*temp0 - temp0 = veff(2)*clv_d(n) - dcvfy_d_diff = temp0*(cdv_clv*veffmag_diff+veffmag* + dcvfx_d = veffmag*cdv_clv*temp1 + temp1 = veff(2)*clv_d(n) + dcvfy_d_diff = temp1*(cdv_clv*veffmag_diff+veffmag* + cdv_clv_diff) + veffmag*cdv_clv*(clv_d(n)*veff_diff(2)+ + veff(2)*clv_d_diff(n)) - dcvfy_d = veffmag*cdv_clv*temp0 - temp0 = veff(3)*clv_d(n) - dcvfz_d_diff = temp0*(cdv_clv*veffmag_diff+veffmag* + dcvfy_d = veffmag*cdv_clv*temp1 + temp1 = veff(3)*clv_d(n) + dcvfz_d_diff = temp1*(cdv_clv*veffmag_diff+veffmag* + cdv_clv_diff) + veffmag*cdv_clv*(clv_d(n)*veff_diff(3)+ + veff(3)*clv_d_diff(n)) - dcvfz_d = veffmag*cdv_clv*temp0 + dcvfz_d = veffmag*cdv_clv*temp1 C cfx_d_diff(n) = cfx_d_diff(n) + dcvfx_d_diff cfx_d(n) = cfx_d(n) + dcvfx_d @@ -1850,8 +1905,11 @@ SUBROUTINE SFFORC_D() C C...Store strip X,Y,Z body axes forces and moments C referred to c/4 and strip area and chord + cf_lstrp_diff(1, j) = cfx_diff cf_lstrp(1, j) = cfx + cf_lstrp_diff(2, j) = cfy_diff cf_lstrp(2, j) = cfy + cf_lstrp_diff(3, j) = cfz_diff cf_lstrp(3, j) = cfz cm_lstrp(1, j) = cmx cm_lstrp(2, j) = cmy @@ -1936,53 +1994,53 @@ SUBROUTINE SFFORC_D() r(3) = rc4(3) - xyzref(3) C... Strip moments in body axes about the case moment reference point XYZREF C normalized by strip area and chord - temp0 = (cfz*r(2)-cfy*r(3))/cr + temp1 = (cfz*r(2)-cfy*r(3))/cr cmstrp_diff(1, j) = cmx_diff + (r(2)*cfz_diff+cfz*r_diff(2)-r(3) - + *cfy_diff-cfy*r_diff(3)-temp0*cr_diff)/cr - cmstrp(1, j) = cmx + temp0 - temp0 = (cfx*r(3)-cfz*r(1))/cr + + *cfy_diff-cfy*r_diff(3)-temp1*cr_diff)/cr + cmstrp(1, j) = cmx + temp1 + temp1 = (cfx*r(3)-cfz*r(1))/cr cmstrp_diff(2, j) = cmy_diff + (r(3)*cfx_diff+cfx*r_diff(3)-r(1) - + *cfz_diff-cfz*r_diff(1)-temp0*cr_diff)/cr - cmstrp(2, j) = cmy + temp0 - temp0 = (cfy*r(1)-cfx*r(2))/cr + + *cfz_diff-cfz*r_diff(1)-temp1*cr_diff)/cr + cmstrp(2, j) = cmy + temp1 + temp1 = (cfy*r(1)-cfx*r(2))/cr cmstrp_diff(3, j) = cmz_diff + (r(1)*cfy_diff+cfy*r_diff(1)-r(2) - + *cfx_diff-cfx*r_diff(2)-temp0*cr_diff)/cr - cmstrp(3, j) = cmz + temp0 + + *cfx_diff-cfx*r_diff(2)-temp1*cr_diff)/cr + cmstrp(3, j) = cmz + temp1 C DO n=1,numax - temp0 = (cfz_u(n)*r(2)-cfy_u(n)*r(3))/cr + temp1 = (cfz_u(n)*r(2)-cfy_u(n)*r(3))/cr cmst_u_diff(1, j, n) = cmx_u_diff(n) + (r(2)*cfz_u_diff(n)+ + cfz_u(n)*r_diff(2)-r(3)*cfy_u_diff(n)-cfy_u(n)*r_diff(3)- - + temp0*cr_diff)/cr - cmst_u(1, j, n) = cmx_u(n) + temp0 - temp0 = (cfx_u(n)*r(3)-cfz_u(n)*r(1))/cr + + temp1*cr_diff)/cr + cmst_u(1, j, n) = cmx_u(n) + temp1 + temp1 = (cfx_u(n)*r(3)-cfz_u(n)*r(1))/cr cmst_u_diff(2, j, n) = cmy_u_diff(n) + (r(3)*cfx_u_diff(n)+ + cfx_u(n)*r_diff(3)-r(1)*cfz_u_diff(n)-cfz_u(n)*r_diff(1)- - + temp0*cr_diff)/cr - cmst_u(2, j, n) = cmy_u(n) + temp0 - temp0 = (cfy_u(n)*r(1)-cfx_u(n)*r(2))/cr + + temp1*cr_diff)/cr + cmst_u(2, j, n) = cmy_u(n) + temp1 + temp1 = (cfy_u(n)*r(1)-cfx_u(n)*r(2))/cr cmst_u_diff(3, j, n) = cmz_u_diff(n) + (r(1)*cfy_u_diff(n)+ + cfy_u(n)*r_diff(1)-r(2)*cfx_u_diff(n)-cfx_u(n)*r_diff(2)- - + temp0*cr_diff)/cr - cmst_u(3, j, n) = cmz_u(n) + temp0 + + temp1*cr_diff)/cr + cmst_u(3, j, n) = cmz_u(n) + temp1 ENDDO C DO n=1,ncontrol - temp0 = (cfz_d(n)*r(2)-cfy_d(n)*r(3))/cr + temp1 = (cfz_d(n)*r(2)-cfy_d(n)*r(3))/cr cmst_d_diff(1, j, n) = cmx_d_diff(n) + (r(2)*cfz_d_diff(n)+ + cfz_d(n)*r_diff(2)-r(3)*cfy_d_diff(n)-cfy_d(n)*r_diff(3)- - + temp0*cr_diff)/cr - cmst_d(1, j, n) = cmx_d(n) + temp0 - temp0 = (cfx_d(n)*r(3)-cfz_d(n)*r(1))/cr + + temp1*cr_diff)/cr + cmst_d(1, j, n) = cmx_d(n) + temp1 + temp1 = (cfx_d(n)*r(3)-cfz_d(n)*r(1))/cr cmst_d_diff(2, j, n) = cmy_d_diff(n) + (r(3)*cfx_d_diff(n)+ + cfx_d(n)*r_diff(3)-r(1)*cfz_d_diff(n)-cfz_d(n)*r_diff(1)- - + temp0*cr_diff)/cr - cmst_d(2, j, n) = cmy_d(n) + temp0 - temp0 = (cfy_d(n)*r(1)-cfx_d(n)*r(2))/cr + + temp1*cr_diff)/cr + cmst_d(2, j, n) = cmy_d(n) + temp1 + temp1 = (cfy_d(n)*r(1)-cfx_d(n)*r(2))/cr cmst_d_diff(3, j, n) = cmz_d_diff(n) + (r(1)*cfy_d_diff(n)+ + cfy_d(n)*r_diff(1)-r(2)*cfx_d_diff(n)-cfx_d(n)*r_diff(2)- - + temp0*cr_diff)/cr - cmst_d(3, j, n) = cmz_d(n) + temp0 + + temp1*cr_diff)/cr + cmst_d(3, j, n) = cmz_d(n) + temp1 ENDDO C DO n=1,ndesign @@ -1992,8 +2050,16 @@ SUBROUTINE SFFORC_D() ENDDO C C...Components of X,Y,Z forces in local strip axes + cl_lstrp_diff(j) = cfx*ulift_diff(1) + ulift(1)*cfx_diff + cfy* + + ulift_diff(2) + ulift(2)*cfy_diff + cfz*ulift_diff(3) + ulift( + + 3)*cfz_diff cl_lstrp(j) = ulift(1)*cfx + ulift(2)*cfy + ulift(3)*cfz + cd_lstrp_diff(j) = cfx*udrag_diff(1) + udrag(1)*cfx_diff + cfy* + + udrag_diff(2) + udrag(2)*cfy_diff + cfz*udrag_diff(3) + udrag( + + 3)*cfz_diff cd_lstrp(j) = udrag(1)*cfx + udrag(2)*cfy + udrag(3)*cfz + cmc4_lstrp_diff(j) = cmy*ensz_diff(j) + ensz(j)*cmy_diff - cmz* + + ensy_diff(j) - ensy(j)*cmz_diff cmc4_lstrp(j) = ensz(j)*cmy - ensy(j)*cmz C C...Axial/normal forces and lift/drag in plane normal to dihedral of strip @@ -2032,19 +2098,34 @@ SUBROUTINE SFFORC_D() END IF C C------ spanwise and perpendicular velocity components + vspan_diff = ess(1, j)*veff_diff(1) + veff(1)*ess_diff(1, j) + + + ess(2, j)*veff_diff(2) + veff(2)*ess_diff(2, j) + ess(3, j)* + + veff_diff(3) + veff(3)*ess_diff(3, j) vspan = veff(1)*ess(1, j) + veff(2)*ess(2, j) + veff(3)*ess(3, j + ) + vperp_diff(1) = veff_diff(1) - vspan*ess_diff(1, j) - ess(1, j)* + + vspan_diff vperp(1) = veff(1) - ess(1, j)*vspan + vperp_diff(2) = veff_diff(2) - vspan*ess_diff(2, j) - ess(2, j)* + + vspan_diff vperp(2) = veff(2) - ess(2, j)*vspan + vperp_diff(3) = veff_diff(3) - vspan*ess_diff(3, j) - ess(3, j)* + + vspan_diff vperp(3) = veff(3) - ess(3, j)*vspan C + vpsq_diff = 2*vperp(1)*vperp_diff(1) + 2*vperp(2)*vperp_diff(2) + + + 2*vperp(3)*vperp_diff(3) vpsq = vperp(1)**2 + vperp(2)**2 + vperp(3)**2 IF (vpsq .EQ. 0.0) THEN vpsqi = 1.0 + vpsqi_diff = 0.D0 ELSE + vpsqi_diff = -(vpsq_diff/vpsq**2) vpsqi = 1.0/vpsq END IF Ccc CLT_LSTRP(J) = CN_LSTRP(J) * VPSQI + clt_lstrp_diff(j) = vpsqi*cl_lstrp_diff(j) + cl_lstrp(j)* + + vpsqi_diff clt_lstrp(j) = cl_lstrp(j)*vpsqi cla_lstrp(j) = cl_lstrp(j)*vsqi C @@ -2055,21 +2136,51 @@ SUBROUTINE SFFORC_D() r(2) = rc4(2) - rle(2, j) r_diff(3) = rc4_diff(3) - rle_diff(3, j) r(3) = rc4(3) - rle(3, j) + delx_diff = rle2_diff(1, j) - rle1_diff(1, j) delx = rle2(1, j) - rle1(1, j) + dely_diff = rle2_diff(2, j) - rle1_diff(2, j) dely = rle2(2, j) - rle1(2, j) + delz_diff = rle2_diff(3, j) - rle1_diff(3, j) delz = rle2(3, j) - rle1(3, j) C IF (imags(lssurf(j)) .LT. 0) THEN + delx_diff = -delx_diff delx = -delx + dely_diff = -dely_diff dely = -dely + delz_diff = -delz_diff delz = -delz END IF + arg1_diff = 2*delx*delx_diff + 2*dely*dely_diff + 2*delz* + + delz_diff arg1 = delx**2 + dely**2 + delz**2 - dmag = SQRT(arg1) + temp1 = SQRT(arg1) + IF (arg1 .EQ. 0.D0) THEN + dmag_diff = 0.D0 + ELSE + dmag_diff = arg1_diff/(2.0*temp1) + END IF + dmag = temp1 + cmle_lstrp_diff(j) = 0.D0 cmle_lstrp(j) = 0.0 - IF (dmag .NE. 0.0) cmle_lstrp(j) = delx/dmag*(cmx+(cfz*r(2)-cfy* - + r(3))/cr) + dely/dmag*(cmy+(cfx*r(3)-cfz*r(1))/cr) + delz/ - + dmag*(cmz+(cfy*r(1)-cfx*r(2))/cr) + IF (dmag .NE. 0.0) THEN + temp1 = delx/dmag + temp = (cfz*r(2)-cfy*r(3))/cr + temp2 = dely/dmag + temp3 = (cfx*r(3)-cfz*r(1))/cr + temp4 = delz/dmag + temp5 = (cfy*r(1)-cfx*r(2))/cr + cmle_lstrp_diff(j) = temp1*(cmx_diff+(r(2)*cfz_diff+cfz*r_diff + + (2)-r(3)*cfy_diff-cfy*r_diff(3)-temp*cr_diff)/cr) + (cmx+ + + temp)*(delx_diff-temp1*dmag_diff)/dmag + temp2*(cmy_diff+(r( + + 3)*cfx_diff+cfx*r_diff(3)-r(1)*cfz_diff-cfz*r_diff(1)-temp3* + + cr_diff)/cr) + (cmy+temp3)*(dely_diff-temp2*dmag_diff)/dmag + + + temp4*(cmz_diff+(r(1)*cfy_diff+cfy*r_diff(1)-r(2)*cfx_diff + + -cfx*r_diff(2)-temp5*cr_diff)/cr) + (cmz+temp5)*(delz_diff- + + temp4*dmag_diff)/dmag + cmle_lstrp(j) = (cmx+temp)*temp1 + (cmy+temp3)*temp2 + (cmz+ + + temp5)*temp4 + END IF ENDDO cdtot_diff = 0.D0 cytot_diff = 0.D0 @@ -2298,145 +2409,145 @@ SUBROUTINE SFFORC_D() enave(3) = enave(3) + sr*ensz(j) C C--- Surface lift and drag referenced to case SREF, CREF, BREF - temp1 = sr/sref - cdsurf_diff(is) = cdsurf_diff(is) + temp1*cdstrp_diff(j) + - + cdstrp(j)*(sr_diff-temp1*sref_diff)/sref - cdsurf(is) = cdsurf(is) + cdstrp(j)*temp1 - temp1 = sr/sref - cysurf_diff(is) = cysurf_diff(is) + temp1*cystrp_diff(j) + - + cystrp(j)*(sr_diff-temp1*sref_diff)/sref - cysurf(is) = cysurf(is) + cystrp(j)*temp1 - temp1 = sr/sref - clsurf_diff(is) = clsurf_diff(is) + temp1*clstrp_diff(j) + - + clstrp(j)*(sr_diff-temp1*sref_diff)/sref - clsurf(is) = clsurf(is) + clstrp(j)*temp1 + temp0 = sr/sref + cdsurf_diff(is) = cdsurf_diff(is) + temp0*cdstrp_diff(j) + + + cdstrp(j)*(sr_diff-temp0*sref_diff)/sref + cdsurf(is) = cdsurf(is) + cdstrp(j)*temp0 + temp0 = sr/sref + cysurf_diff(is) = cysurf_diff(is) + temp0*cystrp_diff(j) + + + cystrp(j)*(sr_diff-temp0*sref_diff)/sref + cysurf(is) = cysurf(is) + cystrp(j)*temp0 + temp0 = sr/sref + clsurf_diff(is) = clsurf_diff(is) + temp0*clstrp_diff(j) + + + clstrp(j)*(sr_diff-temp0*sref_diff)/sref + clsurf(is) = clsurf(is) + clstrp(j)*temp0 C--- Surface body axes forces referenced to case SREF, CREF, BREF - temp1 = sr/sref - cfsurf_diff(1, is) = cfsurf_diff(1, is) + temp1*cfstrp_diff(1 - + , j) + cfstrp(1, j)*(sr_diff-temp1*sref_diff)/sref - cfsurf(1, is) = cfsurf(1, is) + cfstrp(1, j)*temp1 - temp1 = sr/sref - cfsurf_diff(2, is) = cfsurf_diff(2, is) + temp1*cfstrp_diff(2 - + , j) + cfstrp(2, j)*(sr_diff-temp1*sref_diff)/sref - cfsurf(2, is) = cfsurf(2, is) + cfstrp(2, j)*temp1 - temp1 = sr/sref - cfsurf_diff(3, is) = cfsurf_diff(3, is) + temp1*cfstrp_diff(3 - + , j) + cfstrp(3, j)*(sr_diff-temp1*sref_diff)/sref - cfsurf(3, is) = cfsurf(3, is) + cfstrp(3, j)*temp1 + temp0 = sr/sref + cfsurf_diff(1, is) = cfsurf_diff(1, is) + temp0*cfstrp_diff(1 + + , j) + cfstrp(1, j)*(sr_diff-temp0*sref_diff)/sref + cfsurf(1, is) = cfsurf(1, is) + cfstrp(1, j)*temp0 + temp0 = sr/sref + cfsurf_diff(2, is) = cfsurf_diff(2, is) + temp0*cfstrp_diff(2 + + , j) + cfstrp(2, j)*(sr_diff-temp0*sref_diff)/sref + cfsurf(2, is) = cfsurf(2, is) + cfstrp(2, j)*temp0 + temp0 = sr/sref + cfsurf_diff(3, is) = cfsurf_diff(3, is) + temp0*cfstrp_diff(3 + + , j) + cfstrp(3, j)*(sr_diff-temp0*sref_diff)/sref + cfsurf(3, is) = cfsurf(3, is) + cfstrp(3, j)*temp0 C--- Surface body axes moments referenced to case SREF, CREF, BREF about XYZREF - temp1 = cmstrp(1, j)*sr*cr/(sref*bref) + temp0 = cmstrp(1, j)*sr*cr/(sref*bref) cmsurf_diff(1, is) = cmsurf_diff(1, is) + (sr*cr*cmstrp_diff(1 - + , j)+cmstrp(1, j)*(cr*sr_diff+sr*cr_diff)-temp1*(bref* + + , j)+cmstrp(1, j)*(cr*sr_diff+sr*cr_diff)-temp0*(bref* + sref_diff+sref*bref_diff))/(sref*bref) - cmsurf(1, is) = cmsurf(1, is) + temp1 - temp1 = cmstrp(2, j)*sr*cr/(sref*cref) + cmsurf(1, is) = cmsurf(1, is) + temp0 + temp0 = cmstrp(2, j)*sr*cr/(sref*cref) cmsurf_diff(2, is) = cmsurf_diff(2, is) + (sr*cr*cmstrp_diff(2 - + , j)+cmstrp(2, j)*(cr*sr_diff+sr*cr_diff)-temp1*(cref* + + , j)+cmstrp(2, j)*(cr*sr_diff+sr*cr_diff)-temp0*(cref* + sref_diff+sref*cref_diff))/(sref*cref) - cmsurf(2, is) = cmsurf(2, is) + temp1 - temp1 = cmstrp(3, j)*sr*cr/(sref*bref) + cmsurf(2, is) = cmsurf(2, is) + temp0 + temp0 = cmstrp(3, j)*sr*cr/(sref*bref) cmsurf_diff(3, is) = cmsurf_diff(3, is) + (sr*cr*cmstrp_diff(3 - + , j)+cmstrp(3, j)*(cr*sr_diff+sr*cr_diff)-temp1*(bref* + + , j)+cmstrp(3, j)*(cr*sr_diff+sr*cr_diff)-temp0*(bref* + sref_diff+sref*bref_diff))/(sref*bref) - cmsurf(3, is) = cmsurf(3, is) + temp1 + cmsurf(3, is) = cmsurf(3, is) + temp0 C C--- Bug fix, HHY/S.Allmaras C--- Surface viscous drag referenced to case SREF, CREF, BREF - temp1 = sr/sref - cdvsurf_diff(is) = cdvsurf_diff(is) + temp1*cdv_lstrp_diff(j) - + + cdv_lstrp(j)*(sr_diff-temp1*sref_diff)/sref - cdvsurf(is) = cdvsurf(is) + cdv_lstrp(j)*temp1 -C - temp1 = sr/sref - cds_a_diff(is) = cds_a_diff(is) + temp1*cdst_a_diff(j) + - + cdst_a(j)*(sr_diff-temp1*sref_diff)/sref - cds_a(is) = cds_a(is) + cdst_a(j)*temp1 + temp0 = sr/sref + cdvsurf_diff(is) = cdvsurf_diff(is) + temp0*cdv_lstrp_diff(j) + + + cdv_lstrp(j)*(sr_diff-temp0*sref_diff)/sref + cdvsurf(is) = cdvsurf(is) + cdv_lstrp(j)*temp0 +C + temp0 = sr/sref + cds_a_diff(is) = cds_a_diff(is) + temp0*cdst_a_diff(j) + + + cdst_a(j)*(sr_diff-temp0*sref_diff)/sref + cds_a(is) = cds_a(is) + cdst_a(j)*temp0 cys_a(is) = cys_a(is) + cyst_a(j)*sr/sref - temp1 = sr/sref - cls_a_diff(is) = cls_a_diff(is) + temp1*clst_a_diff(j) + - + clst_a(j)*(sr_diff-temp1*sref_diff)/sref - cls_a(is) = cls_a(is) + clst_a(j)*temp1 + temp0 = sr/sref + cls_a_diff(is) = cls_a_diff(is) + temp0*clst_a_diff(j) + + + clst_a(j)*(sr_diff-temp0*sref_diff)/sref + cls_a(is) = cls_a(is) + clst_a(j)*temp0 C DO n=1,numax - temp1 = sr/sref - cds_u_diff(is, n) = cds_u_diff(is, n) + temp1*cdst_u_diff(j - + , n) + cdst_u(j, n)*(sr_diff-temp1*sref_diff)/sref - cds_u(is, n) = cds_u(is, n) + cdst_u(j, n)*temp1 - temp1 = sr/sref - cys_u_diff(is, n) = cys_u_diff(is, n) + temp1*cyst_u_diff(j - + , n) + cyst_u(j, n)*(sr_diff-temp1*sref_diff)/sref - cys_u(is, n) = cys_u(is, n) + cyst_u(j, n)*temp1 - temp1 = sr/sref - cls_u_diff(is, n) = cls_u_diff(is, n) + temp1*clst_u_diff(j - + , n) + clst_u(j, n)*(sr_diff-temp1*sref_diff)/sref - cls_u(is, n) = cls_u(is, n) + clst_u(j, n)*temp1 + temp0 = sr/sref + cds_u_diff(is, n) = cds_u_diff(is, n) + temp0*cdst_u_diff(j + + , n) + cdst_u(j, n)*(sr_diff-temp0*sref_diff)/sref + cds_u(is, n) = cds_u(is, n) + cdst_u(j, n)*temp0 + temp0 = sr/sref + cys_u_diff(is, n) = cys_u_diff(is, n) + temp0*cyst_u_diff(j + + , n) + cyst_u(j, n)*(sr_diff-temp0*sref_diff)/sref + cys_u(is, n) = cys_u(is, n) + cyst_u(j, n)*temp0 + temp0 = sr/sref + cls_u_diff(is, n) = cls_u_diff(is, n) + temp0*clst_u_diff(j + + , n) + clst_u(j, n)*(sr_diff-temp0*sref_diff)/sref + cls_u(is, n) = cls_u(is, n) + clst_u(j, n)*temp0 C DO l=1,3 - temp1 = sr/sref - cfs_u_diff(l, is, n) = cfs_u_diff(l, is, n) + temp1* - + cfst_u_diff(l, j, n) + cfst_u(l, j, n)*(sr_diff-temp1* + temp0 = sr/sref + cfs_u_diff(l, is, n) = cfs_u_diff(l, is, n) + temp0* + + cfst_u_diff(l, j, n) + cfst_u(l, j, n)*(sr_diff-temp0* + sref_diff)/sref - cfs_u(l, is, n) = cfs_u(l, is, n) + cfst_u(l, j, n)*temp1 + cfs_u(l, is, n) = cfs_u(l, is, n) + cfst_u(l, j, n)*temp0 ENDDO - temp1 = cmst_u(1, j, n)*sr*cr/(sref*bref) + temp0 = cmst_u(1, j, n)*sr*cr/(sref*bref) cms_u_diff(1, is, n) = cms_u_diff(1, is, n) + (sr*cr* + cmst_u_diff(1, j, n)+cmst_u(1, j, n)*(cr*sr_diff+sr* - + cr_diff)-temp1*(bref*sref_diff+sref*bref_diff))/(sref*bref + + cr_diff)-temp0*(bref*sref_diff+sref*bref_diff))/(sref*bref + ) - cms_u(1, is, n) = cms_u(1, is, n) + temp1 - temp1 = cmst_u(2, j, n)*sr*cr/(sref*cref) + cms_u(1, is, n) = cms_u(1, is, n) + temp0 + temp0 = cmst_u(2, j, n)*sr*cr/(sref*cref) cms_u_diff(2, is, n) = cms_u_diff(2, is, n) + (sr*cr* + cmst_u_diff(2, j, n)+cmst_u(2, j, n)*(cr*sr_diff+sr* - + cr_diff)-temp1*(cref*sref_diff+sref*cref_diff))/(sref*cref + + cr_diff)-temp0*(cref*sref_diff+sref*cref_diff))/(sref*cref + ) - cms_u(2, is, n) = cms_u(2, is, n) + temp1 - temp1 = cmst_u(3, j, n)*sr*cr/(sref*bref) + cms_u(2, is, n) = cms_u(2, is, n) + temp0 + temp0 = cmst_u(3, j, n)*sr*cr/(sref*bref) cms_u_diff(3, is, n) = cms_u_diff(3, is, n) + (sr*cr* + cmst_u_diff(3, j, n)+cmst_u(3, j, n)*(cr*sr_diff+sr* - + cr_diff)-temp1*(bref*sref_diff+sref*bref_diff))/(sref*bref + + cr_diff)-temp0*(bref*sref_diff+sref*bref_diff))/(sref*bref + ) - cms_u(3, is, n) = cms_u(3, is, n) + temp1 + cms_u(3, is, n) = cms_u(3, is, n) + temp0 ENDDO C DO n=1,ncontrol - temp1 = sr/sref - cds_d_diff(is, n) = cds_d_diff(is, n) + temp1*cdst_d_diff(j - + , n) + cdst_d(j, n)*(sr_diff-temp1*sref_diff)/sref - cds_d(is, n) = cds_d(is, n) + cdst_d(j, n)*temp1 - temp1 = sr/sref - cys_d_diff(is, n) = cys_d_diff(is, n) + temp1*cyst_d_diff(j - + , n) + cyst_d(j, n)*(sr_diff-temp1*sref_diff)/sref - cys_d(is, n) = cys_d(is, n) + cyst_d(j, n)*temp1 - temp1 = sr/sref - cls_d_diff(is, n) = cls_d_diff(is, n) + temp1*clst_d_diff(j - + , n) + clst_d(j, n)*(sr_diff-temp1*sref_diff)/sref - cls_d(is, n) = cls_d(is, n) + clst_d(j, n)*temp1 + temp0 = sr/sref + cds_d_diff(is, n) = cds_d_diff(is, n) + temp0*cdst_d_diff(j + + , n) + cdst_d(j, n)*(sr_diff-temp0*sref_diff)/sref + cds_d(is, n) = cds_d(is, n) + cdst_d(j, n)*temp0 + temp0 = sr/sref + cys_d_diff(is, n) = cys_d_diff(is, n) + temp0*cyst_d_diff(j + + , n) + cyst_d(j, n)*(sr_diff-temp0*sref_diff)/sref + cys_d(is, n) = cys_d(is, n) + cyst_d(j, n)*temp0 + temp0 = sr/sref + cls_d_diff(is, n) = cls_d_diff(is, n) + temp0*clst_d_diff(j + + , n) + clst_d(j, n)*(sr_diff-temp0*sref_diff)/sref + cls_d(is, n) = cls_d(is, n) + clst_d(j, n)*temp0 C DO l=1,3 - temp1 = sr/sref - cfs_d_diff(l, is, n) = cfs_d_diff(l, is, n) + temp1* - + cfst_d_diff(l, j, n) + cfst_d(l, j, n)*(sr_diff-temp1* + temp0 = sr/sref + cfs_d_diff(l, is, n) = cfs_d_diff(l, is, n) + temp0* + + cfst_d_diff(l, j, n) + cfst_d(l, j, n)*(sr_diff-temp0* + sref_diff)/sref - cfs_d(l, is, n) = cfs_d(l, is, n) + cfst_d(l, j, n)*temp1 + cfs_d(l, is, n) = cfs_d(l, is, n) + cfst_d(l, j, n)*temp0 ENDDO - temp1 = cmst_d(1, j, n)*sr*cr/(sref*bref) + temp0 = cmst_d(1, j, n)*sr*cr/(sref*bref) cms_d_diff(1, is, n) = cms_d_diff(1, is, n) + (sr*cr* + cmst_d_diff(1, j, n)+cmst_d(1, j, n)*(cr*sr_diff+sr* - + cr_diff)-temp1*(bref*sref_diff+sref*bref_diff))/(sref*bref + + cr_diff)-temp0*(bref*sref_diff+sref*bref_diff))/(sref*bref + ) - cms_d(1, is, n) = cms_d(1, is, n) + temp1 - temp1 = cmst_d(2, j, n)*sr*cr/(sref*cref) + cms_d(1, is, n) = cms_d(1, is, n) + temp0 + temp0 = cmst_d(2, j, n)*sr*cr/(sref*cref) cms_d_diff(2, is, n) = cms_d_diff(2, is, n) + (sr*cr* + cmst_d_diff(2, j, n)+cmst_d(2, j, n)*(cr*sr_diff+sr* - + cr_diff)-temp1*(cref*sref_diff+sref*cref_diff))/(sref*cref + + cr_diff)-temp0*(cref*sref_diff+sref*cref_diff))/(sref*cref + ) - cms_d(2, is, n) = cms_d(2, is, n) + temp1 - temp1 = cmst_d(3, j, n)*sr*cr/(sref*bref) + cms_d(2, is, n) = cms_d(2, is, n) + temp0 + temp0 = cmst_d(3, j, n)*sr*cr/(sref*bref) cms_d_diff(3, is, n) = cms_d_diff(3, is, n) + (sr*cr* + cmst_d_diff(3, j, n)+cmst_d(3, j, n)*(cr*sr_diff+sr* - + cr_diff)-temp1*(bref*sref_diff+sref*bref_diff))/(sref*bref + + cr_diff)-temp0*(bref*sref_diff+sref*bref_diff))/(sref*bref + ) - cms_d(3, is, n) = cms_d(3, is, n) + temp1 + cms_d(3, is, n) = cms_d(3, is, n) + temp0 ENDDO C DO n=1,ndesign diff --git a/src/ad_src/forward_ad_src/atpforc_d.f b/src/ad_src/forward_ad_src/atpforc_d.f index 54d74141..a98e6cea 100644 --- a/src/ad_src/forward_ad_src/atpforc_d.f +++ b/src/ad_src/forward_ad_src/atpforc_d.f @@ -2,7 +2,7 @@ C Tapenade 3.16 (develop) - 15 Jan 2021 14:26 C C Differentiation of tpforc in forward (tangent) mode (with options i4 dr8 r8): -C variations of useful results: clff cyff cdff spanef +C variations of useful results: clff cyff cdff spanef dwwake C with respect to varying inputs: sref bref chord rv1 rv2 rc C gam C*********************************************************************** @@ -30,6 +30,7 @@ SUBROUTINE TPFORC_D() INCLUDE 'AVL_ad_seeds.inc' C REAL ny, nz + REAL ny_diff, nz_diff REAL vy_u(numax), vz_u(numax), vy_d(ndmax), vz_d(ndmax), vy_g( + ngmax), vz_g(ngmax) REAL p(3, 3), p_m(3, 3), p_a(3, 3), p_b(3, 3) @@ -55,6 +56,7 @@ SUBROUTINE TPFORC_D() REAL dzt REAL dzt_diff REAL dst + REAL dst_diff INTRINSIC SQRT REAL ycntr REAL ycntr_diff @@ -199,6 +201,9 @@ SUBROUTINE TPFORC_D() clff_diff = 0.D0 cyff_diff = 0.D0 cdff_diff = 0.D0 + DO ii1=1,NSTRIP + dwwake_diff(ii1) = 0.D0 + ENDDO C C...Find the normal velocity across each strip at the projected control C point location @@ -208,10 +213,19 @@ SUBROUTINE TPFORC_D() dyt = rt2(2, jc) - rt1(2, jc) dzt_diff = rt2_diff(3, jc) - rt1_diff(3, jc) dzt = rt2(3, jc) - rt1(3, jc) + arg1_diff = 2*dyt*dyt_diff + 2*dzt*dzt_diff arg1 = dyt*dyt + dzt*dzt - dst = SQRT(arg1) + temp = SQRT(arg1) + IF (arg1 .EQ. 0.D0) THEN + dst_diff = 0.D0 + ELSE + dst_diff = arg1_diff/(2.0*temp) + END IF + dst = temp C + ny_diff = -((dzt_diff-dzt*dst_diff/dst)/dst) ny = -(dzt/dst) + nz_diff = (dyt_diff-dyt*dst_diff/dst)/dst nz = dyt/dst ycntr_diff = rtc_diff(2, jc) ycntr = rtc(2, jc) @@ -450,6 +464,8 @@ SUBROUTINE TPFORC_D() ENDDO C C + dwwake_diff(jc) = -(vy*ny_diff) - ny*vy_diff - vz*nz_diff - nz* + + vz_diff dwwake(jc) = -(ny*vy+nz*vz) C C...Trefftz-plane drag is kinetic energy in crossflow diff --git a/src/ad_src/reverse_ad_src/aero_b.f b/src/ad_src/reverse_ad_src/aero_b.f index d9f76a25..46a3bc51 100644 --- a/src/ad_src/reverse_ad_src/aero_b.f +++ b/src/ad_src/reverse_ad_src/aero_b.f @@ -14,7 +14,9 @@ C cytot_rx crtot_rx cmtot_rx cntot_rx cdtot_ry cltot_ry C cytot_ry crtot_ry cmtot_ry cntot_ry cdtot_rz cltot_rz C cytot_rz crtot_rz cmtot_rz cntot_rz xnp sm bb -C rr +C rr cdstrp clstrp cfstrp cmstrp cf_lstrp cd_lstrp +C cl_lstrp cdv_lstrp clt_lstrp cmc4_lstrp cmle_lstrp +C cnc dwwake C with respect to varying inputs: alfa vinf vinf_a vinf_b wrot C sref cref bref xyzref mach cdref clff cyff cdff C spanef cdtot cytot cltot cdtot_d cytot_d cltot_d @@ -28,9 +30,12 @@ C cytot_rx crtot_rx cmtot_rx cntot_rx cdtot_ry cltot_ry C cytot_ry crtot_ry cmtot_ry cntot_ry cdtot_rz cltot_rz C cytot_rz crtot_rz cmtot_rz cntot_rz xnp sm bb -C rr rle chord rle1 chord1 rle2 chord2 wstrip ensy -C ensz xsref ysref zsref rv1 rv2 rv rc gam gam_u -C gam_d src src_u vv vv_u vv_d wv wv_u wv_d +C rr rle chord rle1 chord1 rle2 chord2 wstrip ess +C ensy ensz xsref ysref zsref cdstrp clstrp cfstrp +C cmstrp cf_lstrp cd_lstrp cl_lstrp cdv_lstrp clt_lstrp +C cmc4_lstrp cmle_lstrp cnc dwwake rv1 rv2 rv rc +C gam gam_u gam_d src src_u vv vv_u vv_d wv wv_u +C wv_d C RW status of diff variables: alfa:out vinf:out vinf_a:out vinf_b:out C wrot:out sref:out cref:out bref:out xyzref:out C mach:out cdref:out clff:in-zero cyff:in-zero cdff:in-zero @@ -56,11 +61,14 @@ C cytot_rz:in-zero crtot_rz:in-zero cmtot_rz:in-zero C cntot_rz:in-zero xnp:in-out sm:in-out bb:in-out C rr:in-out rle:out chord:out rle1:out chord1:out -C rle2:out chord2:out wstrip:out ensy:out ensz:out -C xsref:out ysref:out zsref:out rv1:out rv2:out -C rv:out rc:out gam:out gam_u:out gam_d:out src:out -C src_u:out vv:out vv_u:out vv_d:out wv:out wv_u:out -C wv_d:out +C rle2:out chord2:out wstrip:out ess:out ensy:out +C ensz:out xsref:out ysref:out zsref:out cdstrp:in-out +C clstrp:in-out cfstrp:in-out cmstrp:in-out cf_lstrp:in-out +C cd_lstrp:in-out cl_lstrp:in-out cdv_lstrp:in-out +C clt_lstrp:in-out cmc4_lstrp:in-out cmle_lstrp:in-out +C cnc:in-out dwwake:in-out rv1:out rv2:out rv:out +C rc:out gam:out gam_u:out gam_d:out src:out src_u:out +C vv:out vv_u:out vv_d:out wv:out wv_u:out wv_d:out C*********************************************************************** C Module: aero.f C @@ -431,12 +439,16 @@ SUBROUTINE AERO_B() C xyzref cdtot cytot cltot cdtot_a cltot_a cdtot_u C cytot_u cltot_u cdtot_d cytot_d cltot_d cftot C cftot_u cftot_d cmtot cmtot_u cmtot_d cdvtot chord -C rv1 rv2 gam +C cdstrp clstrp cfstrp cmstrp cf_lstrp cd_lstrp +C cl_lstrp cdv_lstrp clt_lstrp cmc4_lstrp cmle_lstrp +C cnc rv1 rv2 gam C with respect to varying inputs: alfa vinf wrot sref cref bref C xyzref cdtot_d cytot_d cltot_d cftot cftot_d cmtot C cmtot_d rle chord rle1 chord1 rle2 chord2 wstrip -C ensy ensz xsref ysref zsref rv1 rv2 rv gam gam_u -C gam_d vv vv_u vv_d wv wv_u wv_d +C ess ensy ensz xsref ysref zsref cdstrp clstrp +C cfstrp cmstrp cf_lstrp cd_lstrp cl_lstrp cdv_lstrp +C clt_lstrp cmc4_lstrp cmle_lstrp cnc rv1 rv2 rv +C gam gam_u gam_d vv vv_u vv_d wv wv_u wv_d C AERO C C @@ -454,6 +466,7 @@ SUBROUTINE SFFORC_B() REAL vrot(3), vrot_u(3), wrot_u(3) REAL vrot_diff(3), vrot_u_diff(3), wrot_u_diff(3) REAL vperp(3) + REAL vperp_diff(3) REAL g(3), r(3), rh(3), mh(3) REAL g_diff(3), r_diff(3) REAL f(3), f_u(3, 6), f_d(3, ndmax), f_g(3, ngmax) @@ -568,12 +581,19 @@ SUBROUTINE SFFORC_B() REAL vsq REAL vsqi REAL vspan + REAL vspan_diff REAL vpsq + REAL vpsq_diff REAL vpsqi + REAL vpsqi_diff REAL delx + REAL delx_diff REAL dely + REAL dely_diff REAL delz + REAL delz_diff REAL dmag + REAL dmag_diff INTEGER is INTEGER jj REAL dcm @@ -582,13 +602,22 @@ SUBROUTINE SFFORC_B() REAL result1 REAL result1_diff REAL result2 + REAL temp REAL temp_diff INTEGER ii1 - INTEGER branch + REAL temp0 REAL temp_diff0 + REAL temp1 + REAL temp2 REAL temp_diff1 - REAL(kind=8) temp_diff2 - REAL(kind=avl_real) temp_diff3 + REAL temp_diff2 + REAL temp3 + REAL temp4 + REAL temp_diff3 + REAL temp_diff4 + REAL(kind=8) temp_diff5 + INTEGER branch + REAL(kind=avl_real) temp_diff6 INTEGER ii2 INTEGER ii3 DATA icrs /2, 3, 1/ @@ -619,6 +648,7 @@ SUBROUTINE SFFORC_B() CALL PUSHREAL8ARRAY(f_d, 3*ndmax) CALL PUSHREAL8ARRAY(ulift, 3) CALL PUSHREAL8ARRAY(rc4, 3) + CALL PUSHREAL8ARRAY(vperp, 3) CALL PUSHREAL8ARRAY(veffmag_u, numax) CALL PUSHREAL8ARRAY(fgam, 3) CALL PUSHREAL8ARRAY(f_u, 3*6) @@ -1639,15 +1669,9 @@ SUBROUTINE SFFORC_B() DO ii1=1,NSTRIP wstrip_diff(ii1) = 0.D0 ENDDO - DO ii1=1,NSTRIP - cdstrp_diff(ii1) = 0.D0 - ENDDO DO ii1=1,NSTRIP cystrp_diff(ii1) = 0.D0 ENDDO - DO ii1=1,NSTRIP - clstrp_diff(ii1) = 0.D0 - ENDDO DO ii1=1,NSTRIP cdst_a_diff(ii1) = 0.D0 ENDDO @@ -1684,11 +1708,6 @@ SUBROUTINE SFFORC_B() clst_d_diff(ii2, ii1) = 0.D0 ENDDO ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - cfstrp_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO DO ii1=1,numax DO ii2=1,NSTRIP DO ii3=1,3 @@ -1703,11 +1722,6 @@ SUBROUTINE SFFORC_B() ENDDO ENDDO ENDDO - DO ii1=1,NSTRIP - DO ii2=1,3 - cmstrp_diff(ii2, ii1) = 0.D0 - ENDDO - ENDDO DO ii1=1,numax DO ii2=1,NSTRIP DO ii3=1,3 @@ -1722,9 +1736,6 @@ SUBROUTINE SFFORC_B() ENDDO ENDDO ENDDO - DO ii1=1,NSTRIP - cdv_lstrp_diff(ii1) = 0.D0 - ENDDO DO is=nsurf,1,-1 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN @@ -1784,169 +1795,169 @@ SUBROUTINE SFFORC_B() sr_diff = 0.D0 cr_diff = 0.D0 DO n=ncontrol,1,-1 - temp_diff2 = cms_d_diff(3, is, n)/(sref*bref) + temp_diff5 = cms_d_diff(3, is, n)/(sref*bref) cmst_d_diff(3, j, n) = cmst_d_diff(3, j, n) + sr*cr* - + temp_diff2 - temp_diff1 = cmst_d(3, j, n)*temp_diff2 - temp_diff3 = -(cmst_d(3, j, n)*sr*cr*temp_diff2/(sref*bref)) - sref_diff = sref_diff + bref*temp_diff3 - bref_diff = bref_diff + sref*temp_diff3 - sr_diff = sr_diff + cr*temp_diff1 - cr_diff = cr_diff + sr*temp_diff1 - temp_diff2 = cms_d_diff(2, is, n)/(sref*cref) + + temp_diff5 + temp_diff3 = cmst_d(3, j, n)*temp_diff5 + temp_diff6 = -(cmst_d(3, j, n)*sr*cr*temp_diff5/(sref*bref)) + sref_diff = sref_diff + bref*temp_diff6 + bref_diff = bref_diff + sref*temp_diff6 + sr_diff = sr_diff + cr*temp_diff3 + cr_diff = cr_diff + sr*temp_diff3 + temp_diff5 = cms_d_diff(2, is, n)/(sref*cref) cmst_d_diff(2, j, n) = cmst_d_diff(2, j, n) + sr*cr* - + temp_diff2 - temp_diff1 = cmst_d(2, j, n)*temp_diff2 - temp_diff3 = -(cmst_d(2, j, n)*sr*cr*temp_diff2/(sref*cref)) - sref_diff = sref_diff + cref*temp_diff3 - cref_diff = cref_diff + sref*temp_diff3 - sr_diff = sr_diff + cr*temp_diff1 - cr_diff = cr_diff + sr*temp_diff1 - temp_diff2 = cms_d_diff(1, is, n)/(sref*bref) + + temp_diff5 + temp_diff3 = cmst_d(2, j, n)*temp_diff5 + temp_diff6 = -(cmst_d(2, j, n)*sr*cr*temp_diff5/(sref*cref)) + sref_diff = sref_diff + cref*temp_diff6 + cref_diff = cref_diff + sref*temp_diff6 + sr_diff = sr_diff + cr*temp_diff3 + cr_diff = cr_diff + sr*temp_diff3 + temp_diff5 = cms_d_diff(1, is, n)/(sref*bref) cmst_d_diff(1, j, n) = cmst_d_diff(1, j, n) + sr*cr* - + temp_diff2 - temp_diff1 = cmst_d(1, j, n)*temp_diff2 - temp_diff3 = -(cmst_d(1, j, n)*sr*cr*temp_diff2/(sref*bref)) - sref_diff = sref_diff + bref*temp_diff3 - bref_diff = bref_diff + sref*temp_diff3 - sr_diff = sr_diff + cr*temp_diff1 - cr_diff = cr_diff + sr*temp_diff1 + + temp_diff5 + temp_diff3 = cmst_d(1, j, n)*temp_diff5 + temp_diff6 = -(cmst_d(1, j, n)*sr*cr*temp_diff5/(sref*bref)) + sref_diff = sref_diff + bref*temp_diff6 + bref_diff = bref_diff + sref*temp_diff6 + sr_diff = sr_diff + cr*temp_diff3 + cr_diff = cr_diff + sr*temp_diff3 DO l=3,1,-1 cfst_d_diff(l, j, n) = cfst_d_diff(l, j, n) + sr* + cfs_d_diff(l, is, n)/sref - temp_diff2 = cfst_d(l, j, n)*cfs_d_diff(l, is, n)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cfst_d(l, j, n)*cfs_d_diff(l, is, n)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref ENDDO clst_d_diff(j, n) = clst_d_diff(j, n) + sr*cls_d_diff(is, n) + /sref - temp_diff2 = clst_d(j, n)*cls_d_diff(is, n)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = clst_d(j, n)*cls_d_diff(is, n)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cyst_d_diff(j, n) = cyst_d_diff(j, n) + sr*cys_d_diff(is, n) + /sref - temp_diff2 = cyst_d(j, n)*cys_d_diff(is, n)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cyst_d(j, n)*cys_d_diff(is, n)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cdst_d_diff(j, n) = cdst_d_diff(j, n) + sr*cds_d_diff(is, n) + /sref - temp_diff2 = cdst_d(j, n)*cds_d_diff(is, n)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cdst_d(j, n)*cds_d_diff(is, n)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref ENDDO DO n=numax,1,-1 - temp_diff2 = cms_u_diff(3, is, n)/(sref*bref) + temp_diff5 = cms_u_diff(3, is, n)/(sref*bref) cmst_u_diff(3, j, n) = cmst_u_diff(3, j, n) + sr*cr* - + temp_diff2 - temp_diff1 = cmst_u(3, j, n)*temp_diff2 - temp_diff3 = -(cmst_u(3, j, n)*sr*cr*temp_diff2/(sref*bref)) - sref_diff = sref_diff + bref*temp_diff3 - bref_diff = bref_diff + sref*temp_diff3 - sr_diff = sr_diff + cr*temp_diff1 - cr_diff = cr_diff + sr*temp_diff1 - temp_diff2 = cms_u_diff(2, is, n)/(sref*cref) + + temp_diff5 + temp_diff3 = cmst_u(3, j, n)*temp_diff5 + temp_diff6 = -(cmst_u(3, j, n)*sr*cr*temp_diff5/(sref*bref)) + sref_diff = sref_diff + bref*temp_diff6 + bref_diff = bref_diff + sref*temp_diff6 + sr_diff = sr_diff + cr*temp_diff3 + cr_diff = cr_diff + sr*temp_diff3 + temp_diff5 = cms_u_diff(2, is, n)/(sref*cref) cmst_u_diff(2, j, n) = cmst_u_diff(2, j, n) + sr*cr* - + temp_diff2 - temp_diff1 = cmst_u(2, j, n)*temp_diff2 - temp_diff3 = -(cmst_u(2, j, n)*sr*cr*temp_diff2/(sref*cref)) - sref_diff = sref_diff + cref*temp_diff3 - cref_diff = cref_diff + sref*temp_diff3 - sr_diff = sr_diff + cr*temp_diff1 - cr_diff = cr_diff + sr*temp_diff1 - temp_diff2 = cms_u_diff(1, is, n)/(sref*bref) + + temp_diff5 + temp_diff3 = cmst_u(2, j, n)*temp_diff5 + temp_diff6 = -(cmst_u(2, j, n)*sr*cr*temp_diff5/(sref*cref)) + sref_diff = sref_diff + cref*temp_diff6 + cref_diff = cref_diff + sref*temp_diff6 + sr_diff = sr_diff + cr*temp_diff3 + cr_diff = cr_diff + sr*temp_diff3 + temp_diff5 = cms_u_diff(1, is, n)/(sref*bref) cmst_u_diff(1, j, n) = cmst_u_diff(1, j, n) + sr*cr* - + temp_diff2 - temp_diff1 = cmst_u(1, j, n)*temp_diff2 - temp_diff3 = -(cmst_u(1, j, n)*sr*cr*temp_diff2/(sref*bref)) - sref_diff = sref_diff + bref*temp_diff3 - bref_diff = bref_diff + sref*temp_diff3 - sr_diff = sr_diff + cr*temp_diff1 - cr_diff = cr_diff + sr*temp_diff1 + + temp_diff5 + temp_diff3 = cmst_u(1, j, n)*temp_diff5 + temp_diff6 = -(cmst_u(1, j, n)*sr*cr*temp_diff5/(sref*bref)) + sref_diff = sref_diff + bref*temp_diff6 + bref_diff = bref_diff + sref*temp_diff6 + sr_diff = sr_diff + cr*temp_diff3 + cr_diff = cr_diff + sr*temp_diff3 DO l=3,1,-1 cfst_u_diff(l, j, n) = cfst_u_diff(l, j, n) + sr* + cfs_u_diff(l, is, n)/sref - temp_diff2 = cfst_u(l, j, n)*cfs_u_diff(l, is, n)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cfst_u(l, j, n)*cfs_u_diff(l, is, n)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref ENDDO clst_u_diff(j, n) = clst_u_diff(j, n) + sr*cls_u_diff(is, n) + /sref - temp_diff2 = clst_u(j, n)*cls_u_diff(is, n)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = clst_u(j, n)*cls_u_diff(is, n)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cyst_u_diff(j, n) = cyst_u_diff(j, n) + sr*cys_u_diff(is, n) + /sref - temp_diff2 = cyst_u(j, n)*cys_u_diff(is, n)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cyst_u(j, n)*cys_u_diff(is, n)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cdst_u_diff(j, n) = cdst_u_diff(j, n) + sr*cds_u_diff(is, n) + /sref - temp_diff2 = cdst_u(j, n)*cds_u_diff(is, n)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cdst_u(j, n)*cds_u_diff(is, n)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref ENDDO clst_a_diff(j) = clst_a_diff(j) + sr*cls_a_diff(is)/sref - temp_diff2 = clst_a(j)*cls_a_diff(is)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = clst_a(j)*cls_a_diff(is)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cdst_a_diff(j) = cdst_a_diff(j) + sr*cds_a_diff(is)/sref - temp_diff2 = cdst_a(j)*cds_a_diff(is)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cdst_a(j)*cds_a_diff(is)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cdv_lstrp_diff(j) = cdv_lstrp_diff(j) + sr*cdvsurf_diff(is)/ + sref - temp_diff2 = cdv_lstrp(j)*cdvsurf_diff(is)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref - temp_diff2 = cmsurf_diff(3, is)/(sref*bref) - cmstrp_diff(3, j) = cmstrp_diff(3, j) + sr*cr*temp_diff2 - temp_diff1 = cmstrp(3, j)*temp_diff2 - temp_diff3 = -(cmstrp(3, j)*sr*cr*temp_diff2/(sref*bref)) - sref_diff = sref_diff + bref*temp_diff3 - bref_diff = bref_diff + sref*temp_diff3 - sr_diff = sr_diff + cr*temp_diff1 - cr_diff = cr_diff + sr*temp_diff1 - temp_diff2 = cmsurf_diff(2, is)/(sref*cref) - cmstrp_diff(2, j) = cmstrp_diff(2, j) + sr*cr*temp_diff2 - temp_diff1 = cmstrp(2, j)*temp_diff2 - temp_diff3 = -(cmstrp(2, j)*sr*cr*temp_diff2/(sref*cref)) - sref_diff = sref_diff + cref*temp_diff3 - cref_diff = cref_diff + sref*temp_diff3 - sr_diff = sr_diff + cr*temp_diff1 - cr_diff = cr_diff + sr*temp_diff1 - temp_diff2 = cmsurf_diff(1, is)/(sref*bref) - cmstrp_diff(1, j) = cmstrp_diff(1, j) + sr*cr*temp_diff2 - temp_diff1 = cmstrp(1, j)*temp_diff2 - temp_diff3 = -(cmstrp(1, j)*sr*cr*temp_diff2/(sref*bref)) - bref_diff = bref_diff + sref*temp_diff3 - cr_diff = cr_diff + sr*temp_diff1 + temp_diff5 = cdv_lstrp(j)*cdvsurf_diff(is)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref + temp_diff5 = cmsurf_diff(3, is)/(sref*bref) + cmstrp_diff(3, j) = cmstrp_diff(3, j) + sr*cr*temp_diff5 + temp_diff3 = cmstrp(3, j)*temp_diff5 + temp_diff6 = -(cmstrp(3, j)*sr*cr*temp_diff5/(sref*bref)) + sref_diff = sref_diff + bref*temp_diff6 + bref_diff = bref_diff + sref*temp_diff6 + sr_diff = sr_diff + cr*temp_diff3 + cr_diff = cr_diff + sr*temp_diff3 + temp_diff5 = cmsurf_diff(2, is)/(sref*cref) + cmstrp_diff(2, j) = cmstrp_diff(2, j) + sr*cr*temp_diff5 + temp_diff3 = cmstrp(2, j)*temp_diff5 + temp_diff6 = -(cmstrp(2, j)*sr*cr*temp_diff5/(sref*cref)) + sref_diff = sref_diff + cref*temp_diff6 + cref_diff = cref_diff + sref*temp_diff6 + sr_diff = sr_diff + cr*temp_diff3 + cr_diff = cr_diff + sr*temp_diff3 + temp_diff5 = cmsurf_diff(1, is)/(sref*bref) + cmstrp_diff(1, j) = cmstrp_diff(1, j) + sr*cr*temp_diff5 + temp_diff3 = cmstrp(1, j)*temp_diff5 + temp_diff6 = -(cmstrp(1, j)*sr*cr*temp_diff5/(sref*bref)) + bref_diff = bref_diff + sref*temp_diff6 + cr_diff = cr_diff + sr*temp_diff3 cfstrp_diff(3, j) = cfstrp_diff(3, j) + sr*cfsurf_diff(3, is)/ + sref - temp_diff2 = cfstrp(3, j)*cfsurf_diff(3, is)/sref - sref_diff = sref_diff + bref*temp_diff3 - sr*temp_diff2/sref - sr_diff = sr_diff + cr*temp_diff1 + temp_diff2 + temp_diff5 = cfstrp(3, j)*cfsurf_diff(3, is)/sref + sref_diff = sref_diff + bref*temp_diff6 - sr*temp_diff5/sref + sr_diff = sr_diff + cr*temp_diff3 + temp_diff5 cfstrp_diff(2, j) = cfstrp_diff(2, j) + sr*cfsurf_diff(2, is)/ + sref - temp_diff2 = cfstrp(2, j)*cfsurf_diff(2, is)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cfstrp(2, j)*cfsurf_diff(2, is)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cfstrp_diff(1, j) = cfstrp_diff(1, j) + sr*cfsurf_diff(1, is)/ + sref - temp_diff2 = cfstrp(1, j)*cfsurf_diff(1, is)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cfstrp(1, j)*cfsurf_diff(1, is)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref clstrp_diff(j) = clstrp_diff(j) + sr*clsurf_diff(is)/sref - temp_diff2 = clstrp(j)*clsurf_diff(is)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = clstrp(j)*clsurf_diff(is)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cystrp_diff(j) = cystrp_diff(j) + sr*cysurf_diff(is)/sref - temp_diff2 = cystrp(j)*cysurf_diff(is)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cystrp(j)*cysurf_diff(is)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref cdstrp_diff(j) = cdstrp_diff(j) + sr*cdsurf_diff(is)/sref - temp_diff2 = cdstrp(j)*cdsurf_diff(is)/sref - sr_diff = sr_diff + temp_diff2 - sref_diff = sref_diff - sr*temp_diff2/sref + temp_diff5 = cdstrp(j)*cdsurf_diff(is)/sref + sr_diff = sr_diff + temp_diff5 + sref_diff = sref_diff - sr*temp_diff5/sref chord_diff(j) = chord_diff(j) + cr_diff + wstrip(j)*sr_diff wstrip_diff(j) = wstrip_diff(j) + chord(j)*sr_diff ENDDO @@ -2004,6 +2015,11 @@ SUBROUTINE SFFORC_B() DO ii1=1,NSTRIP chord2_diff(ii1) = 0.D0 ENDDO + DO ii1=1,NSTRIP + DO ii2=1,3 + ess_diff(ii2, ii1) = 0.D0 + ENDDO + ENDDO DO ii1=1,NSTRIP ensy_diff(ii1) = 0.D0 ENDDO @@ -2157,6 +2173,9 @@ SUBROUTINE SFFORC_B() rc4_diff(ii1) = 0.D0 ENDDO sina_diff = 0.D0 + DO ii1=1,3 + vperp_diff(ii1) = 0.D0 + ENDDO DO ii1=1,numax veffmag_u_diff(ii1) = 0.D0 ENDDO @@ -2201,6 +2220,7 @@ SUBROUTINE SFFORC_B() CALL POPREAL8ARRAY(f_u, 3*6) CALL POPREAL8ARRAY(fgam, 3) CALL POPREAL8ARRAY(veffmag_u, numax) + CALL POPREAL8ARRAY(vperp, 3) CALL POPREAL8ARRAY(rc4, 3) CALL POPREAL8ARRAY(ulift, 3) CALL POPREAL8ARRAY(f_d, 3*ndmax) @@ -2305,6 +2325,9 @@ SUBROUTINE SFFORC_B() cfx = 0. cfy = 0. cfz = 0. + cmx = 0. + cmy = 0. + cmz = 0. C DO n=1,numax cfx_u(n) = 0. @@ -3001,6 +3024,7 @@ SUBROUTINE SFFORC_B() C normalized by strip area and chord C C...Components of X,Y,Z forces in local strip axes + cl_lstrp(j) = ulift(1)*cfx + ulift(2)*cfy + ulift(3)*cfz C C...Axial/normal forces and lift/drag in plane normal to dihedral of strip C...CN,CA forces are rotated to be in and normal to strip incidence @@ -3016,22 +3040,166 @@ SUBROUTINE SFFORC_B() C print *,"WROT ",WROT C C------ set total effective velocity = freestream + rotation + CALL CROSS(rrot, wrot, vrot) + CALL PUSHREAL8(veff(1)) + veff(1) = vinf(1) + vrot(1) + CALL PUSHREAL8(veff(2)) + veff(2) = vinf(2) + vrot(2) + CALL PUSHREAL8(veff(3)) + veff(3) = vinf(3) + vrot(3) +C +C +C------ spanwise and perpendicular velocity components + vspan = veff(1)*ess(1, j) + veff(2)*ess(2, j) + veff(3)*ess(3, j + + ) + vperp(1) = veff(1) - ess(1, j)*vspan + vperp(2) = veff(2) - ess(2, j)*vspan + vperp(3) = veff(3) - ess(3, j)*vspan C + vpsq = vperp(1)**2 + vperp(2)**2 + vperp(3)**2 + IF (vpsq .EQ. 0.0) THEN + vpsqi = 1.0 + CALL PUSHCONTROL1B(0) + ELSE + vpsqi = 1.0/vpsq + CALL PUSHCONTROL1B(1) + END IF +Ccc CLT_LSTRP(J) = CN_LSTRP(J) * VPSQI +C +C--- Moment about strip LE midpoint in direction of LE segment + CALL PUSHREAL8(r(1)) + r(1) = rc4(1) - rle(1, j) + CALL PUSHREAL8(r(2)) + r(2) = rc4(2) - rle(2, j) + CALL PUSHREAL8(r(3)) + r(3) = rc4(3) - rle(3, j) + delx = rle2(1, j) - rle1(1, j) + dely = rle2(2, j) - rle1(2, j) + delz = rle2(3, j) - rle1(3, j) +C + IF (imags(lssurf(j)) .LT. 0) THEN + delx = -delx + dely = -dely + delz = -delz + CALL PUSHCONTROL1B(0) + ELSE + CALL PUSHCONTROL1B(1) + END IF + dmag = SQRT(delx**2 + dely**2 + delz**2) + IF (dmag .NE. 0.0) THEN + temp = delx/dmag + temp0 = (cfz*r(2)-cfy*r(3))/cr + temp1 = dely/dmag + temp2 = (cfx*r(3)-cfz*r(1))/cr + temp3 = delz/dmag + temp4 = (cfy*r(1)-cfx*r(2))/cr + cmx_diff = temp*cmle_lstrp_diff(j) + temp_diff0 = temp*cmle_lstrp_diff(j)/cr + temp_diff = (cmx+temp0)*cmle_lstrp_diff(j)/dmag + cmy_diff = temp1*cmle_lstrp_diff(j) + temp_diff1 = temp1*cmle_lstrp_diff(j)/cr + temp_diff2 = (cmy+temp2)*cmle_lstrp_diff(j)/dmag + cmz_diff = temp3*cmle_lstrp_diff(j) + temp_diff3 = temp3*cmle_lstrp_diff(j)/cr + temp_diff4 = (cmz+temp4)*cmle_lstrp_diff(j)/dmag + cmle_lstrp_diff(j) = 0.D0 + delz_diff = temp_diff4 + dmag_diff = -(temp3*temp_diff4) - temp1*temp_diff2 - temp* + + temp_diff + cfy_diff = r(1)*temp_diff3 - r(3)*temp_diff0 + r_diff(1) = r_diff(1) + cfy*temp_diff3 - cfz*temp_diff1 + cfx_diff = r(3)*temp_diff1 - r(2)*temp_diff3 + r_diff(2) = r_diff(2) + cfz*temp_diff0 - cfx*temp_diff3 + cr_diff = -(temp4*temp_diff3) - temp2*temp_diff1 - temp0* + + temp_diff0 + dely_diff = temp_diff2 + r_diff(3) = r_diff(3) + cfx*temp_diff1 - cfy*temp_diff0 + cfz_diff = r(2)*temp_diff0 - r(1)*temp_diff1 + delx_diff = temp_diff + ELSE + cfx_diff = 0.D0 + cfy_diff = 0.D0 + cfz_diff = 0.D0 + delx_diff = 0.D0 + dely_diff = 0.D0 + delz_diff = 0.D0 + cmx_diff = 0.D0 + cmy_diff = 0.D0 + cmz_diff = 0.D0 + dmag_diff = 0.D0 + cr_diff = 0.D0 + END IF + cmle_lstrp_diff(j) = 0.D0 + IF (delx**2 + dely**2 + delz**2 .EQ. 0.D0) THEN + temp_diff = 0.D0 + ELSE + temp_diff = dmag_diff/(2.0*SQRT(delx**2+dely**2+delz**2)) + END IF + delx_diff = delx_diff + 2*delx*temp_diff + dely_diff = dely_diff + 2*dely*temp_diff + delz_diff = delz_diff + 2*delz*temp_diff + CALL POPCONTROL1B(branch) + IF (branch .EQ. 0) THEN + delz_diff = -delz_diff + dely_diff = -dely_diff + delx_diff = -delx_diff + END IF + rle2_diff(3, j) = rle2_diff(3, j) + delz_diff + rle1_diff(3, j) = rle1_diff(3, j) - delz_diff + rle2_diff(2, j) = rle2_diff(2, j) + dely_diff + rle1_diff(2, j) = rle1_diff(2, j) - dely_diff + rle2_diff(1, j) = rle2_diff(1, j) + delx_diff + rle1_diff(1, j) = rle1_diff(1, j) - delx_diff + CALL POPREAL8(r(3)) rc4_diff(3) = rc4_diff(3) + r_diff(3) rle_diff(3, j) = rle_diff(3, j) - r_diff(3) r_diff(3) = 0.D0 + CALL POPREAL8(r(2)) rc4_diff(2) = rc4_diff(2) + r_diff(2) rle_diff(2, j) = rle_diff(2, j) - r_diff(2) r_diff(2) = 0.D0 + CALL POPREAL8(r(1)) rc4_diff(1) = rc4_diff(1) + r_diff(1) rle_diff(1, j) = rle_diff(1, j) - r_diff(1) r_diff(1) = 0.D0 + cl_lstrp_diff(j) = cl_lstrp_diff(j) + vpsqi*clt_lstrp_diff(j) + vpsqi_diff = cl_lstrp(j)*clt_lstrp_diff(j) + clt_lstrp_diff(j) = 0.D0 + CALL POPCONTROL1B(branch) + IF (branch .EQ. 0) THEN + vpsq_diff = 0.D0 + ELSE + vpsq_diff = -(vpsqi_diff/vpsq**2) + END IF + vperp_diff(1) = vperp_diff(1) + 2*vperp(1)*vpsq_diff + vperp_diff(2) = vperp_diff(2) + 2*vperp(2)*vpsq_diff + vperp_diff(3) = vperp_diff(3) + 2*vperp(3)*vpsq_diff + vspan_diff = -(ess(3, j)*vperp_diff(3)) - ess(2, j)*vperp_diff(2 + + ) - ess(1, j)*vperp_diff(1) + veff_diff(3) = veff_diff(3) + vperp_diff(3) + ess(3, j)* + + vspan_diff + ess_diff(3, j) = ess_diff(3, j) + veff(3)*vspan_diff - vspan* + + vperp_diff(3) + vperp_diff(3) = 0.D0 + veff_diff(2) = veff_diff(2) + vperp_diff(2) + ess(2, j)* + + vspan_diff + ess_diff(2, j) = ess_diff(2, j) + veff(2)*vspan_diff - vspan* + + vperp_diff(2) + vperp_diff(2) = 0.D0 + veff_diff(1) = veff_diff(1) + vperp_diff(1) + ess(1, j)* + + vspan_diff + ess_diff(1, j) = ess_diff(1, j) + veff(1)*vspan_diff - vspan* + + vperp_diff(1) + vperp_diff(1) = 0.D0 + CALL POPREAL8(veff(3)) vinf_diff(3) = vinf_diff(3) + veff_diff(3) vrot_diff(3) = vrot_diff(3) + veff_diff(3) veff_diff(3) = 0.D0 + CALL POPREAL8(veff(2)) vinf_diff(2) = vinf_diff(2) + veff_diff(2) vrot_diff(2) = vrot_diff(2) + veff_diff(2) veff_diff(2) = 0.D0 + CALL POPREAL8(veff(1)) vinf_diff(1) = vinf_diff(1) + veff_diff(1) vrot_diff(1) = vrot_diff(1) + veff_diff(1) veff_diff(1) = 0.D0 @@ -3048,84 +3216,102 @@ SUBROUTINE SFFORC_B() xsref_diff(j) = xsref_diff(j) + rrot_diff(1) xyzref_diff(1) = xyzref_diff(1) - rrot_diff(1) rrot_diff(1) = 0.D0 - cr_diff = 0.D0 + ensz_diff(j) = ensz_diff(j) + cmy*cmc4_lstrp_diff(j) + cmy_diff = cmy_diff + ensz(j)*cmc4_lstrp_diff(j) + ensy_diff(j) = ensy_diff(j) - cmz*cmc4_lstrp_diff(j) + cmz_diff = cmz_diff - ensy(j)*cmc4_lstrp_diff(j) + cmc4_lstrp_diff(j) = 0.D0 + udrag_diff(1) = udrag_diff(1) + cfx*cd_lstrp_diff(j) + cfx_diff = cfx_diff + udrag(1)*cd_lstrp_diff(j) + ulift(1)* + + cl_lstrp_diff(j) + udrag_diff(2) = udrag_diff(2) + cfy*cd_lstrp_diff(j) + cfy_diff = cfy_diff + udrag(2)*cd_lstrp_diff(j) + ulift(2)* + + cl_lstrp_diff(j) + udrag_diff(3) = udrag_diff(3) + cfz*cd_lstrp_diff(j) + cfz_diff = cfz_diff + udrag(3)*cd_lstrp_diff(j) + ulift(3)* + + cl_lstrp_diff(j) + cd_lstrp_diff(j) = 0.D0 + ulift_diff(1) = ulift_diff(1) + cfx*cl_lstrp_diff(j) + ulift_diff(2) = ulift_diff(2) + cfy*cl_lstrp_diff(j) + ulift_diff(3) = ulift_diff(3) + cfz*cl_lstrp_diff(j) + cl_lstrp_diff(j) = 0.D0 C$BWD-OF II-LOOP DO n=1,ncontrol cmz_d_diff(n) = cmz_d_diff(n) + cmst_d_diff(3, j, n) - temp_diff1 = cmst_d_diff(3, j, n)/cr + temp_diff3 = cmst_d_diff(3, j, n)/cr cmst_d_diff(3, j, n) = 0.D0 - cfy_d_diff(n) = cfy_d_diff(n) + r(1)*temp_diff1 - r_diff(1) = r_diff(1) + cfy_d(n)*temp_diff1 - cfx_d_diff(n) = cfx_d_diff(n) - r(2)*temp_diff1 - r_diff(2) = r_diff(2) - cfx_d(n)*temp_diff1 - cr_diff = cr_diff - (cfy_d(n)*r(1)-cfx_d(n)*r(2))*temp_diff1/ + cfy_d_diff(n) = cfy_d_diff(n) + r(1)*temp_diff3 + r_diff(1) = r_diff(1) + cfy_d(n)*temp_diff3 + cfx_d_diff(n) = cfx_d_diff(n) - r(2)*temp_diff3 + r_diff(2) = r_diff(2) - cfx_d(n)*temp_diff3 + cr_diff = cr_diff - (cfy_d(n)*r(1)-cfx_d(n)*r(2))*temp_diff3/ + cr cmy_d_diff(n) = cmy_d_diff(n) + cmst_d_diff(2, j, n) - temp_diff1 = cmst_d_diff(2, j, n)/cr + temp_diff3 = cmst_d_diff(2, j, n)/cr cmst_d_diff(2, j, n) = 0.D0 - cfx_d_diff(n) = cfx_d_diff(n) + r(3)*temp_diff1 - r_diff(3) = r_diff(3) + cfx_d(n)*temp_diff1 - cfz_d_diff(n) = cfz_d_diff(n) - r(1)*temp_diff1 - r_diff(1) = r_diff(1) - cfz_d(n)*temp_diff1 - cr_diff = cr_diff - (cfx_d(n)*r(3)-cfz_d(n)*r(1))*temp_diff1/ + cfx_d_diff(n) = cfx_d_diff(n) + r(3)*temp_diff3 + r_diff(3) = r_diff(3) + cfx_d(n)*temp_diff3 + cfz_d_diff(n) = cfz_d_diff(n) - r(1)*temp_diff3 + r_diff(1) = r_diff(1) - cfz_d(n)*temp_diff3 + cr_diff = cr_diff - (cfx_d(n)*r(3)-cfz_d(n)*r(1))*temp_diff3/ + cr cmx_d_diff(n) = cmx_d_diff(n) + cmst_d_diff(1, j, n) - temp_diff1 = cmst_d_diff(1, j, n)/cr + temp_diff3 = cmst_d_diff(1, j, n)/cr cmst_d_diff(1, j, n) = 0.D0 - cfz_d_diff(n) = cfz_d_diff(n) + r(2)*temp_diff1 - r_diff(2) = r_diff(2) + cfz_d(n)*temp_diff1 - cfy_d_diff(n) = cfy_d_diff(n) - r(3)*temp_diff1 - r_diff(3) = r_diff(3) - cfy_d(n)*temp_diff1 - cr_diff = cr_diff - (cfz_d(n)*r(2)-cfy_d(n)*r(3))*temp_diff1/ + cfz_d_diff(n) = cfz_d_diff(n) + r(2)*temp_diff3 + r_diff(2) = r_diff(2) + cfz_d(n)*temp_diff3 + cfy_d_diff(n) = cfy_d_diff(n) - r(3)*temp_diff3 + r_diff(3) = r_diff(3) - cfy_d(n)*temp_diff3 + cr_diff = cr_diff - (cfz_d(n)*r(2)-cfy_d(n)*r(3))*temp_diff3/ + cr ENDDO C$BWD-OF II-LOOP DO n=1,numax cmz_u_diff(n) = cmz_u_diff(n) + cmst_u_diff(3, j, n) - temp_diff1 = cmst_u_diff(3, j, n)/cr + temp_diff3 = cmst_u_diff(3, j, n)/cr cmst_u_diff(3, j, n) = 0.D0 - cfy_u_diff(n) = cfy_u_diff(n) + r(1)*temp_diff1 - r_diff(1) = r_diff(1) + cfy_u(n)*temp_diff1 - cfx_u_diff(n) = cfx_u_diff(n) - r(2)*temp_diff1 - r_diff(2) = r_diff(2) - cfx_u(n)*temp_diff1 - cr_diff = cr_diff - (cfy_u(n)*r(1)-cfx_u(n)*r(2))*temp_diff1/ + cfy_u_diff(n) = cfy_u_diff(n) + r(1)*temp_diff3 + r_diff(1) = r_diff(1) + cfy_u(n)*temp_diff3 + cfx_u_diff(n) = cfx_u_diff(n) - r(2)*temp_diff3 + r_diff(2) = r_diff(2) - cfx_u(n)*temp_diff3 + cr_diff = cr_diff - (cfy_u(n)*r(1)-cfx_u(n)*r(2))*temp_diff3/ + cr cmy_u_diff(n) = cmy_u_diff(n) + cmst_u_diff(2, j, n) - temp_diff1 = cmst_u_diff(2, j, n)/cr + temp_diff3 = cmst_u_diff(2, j, n)/cr cmst_u_diff(2, j, n) = 0.D0 - cfx_u_diff(n) = cfx_u_diff(n) + r(3)*temp_diff1 - r_diff(3) = r_diff(3) + cfx_u(n)*temp_diff1 - cfz_u_diff(n) = cfz_u_diff(n) - r(1)*temp_diff1 - r_diff(1) = r_diff(1) - cfz_u(n)*temp_diff1 - cr_diff = cr_diff - (cfx_u(n)*r(3)-cfz_u(n)*r(1))*temp_diff1/ + cfx_u_diff(n) = cfx_u_diff(n) + r(3)*temp_diff3 + r_diff(3) = r_diff(3) + cfx_u(n)*temp_diff3 + cfz_u_diff(n) = cfz_u_diff(n) - r(1)*temp_diff3 + r_diff(1) = r_diff(1) - cfz_u(n)*temp_diff3 + cr_diff = cr_diff - (cfx_u(n)*r(3)-cfz_u(n)*r(1))*temp_diff3/ + cr cmx_u_diff(n) = cmx_u_diff(n) + cmst_u_diff(1, j, n) - temp_diff1 = cmst_u_diff(1, j, n)/cr + temp_diff3 = cmst_u_diff(1, j, n)/cr cmst_u_diff(1, j, n) = 0.D0 - cfz_u_diff(n) = cfz_u_diff(n) + r(2)*temp_diff1 - r_diff(2) = r_diff(2) + cfz_u(n)*temp_diff1 - cfy_u_diff(n) = cfy_u_diff(n) - r(3)*temp_diff1 - r_diff(3) = r_diff(3) - cfy_u(n)*temp_diff1 - cr_diff = cr_diff - (cfz_u(n)*r(2)-cfy_u(n)*r(3))*temp_diff1/ + cfz_u_diff(n) = cfz_u_diff(n) + r(2)*temp_diff3 + r_diff(2) = r_diff(2) + cfz_u(n)*temp_diff3 + cfy_u_diff(n) = cfy_u_diff(n) - r(3)*temp_diff3 + r_diff(3) = r_diff(3) - cfy_u(n)*temp_diff3 + cr_diff = cr_diff - (cfz_u(n)*r(2)-cfy_u(n)*r(3))*temp_diff3/ + cr ENDDO - cmz_diff = cmstrp_diff(3, j) + cmz_diff = cmz_diff + cmstrp_diff(3, j) temp_diff = cmstrp_diff(3, j)/cr cmstrp_diff(3, j) = 0.D0 - cfy_diff = r(1)*temp_diff + cfy_diff = cfy_diff + r(1)*temp_diff r_diff(1) = r_diff(1) + cfy*temp_diff - cfx_diff = -(r(2)*temp_diff) + cfx_diff = cfx_diff - r(2)*temp_diff r_diff(2) = r_diff(2) - cfx*temp_diff cr_diff = cr_diff - (cfy*r(1)-cfx*r(2))*temp_diff/cr - cmy_diff = cmstrp_diff(2, j) + cmy_diff = cmy_diff + cmstrp_diff(2, j) temp_diff = cmstrp_diff(2, j)/cr cmstrp_diff(2, j) = 0.D0 cfx_diff = cfx_diff + r(3)*temp_diff r_diff(3) = r_diff(3) + cfx*temp_diff - cfz_diff = -(r(1)*temp_diff) + cfz_diff = cfz_diff - r(1)*temp_diff r_diff(1) = r_diff(1) - cfz*temp_diff cr_diff = cr_diff - (cfx*r(3)-cfz*r(1))*temp_diff/cr - cmx_diff = cmstrp_diff(1, j) + cmx_diff = cmx_diff + cmstrp_diff(1, j) temp_diff = cmstrp_diff(1, j)/cr cmstrp_diff(1, j) = 0.D0 cfz_diff = cfz_diff + r(2)*temp_diff @@ -3185,23 +3371,27 @@ SUBROUTINE SFFORC_B() ENDDO cfx_diff = cfx_diff + cosa*cdstrp_diff(j) - cosa*clst_a_diff(j) + - sina*cdst_a_diff(j) - sina*clstrp_diff(j) + cfstrp_diff(1, j - + ) + + ) + cf_lstrp_diff(1, j) cosa_diff = cosa_diff + cfz*cdst_a_diff(j) - cfx*clst_a_diff(j) + + cfz*clstrp_diff(j) + cfx*cdstrp_diff(j) cfz_diff = cfz_diff + cosa*cdst_a_diff(j) - sina*clst_a_diff(j) + + cosa*clstrp_diff(j) + sina*cdstrp_diff(j) + cfstrp_diff(3, j - + ) + + ) + cf_lstrp_diff(3, j) sina_diff = sina_diff + cfz*cdstrp_diff(j) - cfz*clst_a_diff(j) + - cfx*cdst_a_diff(j) - cfx*clstrp_diff(j) clst_a_diff(j) = 0.D0 cdst_a_diff(j) = 0.D0 clstrp_diff(j) = 0.D0 - cfy_diff = cfy_diff + cystrp_diff(j) + cfstrp_diff(2, j) + cfy_diff = cfy_diff + cystrp_diff(j) + cfstrp_diff(2, j) + + + cf_lstrp_diff(2, j) cystrp_diff(j) = 0.D0 cdstrp_diff(j) = 0.D0 cfstrp_diff(3, j) = 0.D0 cfstrp_diff(2, j) = 0.D0 cfstrp_diff(1, j) = 0.D0 + cf_lstrp_diff(3, j) = 0.D0 + cf_lstrp_diff(2, j) = 0.D0 + cf_lstrp_diff(1, j) = 0.D0 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN cdv_clv_diff = 0.D0 @@ -3214,24 +3404,24 @@ SUBROUTINE SFFORC_B() dcvfz_d_diff = cfz_d_diff(n) dcvfy_d_diff = cfy_d_diff(n) dcvfx_d_diff = cfx_d_diff(n) - temp_diff1 = veff(3)*clv_d(n)*dcvfz_d_diff - temp_diff0 = veffmag*cdv_clv*dcvfz_d_diff - veff_diff(3) = veff_diff(3) + clv_d(n)*temp_diff0 - clv_d_diff(n) = clv_d_diff(n) + veff(3)*temp_diff0 - veffmag_diff = veffmag_diff + cdv_clv*temp_diff1 - cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff1 - temp_diff1 = veff(2)*clv_d(n)*dcvfy_d_diff - temp_diff0 = veffmag*cdv_clv*dcvfy_d_diff - veff_diff(2) = veff_diff(2) + clv_d(n)*temp_diff0 - clv_d_diff(n) = clv_d_diff(n) + veff(2)*temp_diff0 - veffmag_diff = veffmag_diff + cdv_clv*temp_diff1 - cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff1 - temp_diff1 = veff(1)*clv_d(n)*dcvfx_d_diff - temp_diff0 = veffmag*cdv_clv*dcvfx_d_diff - veff_diff(1) = veff_diff(1) + clv_d(n)*temp_diff0 - clv_d_diff(n) = clv_d_diff(n) + veff(1)*temp_diff0 - veffmag_diff = veffmag_diff + cdv_clv*temp_diff1 - cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff1 + temp_diff3 = veff(3)*clv_d(n)*dcvfz_d_diff + temp_diff4 = veffmag*cdv_clv*dcvfz_d_diff + veff_diff(3) = veff_diff(3) + clv_d(n)*temp_diff4 + clv_d_diff(n) = clv_d_diff(n) + veff(3)*temp_diff4 + veffmag_diff = veffmag_diff + cdv_clv*temp_diff3 + cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff3 + temp_diff3 = veff(2)*clv_d(n)*dcvfy_d_diff + temp_diff4 = veffmag*cdv_clv*dcvfy_d_diff + veff_diff(2) = veff_diff(2) + clv_d(n)*temp_diff4 + clv_d_diff(n) = clv_d_diff(n) + veff(2)*temp_diff4 + veffmag_diff = veffmag_diff + cdv_clv*temp_diff3 + cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff3 + temp_diff3 = veff(1)*clv_d(n)*dcvfx_d_diff + temp_diff4 = veffmag*cdv_clv*dcvfx_d_diff + veff_diff(1) = veff_diff(1) + clv_d(n)*temp_diff4 + clv_d_diff(n) = clv_d_diff(n) + veff(1)*temp_diff4 + veffmag_diff = veffmag_diff + cdv_clv*temp_diff3 + cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff3 ENDDO cdv_diff = 0.D0 CALL POPREAL8ARRAY(cfz_u, numax) @@ -3242,42 +3432,43 @@ SUBROUTINE SFFORC_B() dcvfz_u_diff = cfz_u_diff(n) dcvfy_u_diff = cfy_u_diff(n) dcvfx_u_diff = cfx_u_diff(n) - temp_diff1 = cdv*dcvfz_u_diff + temp_diff3 = cdv*dcvfz_u_diff cdv_diff = cdv_diff + (veff_u(3, n)*veffmag+veff(3)* + veffmag_u(n))*dcvfz_u_diff + (veff_u(2, n)*veffmag+veff(2) + *veffmag_u(n))*dcvfy_u_diff + (veff_u(1, n)*veffmag+veff(1 + )*veffmag_u(n))*dcvfx_u_diff - temp_diff0 = veff(3)*clv_u(n)*dcvfz_u_diff - temp_diff = veffmag*cdv_clv*dcvfz_u_diff - veff_diff(3) = veff_diff(3) + clv_u(n)*temp_diff + veffmag_u - + (n)*temp_diff1 - clv_u_diff(n) = clv_u_diff(n) + veff(3)*temp_diff - veffmag_diff = veffmag_diff + cdv_clv*temp_diff0 - cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff0 - veff_u_diff(3, n) = veff_u_diff(3, n) + veffmag*temp_diff1 - veffmag_u_diff(n) = veffmag_u_diff(n) + veff(3)*temp_diff1 - temp_diff0 = veff(2)*clv_u(n)*dcvfy_u_diff - veffmag_diff = veffmag_diff + veff_u(3, n)*temp_diff1 + - + cdv_clv*temp_diff0 - temp_diff1 = cdv*dcvfy_u_diff - temp_diff = veffmag*cdv_clv*dcvfy_u_diff - veff_diff(2) = veff_diff(2) + clv_u(n)*temp_diff + veffmag_u - + (n)*temp_diff1 - clv_u_diff(n) = clv_u_diff(n) + veff(2)*temp_diff - cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff0 - veff_u_diff(2, n) = veff_u_diff(2, n) + veffmag*temp_diff1 - temp_diff = cdv*dcvfx_u_diff - veffmag_u_diff(n) = veffmag_u_diff(n) + veff(2)*temp_diff1 + - + veff(1)*temp_diff - temp_diff0 = veff(1)*clv_u(n)*dcvfx_u_diff - veffmag_diff = veffmag_diff + veff_u(2, n)*temp_diff1 + - + cdv_clv*temp_diff0 + veff_u(1, n)*temp_diff + temp_diff4 = veff(3)*clv_u(n)*dcvfz_u_diff + temp_diff1 = veffmag*cdv_clv*dcvfz_u_diff + veff_diff(3) = veff_diff(3) + clv_u(n)*temp_diff1 + + + veffmag_u(n)*temp_diff3 + clv_u_diff(n) = clv_u_diff(n) + veff(3)*temp_diff1 + veffmag_diff = veffmag_diff + cdv_clv*temp_diff4 + veff_u(3 + + , n)*temp_diff3 + cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff4 + veff_u_diff(3, n) = veff_u_diff(3, n) + veffmag*temp_diff3 + veffmag_u_diff(n) = veffmag_u_diff(n) + veff(3)*temp_diff3 + temp_diff3 = cdv*dcvfy_u_diff + temp_diff4 = veff(2)*clv_u(n)*dcvfy_u_diff + temp_diff1 = veffmag*cdv_clv*dcvfy_u_diff + veff_diff(2) = veff_diff(2) + clv_u(n)*temp_diff1 + + + veffmag_u(n)*temp_diff3 + clv_u_diff(n) = clv_u_diff(n) + veff(2)*temp_diff1 + veffmag_diff = veffmag_diff + cdv_clv*temp_diff4 + veff_u(2 + + , n)*temp_diff3 + cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff4 + veff_u_diff(2, n) = veff_u_diff(2, n) + veffmag*temp_diff3 + veffmag_u_diff(n) = veffmag_u_diff(n) + veff(2)*temp_diff3 + temp_diff3 = cdv*dcvfx_u_diff + temp_diff4 = veff(1)*clv_u(n)*dcvfx_u_diff temp_diff1 = veffmag*cdv_clv*dcvfx_u_diff veff_diff(1) = veff_diff(1) + clv_u(n)*temp_diff1 + - + veffmag_u(n)*temp_diff + + veffmag_u(n)*temp_diff3 clv_u_diff(n) = clv_u_diff(n) + veff(1)*temp_diff1 - cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff0 - veff_u_diff(1, n) = veff_u_diff(1, n) + veffmag*temp_diff + veffmag_diff = veffmag_diff + cdv_clv*temp_diff4 + veff_u(1 + + , n)*temp_diff3 + cdv_clv_diff = cdv_clv_diff + veffmag*temp_diff4 + veff_u_diff(1, n) = veff_u_diff(1, n) + veffmag*temp_diff3 + veffmag_u_diff(n) = veffmag_u_diff(n) + veff(1)*temp_diff3 ENDDO udrag_diff(1) = udrag_diff(1) + dcvfx*cdv_lstrp_diff(j) dcvfx_diff = udrag(1)*cdv_lstrp_diff(j) + cfx_diff @@ -3583,30 +3774,30 @@ SUBROUTINE SFFORC_B() dcfz_d = fgam_d(3, n)/sr C C - temp_diff = cmz_d_diff(n)/cr - dcfy_d_diff = r(1)*temp_diff - r_diff(1) = r_diff(1) + dcfy_d*temp_diff - dcfx_d_diff = -(r(2)*temp_diff) - r_diff(2) = r_diff(2) - dcfx_d*temp_diff - cr_diff = cr_diff - (dcfy_d*r(1)-dcfx_d*r(2))*temp_diff/ - + cr - temp_diff = cmy_d_diff(n)/cr - dcfx_d_diff = dcfx_d_diff + r(3)*temp_diff + cfx_d_diff( - + n) - r_diff(3) = r_diff(3) + dcfx_d*temp_diff - dcfz_d_diff = -(r(1)*temp_diff) - r_diff(1) = r_diff(1) - dcfz_d*temp_diff - cr_diff = cr_diff - (dcfx_d*r(3)-dcfz_d*r(1))*temp_diff/ - + cr - temp_diff = cmx_d_diff(n)/cr - dcfz_d_diff = dcfz_d_diff + r(2)*temp_diff + cfz_d_diff( - + n) - r_diff(2) = r_diff(2) + dcfz_d*temp_diff + temp_diff3 = cmz_d_diff(n)/cr + dcfy_d_diff = r(1)*temp_diff3 + r_diff(1) = r_diff(1) + dcfy_d*temp_diff3 + dcfx_d_diff = -(r(2)*temp_diff3) + r_diff(2) = r_diff(2) - dcfx_d*temp_diff3 + cr_diff = cr_diff - (dcfy_d*r(1)-dcfx_d*r(2))*temp_diff3 + + /cr + temp_diff3 = cmy_d_diff(n)/cr + dcfx_d_diff = dcfx_d_diff + r(3)*temp_diff3 + cfx_d_diff + + (n) + r_diff(3) = r_diff(3) + dcfx_d*temp_diff3 + dcfz_d_diff = -(r(1)*temp_diff3) + r_diff(1) = r_diff(1) - dcfz_d*temp_diff3 + cr_diff = cr_diff - (dcfx_d*r(3)-dcfz_d*r(1))*temp_diff3 + + /cr + temp_diff3 = cmx_d_diff(n)/cr + dcfz_d_diff = dcfz_d_diff + r(2)*temp_diff3 + cfz_d_diff + + (n) + r_diff(2) = r_diff(2) + dcfz_d*temp_diff3 dcfy_d_diff = dcfy_d_diff + cfy_d_diff(n) - r(3)* - + temp_diff - r_diff(3) = r_diff(3) - dcfy_d*temp_diff - cr_diff = cr_diff - (dcfz_d*r(2)-dcfy_d*r(3))*temp_diff/ - + cr + + temp_diff3 + r_diff(3) = r_diff(3) - dcfy_d*temp_diff3 + cr_diff = cr_diff - (dcfz_d*r(2)-dcfy_d*r(3))*temp_diff3 + + /cr fgam_d_diff(3, n) = fgam_d_diff(3, n) + dcfz_d_diff/sr sr_diff = sr_diff - fgam_d(3, n)*dcfz_d_diff/sr**2 - + fgam_d(2, n)*dcfy_d_diff/sr**2 - fgam_d(1, n)* @@ -3617,31 +3808,31 @@ SUBROUTINE SFFORC_B() DO n=numax,1,-1 dcfy_u = fgam_u(2, n)/sr dcfx_u = fgam_u(1, n)/sr - temp_diff = cmz_u_diff(n)/cr - dcfy_u_diff = r(1)*temp_diff - r_diff(1) = r_diff(1) + dcfy_u*temp_diff - dcfx_u_diff = -(r(2)*temp_diff) - r_diff(2) = r_diff(2) - dcfx_u*temp_diff - cr_diff = cr_diff - (dcfy_u*r(1)-dcfx_u*r(2))*temp_diff/ - + cr + temp_diff3 = cmz_u_diff(n)/cr + dcfy_u_diff = r(1)*temp_diff3 + r_diff(1) = r_diff(1) + dcfy_u*temp_diff3 + dcfx_u_diff = -(r(2)*temp_diff3) + r_diff(2) = r_diff(2) - dcfx_u*temp_diff3 + cr_diff = cr_diff - (dcfy_u*r(1)-dcfx_u*r(2))*temp_diff3 + + /cr dcfz_u = fgam_u(3, n)/sr - temp_diff = cmy_u_diff(n)/cr - dcfx_u_diff = dcfx_u_diff + r(3)*temp_diff + cfx_u_diff( - + n) - r_diff(3) = r_diff(3) + dcfx_u*temp_diff - dcfz_u_diff = -(r(1)*temp_diff) - r_diff(1) = r_diff(1) - dcfz_u*temp_diff - cr_diff = cr_diff - (dcfx_u*r(3)-dcfz_u*r(1))*temp_diff/ - + cr - temp_diff = cmx_u_diff(n)/cr - dcfz_u_diff = dcfz_u_diff + r(2)*temp_diff + cfz_u_diff( - + n) - r_diff(2) = r_diff(2) + dcfz_u*temp_diff + temp_diff3 = cmy_u_diff(n)/cr + dcfx_u_diff = dcfx_u_diff + r(3)*temp_diff3 + cfx_u_diff + + (n) + r_diff(3) = r_diff(3) + dcfx_u*temp_diff3 + dcfz_u_diff = -(r(1)*temp_diff3) + r_diff(1) = r_diff(1) - dcfz_u*temp_diff3 + cr_diff = cr_diff - (dcfx_u*r(3)-dcfz_u*r(1))*temp_diff3 + + /cr + temp_diff3 = cmx_u_diff(n)/cr + dcfz_u_diff = dcfz_u_diff + r(2)*temp_diff3 + cfz_u_diff + + (n) + r_diff(2) = r_diff(2) + dcfz_u*temp_diff3 dcfy_u_diff = dcfy_u_diff + cfy_u_diff(n) - r(3)* - + temp_diff - r_diff(3) = r_diff(3) - dcfy_u*temp_diff - cr_diff = cr_diff - (dcfz_u*r(2)-dcfy_u*r(3))*temp_diff/ - + cr + + temp_diff3 + r_diff(3) = r_diff(3) - dcfy_u*temp_diff3 + cr_diff = cr_diff - (dcfz_u*r(2)-dcfy_u*r(3))*temp_diff3 + + /cr fgam_u_diff(3, n) = fgam_u_diff(3, n) + dcfz_u_diff/sr sr_diff = sr_diff - fgam_u(3, n)*dcfz_u_diff/sr**2 - + fgam_u(2, n)*dcfy_u_diff/sr**2 - fgam_u(1, n)* @@ -3649,27 +3840,31 @@ SUBROUTINE SFFORC_B() fgam_u_diff(2, n) = fgam_u_diff(2, n) + dcfy_u_diff/sr fgam_u_diff(1, n) = fgam_u_diff(1, n) + dcfx_u_diff/sr ENDDO + temp_diff3 = cmz_diff/cr dcfx = fgam(1)/sr dcfy = fgam(2)/sr - temp_diff = cmz_diff/cr - dcfy_diff = r(1)*temp_diff - r_diff(1) = r_diff(1) + dcfy*temp_diff - dcfx_diff = -(r(2)*temp_diff) - r_diff(2) = r_diff(2) - dcfx*temp_diff - cr_diff = cr_diff - (dcfy*r(1)-dcfx*r(2))*temp_diff/cr dcfz = fgam(3)/sr - temp_diff = cmy_diff/cr - dcfx_diff = dcfx_diff + r(3)*temp_diff + cfx_diff - r_diff(3) = r_diff(3) + dcfx*temp_diff - dcfz_diff = -(r(1)*temp_diff) - r_diff(1) = r_diff(1) - dcfz*temp_diff - cr_diff = cr_diff - (dcfx*r(3)-dcfz*r(1))*temp_diff/cr - temp_diff = cmx_diff/cr - dcfz_diff = dcfz_diff + r(2)*temp_diff + cfz_diff - r_diff(2) = r_diff(2) + dcfz*temp_diff - dcfy_diff = dcfy_diff + cfy_diff - r(3)*temp_diff - r_diff(3) = r_diff(3) - dcfy*temp_diff - cr_diff = cr_diff - (dcfz*r(2)-dcfy*r(3))*temp_diff/cr + cr_diff = cr_diff + (ensy(j)*dcfy+ensz(j)*dcfz)*cnc_diff(j + + ) - (dcfy*r(1)-dcfx*r(2))*temp_diff3/cr + temp_diff5 = cr*cnc_diff(j) + ensy_diff(j) = ensy_diff(j) + dcfy*temp_diff5 + dcfy_diff = ensy(j)*temp_diff5 + r(1)*temp_diff3 + ensz_diff(j) = ensz_diff(j) + dcfz*temp_diff5 + r_diff(1) = r_diff(1) + dcfy*temp_diff3 + dcfx_diff = -(r(2)*temp_diff3) + r_diff(2) = r_diff(2) - dcfx*temp_diff3 + temp_diff3 = cmy_diff/cr + dcfz_diff = ensz(j)*temp_diff5 - r(1)*temp_diff3 + dcfx_diff = dcfx_diff + r(3)*temp_diff3 + cfx_diff + r_diff(3) = r_diff(3) + dcfx*temp_diff3 + r_diff(1) = r_diff(1) - dcfz*temp_diff3 + cr_diff = cr_diff - (dcfx*r(3)-dcfz*r(1))*temp_diff3/cr + temp_diff3 = cmx_diff/cr + dcfz_diff = dcfz_diff + r(2)*temp_diff3 + cfz_diff + r_diff(2) = r_diff(2) + dcfz*temp_diff3 + dcfy_diff = dcfy_diff + cfy_diff - r(3)*temp_diff3 + r_diff(3) = r_diff(3) - dcfy*temp_diff3 + cr_diff = cr_diff - (dcfz*r(2)-dcfy*r(3))*temp_diff3/cr fgam_diff(3) = fgam_diff(3) + dcfz_diff/sr sr_diff = sr_diff - fgam(3)*dcfz_diff/sr**2 - fgam(2)* + dcfy_diff/sr**2 - fgam(1)*dcfx_diff/sr**2 @@ -4022,24 +4217,24 @@ SUBROUTINE SFFORC_B() dcfz_d = fgam_d(3, n)/sr C C - temp_diff = cmz_d_diff(n)/cr - dcfy_d_diff = r(1)*temp_diff - r_diff(1) = r_diff(1) + dcfy_d*temp_diff - dcfx_d_diff = -(r(2)*temp_diff) - r_diff(2) = r_diff(2) - dcfx_d*temp_diff - cr_diff = cr_diff - (dcfy_d*r(1)-dcfx_d*r(2))*temp_diff/cr - temp_diff = cmy_d_diff(n)/cr - dcfx_d_diff = dcfx_d_diff + r(3)*temp_diff + cfx_d_diff(n) - r_diff(3) = r_diff(3) + dcfx_d*temp_diff - dcfz_d_diff = -(r(1)*temp_diff) - r_diff(1) = r_diff(1) - dcfz_d*temp_diff - cr_diff = cr_diff - (dcfx_d*r(3)-dcfz_d*r(1))*temp_diff/cr - temp_diff = cmx_d_diff(n)/cr - dcfz_d_diff = dcfz_d_diff + r(2)*temp_diff + cfz_d_diff(n) - r_diff(2) = r_diff(2) + dcfz_d*temp_diff - dcfy_d_diff = dcfy_d_diff + cfy_d_diff(n) - r(3)*temp_diff - r_diff(3) = r_diff(3) - dcfy_d*temp_diff - cr_diff = cr_diff - (dcfz_d*r(2)-dcfy_d*r(3))*temp_diff/cr + temp_diff3 = cmz_d_diff(n)/cr + dcfy_d_diff = r(1)*temp_diff3 + r_diff(1) = r_diff(1) + dcfy_d*temp_diff3 + dcfx_d_diff = -(r(2)*temp_diff3) + r_diff(2) = r_diff(2) - dcfx_d*temp_diff3 + cr_diff = cr_diff - (dcfy_d*r(1)-dcfx_d*r(2))*temp_diff3/cr + temp_diff3 = cmy_d_diff(n)/cr + dcfx_d_diff = dcfx_d_diff + r(3)*temp_diff3 + cfx_d_diff(n) + r_diff(3) = r_diff(3) + dcfx_d*temp_diff3 + dcfz_d_diff = -(r(1)*temp_diff3) + r_diff(1) = r_diff(1) - dcfz_d*temp_diff3 + cr_diff = cr_diff - (dcfx_d*r(3)-dcfz_d*r(1))*temp_diff3/cr + temp_diff3 = cmx_d_diff(n)/cr + dcfz_d_diff = dcfz_d_diff + r(2)*temp_diff3 + cfz_d_diff(n) + r_diff(2) = r_diff(2) + dcfz_d*temp_diff3 + dcfy_d_diff = dcfy_d_diff + cfy_d_diff(n) - r(3)*temp_diff3 + r_diff(3) = r_diff(3) - dcfy_d*temp_diff3 + cr_diff = cr_diff - (dcfz_d*r(2)-dcfy_d*r(3))*temp_diff3/cr fgam_d_diff(3, n) = fgam_d_diff(3, n) + dcfz_d_diff/sr sr_diff = sr_diff - fgam_d(3, n)*dcfz_d_diff/sr**2 - fgam_d( + 2, n)*dcfy_d_diff/sr**2 - fgam_d(1, n)*dcfx_d_diff/sr**2 @@ -4053,24 +4248,24 @@ SUBROUTINE SFFORC_B() dcfz_u = fgam_u(3, n)/sr C C - temp_diff = cmz_u_diff(n)/cr - dcfy_u_diff = r(1)*temp_diff - r_diff(1) = r_diff(1) + dcfy_u*temp_diff - dcfx_u_diff = -(r(2)*temp_diff) - r_diff(2) = r_diff(2) - dcfx_u*temp_diff - cr_diff = cr_diff - (dcfy_u*r(1)-dcfx_u*r(2))*temp_diff/cr - temp_diff = cmy_u_diff(n)/cr - dcfx_u_diff = dcfx_u_diff + r(3)*temp_diff + cfx_u_diff(n) - r_diff(3) = r_diff(3) + dcfx_u*temp_diff - dcfz_u_diff = -(r(1)*temp_diff) - r_diff(1) = r_diff(1) - dcfz_u*temp_diff - cr_diff = cr_diff - (dcfx_u*r(3)-dcfz_u*r(1))*temp_diff/cr - temp_diff = cmx_u_diff(n)/cr - dcfz_u_diff = dcfz_u_diff + r(2)*temp_diff + cfz_u_diff(n) - r_diff(2) = r_diff(2) + dcfz_u*temp_diff - dcfy_u_diff = dcfy_u_diff + cfy_u_diff(n) - r(3)*temp_diff - r_diff(3) = r_diff(3) - dcfy_u*temp_diff - cr_diff = cr_diff - (dcfz_u*r(2)-dcfy_u*r(3))*temp_diff/cr + temp_diff3 = cmz_u_diff(n)/cr + dcfy_u_diff = r(1)*temp_diff3 + r_diff(1) = r_diff(1) + dcfy_u*temp_diff3 + dcfx_u_diff = -(r(2)*temp_diff3) + r_diff(2) = r_diff(2) - dcfx_u*temp_diff3 + cr_diff = cr_diff - (dcfy_u*r(1)-dcfx_u*r(2))*temp_diff3/cr + temp_diff3 = cmy_u_diff(n)/cr + dcfx_u_diff = dcfx_u_diff + r(3)*temp_diff3 + cfx_u_diff(n) + r_diff(3) = r_diff(3) + dcfx_u*temp_diff3 + dcfz_u_diff = -(r(1)*temp_diff3) + r_diff(1) = r_diff(1) - dcfz_u*temp_diff3 + cr_diff = cr_diff - (dcfx_u*r(3)-dcfz_u*r(1))*temp_diff3/cr + temp_diff3 = cmx_u_diff(n)/cr + dcfz_u_diff = dcfz_u_diff + r(2)*temp_diff3 + cfz_u_diff(n) + r_diff(2) = r_diff(2) + dcfz_u*temp_diff3 + dcfy_u_diff = dcfy_u_diff + cfy_u_diff(n) - r(3)*temp_diff3 + r_diff(3) = r_diff(3) - dcfy_u*temp_diff3 + cr_diff = cr_diff - (dcfz_u*r(2)-dcfy_u*r(3))*temp_diff3/cr fgam_u_diff(3, n) = fgam_u_diff(3, n) + dcfz_u_diff/sr sr_diff = sr_diff - fgam_u(3, n)*dcfz_u_diff/sr**2 - fgam_u( + 2, n)*dcfy_u_diff/sr**2 - fgam_u(1, n)*dcfx_u_diff/sr**2 @@ -4078,24 +4273,28 @@ SUBROUTINE SFFORC_B() fgam_u_diff(1, n) = fgam_u_diff(1, n) + dcfx_u_diff/sr ENDDO CALL POPINTEGER4(n) - temp_diff = cmz_diff/cr - dcfy_diff = r(1)*temp_diff - r_diff(1) = r_diff(1) + dcfy*temp_diff - dcfx_diff = -(r(2)*temp_diff) - r_diff(2) = r_diff(2) - dcfx*temp_diff - cr_diff = cr_diff - (dcfy*r(1)-dcfx*r(2))*temp_diff/cr - temp_diff = cmy_diff/cr - dcfx_diff = dcfx_diff + r(3)*temp_diff + cfx_diff - r_diff(3) = r_diff(3) + dcfx*temp_diff - dcfz_diff = -(r(1)*temp_diff) - r_diff(1) = r_diff(1) - dcfz*temp_diff - cr_diff = cr_diff - (dcfx*r(3)-dcfz*r(1))*temp_diff/cr - temp_diff = cmx_diff/cr - dcfz_diff = dcfz_diff + r(2)*temp_diff + cfz_diff - r_diff(2) = r_diff(2) + dcfz*temp_diff - dcfy_diff = dcfy_diff + cfy_diff - r(3)*temp_diff - r_diff(3) = r_diff(3) - dcfy*temp_diff - cr_diff = cr_diff - (dcfz*r(2)-dcfy*r(3))*temp_diff/cr + temp_diff3 = cmz_diff/cr + cr_diff = cr_diff + (ensy(j)*dcfy+ensz(j)*dcfz)*cnc_diff(j) - + + (dcfy*r(1)-dcfx*r(2))*temp_diff3/cr + temp_diff5 = cr*cnc_diff(j) + ensy_diff(j) = ensy_diff(j) + dcfy*temp_diff5 + dcfy_diff = ensy(j)*temp_diff5 + r(1)*temp_diff3 + ensz_diff(j) = ensz_diff(j) + dcfz*temp_diff5 + r_diff(1) = r_diff(1) + dcfy*temp_diff3 + dcfx_diff = -(r(2)*temp_diff3) + r_diff(2) = r_diff(2) - dcfx*temp_diff3 + temp_diff3 = cmy_diff/cr + dcfz_diff = ensz(j)*temp_diff5 - r(1)*temp_diff3 + dcfx_diff = dcfx_diff + r(3)*temp_diff3 + cfx_diff + r_diff(3) = r_diff(3) + dcfx*temp_diff3 + r_diff(1) = r_diff(1) - dcfz*temp_diff3 + cr_diff = cr_diff - (dcfx*r(3)-dcfz*r(1))*temp_diff3/cr + temp_diff3 = cmx_diff/cr + dcfz_diff = dcfz_diff + r(2)*temp_diff3 + cfz_diff + r_diff(2) = r_diff(2) + dcfz*temp_diff3 + dcfy_diff = dcfy_diff + cfy_diff - r(3)*temp_diff3 + r_diff(3) = r_diff(3) - dcfy*temp_diff3 + cr_diff = cr_diff - (dcfz*r(2)-dcfy*r(3))*temp_diff3/cr fgam_diff(3) = fgam_diff(3) + dcfz_diff/sr sr_diff = sr_diff - fgam(3)*dcfz_diff/sr**2 - fgam(2)* + dcfy_diff/sr**2 - fgam(1)*dcfx_diff/sr**2 @@ -4326,6 +4525,7 @@ SUBROUTINE SFFORC_B() cfy_u_diff(n) = 0.D0 cfx_u_diff(n) = 0.D0 ENDDO + cnc_diff(j) = 0.D0 rle_diff(3, j) = rle_diff(3, j) + rc4_diff(3) rc4_diff(3) = 0.D0 rle_diff(2, j) = rle_diff(2, j) + rc4_diff(2) diff --git a/src/ad_src/reverse_ad_src/atpforc_b.f b/src/ad_src/reverse_ad_src/atpforc_b.f index ca467866..68b3db61 100644 --- a/src/ad_src/reverse_ad_src/atpforc_b.f +++ b/src/ad_src/reverse_ad_src/atpforc_b.f @@ -3,8 +3,9 @@ C C Differentiation of tpforc in reverse (adjoint) mode (with options i4 dr8 r8): C gradient of useful results: bref clff cyff cdff spanef -C with respect to varying inputs: sref bref chord rv1 rv2 rc -C gam +C dwwake +C with respect to varying inputs: sref bref chord dwwake rv1 +C rv2 rc gam C*********************************************************************** C Module: atpforc.f C @@ -30,6 +31,7 @@ SUBROUTINE TPFORC_B() INCLUDE 'AVL_ad_seeds.inc' C REAL ny, nz + REAL ny_diff, nz_diff REAL vy_u(numax), vz_u(numax), vy_d(ndmax), vz_d(ndmax), vy_g( + ngmax), vz_g(ngmax) REAL p(3, 3), p_m(3, 3), p_a(3, 3), p_b(3, 3) @@ -55,6 +57,7 @@ SUBROUTINE TPFORC_B() REAL dzt REAL dzt_diff REAL dst + REAL dst_diff INTRINSIC SQRT REAL ycntr REAL ycntr_diff @@ -89,9 +92,9 @@ SUBROUTINE TPFORC_B() REAL spanef_cy REAL spanef_cd REAL temp + REAL temp_diff REAL temp0 REAL temp1 - REAL temp_diff REAL temp_diff0 REAL temp_diff1 REAL(kind=8) temp2 @@ -150,8 +153,12 @@ SUBROUTINE TPFORC_B() C...Find the normal velocity across each strip at the projected control C point location DO jc=1,nstrip + CALL PUSHREAL8(dyt) dyt = rt2(2, jc) - rt1(2, jc) + CALL PUSHREAL8(dzt) dzt = rt2(3, jc) - rt1(3, jc) + CALL PUSHREAL8(dst) + dst = SQRT(dyt*dyt + dzt*dzt) C ycntr = rtc(2, jc) zcntr = rtc(3, jc) @@ -325,17 +332,17 @@ SUBROUTINE TPFORC_B() dzt = rt2(3, jc) - rt1(3, jc) temp2 = gams(jc)/sref temp_diff2 = (dzt*vy-dyt*vz)*cdff_diff/sref - temp_diff = temp2*cdff_diff - vy_diff = dzt*temp_diff - vz_diff = -(dyt*temp_diff) + temp_diff0 = temp2*cdff_diff + vy_diff = dzt*temp_diff0 + vz_diff = -(dyt*temp_diff0) gams_diff(jc) = gams_diff(jc) + temp_diff2 + dyt*2.0*clff_diff + /sref - dzt*2.0*cyff_diff/sref sref_diff = sref_diff - temp2*temp_diff2 temp_diff2 = -(gams(jc)*2.0*cyff_diff/sref) - dzt_diff = vy*temp_diff + temp_diff2 + dzt_diff = vy*temp_diff0 + temp_diff2 sref_diff = sref_diff - dzt*temp_diff2/sref temp_diff2 = gams(jc)*2.0*clff_diff/sref - dyt_diff = temp_diff2 - vz*temp_diff + dyt_diff = temp_diff2 - vz*temp_diff0 sref_diff = sref_diff - dyt*temp_diff2/sref ELSE dyt_diff = 0.D0 @@ -343,6 +350,13 @@ SUBROUTINE TPFORC_B() vz_diff = 0.D0 dzt_diff = 0.D0 END IF + ny = -(dzt/dst) + nz = dyt/dst + ny_diff = -(vy*dwwake_diff(jc)) + vy_diff = vy_diff - ny*dwwake_diff(jc) + nz_diff = -(vz*dwwake_diff(jc)) + vz_diff = vz_diff - nz*dwwake_diff(jc) + dwwake_diff(jc) = 0.D0 ycntr = rtc(2, jc) zcntr = rtc(3, jc) ycntr_diff = 0.D0 @@ -356,22 +370,22 @@ SUBROUTINE TPFORC_B() dz2 = zcntr - (zoff-rt2(3, jv)) rsq1 = dy1*dy1 + dz1*dz1 rsq2 = dy2*dy2 + dz2*dz2 - temp_diff = hpi*iysym*izsym*vz_diff + temp_diff0 = hpi*iysym*izsym*vz_diff gams_diff(jv) = gams_diff(jv) + (dy2/rsq2-dy1/rsq1)* - + temp_diff - temp_diff0 = gams(jv)*temp_diff - dy2_diff = temp_diff0/rsq2 - rsq2_diff = -(dy2*temp_diff0/rsq2**2) - dy1_diff = -(temp_diff0/rsq1) - rsq1_diff = dy1*temp_diff0/rsq1**2 - temp_diff = hpi*iysym*izsym*vy_diff + + temp_diff0 + temp_diff1 = gams(jv)*temp_diff0 + dy2_diff = temp_diff1/rsq2 + rsq2_diff = -(dy2*temp_diff1/rsq2**2) + dy1_diff = -(temp_diff1/rsq1) + rsq1_diff = dy1*temp_diff1/rsq1**2 + temp_diff0 = hpi*iysym*izsym*vy_diff gams_diff(jv) = gams_diff(jv) + (dz1/rsq1-dz2/rsq2)* - + temp_diff - temp_diff0 = gams(jv)*temp_diff - rsq1_diff = rsq1_diff - dz1*temp_diff0/rsq1**2 - dz1_diff = temp_diff0/rsq1 + 2*dz1*rsq1_diff - rsq2_diff = rsq2_diff + dz2*temp_diff0/rsq2**2 - dz2_diff = 2*dz2*rsq2_diff - temp_diff0/rsq2 + + temp_diff0 + temp_diff1 = gams(jv)*temp_diff0 + rsq1_diff = rsq1_diff - dz1*temp_diff1/rsq1**2 + dz1_diff = temp_diff1/rsq1 + 2*dz1*rsq1_diff + rsq2_diff = rsq2_diff + dz2*temp_diff1/rsq2**2 + dz2_diff = 2*dz2*rsq2_diff - temp_diff1/rsq2 dy2_diff = dy2_diff + 2*dy2*rsq2_diff dy1_diff = dy1_diff + 2*dy1*rsq1_diff zcntr_diff = zcntr_diff + dz2_diff + dz1_diff @@ -389,20 +403,20 @@ SUBROUTINE TPFORC_B() dz2 = zcntr - rt2(3, jv) rsq1 = dy1*dy1 + dz1*dz1 rsq2 = dy2*dy2 + dz2*dz2 - temp_diff = -(hpi*iysym*vz_diff) - gams_diff(jv) = gams_diff(jv) + (dy2/rsq2-dy1/rsq1)*temp_diff - temp_diff0 = gams(jv)*temp_diff - dy2_diff = temp_diff0/rsq2 - rsq2_diff = -(dy2*temp_diff0/rsq2**2) - dy1_diff = -(temp_diff0/rsq1) - rsq1_diff = dy1*temp_diff0/rsq1**2 - temp_diff = -(hpi*iysym*vy_diff) - gams_diff(jv) = gams_diff(jv) + (dz1/rsq1-dz2/rsq2)*temp_diff - temp_diff0 = gams(jv)*temp_diff - rsq1_diff = rsq1_diff - dz1*temp_diff0/rsq1**2 - dz1_diff = temp_diff0/rsq1 + 2*dz1*rsq1_diff - rsq2_diff = rsq2_diff + dz2*temp_diff0/rsq2**2 - dz2_diff = 2*dz2*rsq2_diff - temp_diff0/rsq2 + temp_diff0 = -(hpi*iysym*vz_diff) + gams_diff(jv) = gams_diff(jv) + (dy2/rsq2-dy1/rsq1)*temp_diff0 + temp_diff1 = gams(jv)*temp_diff0 + dy2_diff = temp_diff1/rsq2 + rsq2_diff = -(dy2*temp_diff1/rsq2**2) + dy1_diff = -(temp_diff1/rsq1) + rsq1_diff = dy1*temp_diff1/rsq1**2 + temp_diff0 = -(hpi*iysym*vy_diff) + gams_diff(jv) = gams_diff(jv) + (dz1/rsq1-dz2/rsq2)*temp_diff0 + temp_diff1 = gams(jv)*temp_diff0 + rsq1_diff = rsq1_diff - dz1*temp_diff1/rsq1**2 + dz1_diff = temp_diff1/rsq1 + 2*dz1*rsq1_diff + rsq2_diff = rsq2_diff + dz2*temp_diff1/rsq2**2 + dz2_diff = 2*dz2*rsq2_diff - temp_diff1/rsq2 CALL POPREAL8(rsq2) dy2_diff = dy2_diff + 2*dy2*rsq2_diff CALL POPREAL8(rsq1) @@ -421,22 +435,22 @@ SUBROUTINE TPFORC_B() dz2 = zcntr - (zoff-rt2(3, jv)) rsq1 = dy1*dy1 + dz1*dz1 rsq2 = dy2*dy2 + dz2*dz2 - temp_diff = -(hpi*izsym*vz_diff) + temp_diff0 = -(hpi*izsym*vz_diff) gams_diff(jv) = gams_diff(jv) + (dy2/rsq2-dy1/rsq1)* - + temp_diff - temp_diff0 = gams(jv)*temp_diff - dy2_diff = temp_diff0/rsq2 - rsq2_diff = -(dy2*temp_diff0/rsq2**2) - dy1_diff = -(temp_diff0/rsq1) - rsq1_diff = dy1*temp_diff0/rsq1**2 - temp_diff = -(hpi*izsym*vy_diff) + + temp_diff0 + temp_diff1 = gams(jv)*temp_diff0 + dy2_diff = temp_diff1/rsq2 + rsq2_diff = -(dy2*temp_diff1/rsq2**2) + dy1_diff = -(temp_diff1/rsq1) + rsq1_diff = dy1*temp_diff1/rsq1**2 + temp_diff0 = -(hpi*izsym*vy_diff) gams_diff(jv) = gams_diff(jv) + (dz1/rsq1-dz2/rsq2)* - + temp_diff - temp_diff0 = gams(jv)*temp_diff - rsq1_diff = rsq1_diff - dz1*temp_diff0/rsq1**2 - dz1_diff = temp_diff0/rsq1 + 2*dz1*rsq1_diff - rsq2_diff = rsq2_diff + dz2*temp_diff0/rsq2**2 - dz2_diff = 2*dz2*rsq2_diff - temp_diff0/rsq2 + + temp_diff0 + temp_diff1 = gams(jv)*temp_diff0 + rsq1_diff = rsq1_diff - dz1*temp_diff1/rsq1**2 + dz1_diff = temp_diff1/rsq1 + 2*dz1*rsq1_diff + rsq2_diff = rsq2_diff + dz2*temp_diff1/rsq2**2 + dz2_diff = 2*dz2*rsq2_diff - temp_diff1/rsq2 CALL POPREAL8(rsq2) dy2_diff = dy2_diff + 2*dy2*rsq2_diff CALL POPREAL8(rsq1) @@ -454,38 +468,38 @@ SUBROUTINE TPFORC_B() dy2 = ycntr - rt2(2, jv) gams_diff(jv) = gams_diff(jv) + (dy2/rsq2-dy1/rsq1)*hpi* + vz_diff + (dz1/rsq1-dz2/rsq2)*hpi*vy_diff - temp_diff = gams(jv)*hpi*vz_diff - dy2_diff = temp_diff/rsq2 - rsq2_diff = -(dy2*temp_diff/rsq2**2) - dy1_diff = -(temp_diff/rsq1) - rsq1_diff = dy1*temp_diff/rsq1**2 - temp_diff = gams(jv)*hpi*vy_diff - dz1_diff = temp_diff/rsq1 - rsq1_diff = rsq1_diff - dz1*temp_diff/rsq1**2 - dz2_diff = -(temp_diff/rsq2) - rsq2_diff = rsq2_diff + dz2*temp_diff/rsq2**2 + temp_diff0 = gams(jv)*hpi*vz_diff + dy2_diff = temp_diff0/rsq2 + rsq2_diff = -(dy2*temp_diff0/rsq2**2) + dy1_diff = -(temp_diff0/rsq1) + rsq1_diff = dy1*temp_diff0/rsq1**2 + temp_diff0 = gams(jv)*hpi*vy_diff + dz1_diff = temp_diff0/rsq1 + rsq1_diff = rsq1_diff - dz1*temp_diff0/rsq1**2 + dz2_diff = -(temp_diff0/rsq2) + rsq2_diff = rsq2_diff + dz2*temp_diff0/rsq2**2 CALL POPREAL8(rsq2) temp1 = dy2*dy2 + dz2*dz2 IF (temp1**2 + rcore**4 .EQ. 0.D0) THEN - temp_diff0 = 0.D0 + temp_diff1 = 0.D0 ELSE - temp_diff0 = rsq2_diff/(2.0*SQRT(temp1**2+rcore**4)) + temp_diff1 = rsq2_diff/(2.0*SQRT(temp1**2+rcore**4)) END IF - temp_diff = 2*temp1*temp_diff0 - rcore_diff = 4*rcore**3*temp_diff0 - dy2_diff = dy2_diff + 2*dy2*temp_diff - dz2_diff = dz2_diff + 2*dz2*temp_diff + temp_diff0 = 2*temp1*temp_diff1 + rcore_diff = 4*rcore**3*temp_diff1 + dy2_diff = dy2_diff + 2*dy2*temp_diff0 + dz2_diff = dz2_diff + 2*dz2*temp_diff0 CALL POPREAL8(rsq1) temp1 = dy1*dy1 + dz1*dz1 IF (temp1**2 + rcore**4 .EQ. 0.D0) THEN - temp_diff0 = 0.D0 + temp_diff1 = 0.D0 ELSE - temp_diff0 = rsq1_diff/(2.0*SQRT(temp1**2+rcore**4)) + temp_diff1 = rsq1_diff/(2.0*SQRT(temp1**2+rcore**4)) END IF - temp_diff = 2*temp1*temp_diff0 - rcore_diff = rcore_diff + 4*rcore**3*temp_diff0 - dy1_diff = dy1_diff + 2*dy1*temp_diff - dz1_diff = dz1_diff + 2*dz1*temp_diff + temp_diff0 = 2*temp1*temp_diff1 + rcore_diff = rcore_diff + 4*rcore**3*temp_diff1 + dy1_diff = dy1_diff + 2*dy1*temp_diff0 + dz1_diff = dz1_diff + 2*dz1*temp_diff0 zcntr_diff = zcntr_diff + dz2_diff + dz1_diff rt2_diff(3, jv) = rt2_diff(3, jv) - dz2_diff rt1_diff(3, jv) = rt1_diff(3, jv) - dz1_diff @@ -507,23 +521,34 @@ SUBROUTINE TPFORC_B() temp = rt2(3, jv) - rt1(3, jv) temp0 = rt2(2, jv) - rt1(2, jv) IF (temp0**2 + temp**2 .EQ. 0.D0) THEN - temp_diff = 0.D0 + temp_diff0 = 0.D0 ELSE - temp_diff = dsyz_diff/(2.0*SQRT(temp0**2+temp**2)) + temp_diff0 = dsyz_diff/(2.0*SQRT(temp0**2+temp**2)) END IF - temp_diff0 = 2*temp0*temp_diff - temp_diff1 = 2*temp*temp_diff - rt2_diff(3, jv) = rt2_diff(3, jv) + temp_diff1 - rt1_diff(3, jv) = rt1_diff(3, jv) - temp_diff1 - rt2_diff(2, jv) = rt2_diff(2, jv) + temp_diff0 - rt1_diff(2, jv) = rt1_diff(2, jv) - temp_diff0 + temp_diff1 = 2*temp0*temp_diff0 + temp_diff = 2*temp*temp_diff0 + rt2_diff(3, jv) = rt2_diff(3, jv) + temp_diff + rt1_diff(3, jv) = rt1_diff(3, jv) - temp_diff + rt2_diff(2, jv) = rt2_diff(2, jv) + temp_diff1 + rt1_diff(2, jv) = rt1_diff(2, jv) - temp_diff1 ENDDO + dst_diff = dzt*ny_diff/dst**2 - dyt*nz_diff/dst**2 + IF (dyt**2 + dzt**2 .EQ. 0.D0) THEN + temp_diff = 0.D0 + ELSE + temp_diff = dst_diff/(2.0*SQRT(dyt**2+dzt**2)) + END IF CALL POPREAL8(vz) CALL POPREAL8(vy) rtc_diff(3, jc) = rtc_diff(3, jc) + zcntr_diff rtc_diff(2, jc) = rtc_diff(2, jc) + ycntr_diff + dyt_diff = dyt_diff + nz_diff/dst + 2*dyt*temp_diff + dzt_diff = dzt_diff + 2*dzt*temp_diff - ny_diff/dst + CALL POPREAL8(dst) + CALL POPREAL8(dzt) rt2_diff(3, jc) = rt2_diff(3, jc) + dzt_diff rt1_diff(3, jc) = rt1_diff(3, jc) - dzt_diff + CALL POPREAL8(dyt) rt2_diff(2, jc) = rt2_diff(2, jc) + dyt_diff rt1_diff(2, jc) = rt1_diff(2, jc) - dyt_diff ENDDO diff --git a/tests/test_body_axis_derivs_partial_derivs.py b/tests/test_body_axis_derivs_partial_derivs.py index 73594027..a471e97e 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 f1a09393..bb82ec26 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 2a589e18..8c3ed8bc 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") + 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-7, + 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 d1527ca3..4869b845 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 00000000..1e802d6d --- /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 e404a845..4491f1c7 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 5d562ada..e7b87c3b 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 df1bec44..00000000 Binary files a/tests/wing_mesh.npy and /dev/null differ