fix derivs for geoms with bodies

docs/installation.md

@@ -6,7 +6,7 @@ The recommended method to install `optvl` is via pip:
 ```shell
 pip install optvl
 ```
-This package is packaged with an OpenBLAS linear solver for quicker analysis.
+This python package is comes with an OpenBLAS linear solver for quicker analysis.
 
 ### Supported Platforms
 Currently, `optvl` supports Linux, macOS (Apple Silicon and Intel), and Windows!

examples/run_cl_alpha_test.py

@@ -0,0 +1,58 @@
+from optvl import OVLSolver
+import numpy as np
+
+
+def finite_diff(ovl, var, step=1e-6):
+    forces = ovl.get_total_forces()
+    val = ovl.get_variable(var)
+    ovl.set_variable(var, val + step)
+    ovl.execute_run()
+    dforces = ovl.get_total_forces()
+
+    for key in dforces:
+        dforces[key] = (dforces[key] - forces[key]) / step
+
+    return dforces
+
+
+for geom in ["../geom_files/rectangle.avl", "../geom_files/aircraft.avl", "../geom_files/supra.avl"]:
+    ovl = OVLSolver(geo_file=geom, debug=False)
+
+    alpha = 10.0
+    ovl.set_variable('alpha', alpha)
+    ovl.execute_run()
+
+    forces = ovl.get_total_forces()
+    R = np.array([[np.sin(np.pi/180*alpha), 0, -1*np.cos(np.pi/180*alpha)]])
+
+    CF = np.array([[forces['CX'], forces['CY'], forces['CZ']]]).T
+
+    # Compute dCL/dalpha using chain rule:
+    # CL(alpha) = R(alpha) @ CF(alpha)
+    # dCL/dalpha = dR/dalpha @ CF + R @ dCF/dalpha
+
+    dR_dalpha = np.array([[np.cos(np.pi/180*alpha)*np.pi/180, 0, np.sin(np.pi/180*alpha)*np.pi/180]])
+
+    sens = ovl.execute_run_sensitivities(['CL', 'CX', 'CY', 'CZ'])
+    dCF_dalpha = np.array([[sens['CX']['alpha'], sens['CY']['alpha'], sens['CZ']['alpha']]]).T
+
+    body_derivs = ovl.get_body_axis_derivs()
+    dCF_dV = np.zeros((3, 6))
+    for idx_vel, vel_comp in enumerate(["u", "v", "w", "p", "q", "r"]):
+        for idx_f, f_comp in enumerate(["X", "Y", "Z"]):
+            dCF_dV[idx_f, idx_vel] = body_derivs[f"dC{f_comp}/d{vel_comp}"]
+
+    dV_dalpha = np.array([[-np.sin(np.pi/180*alpha)*np.pi/180, 0, np.cos(np.pi/180*alpha)*np.pi/180, 0, 0, 0]]).T
+
+    stab_derivs = ovl.get_stab_derivs()
+    dforces = finite_diff(ovl, 'alpha', step=1e-8)
+    dCL_dalpha_avl = stab_derivs['dCL/dalpha'] * np.pi/180  # stab derivs use radians, others use degrees
+    dCL_dalpha_p_ad = (dR_dalpha @ CF + R @ dCF_dalpha)[0, 0]
+    dCL_dalpha_analytic = (dR_dalpha @ CF + R @ dCF_dV @ dV_dalpha)[0, 0]
+
+    print(f'--- comparison for {geom} ----')
+    print('AVL    ', dCL_dalpha_avl)
+    print('full AD', sens['CL']['alpha'], (sens['CL']['alpha'] - dCL_dalpha_avl) / dCL_dalpha_avl)
+    print('part AD', dCL_dalpha_p_ad, (dCL_dalpha_p_ad - dCL_dalpha_avl) / dCL_dalpha_avl)
+    print('hand   ', dCL_dalpha_analytic, (dCL_dalpha_analytic - dCL_dalpha_avl) / dCL_dalpha_avl)
+    print('FD     ', dforces['CL'], (dforces['CL'] - dCL_dalpha_avl) / dCL_dalpha_avl)

