diff --git a/geom_files/make_pkls/make_aircraft_L1.py b/geom_files/make_pkls/make_aircraft_L1.py new file mode 100644 index 0000000..3c0e145 --- /dev/null +++ b/geom_files/make_pkls/make_aircraft_L1.py @@ -0,0 +1,130 @@ +# ============================================================================= +# Extension modules +# ============================================================================= +from optvl import OVLSolver + +# ============================================================================= +# Standard Python Modules +# ============================================================================= +import os + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np +import pickle +from openaerostruct.meshing.mesh_generator import generate_mesh +from openaerostruct.geometry.utils import taper + + +mesh_dict = { + "num_y": 5, + "num_x": 3, + "wing_type": "rect", + "symmetry": True, + "span": 2.0, + "root_chord": 0.4, + "span_cos_spacing": 0.0, + "chord_cos_spacing": 1.0, + } +mesh_wing = generate_mesh(mesh_dict) + +taper(mesh_wing,0.35,True,ref_axis_pos=0.0) + + +mesh_dict = { + "num_y": 3, + "num_x": 11, + "wing_type": "rect", + "symmetry": True, + "span": 1.0, + "root_chord": 0.325, + "span_cos_spacing": 0.0, + "chord_cos_spacing": 1.0, + } +mesh_tail = generate_mesh(mesh_dict) + +mesh_wing[:,:,1] = -mesh_wing[:,:,1] +mesh_wing = mesh_wing[:,::-1,:] + +mesh_tail[:,:,1] = -mesh_tail[:,:,1] +mesh_tail = mesh_tail[:,::-1,:] + +mesh_tail[:,:,0] += 1.5 + + +surf = { + "Wing": { + # General + # "component": np.int32(0), # logical surface component index (for grouping interacting surfaces, see AVL manual) + # "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh_wing), # (nx,ny,3) numpy array containing mesh coordinates + + }, + "Horizontal Tail": { + # General + # "component": np.int32(1), # logical surface component index (for grouping interacting surfaces, see AVL manual) + # "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh_tail), # (nx,ny,3) numpy array containing mesh coordinates + # Control Surface Specification + "control_assignments": { + "Elevator" : {"assignment":np.arange(0,mesh_tail.shape[1]), + "xhinged": 0.0, # x/c location of hinge + "vhinged": np.array([0,1,0]), # vector giving hinge axis about which surface rotates + "gaind": -1.0, # control surface gain + "refld": 1.0 # control surface reflection, sign of deflection for duplicated surface + } + }, + + } +} + + +input_dict = { + "title": "MACH MDAO AVL", + "mach": np.float64(0.0), + "iysym": np.int32(0), + "izsym": np.int32(0), + "zsym": np.float64(0.0), + "Sref": np.float64(1.13), + "Cref": np.float64(0.4), + "Bref": np.float64(1.0), + "XYZref": np.array([0.0, 0.0, 0.0],dtype=np.float64), + "CDp": np.float64(0.0), + "surfaces": surf, + # "bodies": fuselage, + # Global Control and DV info + "dname": ["Elevator"], # Name of control input for each corresonding index +} + +# For verification +# base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder +# geom_dir = os.path.join(base_dir, '..', 'geom_files') +# rect_file = os.path.join(geom_dir, 'aircraft_L1.avl') + + +# solver = OVLSolver(input_dict=input_dict,debug=True) +# solver = OVLSolver(geo_file=rect_file,debug=True) + +# solver.set_variable("alpha", 25.0) +# solver.set_variable("beta", 5.0) +# solver.execute_run() + +with open("aircraft_L1.pkl", 'wb') as f: + pickle.dump(input_dict, f) diff --git a/geom_files/make_pkls/make_rect_with_body.py b/geom_files/make_pkls/make_rect_with_body.py new file mode 100644 index 0000000..11bc622 --- /dev/null +++ b/geom_files/make_pkls/make_rect_with_body.py @@ -0,0 +1,103 @@ +# ============================================================================= +# Extension modules +# ============================================================================= +from optvl import OVLSolver + +# ============================================================================= +# Standard Python Modules +# ============================================================================= +import os + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np +import pickle + + +mesh = np.zeros((2,2,3)) +mesh[1,:,0] = 1.0 +mesh[:,1,1] = 1.0 +mesh[:,:,2] = 1.0 + + +surf = { + "Wing": { + # General + # "component": np.int32(1), # logical surface component index (for grouping interacting surfaces, see AVL manual) + # "yduplicate": np.float64(0.0), # surface is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0], dtype=np.float64 + ), # scaling factors applied to all x,y,z coordinates (chords arealso scaled by Xscale) + "translate": np.array( + [0.0, 0.0, 0.0], dtype=np.float64 + ), # offset added on to all X,Y,Z values in this surface + # Geometry: Mesh + "mesh": np.float64(mesh), # (nx,ny,3) numpy array containing mesh coordinates + # Control Surface Specification + "control_assignments": { + "Elevator" : {"assignment":np.arange(0,mesh.shape[1]), + "xhinged": 0.5, # x/c location of hinge + "vhinged": np.array([0,1,0]), # vector giving hinge axis about which surface rotates + "gaind": -1.0, # control surface gain + "refld": 1.0 # control surface reflection, sign of deflection for duplicated surface + } + }, + + } +} + +fuselage = {"Fuse pod": { + # General + # 'yduplicate': np.float64(0), # body is duplicated over the ysymm plane + # Geometry + "scale": np.array( + [1.0, 1.0, 1.0] + ), # scaling factors applied to all x,y,z coordinates (chords areal so scaled by Xscale) + "translate": np.array([0.0, 0.0, 0.0]), # offset added on to all X,Y,Z values in this surface + # Discretization + "nvb": np.int32(4), # number of source-line nodes + "bspace": np.float64(2.0), # lengthwise node spacing parameter + "bfile": "../geom_files/fuseSimple.dat", # body oml file name +} +} + +input_dict = { + "title": "MACH MDAO AVL", + "mach": np.float64(0.0), + "iysym": np.int32(0), + "izsym": np.int32(0), + "zsym": np.float64(0.0), + "Sref": np.float64(1.123), + "Cref": np.float64(0.25), + "Bref": np.float64(6.01), + "XYZref": np.array([0.0, 0, 0],dtype=np.float64), + "CDp": np.float64(0.0), + "surfaces": surf, + "bodies": fuselage, + # Global Control and DV info + "dname": ["Elevator"], # Name of control input for each corresonding index +} + +# For verification +# base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder +# geom_dir = os.path.join(base_dir, '..', 'geom_files') +# rect_file = os.path.join(geom_dir, 'rect_with_body.avl') + + +# solver = OVLSolver(input_dict=input_dict,debug=True) +# solver = OVLSolver(geo_file=rect_file,debug=True) + +# solver.set_variable("alpha", 25.0) +# solver.set_variable("beta", 5.0) +# solver.execute_run() + +# solver.plot_geom() + +with open("rect_with_body.pkl", 'wb') as f: + pickle.dump(input_dict, f) + + + + diff --git a/geom_files/rect_out.avl b/geom_files/rect_out.avl new file mode 100644 index 0000000..1cbd4d4 --- /dev/null +++ b/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 diff --git a/optvl/om_wrapper.py b/optvl/om_wrapper.py index 9a43395..93d7a51 100644 --- a/optvl/om_wrapper.py +++ b/optvl/om_wrapper.py @@ -24,7 +24,8 @@ class OVLGroup(om.Group): """ 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) @@ -40,6 +41,7 @@ def initialize(self): def setup(self): geom_file = self.options["geom_file"] + input_dict = self.options["input_dict"] mass_file = self.options["mass_file"] input_param_vals = self.options["input_param_vals"] @@ -50,7 +52,7 @@ def setup(self): 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", @@ -91,6 +93,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 +110,43 @@ 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"): + # 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) + if add_as == "inputs": + self.add_input(mesh_key, val=mesh, units="m", tags="geom_mesh") + elif add_as == "outputs": + self.add_output(mesh_key, val=mesh, units="m", tags="geom_mesh") - 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_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 +172,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 +193,16 @@ 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] = {} + # no scalars in meshes + surf_mesh_data[surf][param] = inputs[input_var] + return surf_data, surf_mesh_data def om_surf_dict_to_input(surf_dict): @@ -241,16 +274,23 @@ 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") 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) + self.ovl.set_mesh(idx_surf,surf_mesh_data[surf]["mesh"]) gam_arr = outputs["gamma"] gam_d_arr = outputs["gamma_d"] @@ -276,11 +316,17 @@ 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) 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) + self.ovl.set_mesh(idx_surf,surf_mesh_data[surf]["mesh"]) # def_dict = self.ovl.get_control_deflections() print("executing ovl run") @@ -313,7 +359,7 @@ 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) param_seeds = {} for param in self.ovl.param_idx_dict: @@ -325,8 +371,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,9 +386,9 @@ 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 + res_seeds=res_seeds, res_d_seeds=res_d_seeds, res_u_seeds=res_u_seeds,reshape_mesh_seeds = True ) ) @@ -356,10 +402,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: @@ -426,6 +475,7 @@ def setup(self): self.control_names = add_ovl_controls_as_inputs(self, self.ovl) add_ovl_geom_vars(self, self.ovl, add_as="inputs", include_airfoil_geom=input_airfoil_geom) + add_ovl_mesh_vars(self, self.ovl, add_as="inputs") # add the outputs for func_key in self.ovl.case_var_to_fort_var: @@ -455,11 +505,17 @@ 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) 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) + self.ovl.set_mesh(idx_surf,surf_mesh_data[surf]["mesh"]) gam_arr = inputs["gamma"] gam_d_arr = inputs["gamma_d"] @@ -538,11 +594,12 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): if ref in d_inputs: ref_seeds[ref] = d_inputs[ref] - geom_seeds = self.om_input_to_surf_dict(self, d_inputs) + geom_seeds, mesh_seeds = self.om_input_to_surf_dict(self, d_inputs) func_seeds, _, csd_seeds, stab_derivs_seeds, body_axis_seeds, _, _ = self.ovl._execute_jac_vec_prod_fwd( con_seeds=con_seeds, geom_seeds=geom_seeds, + mesh_seeds=mesh_seeds, gamma_seeds=gamma_seeds, gamma_d_seeds=gamma_d_seeds, gamma_u_seeds=gamma_u_seeds, @@ -610,12 +667,13 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): # print(var_name, body_axis_seeds[func_key]) print(f" running rev mode derivs for {func_key}") - con_seeds, geom_seeds, mesh_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = ( + con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds = ( self.ovl._execute_jac_vec_prod_rev( func_seeds=func_seeds, consurf_derivs_seeds=csd_seeds, stab_derivs_seeds=stab_derivs_seeds, body_axis_derivs_seeds=body_axis_seeds, + reshape_mesh_seeds = True, ) ) @@ -629,10 +687,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: @@ -731,9 +792,10 @@ def setup(self): 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) + add_ovl_mesh_vars(self, avl, add_as="outputs") - if mesh_output: - add_ovl_mesh_out_as_output(self,avl) + # if mesh_output: + # add_ovl_mesh_out_as_output(self,avl) class Differencer(om.ExplicitComponent): def setup(self): diff --git a/optvl/optvl_class.py b/optvl/optvl_class.py index 8ab2a53..bbd9d12 100644 --- a/optvl/optvl_class.py +++ b/optvl/optvl_class.py @@ -403,10 +403,11 @@ def __init__( deriv_key = self._get_deriv_key(var, func) self.case_body_derivs_to_fort_var[deriv_key] = ["CASE_R", f"{func_to_prefix[func]}TOT_U_BA", idx_var] - # In the case where we used a file then we have to initialize these before _init_map_data so ad seeds work correctly + # In the case where we used a file then we have to initialize these before _init_map_data so ad seeds work correclty if not input_dict: self.mesh_idx_first = np.zeros(self.get_num_surfaces(),dtype=np.int32) self.y_offsets = np.zeros(self.get_num_surfaces(),dtype=np.float64) + self.point_sets = {} # the case parameters are stored in a 1d array, # these indices correspond to the position of each parameter in that arra @@ -415,6 +416,9 @@ def __init__( # set the default solver tolerance self.set_avl_fort_arr("CASE_R", "EXEC_TOL", 2e-5) + # init the DVGeo variable + self.DVGeo = None + if timing: print(f"AVL init took {time.time() - start_time} seconds") @@ -788,6 +792,9 @@ def check_type(key, avl_vars, given_val, cast_type=True): # a duplicated mesh as its normally only passed as a dummy argument into the SDUPL subroutine and not stored in the fortran layer. self.y_offsets = np.zeros(self.get_num_surfaces(),dtype=np.float64) + # Class dictionary to store the point set names if a DVGeo object is used + self.point_sets = {} + # Load surfaces if num_surfs > 0: surf_names = list(input_dict["surfaces"].keys()) @@ -1058,8 +1065,10 @@ def check_type(key, avl_vars, given_val, cast_type=True): surf_dict["flatten mesh"] = True self.set_mesh(idx_surf, surf_dict["mesh"],flatten=surf_dict["flatten mesh"],update_nvs=True,update_nvc=True) # set_mesh handles the Fortran indexing and ordering self.avl.makesurf_mesh(idx_surf + 1) #+1 for Fortran indexing + self.point_sets[idx_surf] = "optvl_%s_coords" % surf_name # store the pointset name for DVGeo later else: self.avl.makesurf(idx_surf + 1) # +1 to convert to 1 based indexing + self.point_sets[idx_surf] = None # No pointset available if no mesh if "yduplicate" in surf_dict.keys(): self.avl.sdupl(idx_surf + 1, surf_dict["yduplicate"], "YDUP") @@ -3149,6 +3158,86 @@ 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. + + Examples + -------- + >>> CFDsolver = (comm=comm, options=CFDoptions) + >>> CFDsolver.setDVGeo(DVGeo) + """ + + self.DVGeo = DVGeo + + # save the common kwargs dict. default is empty + if pointSetKwargs is None: + self.pointSetKwargs = {} + else: + self.pointSetKwargs = pointSetKwargs + + # save if we have customPointSetFamilies. this default is not mutable so we can just set it as is. + self.customPointSetFamilies = customPointSetFamilies + + def update_DVGeo(self): + """If DVGeo is present this function will embed all the meshes into it if needed and then perform an update + of the pointset and set the meshes back into AVL. + """ + # Check if we have an DVGeo object to deal with: + if self.DVGeo is not None: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a mesh, skip it + + # Embed the points if they haven't been already + if point_set_name not in self.DVGeo.points: + mesh = self.get_mesh(idx_surf) + coords0 = mesh.transpose((1,0,2)).reshape((mesh.shape[0]*mesh.shape[1],3)) + self.DVGeo.addPointSet(coords0, point_set_name, **self.pointSetKwargs) + # print(f"Embedeed point set {point_set_name}!") + + # Check if our point-set is up to date and update the mesh accordingly + if not self.DVGeo.pointSetUpToDate(point_set_name): + coords = self.DVGeo.update(point_set_name) + mesh_old = self.get_mesh(idx_surf) + mesh_new = copy.deepcopy(coords.reshape((mesh_old.shape[1],mesh_old.shape[0],3)).transpose((1,0,2))) + self.set_mesh(idx_surf,mesh_new) + + # Now update all of our surfaces in AVL + self.avl.update_surfaces() # region --- Utility functions def get_num_surfaces(self) -> int: @@ -3828,6 +3917,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 +3932,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 +3966,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 +3993,38 @@ def _execute_jac_vec_prod_fwd( # self.clear_ad_seeds() self.set_variable_ad_seeds(con_seeds) self.set_geom_ad_seeds(geom_seeds) - self.set_mesh_ad_seeds(mesh_seeds) + # self.set_mesh_ad_seeds(mesh_seeds) self.set_gamma_ad_seeds(gamma_seeds) self.set_gamma_d_ad_seeds(gamma_d_seeds) self.set_gamma_u_ad_seeds(gamma_u_seeds) self.set_parameter_ad_seeds(param_seeds) self.set_reference_ad_seeds(ref_seeds) + # Since DVGeo seeds operate entirely within the python layer we set them here + if self.DVGeo is not None and dvgeo_seeds is not None: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a pointset, skip it + + # If no mesh seed was provided for the surface we will need to start it at zero + if surface not in mesh_seeds.keys(): + nx = self.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + mesh_seeds[surface] = {} + mesh_seeds[surface]["mesh"] = np.zeros((nx*ny,3)) + + # Loop over the design variables and accumulate the sensitivity product into the mesh_seeds + mesh_seeds[surface]["mesh"] += self.DVGeo.totalSensitivityProd(dvgeo_seeds[surface], point_set_name).reshape( + mesh_seeds[surface]["mesh"].shape + ) + + self.set_mesh_ad_seeds(mesh_seeds) + self.avl.update_surfaces_d() self.avl.get_res_d() self.avl.velsum_d() @@ -3942,6 +4061,37 @@ def _execute_jac_vec_prod_fwd( self.set_parameter_ad_seeds(param_seeds, mode="FD", scale=step) self.set_reference_ad_seeds(ref_seeds, mode="FD", scale=step) + + # Since DVGeo operates entirely within the python layer we have have to do this + if self.DVGeo is not None and dvgeo_seeds is not None: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a pointset, skip it + + # Apply the FD step to the DV seeds + for dv in dvgeo_seeds[surface].keys(): + # current design var values + currentDV = self.DVGeo.getValues()[dv] + + # Set the updated DVs + self.DVGeo.setDesignVars({dv: currentDV + dvgeo_seeds[surface][dv]*step}) + + # get mesh size + nx = self.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + + # Compute the perturbed mesh values + coords = self.DVGeo.update(point_set_name) + mesh_pertub = copy.deepcopy(coords.reshape((ny,nx,3)).transpose((1,0,2))) + + self.set_mesh(idx_surf,mesh_pertub) + + # propogate the seeds through without resolving self.avl.update_surfaces() self.avl.get_res() @@ -3966,6 +4116,36 @@ def _execute_jac_vec_prod_fwd( self.set_parameter_ad_seeds(param_seeds, mode="FD", scale=-1 * step) self.set_reference_ad_seeds(ref_seeds, mode="FD", scale=-1 * step) + # Set the mesh seeds back + if self.DVGeo is not None and dvgeo_seeds is not None: + # Loop over all surfaces + for surface in self.unique_surface_names: + # Get the pointset name + idx_surf = self.get_surface_index(surf_name=surface) + if idx_surf in self.point_sets.keys(): + point_set_name = self.point_sets[idx_surf] + else: + continue # This surface doesn't have a pointset, skip it + + # Restore the orignal dvgeo seeds + for dv in dvgeo_seeds[surface].keys(): + # current design var values + currentDV = self.DVGeo.getValues()[dv] + + # Set the updated DVs + self.DVGeo.setDesignVars({dv: currentDV - dvgeo_seeds[surface][dv]*step}) + + # get mesh size + nx = self.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + ny = self.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + + # Compute the original mesh values + coords = self.DVGeo.update(point_set_name) + mesh_orig = copy.deepcopy(coords.reshape((ny,nx,3)).transpose((1,0,2))) + + self.set_mesh(idx_surf,mesh_orig) + + self.avl.update_surfaces() self.avl.get_res() self.avl.velsum() @@ -4024,6 +4204,7 @@ def _execute_jac_vec_prod_rev( body_axis_derivs_seeds: Optional[Dict[str, float]] = None, res_d_seeds: Optional[np.ndarray] = None, res_u_seeds: Optional[np.ndarray] = None, + reshape_mesh_seeds: Optional[bool] = False, print_timings: Optional[bool] = False, ) -> Tuple[ Dict[str, float], @@ -4044,6 +4225,7 @@ def _execute_jac_vec_prod_rev( body_axis_derivs_seeds: body axis derivatives AD seeds res_d_seeds: dResidual/d(Controls Deflection) AD seeds res_u_seeds: dResidual/d(flight condition) AD seeds + reshape_mesh_seeds: Reshape the mesh_seeds output to match the shape of the input mesh array print_timings: flag to show timing data Returns: @@ -4130,6 +4312,43 @@ 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) + ) + + if reshape_mesh_seeds: + for surface in self.unique_surface_names: + idx_surf = self.get_surface_index(surf_name=surface) + if self.avl.SURF_MESH_L.LSURFMSH[idx_surf]: + mesh_seed = mesh_seeds[surface]["mesh"] + # nx = self.avl.SURF_GEOM_I.NVC[idx_surf] + 1 + # ny = self.avl.SURF_GEOM_I.NVS[idx_surf] + 1 + mesh_cur = self.get_mesh(idx_surf) + # copy.deepcopy(coords.reshape((mesh_old.shape[1],mesh_old.shape[0],3)).transpose((1,0,2))) + mesh_seed = np.reshape(mesh_seed,(mesh_cur.shape[1],mesh_cur.shape[0],3)) + mesh_seed = np.transpose(mesh_seed,(1,0,2)) + mesh_seeds[surface]["mesh"] = mesh_seed + # mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seed.reshape((mesh_cur.shape[1],mesh_cur.shape[0],3)).transpose((1,0,2))) + # mesh_seeds[surface]["mesh"] = copy.deepcopy(mesh_seed.reshape((nx,ny,3))) + self.set_function_ad_seeds(func_seeds, scale=0.0) self.set_residual_ad_seeds(res_seeds, scale=0.0) self.set_residual_d_ad_seeds(res_d_seeds, scale=0.0) @@ -4144,7 +4363,7 @@ def _execute_jac_vec_prod_rev( if print_timings: print(f" Total Time: {time.time() - time_start}") - return con_seeds, geom_seeds, mesh_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds + return con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds def execute_run_sensitivities( self, @@ -4178,7 +4397,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 +4414,7 @@ def execute_run_sensitivities( # get the resulting adjoint vector (dfunc/dRes) from fortran dfdR = self.get_residual_ad_seeds() # self.clear_ad_seeds() - con_seeds, geom_seeds, mesh_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( + con_seeds, geom_seeds, mesh_seeds, dvgeo_seeds, _, _, _, param_seeds, ref_seeds = self._execute_jac_vec_prod_rev( func_seeds={func: 1.0}, res_seeds=dfdR ) if print_timings: @@ -4205,7 +4424,7 @@ def execute_run_sensitivities( sens[func].update(con_seeds) # I don't know if it's worth combining geom_seeds and mesh_seeds into one just to make this one part less nasty for key in geom_seeds: - sens[func][key] = geom_seeds[key] | mesh_seeds[key] + sens[func][key] = geom_seeds[key] | mesh_seeds[key] | dvgeo_seeds[key] # sens[func].update(geom_seeds) # sens[func].update(mesh_seeds) sens[func].update(param_seeds) @@ -4222,7 +4441,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 +4461,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 +4470,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 +4487,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 +4507,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 +4516,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 +4535,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 +4555,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 +4564,10 @@ def execute_run_sensitivities( time_last = time.time() sens[func_key].update(con_seeds) - for surf_key in geom_seeds: - sens[func_key][surf_key] = geom_seeds[surf_key] | mesh_seeds[surf_key] + for key in geom_seeds: + sens[func_key][key] = geom_seeds[key] | mesh_seeds[key] | dvgeo_seeds[key] + # sens[func_key].update(geom_seeds) + # sens[func_key].update(mesh_seeds) sens[func_key].update(param_seeds) sens[func_key].update(ref_seeds) diff --git a/optvl/utils/ffd_utils.py b/optvl/utils/ffd_utils.py new file mode 100644 index 0000000..e3cdc7e --- /dev/null +++ b/optvl/utils/ffd_utils.py @@ -0,0 +1,76 @@ +# ============================================================================= +# Standard Python Modules +# ============================================================================= +import os + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np + +# ============================================================================= +# Local modules +# ============================================================================= + + +def write_FFD_file(mesh, mx, my, filename="ffd", cushion=0.05): + """Utility functions that generates a box FFD around a wing mesh. + Based on the utility function in OpenAeroStruct. + + Args: + mesh: mesh array for the surface to be embedded in the FFD np.ndarray (nx,ny,3) + mx: number of control points in the x-direction + my: number of control points in the y-direction + filename: name of the ffd file to write without extension.Defaults to "ffd.dat". + cushion: Amount of cushion to apply from the LE/TE and tips. Defaults to 0.05. + + Returns: + filename of the written FFD file as a str + """ + nx, ny = mesh.shape[:2] + + half_ffd = np.zeros((mx, my, 3)) + + LE = mesh[0, :, :] + TE = mesh[-1, :, :] + + half_ffd[0, :, 0] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), LE[:, 0]) + half_ffd[0, :, 1] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), LE[:, 1]) + half_ffd[0, :, 2] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), LE[:, 2]) + + half_ffd[-1, :, 0] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), TE[:, 0]) + half_ffd[-1, :, 1] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), TE[:, 1]) + half_ffd[-1, :, 2] = np.interp(np.linspace(0, 1, my), np.linspace(0, 1, ny), TE[:, 2]) + + for i in range(my): + half_ffd[:, i, 0] = np.linspace(half_ffd[0, i, 0], half_ffd[-1, i, 0], mx) + half_ffd[:, i, 1] = np.linspace(half_ffd[0, i, 1], half_ffd[-1, i, 1], mx) + half_ffd[:, i, 2] = np.linspace(half_ffd[0, i, 2], half_ffd[-1, i, 2], mx) + + half_ffd[0, :, 0] -= cushion + half_ffd[-1, :, 0] += cushion + half_ffd[:, 0, 1] -= cushion + half_ffd[:, -1, 1] += cushion + + bottom_ffd = half_ffd.copy() + bottom_ffd[:, :, 2] -= cushion + + top_ffd = half_ffd.copy() + top_ffd[:, :, 2] += cushion + + ffd = np.vstack((bottom_ffd, top_ffd)) + + filename = filename + ".xyz" + + with open(filename, "w") as f: + f.write("1\n") + f.write("{} {} {}\n".format(mx, 2, my)) + x = np.array_str(ffd[:, :, 0].flatten(order="F"))[1:-1] + "\n" + y = np.array_str(ffd[:, :, 1].flatten(order="F"))[1:-1] + "\n" + z = np.array_str(ffd[:, :, 2].flatten(order="F"))[1:-1] + "\n" + + f.write(x) + f.write(y) + f.write(z) + + return filename \ No newline at end of file diff --git a/tests/test_body_axis_derivs_partial_derivs.py b/tests/test_body_axis_derivs_partial_derivs.py index 7359402..a471e97 100644 --- a/tests/test_body_axis_derivs_partial_derivs.py +++ b/tests/test_body_axis_derivs_partial_derivs.py @@ -196,7 +196,7 @@ def test_rev_gamma_u(self): # for var_key in bd_d_fwd[deriv_func]: bd_d_rev = {deriv_func: 1.0} - gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(body_axis_derivs_seeds=bd_d_rev)[5] + gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(body_axis_derivs_seeds=bd_d_rev)[6] rev_sum = np.sum(gamma_u_seeds_rev * gamma_u_seeds_fwd) @@ -247,7 +247,7 @@ def test_rev_ref(self): for deriv_func, var_dict in self.ovl_solver.case_body_derivs_to_fort_var.items(): body_axis_deriv_seeds_rev[deriv_func] = np.random.rand(1)[0] - ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(body_axis_derivs_seeds=body_axis_deriv_seeds_rev)[7] + ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(body_axis_derivs_seeds=body_axis_deriv_seeds_rev)[8] self.ovl_solver.clear_ad_seeds_fast() diff --git a/tests/test_consurf_partial_derivs.py b/tests/test_consurf_partial_derivs.py index f1a0939..bb82ec2 100644 --- a/tests/test_consurf_partial_derivs.py +++ b/tests/test_consurf_partial_derivs.py @@ -168,7 +168,7 @@ def test_rev_gamma_d(self): res_d_seeds_fwd = self.ovl_solver._execute_jac_vec_prod_fwd(gamma_d_seeds=gamma_d_seeds_fwd)[5] self.ovl_solver.clear_ad_seeds_fast() - gamma_d_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_d_seeds=res_d_seeds_rev)[4] + gamma_d_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_d_seeds=res_d_seeds_rev)[5] gamma_sum = np.sum(gamma_d_seeds_rev * gamma_d_seeds_fwd) res_sum = np.sum(res_d_seeds_rev * res_d_seeds_fwd) @@ -359,8 +359,7 @@ def test_rev_gamma_d(self): for deriv_func in cs_d_fwd: cs_d_rev = {deriv_func: 1.0} - gamma_d_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(consurf_derivs_seeds=cs_d_rev)[4] - + gamma_d_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(consurf_derivs_seeds=cs_d_rev)[5] rev_sum = np.sum(gamma_d_seeds_rev * gamma_d_seeds_fwd) fwd_sum = np.sum(cs_d_fwd[deriv_func]) diff --git a/tests/test_partial_derivs.py b/tests/test_partial_derivs.py index d1527ca..4869b84 100644 --- a/tests/test_partial_derivs.py +++ b/tests/test_partial_derivs.py @@ -16,6 +16,7 @@ import numpy as np + base_dir = os.path.dirname(os.path.abspath(__file__)) # Path to current folder geom_dir = os.path.join(base_dir, '..', 'geom_files') @@ -215,7 +216,7 @@ def test_rev_gamma(self): self.ovl_solver.clear_ad_seeds_fast() for func_key in self.ovl_solver.case_var_to_fort_var: - gamma_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[3] + gamma_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[4] rev_sum = np.sum(gamma_seeds_rev * gamma_seeds_fwd) fwd_sum = np.sum(func_seeds_fwd[func_key]) @@ -262,7 +263,7 @@ def test_rev_param(self): self.ovl_solver.clear_ad_seeds_fast() for func_key in self.ovl_solver.case_var_to_fort_var: - param_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[6] + param_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[7] # print(f"{func_key} wrt {param_key}", "fwd ", func_seeds_fwd[func_key], "rev", param_seeds_rev[param_key]) tol = 1e-14 @@ -316,7 +317,7 @@ def test_rev_ref(self): self.ovl_solver.clear_ad_seeds_fast() for func_key in self.ovl_solver.case_var_to_fort_var: - ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[7] + ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(func_seeds={func_key: 1.0})[8] # print(f"{func_key} wrt {ref_key}", "fwd ", func_seeds_fwd[func_key], "rev", ref_seeds_rev[ref_key]) tol = 1e-14 @@ -465,7 +466,7 @@ def test_fwd_param(self): def test_rev_param(self): num_res = self.ovl_solver.get_mesh_size() res_seeds_rev = np.random.rand(num_res) - param_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_seeds=res_seeds_rev)[6] + param_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_seeds=res_seeds_rev)[7] self.ovl_solver.clear_ad_seeds_fast() @@ -498,7 +499,7 @@ def test_fwd_ref(self): def test_rev_ref(self): num_res = self.ovl_solver.get_mesh_size() res_seeds_rev = np.random.rand(num_res) - ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_seeds=res_seeds_rev)[7] + ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_seeds=res_seeds_rev)[8] self.ovl_solver.clear_ad_seeds_fast() diff --git a/tests/test_pygeo.py b/tests/test_pygeo.py new file mode 100644 index 0000000..eeb2e4c --- /dev/null +++ b/tests/test_pygeo.py @@ -0,0 +1,463 @@ +# ============================================================================= +# 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 + +# --------------------------------------------------------------------------- +# 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) + +from wing_mesh import mesh # (nx, ny, 3) numpy array + +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 + + # NOTE: Current having trouble with global DVs in tests + nrefaxpts = self.ovl_solver.DVGeo.addRefAxis("c4", xFraction=0.25, alignIndex="k") + + def twist(val, geo): + for i in range(nrefaxpts): + geo.rot_y["c4"].coef[i] = val[i] + + def sweep(val, geo): + # the extractCoef method gets the unperturbed ref axis control points + C = geo.extractCoef("c4") + C_orig = C.copy() + # we will sweep the wing about the first point in the ref axis + sweep_ref_pt = C_orig[0, :] + + theta = -val[0] * np.pi / 180 + # rot_mtx = np.array([[np.cos(theta), 0.0, -np.sin(theta)], [0.0, 1.0, 0.0], [np.sin(theta), 0.0, np.cos(theta)]]) + rot_mtx = np.array([[np.cos(theta), -np.sin(theta), 0.0], [np.sin(theta), np.cos(theta), 0.0], [0.0, 0.0, 1.0]]) + + # modify the control points of the ref axis + # by applying a rotation about the first point in the x-z plane + for i in range(nrefaxpts): + # get the vector from each ref axis point to the wing root + vec = C[i, :] - sweep_ref_pt + # need to now rotate this by the sweep angle and add back the wing root loc + C[i, :] = sweep_ref_pt + rot_mtx @ vec + # use the restoreCoef method to put the control points back in the right place + geo.restoreCoef(C, "c4") + + self.ovl_solver.DVGeo.addGlobalDV("twist", func=twist, value=np.zeros(nrefaxpts), lower=-10, upper=10, scale=0.05) + self.ovl_solver.DVGeo.addGlobalDV("sweep", func=sweep, value=0.0, lower=0, upper=45, scale=0.05) + self.ovl_solver.DVGeo.addLocalDV("shape", lower=-0.25, upper=0.25, axis="z", scale=1.0) + + # Execute the solver + self.ovl_solver.update_DVGeo() + self.ovl_solver.set_variable("alpha", 5.0) + self.ovl_solver.set_variable("beta", 0.0) + self.ovl_solver.set_parameter("Mach", 0.0) + self.ovl_solver.execute_run() + + # compare the analytical gradients with finite difference for dvgeo + surf_key = list(self.ovl_solver.surf_mesh_to_fort_var.keys())[0] + dvgeo_vars = self.ovl_solver.DVGeo.getVarNames() + cs_names = self.ovl_solver.get_control_names() + + consurf_vars = [] + for func_key in self.ovl_solver.case_derivs_to_fort_var: + consurf_vars.append(self.ovl_solver._get_deriv_key(cs_names[0], func_key)) + + func_vars = self.ovl_solver.case_var_to_fort_var + stab_derivs = self.ovl_solver.case_stab_derivs_to_fort_var + body_axis_derivs = self.ovl_solver.case_body_derivs_to_fort_var + + sens = self.ovl_solver.execute_run_sensitivities( + func_vars, + consurf_derivs=consurf_vars, + stab_derivs=stab_derivs, + body_axis_derivs=body_axis_derivs, + print_timings=False, + ) + + # for con_key in self.ovl_solver.con_var_to_fort_var: + sens_FD = {} + for surf_key in self.ovl_solver.surf_geom_to_fort_var: + sens_FD[surf_key] = {} + for dvgeo_var_key in dvgeo_vars: + # arr = self.ovl_solver.get_mesh(self.ovl_solver.get_surface_index(surf_key)).reshape(-1,3) + arr = self.ovl_solver.DVGeo.getValues()[dvgeo_var_key] + np.random.seed(arr.size) + rand_arr = np.random.rand(*arr.shape) + rand_arr /= np.linalg.norm(rand_arr) + + func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif({surf_key: {dvgeo_var_key: rand_arr}}, step=1.0e-7) + + for func_key in func_vars: + dvgeo_var_dot = np.sum(sens[func_key][surf_key][dvgeo_var_key] * rand_arr) + func_dot = func_seeds[func_key] + + rel_err = np.abs(dvgeo_var_dot - func_dot) / np.abs(func_dot + 1e-20) + + # print( + # f"{func_key:5} wrt {surf_key}:{dvgeo_var_key:10} | AD:{dvgeo_var_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" + # ) + tol = 1e-7 + if np.abs(dvgeo_var_dot) < tol or np.abs(func_dot) < tol: + # If either value is basically zero, use an absolute tolerance + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + atol=1e-4, + err_msg=f"{func_key:5} wrt {surf_key}:{dvgeo_var_key:10}", + ) + else: + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + rtol=5e-3, + err_msg=f"{func_key:5} wrt {surf_key}:{dvgeo_var_key:10}", + ) + + for func_key in consurf_vars: + # for cs_key in consurf_vars[func_key]: + dvgeo_var_dot = np.sum(sens[func_key][surf_key][dvgeo_var_key] * rand_arr) + func_dot = consurf_deriv_seeds[func_key] + + # rel_err = np.abs(dvgeo_var_dot - func_dot) / np.abs(func_dot + 1e-20) + # print( + # f"{func_key} wrt {surf_key}:{mesh_key:10} | AD:{dvgeo_var_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" + # ) + + tol = 1e-8 + if np.abs(dvgeo_var_dot) < tol or np.abs(func_dot) < tol: + # If either value is basically zero, use an absolute tolerance + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + atol=1e-4, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + else: + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + rtol=6e-3, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + + for func_key in stab_derivs_seeds: + if func_key == "spiral parameter" or func_key == "lateral parameter": + continue # These derivatives blow up in different way for both adjoint and FD + dvgeo_var_dot = np.sum(sens[func_key][surf_key][dvgeo_var_key] * rand_arr) + func_dot = stab_derivs_seeds[func_key] + + rel_err = np.abs(dvgeo_var_dot - func_dot) / np.abs(func_dot + 1e-20) + + # print( + # f"{func_key} wrt {surf_key}:{dvgeo_var_key:10} | AD:{dvgeo_var_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" + # ) + + tol = 5e-7 + if np.abs(dvgeo_var_dot) < tol or np.abs(func_dot) < tol: + # If either value is basically zero, use an absolute tolerance + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + atol=5e-9, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + else: + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + rtol=6e-3, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + + for func_key in body_axis_derivs_seeds: + dvgeo_var_dot = np.sum(sens[func_key][surf_key][dvgeo_var_key] * rand_arr) + func_dot = body_axis_derivs_seeds[func_key] + + rel_err = np.abs(dvgeo_var_dot - func_dot) / np.abs(func_dot + 1e-20) + + # print( + # f"{func_key} wrt {surf_key}:{dvgeo_var_key:10} | AD:{dvgeo_var_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" + # ) + + tol = 1e-6 + if np.abs(dvgeo_var_dot) < tol or np.abs(func_dot) < tol: + # If either value is basically zero, use an absolute tolerance + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + atol=5e-8, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + else: + np.testing.assert_allclose( + dvgeo_var_dot, + func_dot, + rtol=6e-3, + err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", + ) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/tests/test_stab_derivs_partial_derivs.py b/tests/test_stab_derivs_partial_derivs.py index e404a84..4491f1c 100644 --- a/tests/test_stab_derivs_partial_derivs.py +++ b/tests/test_stab_derivs_partial_derivs.py @@ -165,7 +165,7 @@ def test_rev_gamma_u(self): res_u_seeds_fwd = self.ovl_solver._execute_jac_vec_prod_fwd(gamma_u_seeds=gamma_u_seeds_fwd)[6] self.ovl_solver.clear_ad_seeds_fast() - gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_u_seeds=res_u_seeds_rev)[5] + gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(res_u_seeds=res_u_seeds_rev)[6] gamma_sum = np.sum(gamma_u_seeds_rev * gamma_u_seeds_fwd) res_sum = np.sum(res_u_seeds_rev * res_u_seeds_fwd) @@ -378,7 +378,7 @@ def test_rev_gamma_u(self): # for var_key in sd_d_fwd[deriv_func]: sd_d_rev = {deriv_func: 1.0} - gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(stab_derivs_seeds=sd_d_rev)[5] + gamma_u_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(stab_derivs_seeds=sd_d_rev)[6] rev_sum = np.sum(gamma_u_seeds_rev * gamma_u_seeds_fwd) @@ -425,7 +425,7 @@ def test_rev_ref(self): for deriv_func, var_dict in self.ovl_solver.case_stab_derivs_to_fort_var.items(): stab_deriv_seeds_rev[deriv_func] = np.random.rand(1)[0] - ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(stab_derivs_seeds=stab_deriv_seeds_rev)[7] + ref_seeds_rev = self.ovl_solver._execute_jac_vec_prod_rev(stab_derivs_seeds=stab_deriv_seeds_rev)[8] self.ovl_solver.clear_ad_seeds_fast() diff --git a/tests/test_total_derivs.py b/tests/test_total_derivs.py index 5d562ad..e7b87c3 100644 --- a/tests/test_total_derivs.py +++ b/tests/test_total_derivs.py @@ -211,6 +211,7 @@ def test_geom(self): surf_key = list(self.ovl_solver.surf_geom_to_fort_var.keys())[0] geom_vars = self.ovl_solver.surf_geom_to_fort_var[surf_key] + # geom_vars += self.ovl_solver.surf_mesh_to_fort_var[surf_key] cs_names = self.ovl_solver.get_control_names() consurf_vars = [] diff --git a/tests/wing_mesh.npy b/tests/wing_mesh.npy deleted file mode 100644 index df1bec4..0000000 Binary files a/tests/wing_mesh.npy and /dev/null differ