diff --git a/geom_files/rect_out.avl b/geom_files/rect_out.avl deleted file mode 100644 index 420353ef..00000000 --- a/geom_files/rect_out.avl +++ /dev/null @@ -1,41 +0,0 @@ -# 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/optvl/optvl_class.py b/optvl/optvl_class.py index c4260525..ee72b1b2 100644 --- a/optvl/optvl_class.py +++ b/optvl/optvl_class.py @@ -216,6 +216,11 @@ class OVLSolver(object): "XYZref": ["CASE_R","XYZREF0"], "CDp": ["CASE_R", "CDREF0"], } + + viz_aifoil_data_list = [ + "casec", + "tasec" + ] # fmt: on @@ -325,11 +330,13 @@ def __init__( } # control surfaces added in __init__ - # TODO: the keys of this dict aren't used - self.con_var_to_fort_var = { - "alpha": ["CASE_R", "ALFA"], - "beta": ["CASE_R", "BETA"], - } + self.con_var_list = [ + "alpha", + "beta", + "roll rate", + "pitch rate", + "yaw rate", + ] control_names = self.get_control_names() self.dindex_to_con_surf = OrderedDict() @@ -342,7 +349,7 @@ def __init__( idx_control_start = np.max([x for x in self.conval_idx_dict.values()]) + 1 for idx_c_var, c_name in enumerate(control_names): self.conval_idx_dict[c_name] = idx_control_start + idx_c_var - self.con_var_to_fort_var[c_name] = ["CASE_R", "DELCON"] + self.con_var_list.append(c_name) var_to_suffix = { "alpha": "AL", @@ -949,6 +956,7 @@ def check_type(key, avl_vars, given_val, cast_type=True): for key in self.surf_section_geom_to_fort_var[surf_name]: avl_vars_secs = self.surf_section_geom_to_fort_var[surf_name][key] + avl_vars = (avl_vars_secs[0], avl_vars_secs[1], avl_vars_secs[2][j]) if key not in surf_dict: @@ -3416,11 +3424,39 @@ def _get_deriv_key(self, var: str, func: str) -> str: # --------------------------- # --- Derivative routines --- # --------------------------- + + # --- utils --- + def zero_seed_dict(self, seed_dict: Dict) -> Dict: + for key, value in seed_dict.items(): + seed_dict[key] = self._zero_value(value) + return seed_dict + + def _zero_value(self, value): + if isinstance(value, dict): + for k, v in value.items(): + value[k] = self._zero_value(v) + return value + elif isinstance(value, list): + return [self._zero_value(v) for v in value] + elif isinstance(value, np.ndarray): + value[...] = 0 + return value + else: + return 0 + + def deep_update(self, d, other): + for k, v in other.items(): + if isinstance(v, dict) and isinstance(d.get(k), dict): + self.deep_update(d[k], v) + else: + d[k] = v + return d + # --- input ad seeds --- def get_variable_ad_seeds(self) -> Dict[str, float]: var_seeds = {} - for con in self.con_var_to_fort_var: + for con in self.con_var_list: idx_con = self.conval_idx_dict[con] blk = "CASE_R" + self.ad_suffix var = "CONVAL" + self.ad_suffix @@ -3450,7 +3486,7 @@ def set_variable_ad_seeds(self, con_seeds: Dict[str, Dict[str, float]], mode: st self.set_avl_fort_arr(blk, var, val, slicer=slicer) elif mode == "FD": - # # reverse lookup in the con_var_to_fort_var dict + # # reverse lookup in the con_var_list dict if con in self.con_surf_to_dindex: val = self.get_control_deflection(con) else: @@ -3544,23 +3580,56 @@ def get_geom_ad_seeds(self) -> Dict[str, Dict[str, float]]: var += self.ad_suffix geom_seeds[surf_key][geom_key] = copy.deepcopy(self.get_avl_fort_arr(blk, var, slicer=slicer)) + + for geom_key in self.surf_section_geom_to_fort_var[surf_key]: + if geom_key in self.viz_aifoil_data_list: + continue + + blk, var, sec_slicers = self.surf_section_geom_to_fort_var[surf_key][geom_key] + blk += self.ad_suffix + var += self.ad_suffix + + geom_seeds[surf_key][geom_key] = [] + for sec_slice in (sec_slicers): + geom_seeds[surf_key][geom_key].append(copy.deepcopy(self.get_avl_fort_arr(blk, var, slicer=sec_slice))) return geom_seeds def set_geom_ad_seeds(self, geom_seeds: Dict[str, float], mode: str = "AD", scale=1.0) -> None: for surf_key in geom_seeds: for geom_key in geom_seeds[surf_key]: - blk, var, slicer = self.surf_geom_to_fort_var[surf_key][geom_key] + if geom_key in self.surf_geom_to_fort_var[surf_key]: + blk, var, slicer = self.surf_geom_to_fort_var[surf_key][geom_key] + + if mode == "AD": + blk += self.ad_suffix + var += self.ad_suffix + val = geom_seeds[surf_key][geom_key] * scale + elif mode == "FD": + val = self.get_avl_fort_arr(blk, var, slicer=slicer) + val += geom_seeds[surf_key][geom_key] * scale + + self.set_avl_fort_arr(blk, var, val, slicer=slicer) + + elif geom_key in self.viz_aifoil_data_list: + raise ValueError(f"Can not set key {geom_key}. it is only used for viz") + elif geom_key in self.surf_section_geom_to_fort_var[surf_key]: + + blk, var, sec_slicers = self.surf_section_geom_to_fort_var[surf_key][geom_key] - if mode == "AD": - blk += self.ad_suffix - var += self.ad_suffix - val = geom_seeds[surf_key][geom_key] * scale - elif mode == "FD": - val = self.get_avl_fort_arr(blk, var, slicer=slicer) - val += geom_seeds[surf_key][geom_key] * scale - # print(blk, var, val, slicer) - self.set_avl_fort_arr(blk, var, val, slicer=slicer) + if mode == "AD": + blk += self.ad_suffix + var += self.ad_suffix + + for idx_slice, sec_slice in enumerate(sec_slicers): + if mode == "AD": + val = geom_seeds[surf_key][geom_key][idx_slice] * scale + elif mode == "FD": + val = self.get_avl_fort_arr(blk, var, slicer=sec_slice) + val += geom_seeds[surf_key][geom_key][idx_slice] * scale + + self.set_avl_fort_arr(blk, var, val, slicer=sec_slice) + def get_mesh_ad_seeds(self) -> Dict[str, Dict[str, float]]: mesh_seeds = {} @@ -3940,15 +4009,41 @@ def _execute_jac_vec_prod_fwd( mesh_size = self.get_mesh_size() num_control_surfs = self.get_num_control_surfs() - - if con_seeds is None: - con_seeds = {} - - if geom_seeds is None: - geom_seeds = {} - - if mesh_seeds is None: - mesh_seeds = {} + + # start from zero'd con and geom seeds + + # --- con seeds --- + con_seeds_full = self.get_variable_ad_seeds() + self.zero_seed_dict(con_seeds_full) + if con_seeds is not None: + self.deep_update(con_seeds_full, con_seeds) + + # --- geom --- + geom_seeds_full = self.get_geom_ad_seeds() + self.zero_seed_dict(geom_seeds_full) + if geom_seeds is not None: + self.deep_update(geom_seeds_full, geom_seeds) + + + # --- mesh --- + mesh_seeds_full = self.get_mesh_ad_seeds() + self.zero_seed_dict(mesh_seeds_full) + if mesh_seeds is not None: + self.deep_update(mesh_seeds_full, mesh_seeds) + + # --- param --- + param_seeds_full = self.get_parameter_ad_seeds() + self.zero_seed_dict(param_seeds_full) + if param_seeds is not None: + self.deep_update(param_seeds_full, param_seeds) + + # --- ref --- + ref_seeds_full = self.get_reference_ad_seeds() + self.zero_seed_dict(ref_seeds_full) + if ref_seeds is not None: + self.deep_update(ref_seeds_full, ref_seeds) + + # The arrays can be set to zero if they don't exsist if self.DVGeo is None: dvgeo_seeds = {} @@ -3962,27 +4057,25 @@ def _execute_jac_vec_prod_fwd( if gamma_u_seeds is None: gamma_u_seeds = np.zeros((self.NUMAX, mesh_size)) - if param_seeds is None: - param_seeds = {} - - if ref_seeds is None: - ref_seeds = {} res_slice = (slice(0, mesh_size),) res_d_slice = (slice(0, num_control_surfs), slice(0, mesh_size)) res_u_slice = (slice(0, self.NUMAX), slice(0, mesh_size)) - + if mode == "AD": - # set derivative seeds + # The following body seeds are in the process of being supported + # because they are not explicitly passed, zero them out for now + body_seeds = np.zeros((self.NLMAX, 3)) + self.set_avl_fort_arr("VRTX_R_DIFF", "RL_DIFF", body_seeds) + # 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_variable_ad_seeds(con_seeds_full) + self.set_geom_ad_seeds(geom_seeds_full) 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) + self.set_parameter_ad_seeds(param_seeds_full) + self.set_reference_ad_seeds(ref_seeds_full) # 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: @@ -3996,18 +4089,18 @@ def _execute_jac_vec_prod_fwd( 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(): + if surface not in mesh_seeds_full.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)) + mesh_seeds_full[surface] = {} + mesh_seeds_full[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 + # Loop over the design variables and accumulate the sensitivity product into the mesh_seeds_full + mesh_seeds_full[surface]["mesh"] += self.DVGeo.totalSensitivityProd(dvgeo_seeds[surface], point_set_name).reshape( + mesh_seeds_full[surface]["mesh"].shape ) - self.set_mesh_ad_seeds(mesh_seeds) + self.set_mesh_ad_seeds(mesh_seeds_full) self.avl.update_surfaces_d() self.avl.get_res_d() @@ -4023,27 +4116,29 @@ def _execute_jac_vec_prod_fwd( res_d_seeds = self.get_residual_d_ad_seeds() res_u_seeds = self.get_residual_u_ad_seeds() - self.set_variable_ad_seeds(con_seeds, scale=0.0) - self.set_geom_ad_seeds(geom_seeds, scale=0.0) - self.set_mesh_ad_seeds(mesh_seeds, scale=0.0) + self.set_variable_ad_seeds(con_seeds_full, scale=0.0) + self.set_geom_ad_seeds(geom_seeds_full, scale=0.0) + self.set_mesh_ad_seeds(mesh_seeds_full, scale=0.0) self.set_gamma_ad_seeds(gamma_seeds, scale=0.0) self.set_gamma_d_ad_seeds(gamma_d_seeds, scale=0.0) self.set_gamma_u_ad_seeds(gamma_u_seeds, scale=0.0) - self.set_parameter_ad_seeds(param_seeds, scale=0.0) - self.set_reference_ad_seeds(ref_seeds, scale=0.0) + self.set_parameter_ad_seeds(param_seeds_full, scale=0.0) + self.set_reference_ad_seeds(ref_seeds_full, scale=0.0) - # TODO: remove?? + # TODO: revome this line and see if any of the tests break + # one should not have to zero the gam seed since they are set directly, right? + # also the above line should do the zeroing right, set_gamma_ad_seeds self.set_avl_fort_arr("VRTX_R_DIFF", "GAM_DIFF", gamma_seeds * 0.0, slicer=res_slice) if mode == "FD": - self.set_variable_ad_seeds(con_seeds, mode="FD", scale=step) - self.set_geom_ad_seeds(geom_seeds, mode="FD", scale=step) - self.set_mesh_ad_seeds(mesh_seeds, mode="FD", scale=step) + self.set_variable_ad_seeds(con_seeds_full, mode="FD", scale=step) + self.set_geom_ad_seeds(geom_seeds_full, mode="FD", scale=step) + self.set_mesh_ad_seeds(mesh_seeds_full, mode="FD", scale=step) self.set_gamma_ad_seeds(gamma_seeds, mode="FD", scale=step) self.set_gamma_d_ad_seeds(gamma_d_seeds, mode="FD", scale=step) self.set_gamma_u_ad_seeds(gamma_u_seeds, mode="FD", scale=step) - self.set_parameter_ad_seeds(param_seeds, mode="FD", scale=step) - self.set_reference_ad_seeds(ref_seeds, mode="FD", scale=step) + self.set_parameter_ad_seeds(param_seeds_full, mode="FD", scale=step) + self.set_reference_ad_seeds(ref_seeds_full, mode="FD", scale=step) # Since DVGeo operates entirely within the python layer we have have to do this @@ -4091,14 +4186,14 @@ def _execute_jac_vec_prod_fwd( res_d_peturbed = copy.deepcopy(self.get_avl_fort_arr("VRTX_R", "RES_D", slicer=res_d_slice)) res_u_peturbed = copy.deepcopy(self.get_avl_fort_arr("VRTX_R", "RES_U", slicer=res_u_slice)) - self.set_variable_ad_seeds(con_seeds, mode="FD", scale=-1 * step) - self.set_geom_ad_seeds(geom_seeds, mode="FD", scale=-1 * step) - self.set_mesh_ad_seeds(mesh_seeds, mode="FD", scale=-1 * step) + self.set_variable_ad_seeds(con_seeds_full, mode="FD", scale=-1 * step) + self.set_geom_ad_seeds(geom_seeds_full, mode="FD", scale=-1 * step) + self.set_mesh_ad_seeds(mesh_seeds_full, mode="FD", scale=-1 * step) self.set_gamma_ad_seeds(gamma_seeds, mode="FD", scale=-1 * step) self.set_gamma_d_ad_seeds(gamma_d_seeds, mode="FD", scale=-1 * step) self.set_gamma_u_ad_seeds(gamma_u_seeds, mode="FD", scale=-1 * step) - self.set_parameter_ad_seeds(param_seeds, mode="FD", scale=-1 * step) - self.set_reference_ad_seeds(ref_seeds, mode="FD", scale=-1 * step) + self.set_parameter_ad_seeds(param_seeds_full, mode="FD", scale=-1 * step) + self.set_reference_ad_seeds(ref_seeds_full, mode="FD", scale=-1 * step) # Set the mesh seeds back if self.DVGeo is not None and dvgeo_seeds is not None: @@ -4334,6 +4429,187 @@ def _execute_jac_vec_prod_rev( 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_direct( + self, + con_dvs: Optional[List[str]] = None, + geom_dvs: Optional[List[Tuple[(str, str)]]] = None, + param_dvs: Optional[List[str]] = None, + ref_dvs: Optional[List[str]] = None, + add_stab_derivs: Optional[bool] = False, + add_body_axis_derivs: Optional[bool] = False, + add_consurf_derivs: Optional[bool] = False, + print_timings: Optional[bool] = False, + ) -> Dict[str, Dict[str, float]]: + """Run the sensitivities of the input functionals in adjoint mode + + Returns: + sens: a nested dictionary of sensitivities. The first key is the function and the next keys are for the design variables. + """ + + sens = {} + + if self.get_avl_fort_arr("CASE_L", "LTIMING"): + print_timings = True + + + def make_unit_dicts(surfaces, surf_name, key): + """Yield dicts with each scalar element set to 1.0, others zero.""" + val = surfaces[surf_name][key] + + if np.isscalar(val): + yield {surf_name: {key: 1.0}} + + elif isinstance(val, np.ndarray): + for idx in np.ndindex(val.shape): + out = np.zeros_like(val, dtype=float) + out[idx] = 1.0 + yield {surf_name: {key: out}} + + elif isinstance(val, list): + # list of arrays (like xasec, sasec) + for i, arr in enumerate(val): + for idx in np.ndindex(arr.shape): + out = [np.zeros_like(a, dtype=float) for a in val] + out[i][idx] = 1.0 + yield {surf_name: {key: out}} + + # convenience function to get the correct seeds: + def get_dv_seed_list(): + seeds = {} + + con_seed_list = [] + if con_dvs is not None: + for con in con_dvs: + con_seed_list.append({con:1.0}) + seeds['con_seeds'] = con_seed_list + + geom_seed_list = [] + if geom_dvs is not None: + full_geom_seeds = self.get_geom_ad_seeds() + for dv in geom_dvs: + for dv_dict in make_unit_dicts(full_geom_seeds, dv[0], dv[1]): + geom_seed_list.append(dv_dict) + seeds['geom_seeds'] = geom_seed_list + + param_seed_list = [] + if param_dvs is not None: + for param in param_dvs: + param_seed_list.append({param:1.0}) + seeds['param_seeds'] = param_seed_list + + ref_seed_list = [] + if ref_dvs is not None: + for ref in ref_dvs: + ref_seed_list.append({ref:1.0}) + seeds['ref_seeds'] = ref_seed_list + + # for dv in geom_dvs: + # full_seed = full_geom_seeds[dv[0]][dv[1]] + # if isinstance(full_seed, np.ndarray): + + # for i in range(full_seed.size): + + + return seeds + + + + seeds = get_dv_seed_list() + + + def add_deriv_to_dict(dv_seed, sens_dict, seeds): + + def add_sub_dicts(in_dict, out_dict): + # recurse until we don't reach a dictionary + for key in in_dict: + val = in_dict[key] + + if isinstance(val, dict): + if key not in out_dict: + out_dict[key] = {} + add_sub_dicts(val, out_dict[key]) + elif isinstance(val, np.ndarray): + if key not in out_dict: + out_dict[key] = np.zeros_like(val,dtype=float) + out_dict[key] += val * seeds[func_key] + elif isinstance(val, list): + if key not in out_dict: + out_dict[key] = [np.zeros_like(a,dtype=float) for a in val] + for i, arr in enumerate(val): + out_dict[key][i] += arr *seeds[func_key] + else: # scalar + out_dict[key] = out_dict.get(key, 0.0) + val * seeds[func_key] + + for func_key in seeds: + + if func_key not in sens_dict: + sens_dict[func_key] = {} + + add_sub_dicts(dv_seed, sens_dict[func_key]) + + for dv_type in seeds: + for dv_seed in seeds[dv_type]: + + time_last = time.time() + # compute dR/dX + _, pRpX, _, _, _, pR_dpX, pR_upX = self._execute_jac_vec_prod_fwd( + **{dv_type: dv_seed} + ) + if print_timings: + print(f"Time to get RHS: {time.time() - time_last}") + time_last = time.time() + + # now solve the direct equation + # the RHS is absorbing the negative of the total deriv equation + # DONT forget the negative + self.set_residual_ad_seeds(-1 * pRpX) + + if add_stab_derivs or add_body_axis_derivs: + self.set_residual_u_ad_seeds(-1 * pR_upX) + solve_res_u_drt = True + else: + solve_res_u_drt = False + + if add_consurf_derivs: + self.set_residual_d_ad_seeds(-1 * pR_dpX) + solve_res_d_drt = True + else: + solve_res_d_drt = False + + # solve only the state direct method + self.avl.solve_direct(solve_res_u_drt, solve_res_d_drt) + + if print_timings: + print(f"Time to solve direct: {time.time() - time_last}") + time_last = time.time() + + # get the resulting adjoint vector (dU/dX) from fortran + dUdX = self.get_gamma_ad_seeds() + + # this is harmless even if we aren't solving for extra derivatives + dU_ddX = self.get_gamma_d_ad_seeds() + dU_udX = self.get_gamma_u_ad_seeds() + + # do the last partial derivative to comput the totals with the direct vector + func_seeds, _, consurf_derivs_seeds, stab_derivs_seeds, body_axis_derivs_seeds, _, _, = self._execute_jac_vec_prod_fwd( + gamma_seeds=dUdX, gamma_d_seeds=dU_ddX, gamma_u_seeds=dU_udX, **{dv_type: dv_seed} + ) + + + add_deriv_to_dict(dv_seed, sens, func_seeds) + + if add_consurf_derivs: + add_deriv_to_dict(dv_seed, sens, consurf_derivs_seeds) + + if add_stab_derivs: + add_deriv_to_dict(dv_seed, sens, stab_derivs_seeds) + + if add_body_axis_derivs: + add_deriv_to_dict(dv_seed, sens, body_axis_derivs_seeds) + + return sens + + def execute_run_sensitivities( self, funcs: List[str], @@ -4362,16 +4638,14 @@ def execute_run_sensitivities( # set up and solve the adjoint for each function for func in funcs: sens[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() + + # get the RHS of the adjoint equation (pFpU) _, _, _, _, 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() - # self.clear_ad_seeds() # u solver adjoint equation with RHS self.set_gamma_ad_seeds(-1 * pfpU) solve_gamma_u_adj = False diff --git a/src/aoper.f b/src/aoper.f index 53f238f9..38685ca5 100644 --- a/src/aoper.f +++ b/src/aoper.f @@ -2353,4 +2353,46 @@ subroutine solve_adjoint(solve_stab_deriv_adj, solve_con_surf_adj) enddo endif + end !subroutine solve_adjoint + + subroutine solve_direct(solve_stab_deriv_drt, solve_con_surf_drt) + ! solves the direct equation with the jacobian to compute dx/du + use avl_heap_inc + use avl_heap_diff_inc + include "AVL.INC" + include "AVL_ad_seeds.inc" + integer i + logical :: solve_stab_deriv_drt, solve_con_surf_drt + + CALL SETUP + IF(.NOT.LAIC) THEN + call factor_AIC + ENDIF + + ! assume we have run the other routines to put dr/dx is in gam_diff + do i =1,NVOR + GAM_diff(i) = RES_diff(i) + enddo + + CALL BAKSUB(NVOR,NVOR,AICN_LU,IAPIV,GAM_diff) + + if (solve_con_surf_drt) then + DO IC = 1, NCONTROL + do i =1,NVOR + GAM_D_diff(i,IC) = RES_D_diff(i,IC) + enddo + CALL BAKSUB(NVOR,NVOR,AICN_LU,IAPIV,GAM_D_diff(:,IC)) + enddo + endif + + if (solve_stab_deriv_drt) then + DO IU = 1,6 + do i =1,NVOR + GAM_U_diff(i,IU) = RES_U_diff(i,IU) + enddo + + CALL BAKSUB(NVOR,NVOR,AICN_LU,IAPIV,GAM_U_diff(:,IU)) + enddo + endif + end !subroutine solve_adjoint \ No newline at end of file diff --git a/src/f2py/libavl.pyf b/src/f2py/libavl.pyf index ce5919af..67244940 100644 --- a/src/f2py/libavl.pyf +++ b/src/f2py/libavl.pyf @@ -117,6 +117,10 @@ python module libavl ! in subroutine solve_adjoint(solve_stab_deriv_adj, solve_con_surf_adj) logical :: solve_stab_deriv_adj, solve_con_surf_adj end subroutine solve_adjoint + + subroutine solve_direct(solve_stab_deriv_drt, solve_con_surf_drt) + logical :: solve_stab_deriv_drt, solve_con_surf_drt + end subroutine solve_direct subroutine cpoml(save_file) logical :: save_file diff --git a/tests/test_body_axis_derivs_partial_derivs.py b/tests/test_body_axis_derivs_partial_derivs.py index a471e97e..8d135127 100644 --- a/tests/test_body_axis_derivs_partial_derivs.py +++ b/tests/test_body_axis_derivs_partial_derivs.py @@ -39,7 +39,7 @@ def tearDown(self): print(f"{self.id()} Memory usage: {mb_memory:.2f} MB") def test_fwd_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: bd_d = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[4] bd_d_fd = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0}, mode="FD", step=1e-6)[4] @@ -64,7 +64,7 @@ def test_rev_aero_constraint(self): self.ovl_solver.clear_ad_seeds_fast() - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: body_axis_deriv_seeds_fwd= self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[4] body_axis_deriv_sum = 0.0 diff --git a/tests/test_consurf_partial_derivs.py b/tests/test_consurf_partial_derivs.py index bb82ec26..2b5bd2a1 100644 --- a/tests/test_consurf_partial_derivs.py +++ b/tests/test_consurf_partial_derivs.py @@ -39,7 +39,7 @@ def tearDown(self): print(f"{self.id()} Memory usage: {mb_memory:.2f} MB") def test_fwd_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: res_d_seeds = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[5] res_d_seeds_FD = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0}, mode="FD", step=1e-5)[ @@ -53,7 +53,7 @@ def test_fwd_aero_constraint(self): ) def test_rev_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: res_d_seeds = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[5] num_gamma = self.ovl_solver.get_mesh_size() @@ -197,7 +197,7 @@ def tearDown(self): print(f"{self.id()} Memory usage: {mb_memory:.2f} MB") def test_fwd_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: cs_d = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[2] cs_d_fd = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0}, mode="FD", step=1e-8)[2] @@ -223,7 +223,7 @@ def test_rev_aero_constraint(self): self.ovl_solver.clear_ad_seeds_fast() - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: cs_deriv_seeds_fwd = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[2] cs_deriv_sum = 0.0 diff --git a/tests/test_mesh_totals.py b/tests/test_mesh_totals.py index 62a3f377..aa74b41a 100644 --- a/tests/test_mesh_totals.py +++ b/tests/test_mesh_totals.py @@ -122,10 +122,10 @@ def test_aero_constraint(self): sens_sd = self.ovl_solver.execute_run_sensitivities([], stab_derivs=stab_derivs, print_timings=False) sens_bd = self.ovl_solver.execute_run_sensitivities([], body_axis_derivs=body_axis_derivs, print_timings=False) - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: # for con_key in ['beta']: func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif( - [con_key], {}, {}, {}, {}, step=1.0e-5 + [con_key], {}, {}, {}, {}, step=1.0e-6 ) # for func_key in func_vars: @@ -233,7 +233,7 @@ def test_mesh(self): print_timings=False, ) - # for con_key in self.ovl_solver.con_var_to_fort_var: + # for con_key in self.ovl_solver.con_var_list: sens_FD = {} for surf_key in self.ovl_solver.surf_geom_to_fort_var: sens_FD[surf_key] = {} diff --git a/tests/test_partial_derivs.py b/tests/test_partial_derivs.py index 4869b845..89ffa569 100644 --- a/tests/test_partial_derivs.py +++ b/tests/test_partial_derivs.py @@ -39,7 +39,7 @@ def tearDown(self): print(f"{self.id():80} Memory usage: {mb_memory:.2f} MB") def test_fwd_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: func_seeds = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0}, geom_seeds={})[0] func_seeds_FD = self.ovl_solver._execute_jac_vec_prod_fwd( @@ -66,7 +66,7 @@ def test_fwd_aero_constraint(self): ) def test_rev_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: self.ovl_solver.clear_ad_seeds_fast() func_seeds_fwd = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0}, geom_seeds={})[0] @@ -353,7 +353,7 @@ def tearDown(self): print(f"{self.id():80} Memory usage: {mb_memory:.2f} MB") def test_fwd_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: res_seeds_FD = self.ovl_solver._execute_jac_vec_prod_fwd( con_seeds={con_key: 1.0}, geom_seeds={}, mode="FD", step=1e-8 )[1] @@ -370,7 +370,7 @@ def test_rev_aero_constraint(self): self.ovl_solver.clear_ad_seeds_fast() - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: res_seeds_fwd = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[1] # do dot product diff --git a/tests/test_pygeo.py b/tests/test_pygeo.py index 1e802d6d..05709d27 100644 --- a/tests/test_pygeo.py +++ b/tests/test_pygeo.py @@ -423,7 +423,7 @@ def sweep(val, geo): np.testing.assert_allclose( dvgeo_var_dot, func_dot, - atol=5e-9, + atol=7e-9, err_msg=f"{func_key} wrt {surf_key}:{dvgeo_var_key:10}", ) else: diff --git a/tests/test_stab_derivs_partial_derivs.py b/tests/test_stab_derivs_partial_derivs.py index 4491f1c7..7869a873 100644 --- a/tests/test_stab_derivs_partial_derivs.py +++ b/tests/test_stab_derivs_partial_derivs.py @@ -38,7 +38,7 @@ def tearDown(self): print(f"{self.id()} Memory usage: {mb_memory:.2f} MB") def test_fwd_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: res_u_seeds = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[6] res_u_seeds_FD = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0}, mode="FD", step=1e-5)[ @@ -52,7 +52,7 @@ def test_fwd_aero_constraint(self): ) def test_rev_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: # for con_key in ["beta", "beta"]: res_u_seeds = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[6] @@ -214,7 +214,7 @@ def tearDown(self): print(f"{self.id()} Memory usage: {mb_memory:.2f} MB") def test_fwd_aero_constraint(self): - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: sd_d = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[3] sd_d_fd = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0}, mode="FD", step=1e-6)[3] @@ -249,7 +249,7 @@ def test_rev_aero_constraint(self): self.ovl_solver.clear_ad_seeds_fast() - for con_key in self.ovl_solver.con_var_to_fort_var: + for con_key in self.ovl_solver.con_var_list: stab_deriv_seeds_fwd = self.ovl_solver._execute_jac_vec_prod_fwd(con_seeds={con_key: 1.0})[3] stab_deriv_sum = 0.0 diff --git a/tests/test_total_derivs.py b/tests/test_total_derivs.py index e7b87c3b..ba9f3352 100644 --- a/tests/test_total_derivs.py +++ b/tests/test_total_derivs.py @@ -28,11 +28,11 @@ class TestTotals(unittest.TestCase): # TODO: beta derivatives likely wrong def setUp(self): - self.ovl_solver = OVLSolver(geo_file=geom_file) - self.ovl_solver.set_variable("alpha", 5.0) - self.ovl_solver.set_variable("beta", 0.0) - self.ovl_solver.set_parameter("Mach", 0.8) - self.ovl_solver.execute_run() + self.ovl = OVLSolver(geo_file=geom_file) + self.ovl.set_variable("alpha", 5.0) + self.ovl.set_variable("beta", 0.0) + self.ovl.set_parameter("Mach", 0.8) + self.ovl.execute_run() def tearDown(self): # Get the memory usage of the current process using psutil @@ -45,42 +45,42 @@ def finite_dif(self, con_list, geom_seeds, param_seeds, ref_seeds, step=1e-7): for con in con_list: con_seeds[con] = 1.0 - self.ovl_solver.set_variable_ad_seeds(con_seeds, mode="FD", scale=step) - self.ovl_solver.set_geom_ad_seeds(geom_seeds, mode="FD", scale=step) - self.ovl_solver.set_parameter_ad_seeds(param_seeds, mode="FD", scale=step) - self.ovl_solver.set_reference_ad_seeds(ref_seeds, mode="FD", scale=step) - - 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() - - self.ovl_solver.set_variable_ad_seeds(con_seeds, mode="FD", scale=-1 * step) - self.ovl_solver.set_geom_ad_seeds(geom_seeds, mode="FD", scale=-1 * step) - self.ovl_solver.set_parameter_ad_seeds(param_seeds, mode="FD", scale=-1 * step) - self.ovl_solver.set_reference_ad_seeds(ref_seeds, mode="FD", scale=-1 * step) - - 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() + self.ovl.set_variable_ad_seeds(con_seeds, mode="FD", scale=step) + self.ovl.set_geom_ad_seeds(geom_seeds, mode="FD", scale=step) + self.ovl.set_parameter_ad_seeds(param_seeds, mode="FD", scale=step) + self.ovl.set_reference_ad_seeds(ref_seeds, mode="FD", scale=step) + + self.ovl.avl.update_surfaces() + self.ovl.avl.get_res() + self.ovl.avl.exec_rhs() + self.ovl.avl.get_res() + self.ovl.avl.velsum() + self.ovl.avl.aero() + # self.ovl.execute_run() + coef_data_peturb = self.ovl.get_total_forces() + consurf_derivs_peturb = self.ovl.get_control_stab_derivs() + stab_deriv_derivs_peturb = self.ovl.get_stab_derivs() + body_axis_deriv_petrub = self.ovl.get_body_axis_derivs() + body_forces_peturb = self.ovl.get_body_forces() + + self.ovl.set_variable_ad_seeds(con_seeds, mode="FD", scale=-1 * step) + self.ovl.set_geom_ad_seeds(geom_seeds, mode="FD", scale=-1 * step) + self.ovl.set_parameter_ad_seeds(param_seeds, mode="FD", scale=-1 * step) + self.ovl.set_reference_ad_seeds(ref_seeds, mode="FD", scale=-1 * step) + + self.ovl.avl.update_surfaces() + self.ovl.avl.get_res() + self.ovl.avl.exec_rhs() + self.ovl.avl.get_res() + self.ovl.avl.velsum() + self.ovl.avl.aero() + # self.ovl.execute_run() + + coef_data = self.ovl.get_total_forces() + consurf_derivs = self.ovl.get_control_stab_derivs() + stab_deriv_derivs = self.ovl.get_stab_derivs() + body_axis_deriv = self.ovl.get_body_axis_derivs() + body_forces = self.ovl.get_body_forces() body_func_seeds = {} for body in body_forces: @@ -111,17 +111,17 @@ def finite_dif(self, con_list, geom_seeds, param_seeds, ref_seeds, step=1e-7): def test_aero_constraint(self): # compare the analytical gradients with finite difference for each constraint and function - 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_funcs = self.ovl_solver.execute_run_sensitivities(func_vars) - sens_sd = self.ovl_solver.execute_run_sensitivities([], stab_derivs=stab_derivs, print_timings=False) - sens_bd = self.ovl_solver.execute_run_sensitivities([], body_axis_derivs=body_axis_derivs, print_timings=False) - - for con_key in self.ovl_solver.con_var_to_fort_var: + func_vars = self.ovl.case_var_to_fort_var + stab_derivs = self.ovl.case_stab_derivs_to_fort_var + body_axis_derivs = self.ovl.case_body_derivs_to_fort_var + sens_funcs = self.ovl.execute_run_sensitivities(func_vars) + sens_sd = self.ovl.execute_run_sensitivities([], stab_derivs=stab_derivs, print_timings=False) + sens_bd = self.ovl.execute_run_sensitivities([], body_axis_derivs=body_axis_derivs, print_timings=False) + + for con_key in self.ovl.con_var_list: # for con_key in ['beta']: func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif( - [con_key], {}, {}, {}, step=1.0e-5 + [con_key], {}, {}, {}, step=1.0e-6 ) # for func_key in func_vars: @@ -134,7 +134,7 @@ def test_aero_constraint(self): # print(f"{func_key:5} wrt {con_key:5} | AD:{ad_dot: 5e} FD:{fd_dot: 5e} rel err:{rel_err:.2e}") - tol = 1e-8 + tol = 5e-8 if np.abs(ad_dot) < tol or np.abs(fd_dot) < tol: # If either value is basically zero, use an absolute tolerance np.testing.assert_allclose( @@ -147,7 +147,7 @@ def test_aero_constraint(self): np.testing.assert_allclose( ad_dot, fd_dot, - rtol=5e-5, + rtol=5e-4, err_msg=f"func_key {func_key} w.r.t. {con_key}", ) @@ -161,7 +161,7 @@ def test_aero_constraint(self): # f"{func_key} wrt {con_key} | AD:{ad_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" # ) - tol = 1e-8 + tol = 5e-8 if np.abs(ad_dot) < tol or np.abs(func_dot) < tol: # If either value is basically zero, use an absolute tolerance np.testing.assert_allclose( @@ -188,7 +188,7 @@ def test_aero_constraint(self): # f"{func_key} wrt {con_key} | AD:{ad_dot: 5e} FD:{func_dot: 5e} rel err:{rel_err:.2e}" # ) - tol = 1e-8 + tol = 5e-8 if np.abs(ad_dot) < tol or np.abs(func_dot) < tol: # If either value is basically zero, use an absolute tolerance np.testing.assert_allclose( @@ -209,20 +209,19 @@ def test_geom(self): # compare the analytical gradients with finite difference for each # geometric variable and function - 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() + surf_key = list(self.ovl.surf_geom_to_fort_var.keys())[0] + geom_vars = self.ovl.surf_geom_to_fort_var[surf_key] + cs_names = self.ovl.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)) + for func_key in self.ovl.case_derivs_to_fort_var: + consurf_vars.append(self.ovl._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 + func_vars = self.ovl.case_var_to_fort_var + stab_derivs = self.ovl.case_stab_derivs_to_fort_var + body_axis_derivs = self.ovl.case_body_derivs_to_fort_var - sens = self.ovl_solver.execute_run_sensitivities( + sens = self.ovl.execute_run_sensitivities( func_vars, consurf_derivs=consurf_vars, stab_derivs=stab_derivs, @@ -230,12 +229,12 @@ def test_geom(self): print_timings=False, ) - # for con_key in self.ovl_solver.con_var_to_fort_var: + # for con_key in self.ovl.con_var_list: sens_FD = {} - for surf_key in self.ovl_solver.surf_geom_to_fort_var: + for surf_key in self.ovl.surf_geom_to_fort_var: sens_FD[surf_key] = {} for geom_key in geom_vars: - arr = self.ovl_solver.get_surface_param(surf_key, geom_key) + arr = self.ovl.get_surface_param(surf_key, geom_key) np.random.seed(arr.size) rand_arr = np.random.rand(*arr.shape) rand_arr /= np.linalg.norm(rand_arr) @@ -353,12 +352,12 @@ def test_geom(self): def test_params(self): # compare the analytical gradients with finite difference for each constraint and function - func_vars = self.ovl_solver.case_var_to_fort_var - stab_derivs = self.ovl_solver.case_stab_derivs_to_fort_var + func_vars = self.ovl.case_var_to_fort_var + stab_derivs = self.ovl.case_stab_derivs_to_fort_var - sens = self.ovl_solver.execute_run_sensitivities(func_vars, stab_derivs=stab_derivs) + sens = self.ovl.execute_run_sensitivities(func_vars, stab_derivs=stab_derivs) - for param_key in self.ovl_solver.param_idx_dict: + for param_key in self.ovl.param_idx_dict: func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif( [], {}, {param_key: 1.0}, {}, step=1.0e-6 ) @@ -415,12 +414,12 @@ def test_params(self): def test_ref(self): # compare the analytical gradients with finite difference for each constraint and function - func_vars = self.ovl_solver.case_var_to_fort_var - stab_derivs = self.ovl_solver.case_stab_derivs_to_fort_var + func_vars = self.ovl.case_var_to_fort_var + stab_derivs = self.ovl.case_stab_derivs_to_fort_var - sens = self.ovl_solver.execute_run_sensitivities(func_vars, stab_derivs=stab_derivs) + sens = self.ovl.execute_run_sensitivities(func_vars, stab_derivs=stab_derivs) - for ref_key in self.ovl_solver.ref_var_to_fort_var: + for ref_key in self.ovl.ref_var_to_fort_var: # for con_key in ['beta']: func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif( [], {}, {}, {ref_key: 1.0}, step=1.0e-5 @@ -478,6 +477,48 @@ def test_ref(self): err_msg=f"{func_key} wrt {ref_key}", ) +class TestDirectVsAdjoint(unittest.TestCase): + + def setUp(self): + self.ovl = OVLSolver(geo_file=geom_file) + self.ovl.set_variable("alpha", 5.0) + self.ovl.set_variable("beta", 0.0) + self.ovl.set_parameter("Mach", 0.8) + self.ovl.execute_run() + + 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()} Memory usage: {mb_memory:.2f} MB") + + def test_aero_constraint(self): + # compare the analytical gradients with finite difference for each constraint and function + # func_vars = self.ovl.case_var_to_fort_var + stab_derivs = ["dCL/dalpha"] + # body_axis_derivs = self.ovl.case_body_derivs_to_fort_var + funcs = ["CL", "CD", "Cm"] + con_dvs = ["alpha"] + ref_dvs = ["Sref"] + param_dvs = ["Mach"] + geom_dvs = [("Wing", "scale"),("Wing", "chords")] + + sens_adjoint = self.ovl.execute_run_sensitivities(funcs, stab_derivs=stab_derivs ) + + sens_direct = self.ovl.execute_run_sensitivities_direct(geom_dvs=geom_dvs, con_dvs=con_dvs, ref_dvs=ref_dvs, param_dvs=param_dvs, add_stab_derivs=True) + + for func in funcs + stab_derivs: + for dv in geom_dvs: + # print(f"Adjoint d{func}/d{[dv[0]]} {[dv[1]]} {sens_adjoint[func][dv[0]][dv[1]]}") + # print(f"Direct d{func}/d{[dv[0]]} {[dv[1]]} {sens_direct[func][dv[0]][dv[1]]}") + np.testing.assert_allclose(sens_direct[func][dv[0]][dv[1]], sens_adjoint[func][dv[0]][dv[1]], 1e-14, 1e-14) + + + for dv in con_dvs + ref_dvs + param_dvs: + # print(f"Adjoint d{func}/d{dv} {sens_adjoint[func][dv]}") + # print(f"Direct d{func}/d{dv} {sens_direct[func][dv]}") + np.testing.assert_allclose(sens_direct[func][dv], sens_adjoint[func][dv], 1e-14, 1e-14) + if __name__ == "__main__": unittest.main()