geom_files/aircraft_L1_with_body.avl

@@ -0,0 +1,83 @@
+MACH MDAO AVL
+
+#======================================================
+#------------------- Geometry File --------------------
+#======================================================
+# AVL Conventions
+# SI Used: m, kg, etc
+
+#Mach
+0.0
+#IYsym   IZsym   Zsym
+ 0       0       0
+#Sref    Cref    b_wing
+1.13  0.4  1
+#Xref    Yref    Zref
+0.0 0.0 0.0
+
+
+#======================================================
+#--------------------- Main Wing ----------------------
+#======================================================
+SURFACE
+Wing
+#Nchordwise  Cspace   [Nspan   Sspace]
+     2        1.0     2      0.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.0    0.0    0       0.4     0.0      5      3
+
+SECTION
+#Xle  Yle  Zle  |  Chord   Ainc   Nspan   Sspace
+0.0    1.0    0       0.14    0.0      5      3
+
+
+
+#======================================================
+#------------------- Horizontal Tail ------------------
+#======================================================
+SURFACE
+Horizontal Tail
+#Nchordwise  Cspace   Nspan   Sspace
+10       1.0           1         2 
+
+#------------------TAIL ROOT/ELEVATOR------------------
+SECTION
+#Xle  Yle  Zle  |  Chord   Ainc   Nspan   Sspace
+1.5    0.0    0       0.325     0.0       
+
+CONTROL
+#surface gain xhinge hvec SgnDup
+Elevator -1.00 0.0 0 1 0 1.00
+
+SECTION
+#Xle  Yle  Zle  |  Chord   Ainc   Nspan   Sspace
+1.5    0.5    0       0.325     0.0       
+
+CONTROL
+#surface gain xhinge hvec SgnDup
+Elevator -1.00 0.0 0 1 0 1.00
+
+
+BODY
+Fuse pod
+4   2.0
+#
+TRANSLATE
+# 0.0  0.0  -1.75
+0.0  0.0  0.0
+#
+BFIL
+../geom_files/fuseSimple.dat
+
+
+# -- END OF FILE --

geom_files/fuseSimple.dat

@@ -0,0 +1,6 @@
+_F3J Fuse                       
+1.0 -0.5
+0.5 -0.25
+0.0 -0.5
+0.5 -0.75
+1.0 -0.5

geom_files/rect_with_body.avl

@@ -0,0 +1,59 @@
+MACH MDAO AVL
+#======================================================
+#------------------- Geometry File --------------------
+#======================================================
+# AVL Conventions
+# SI Used: m, kg, etc
+
+#Mach
+0.0
+#IYsym   IZsym   Zsym
+ 0       0       0
+#Sref    Cref    b_wing
+1.123    0.25    6.01
+#Xref    Yref    Zref
+0,0 0.0 0.0
+# CDp
+0.00
+
+#======================================================
+#--------------------- Main Wing ----------------------
+#======================================================
+SURFACE
+Wing
+#Nchordwise  Cspace   [Nspan   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.00    0.0    1.0     1.0     0.0
+CONTROL
+#surface gain xhinge hvec SgnDup
+Elevator -1.00 0.5 0 1 0 1.00
+
+SECTION
+#Xle  Yle  Zle  |  Chord   Ainc   Nspan   Sspace
+0.00    1.0    1.0     1.0     0.0
+CONTROL
+#surface gain xhinge hvec SgnDup
+Elevator -1.00 0.5 0 1 0 1.00
+
+# -- END OF FILE --
+# #=============================================
+BODY
+Fuse pod
+4   2.0
+#
+TRANSLATE
+# 0.0  0.0  -1.75
+0.0  0.0  0.0
+#
+BFIL
+../geom_files/fuseSimple.dat

