Implement Adjoint for custom meshes

_temp/scratch.py

@@ -0,0 +1,102 @@
+import numpy as np
+
+from optvl import OVLSolver
+import pickle
+import os
+
+
+
+mesh = np.load("wing_mesh.npy")
+
+
+with open("wing_mesh.pkl", 'wb') as f:
+    pickle.dump(mesh, f)
+
+
+
+test_mesh = np.zeros((2,2,3))
+test_mesh[1,:,0] = 1.0
+test_mesh[:,1,1] = 1.0
+test_mesh[:,:,2] = 1.0
+
+
+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')
+
+
+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(test_mesh), # (nx,ny,3) numpy array containing mesh coordinates
+        # Control Surface Specification
+        "control_assignments": {
+            "Elevator" : {"assignment":np.arange(0,test_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
+}
+
+
+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)
+
+
+
+

config/defaults/config.MACOS_GFORTRAN.mk

@@ -6,13 +6,11 @@ FF90       = gfortran
 
 
 
-FF90_FLAGS = -fdefault-real-8 -fdefault-double-8 -O2 -fPIC -Wno-align-commons -std=legacy -C -mmacosx-version-min=13.6 
-# FF90_FLAGS = -fdefault-real-8 -fdefault -double-8 -O0 -fPIC -Wno-align-commons -Werror=line-truncation -std=legacy -C -mmacosx-version-min=13.6 -g -fcheck=bounds -finit-real=snan -finit-integer=-999999 -ftrapping-math   -ftrapv
+# FF90_FLAGS = -fdefault-real-8 -fdefault-double-8 -O2 -fPIC -Wno-align-commons -std=legacy -C -mmacosx-version-min=13.6 
+FF90_FLAGS = -fdefault-real-8 -fdefault -double-8 -O0 -fPIC -Wno-align-commons -Werror=line-truncation -std=legacy -C -mmacosx-version-min=13.6 -g -fcheck=bounds -finit-real=snan -finit-integer=-999999 -ftrapping-math   -ftrapv
 
 
 C_FLAGS = -O2 -fPIC -mmacosx-version-min=13.6
-
-F2PY = f2py
 F2PY_FF90 = gfortran
 
 PYTHON = python

geom_files/rect_out.avl

@@ -0,0 +1,41 @@
+# generated using OptVL v1.4.2.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

optvl/om_wrapper.py

@@ -117,7 +117,16 @@ def add_ovl_geom_vars(self, ovl, add_as="inputs", include_airfoil_geom=False):
             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()
+
+    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
 
@@ -331,7 +340,7 @@ 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, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = (
+                con_seeds, geom_seeds, mesh_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
                     )
@@ -601,7 +610,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
                         # print(var_name, body_axis_seeds[func_key])
                         print(f"  running rev mode derivs for {func_key}")
 
-            con_seeds, geom_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = (
+            con_seeds, geom_seeds, mesh_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = (
                 self.ovl._execute_jac_vec_prod_rev(
                     func_seeds=func_seeds,
                     consurf_derivs_seeds=csd_seeds,
@@ -713,14 +722,18 @@ class OVLMeshReader(om.ExplicitComponent):
     def initialize(self):
         self.options.declare("geom_file", types=str)
         self.options.declare("mass_file", default=None)
+        self.options.declare("mesh_output",default=False)
 
     def setup(self):
         geom_file = self.options["geom_file"]
         mass_file = self.options["mass_file"]
+        mesh_output = self.options["mesh_output"]
 
         avl = OVLSolver(geo_file=geom_file, mass_file=mass_file, debug=False)
         add_ovl_geom_vars(self, avl, add_as="outputs", include_airfoil_geom=True)
 
+        if mesh_output:
+            add_ovl_mesh_out_as_output(self,avl)
 
 class Differencer(om.ExplicitComponent):
     def setup(self):