meson.build

@@ -2,7 +2,7 @@ project(
   'optvl',
   'c',
   # don't forget to change the value in pyproject.toml
-  version: '2.2.0',
+  version: '2.2.1',
   license: 'GPL-3.0',
   meson_version: '>= 0.64.0',
   default_options: [

pyproject.toml

@@ -28,7 +28,8 @@ dependencies = [
     "numpy>=1.19",
 ]
 readme = "README.md"
-version = "2.2.0" # this automatically updates __init__.py and setup_deprecated.py
+version = "2.2.1" # this automatically updates __init__.py and setup_deprecated.py
+# don't forget to change the value in meson.build!
 
 
 [tool.cibuildwheel]

src/ad_src/Makefile_tapenade

@@ -28,10 +28,10 @@ PP_FILES = $(addprefix $(PP_DIR)/,$(notdir $(ALL_RES_FILES)))
 fullRoutines = "\
 	update_surfaces(XYZSCAL,XYZTRAN,ADDINC,XYZLES,CHORDS,AINCS,XASEC,SASEC,TASEC,CLCDSEC,CLAF)>(ENC,ENV, DXV, CHORDV,CHORD, CHORD1, CHORD2, RLE,RLE1,RLE2, WSTRIP, RV1,RV2,RV,RC,RS,RL, ENSY, ENSZ, ESS, XSREF,YSREF,ZSREF, ENC_D)\
 	\
-	get_res(GAM, GAM_D, GAM_U, CONVAL, PARVAL, YSYM, ZSYM, ENC, ENV, DXV, CHORDV, RV1, RV2, RV, RC, RS, RL, ENC_D)>(RES, RES_D, RES_U, GAM, GAM_D, GAM_U, VINF, WROT, VINF_A, VINF_B, ALFA, BETA, DELCON, RV1,RV2, RV, RC,DXV, XYZREF, WV_GAM, CDREF, MACH)\
+	get_res(GAM, GAM_D, GAM_U, CONVAL, PARVAL, YSYM, ZSYM, ENC, ENV, DXV, CHORDV, RV1, RV2, RV, RC, RS, RL, ENC_D, XYZREF)>(RES, RES_D, RES_U, GAM, GAM_D, GAM_U, VINF, WROT, VINF_A, VINF_B, ALFA, BETA, DELCON, RV1,RV2, RV, RC,DXV, XYZREF, WV_GAM, CDREF, MACH, SRC, SRC_U)\
 	\
 	VElSUM(GAM, GAM_U, GAM_D, VINF, WROT, WV_GAM)>(VV, VV_U, VV_D, WV, WV_U, WV_D, GAM, GAM_U, GAM_D, VINF, WROT, RV, RV1, RV2) \
-	AERO(GAM, GAM_U, GAM_D, VV, VV_U, VV_D, WV, WV_U, WV_D, RC, RV, RV1, RV2, RLE, XYZREF, CHORD, CHORD1, CHORD2, RLE1, RLE2, WSTRIP, ALFA, MACH, VINF, WROT, VINF_A, VINF_B, ENSZ, ENSY, ESS, SREF, CREF, BREF, XSREF,YSREF,ZSREF,DXV, CDREF)>(\
+	AERO(GAM, GAM_U, GAM_D, VV, VV_U, VV_D, WV, WV_U, WV_D, RC, RV, RV1, RV2, RLE, XYZREF, CHORD, CHORD1, CHORD2, RLE1, RLE2, WSTRIP, ALFA, MACH, VINF, WROT, VINF_A, VINF_B, ENSZ, ENSY, ESS, SREF, CREF, BREF, XSREF,YSREF,ZSREF,DXV, CDREF, SRC, SRC_U)>(\
 	CLTOT, CYTOT, CDTOT, CDVTOT, CDITOT\
 	CLFF, CYFF, CDFF,\
 	CFTOT, \