From 7c11577865e0708830691a6a2aa5954c413db556 Mon Sep 17 00:00:00 2001 From: Sean Engelstad Date: Tue, 22 Oct 2024 15:54:50 -0400 Subject: [PATCH] Compute Panel Dimensions in TACS for Stiffened Panel Optimizations (#339) * fix typo in test bdf utils callbacks * update the nlbgs driver * fix plot_manager with new history file setup * update the plot manager function * update the tacs interface to work with panel length and width constraints * add temp debug statements to tacs interface for ML oneway opt * fix sparse gradients for panel length, width in tacs closed form * update the oneway struct driver to make init transfer default true * fix the sparse gradient for adjacency constraints * remove some of the print statements * remove printouts in funtofem * make sparse opt manager default True * remove plot hist with option * fix linear gradients in pyoptsparse * add new sparse linear constraint settings for efficient memory management * None checks --- funtofem/driver/oneway_struct_driver.py | 1 + funtofem/interface/tacs_interface.py | 339 ++++++++++++------ funtofem/model/composite_function.py | 20 +- funtofem/model/funtofem_model.py | 131 ++++++- funtofem/model/variable.py | 9 + funtofem/optimization/optimization_manager.py | 72 ++-- 6 files changed, 439 insertions(+), 133 deletions(-) diff --git a/funtofem/driver/oneway_struct_driver.py b/funtofem/driver/oneway_struct_driver.py index db4b2490..ca62a3e3 100644 --- a/funtofem/driver/oneway_struct_driver.py +++ b/funtofem/driver/oneway_struct_driver.py @@ -538,6 +538,7 @@ def _transfer_fixed_aero_loads(self): # the mesh for the new structural shape body.transfer_loads(scenario) body.transfer_heat_flux(scenario) + return def solve_forward(self): diff --git a/funtofem/interface/tacs_interface.py b/funtofem/interface/tacs_interface.py index 2e30c2ed..78a6db06 100644 --- a/funtofem/interface/tacs_interface.py +++ b/funtofem/interface/tacs_interface.py @@ -20,7 +20,12 @@ limitations under the License. """ -__all__ = ["TacsInterface", "TacsSteadyInterface"] +__all__ = [ + "TacsInterface", + "TacsSteadyInterface", + "TacsPanelDimensions", + "TacsOutputGenerator", +] from mpi4py import MPI from tacs import pytacs, TACS, functions @@ -122,7 +127,10 @@ class TacsSteadyInterface(SolverInterface): A base class to do coupled steady simulations with TACS """ - PANEL_LENGTH_CONSTR = "length" + LENGTH_VAR = "LENGTH" + LENGTH_CONSTR = "PANEL-LENGTH" + WIDTH_VAR = "WIDTH" + WIDTH_CONSTR = "PANEL-WIDTH" def __init__( self, @@ -138,8 +146,8 @@ def __init__( nprocs=None, relaxation_scheme: AitkenRelaxationTacs = None, debug=False, - panel_length_constraint=None, struct_loads_file=None, + tacs_panel_dimensions=None, ): """ Initialize the TACS implementation of the SolverInterface for the FUNtoFEM @@ -214,6 +222,8 @@ def __init__( if var.analysis_type == "structural": self.struct_variables.append(var) + self.tacs_panel_dimensions = tacs_panel_dimensions + # Set the assembler object - if it exists or not self._initialize_variables( model, @@ -244,12 +254,6 @@ def __init__( # Generate output self.gen_output = gen_output - # create panel length constraints - self.panel_length_constraint = panel_length_constraint - self.panel_length_name = "PanelLengthCon_PanelLength" - - self._eval_panel_length(forward=True, adjoint=True) - # Debug flag self._debug = debug if self.comm.rank != 0: @@ -517,6 +521,19 @@ def set_variables(self, scenario, bodies): The bodies in the model """ + if self.tacs_panel_dimensions is not None: + # compute panel length and width (checks if we need to inside the object) + self.tacs_panel_dimensions.compute_panel_length( + self.assembler, + self.struct_variables, + self.LENGTH_CONSTR, + self.LENGTH_VAR, + ) + self.tacs_panel_dimensions.compute_panel_width( + self.assembler, self.struct_variables, self.WIDTH_CONSTR, self.WIDTH_VAR + ) + + # write F2F variable values into TACS if self.tacs_proc: # Set the design variable values on the processors that # have an instance of TACSAssembler. @@ -532,6 +549,10 @@ def set_variables(self, scenario, bodies): self.assembler.setDesignVars(xvec) + # also copy to the constraints local vectors + if self.tacs_panel_dimensions is not None: + self.tacs_panel_dimensions.setDesignVars(xvec) + return def set_functions(self, scenario, bodies): @@ -569,10 +590,13 @@ def get_functions(self, scenario, bodies): # Evaluate the list of the functions of interest feval = None if self.tacs_proc: + # print(f"evalFunctions", flush=True) feval = self.assembler.evalFunctions(self.scenario_data[scenario].func_list) + # print(f"\tDone with evalFunctions", flush=True) # Broacast the list across all processors - not just structural procs feval = self.comm.bcast(feval, root=0) + # print(f"feval = {feval}", flush=True) # Set the function values on all processors for i, func in enumerate(scenario.functions): @@ -595,12 +619,14 @@ def get_function_gradients(self, scenario, bodies): The bodies in the model """ + # print(f"get function gradients start", flush=True) func_grad = self.scenario_data[scenario].func_grad for ifunc, func in enumerate(scenario.functions): for i, var in enumerate(self.struct_variables): # func.set_gradient_component(var, func_grad[ifunc][i]) func.add_gradient_component(var, func_grad[ifunc][i]) + # print(f"\tdoneget function gradients start", flush=True) return @@ -745,7 +771,9 @@ def iterate(self, scenario, bodies, step): self.res.axpy(-1.0, self.ext_force) # Solve for the update + # print(f"solve linear static analysis", flush=True) self.gmres.solve(self.res, self.update) + # print(f"\tsolved linear static analysis", flush=True) if self.comm.rank == 0 and self.aitken_debug: print(f"TACS iterate step: {step}", flush=True) @@ -812,85 +840,9 @@ def iterate(self, scenario, bodies, step): + scenario.T_ref ) - return fail - - def _eval_panel_length(self, forward=True, adjoint=True): - # whether to do the panel length constraint - _has_panel_length = None - if self.comm.rank == 0: - _has_panel_length = self.panel_length_constraint is not None - _has_panel_length = self.comm.bcast(_has_panel_length, root=0) - - # print(f"has panel length rank {self.comm.rank} = {_has_panel_length}", flush=True) - - # compute the panel length constraint - if _has_panel_length: - if forward: - funcs_dict = None - # print(f"rank {self.comm.rank} enter forward", flush=True) - if self.assembler is not None: - funcs = {} - funcs_dict = {} - ct = 0 - self.panel_length_constraint.evalConstraints(funcs) - # print(f"inside rank 0 check post eval constraints", flush=True) - for func in self.model.composite_functions: - if self.PANEL_LENGTH_CONSTR in func.name: - # assume name of form f"{self.PANEL_LENGTH_CONSTR}-fnum" - func.value = funcs[self.panel_length_name][ct] - ct += 1 - funcs_dict[func.full_name] = func.value - # print(f"inside rank 0 check : funcs dict = {funcs_dict}", flush=True) - - # broadcast the funcs dict to other processors - funcs_dict = self.comm.bcast(funcs_dict, root=0) - - # print(f"rank {self.comm.rank} : forward funcs dict = {funcs_dict}", flush=True) - - for func in self.model.composite_functions: - if func.full_name in list(funcs_dict.keys()): - func.value = funcs_dict[func.full_name] - - # compute the panel length constraint - if adjoint: - grads_dict = None - if self.assembler is not None: - funcSens = {} - if self.comm.rank == 0: - grads_dict = {} - ifunc = 0 - self.panel_length_constraint.evalConstraintsSens(funcSens) - for func in self.model.composite_functions: - if ( - self.PANEL_LENGTH_CONSTR in func.name - and self.comm.rank == 0 - ): - - grads_dict[func.full_name] = {} - - # assume name of form f"{self.PANEL_LENGTH_CONSTR}-fnum" - for ivar, var in enumerate(self.struct_variables): - func.derivatives[var] = funcSens[ - self.panel_length_name - ]["struct"].toarray()[ifunc, ivar] - grads_dict[func.full_name][var.full_name] = ( - func.derivatives[var] - ) - - ifunc += 1 - - # broadcast the funcs dict to other processors - grads_dict = self.comm.bcast(grads_dict, root=0) + # print(f"done with iterate", flush=True) - # print(f"rank {self.comm.rank} : grads dict = {grads_dict}", flush=True) - - for func in self.model.composite_functions: - if func.full_name in list(grads_dict.keys()): - for ivar, var in enumerate(self.struct_variables): - func.derivatives[var] = grads_dict[func.full_name][ - var.full_name - ] - # print(f"rank {self.comm.rank} : d{func.full_name}/d{var.full_name} = {func.derivatives[var]}", flush=True) + return fail def post(self, scenario, bodies): """ @@ -906,8 +858,6 @@ def post(self, scenario, bodies): list of FUNtoFEM bodies """ - self._eval_panel_length(adjoint=False) - # update solution and dv1 state (like _updateAssemblerVars() in pytacs) self.set_variables(scenario, bodies) if self.tacs_proc: @@ -1102,7 +1052,9 @@ def iterate_adjoint(self, scenario, bodies, step): self.assembler.applyBCs(self.res) # Solve structural adjoint equation + # print(f"linear adjoint solve", flush=True) self.gmres.solve(self.res, psi[ifunc]) + # print(f"\tdone withlinear adjoint solve", flush=True) # Aitken adjoint step if self.use_aitken: @@ -1181,6 +1133,8 @@ def iterate_adjoint(self, scenario, bodies, step): self.thermal_index :: ndof ].astype(body.dtype) + # print(f"done with iterate adjoint", flush=True) + return fail def post_adjoint(self, scenario, bodies): @@ -1202,7 +1156,7 @@ def post_adjoint(self, scenario, bodies): list of FUNtoFEM bodies """ - self._eval_panel_length(adjoint=True) + # print(f"begin post adjoint", flush=True) func_grad = [] if self.tacs_proc: @@ -1229,6 +1183,8 @@ def post_adjoint(self, scenario, bodies): # Broadcast the gradients to all processors self.scenario_data[scenario].func_grad = self.comm.bcast(func_grad, root=0) + # print(f"\tdone with post adjoint", flush=True) + return def get_coordinate_derivatives(self, scenario, bodies, step): @@ -1289,6 +1245,8 @@ def create_from_bdf( add_loads=True, # whether it will try to add loads or not relaxation_scheme: AitkenRelaxationTacs = None, struct_loads_file=None, + panel_length_dv_index=None, + panel_width_dv_index=None, ): """ Class method to create a TacsSteadyInterface instance using the pytacs BDF loader @@ -1326,7 +1284,11 @@ def create_from_bdf( assembler = None f5 = None Fvec = None - panel_length_constraint = None + tacs_panel_dimensions: TacsPanelDimensions = TacsPanelDimensions( + comm=comm, + panel_length_dv_index=panel_length_dv_index, + panel_width_dv_index=panel_width_dv_index, + ) if world_rank < nprocs: # Create the assembler class fea_assembler = pytacs.pyTACS(bdf_file, tacs_comm, options=struct_options) @@ -1370,20 +1332,21 @@ def create_from_bdf( Fvec = addLoadsFromBDF(fea_assembler) # Fvec = None - # make the panel length constraint object - has_panel_length_funcs = any( - [ - cls.PANEL_LENGTH_CONSTR in comp_func.name - for comp_func in model.composite_functions - ] - ) - if has_panel_length_funcs: - panel_length_constraint = fea_assembler.createPanelLengthConstraint( - "PanelLengthCon" + if panel_length_dv_index is not None and tacs_panel_dimensions is not None: + tacs_panel_dimensions.panel_length_constr = ( + fea_assembler.createPanelLengthConstraint("PanelLengthCon") + ) + tacs_panel_dimensions.panel_length_constr.addConstraint( + cls.LENGTH_CONSTR, + dvIndex=tacs_panel_dimensions.panel_length_dv_index, + ) + if panel_width_dv_index is not None and tacs_panel_dimensions is not None: + tacs_panel_dimensions.panel_width_constr = ( + fea_assembler.createPanelWidthConstraint("PanelWidthCon") + ) + tacs_panel_dimensions.panel_width_constr.addConstraint( + cls.WIDTH_CONSTR, dvIndex=tacs_panel_dimensions.panel_width_dv_index ) - panel_length_constraint.addConstraint("PanelLength", dvIndex=0) - else: - panel_length_constraint = None # Retrieve the assembler from pyTACS fea_assembler object assembler = fea_assembler.assembler @@ -1392,7 +1355,10 @@ def create_from_bdf( f5 = fea_assembler.outputViewer # Create the output generator - gen_output = TacsOutputGenerator(prefix, f5=f5) + if prefix is None: + gen_output = None + else: + gen_output = TacsOutputGenerator(prefix, f5=f5) # get struct ids for coordinate derivatives and .sens file struct_id = None @@ -1459,9 +1425,9 @@ def create_from_bdf( override_rotx=override_rotx, Fvec=Fvec, debug=debug, - panel_length_constraint=panel_length_constraint, relaxation_scheme=relaxation_scheme, struct_loads_file=struct_loads_file, + tacs_panel_dimensions=tacs_panel_dimensions, ) @@ -1486,3 +1452,164 @@ def __call__(self): def _deallocate(self): """free up memory before delete""" self.f5.__dealloc__() + + +class TacsPanelDimensions: + def __init__(self, comm, panel_length_dv_index: int, panel_width_dv_index: int): + self.comm = comm + self.panel_length_dv_index = panel_length_dv_index + self.panel_width_dv_index = panel_width_dv_index + + self.panel_length_constr = None + self.panel_width_constr = None + + def compute_panel_length( + self, assembler, struct_vars, constr_base_name, var_base_name + ): + if self.panel_length_constr is not None: + length_funcs = None + if assembler is not None: + + # get the panel length from the TACS constraint object + length_funcs = {} + # clear up to date otherwise it might copy the old value + self.panel_length_constr.externalClearUpToDate() + self.panel_length_constr.evalConstraints(length_funcs) + + # assume rank 0 is a TACS proc (this is true as TACS uses rank 0 as root) + length_funcs = self.comm.bcast(length_funcs, root=0) + + # update the panel length and width dimensions into the F2F variables + # these will later be set into the TACS constitutive objects in the self.set_variables() call + length_comp_ct = 0 + first_key = [_ for _ in length_funcs][0] + # constraint values return actual - current (so use this to solve for actual panel length) + for var in struct_vars: + if var_base_name in var.name: + var.value += length_funcs[first_key][length_comp_ct] + length_comp_ct += 1 + return + + def compute_panel_width( + self, assembler, struct_vars, constr_base_name, var_base_name + ): + if self.panel_width_constr is not None: + width_funcs = None + if assembler is not None: + + # get the panel width from the TACS constraint object + width_funcs = {} + self.panel_width_constr.externalClearUpToDate() + self.panel_width_constr.evalConstraints(width_funcs) + + # assume rank 0 is a TACS proc (this is true as TACS uses rank 0 as root) + width_funcs = self.comm.bcast(width_funcs, root=0) + + # update the panel length and width dimensions into the F2F variables + # these will later be set into the TACS constitutive objects in the self.set_variables() call + width_comp_ct = 0 + first_key = [_ for _ in width_funcs][0] + # constraint values return actual - current (so use this to solve for actual panel length) + for var in struct_vars: + if var_base_name in var.name: + var.value += width_funcs[first_key][width_comp_ct] + width_comp_ct += 1 + return + + def setDesignVars(self, xvec): + if self.panel_length_constr is not None: + self.panel_length_constr.setDesignVars(xvec) + if self.panel_width_constr is not None: + self.panel_width_constr.setDesignVars(xvec) + + def _compute_panel_dimension_xpt_sens(self, scenario): + """ + compute panel length and width coordinate derivatives + for stiffened panel constitutive objects and write into the f2f variables + """ + + # TBD - haven't written this yet + # need to multiply the panel length, width DV values by the xpt sens of the constraints + # then add to the struct xpt sens for each function in F2F + func_grad = self.scenario_data[scenario].func_grad + + for ifunc, func in enumerate(scenario.functions): + for i, var in enumerate(self.struct_variables): + pass + + # end of prototype + + grads_dict = None + if self.assembler is not None: + funcSens = {} + if self.comm.rank == 0: + grads_dict = {} + ifunc = 0 + self.panel_length_constraint.evalConstraintsSens(funcSens) + for func in self.model.composite_functions: + if self.PANEL_LENGTH_CONSTR in func.name and self.comm.rank == 0: + + grads_dict[func.full_name] = {} + + # assume name of form f"{self.PANEL_LENGTH_CONSTR}-fnum" + for ivar, var in enumerate(self.struct_variables): + func.derivatives[var] = funcSens[self.panel_length_name][ + "struct" + ].toarray()[ifunc, ivar] + grads_dict[func.full_name][var.full_name] = func.derivatives[ + var + ] + + ifunc += 1 + + # broadcast the funcs dict to other processors + grads_dict = self.comm.bcast(grads_dict, root=0) + + for func in self.model.composite_functions: + if func.full_name in list(grads_dict.keys()): + for ivar, var in enumerate(self.struct_variables): + func.derivatives[var] = grads_dict[func.full_name][var.full_name] + + # # compute the panel length values + # # ----------------------------------------- + # if self.panel_length_dv_index: + # length_funcs = None + # if self.assembler is not None: + + # # get the panel length from the TACS constraint object + # length_funcs = {} + # self.panel_length_constr.evalConstraints(length_funcs) + + # # assume rank 0 is a TACS proc (this is true as TACS uses rank 0 as root) + # length_funcs = self.comm.bcast(length_funcs, root=0) + + # # update the panel length and width dimensions into the F2F variables + # # these will later be set into the TACS constitutive objects in the self.set_variables() call + # length_comp_ct = 0 + # for var in self.struct_variables: + # if self.LENGTH_VAR in var.name: + # var.value = length_funcs[self.LENGTH_CONSTR][length_comp_ct] + # length_comp_ct += 1 + + # # compute the panel width values + # # ----------------------------------------- + # if self.panel_length_dv_index: + # width_funcs = None + # if self.assembler is not None: + + # # get the panel width from the TACS constraint object + # width_funcs = {} + # self.panel_width_constr.evalConstraints(width_funcs) + + # # assume rank 0 is a TACS proc (this is true as TACS uses rank 0 as root) + # width_funcs = self.comm.bcast(width_funcs, root=0) + + # # update the panel length and width dimensions into the F2F variables + # # these will later be set into the TACS constitutive objects in the self.set_variables() call + # width_comp_ct = 0 + # for var in self.struct_variables: + # if self.WIDTH_VAR in var.name: + # var.value = width_funcs[self.WIDTH_CONSTR][width_comp_ct] + # width_comp_ct += 1 + + # return diff --git a/funtofem/model/composite_function.py b/funtofem/model/composite_function.py index 54acfaaa..64b87ffc 100644 --- a/funtofem/model/composite_function.py +++ b/funtofem/model/composite_function.py @@ -68,6 +68,14 @@ def __init__( self.derivatives = {} self.df_dgi = None + @property + def optim_derivatives(self): + optim_derivs = { + var: self.derivatives[var] for var in self.derivatives if not (var.state) + } + # print(f"optim derivatives = {optim_derivs}") + return optim_derivs + @classmethod def external(cls, name, optim=False, plot_name=None): return cls( @@ -94,13 +102,14 @@ def vars_only(self) -> bool: @property def sparse_gradient(self): """used for adjacency constraints, vars only functions""" - np_array = np.array([self.derivatives[var] for var in self.derivatives]) + self.evaluate_gradient() + np_array = np.array([self.derivatives[var] for var in self.optim_derivatives]) # return csr_matrix(np_array, shape=(1,np_array.shape[0])) nvars = np_array.shape[0] cols = np.array( [ ivar - for ivar, var in enumerate(self.derivatives) + for ivar, var in enumerate(self.optim_derivatives) if self.derivatives[var] != 0.0 ] ) @@ -108,7 +117,7 @@ def sparse_gradient(self): vals = np.array( [ self.derivatives[var] - for ivar, var in enumerate(self.derivatives) + for ivar, var in enumerate(self.optim_derivatives) if self.derivatives[var] != 0.0 ] ) @@ -147,6 +156,11 @@ def plot_name(self) -> str: else: return self.full_name + def setup_sparse_gradient(self, f2f_model): + """setup the sparse gradient for adjacency functions""" + self.setup_derivative_dict(f2f_model.get_variables(optim=True)) + return self + def setup_derivative_dict(self, variables): for var in variables: self.derivatives[var] = 0.0 diff --git a/funtofem/model/funtofem_model.py b/funtofem/model/funtofem_model.py index 32cda09a..7e1a95fe 100644 --- a/funtofem/model/funtofem_model.py +++ b/funtofem/model/funtofem_model.py @@ -25,6 +25,7 @@ import numpy as np, os, importlib from .variable import Variable from .function import CompositeFunction +import sys # optional tacs import for caps2tacs tacs_loader = importlib.util.find_spec("tacs") @@ -255,7 +256,7 @@ def _send_flow_variables(self, base): self.flow.set_variables(active_shape_vars, active_aero_vars) return - def get_variables(self, names=None, all=False): + def get_variables(self, names=None, all=False, optim=False): """ Get all the coupled and uncoupled variable objects for the entire model. Coupled variables only appear once. @@ -272,6 +273,9 @@ def get_variables(self, names=None, all=False): var_list: list of variable objects """ + if optim: # return all active and non-state variables + return [var for var in self.get_variables() if not (var.state)] + if names is None: dv = [] for scenario in self.scenarios: @@ -343,7 +347,7 @@ def count_functions(self): return len(self.get_functions()) - def get_functions(self, optim=False, all=False): + def get_functions(self, optim=False, all=False, include_vars_only=True): """ Get all the functions in the model @@ -353,6 +357,9 @@ def get_functions(self, optim=False, all=False): get functions for optimization when True otherwise just analysis functions within drivers all: bool get all functions analysis or composite for unittests + include_vars_only: bool + whether to include composite functions that have variables only or not + changed default to True so it doesn't mess up sparse gradient functionality for now.. Returns ------- @@ -369,8 +376,12 @@ def get_functions(self, optim=False, all=False): functions.extend(scenario.functions) if optim or all: - # for optimization also include composite functions - functions += self.composite_functions + composite_functions = self.composite_functions + if not include_vars_only: + composite_functions = [ + cfunc for cfunc in composite_functions if not (cfunc.vars_only) + ] + functions += composite_functions # filter out only functions with optim True flag, can be set with func.optimize() functions = [func for func in functions if func.optim or all] @@ -408,6 +419,39 @@ def get_function_gradients(self, optim=False, all=False): return gradients + def clear_vars_only_composite_functions(self): + """ + clears all vars only composite functions to save memory on large wing optimizations + after registering them to the optimizer [their sparse gradients depend only on design variables and are not needed after this point.] + In a large HSCT wing optimization case, there were 13k adjacency constraints which took up a considerable amount of space on memory.. + """ + import gc + import sys + + # get vars only composite functions first to delete later + vars_only_cfuncs = [ + cfunc for cfunc in self.composite_functions if cfunc.vars_only + ] + if len(vars_only_cfuncs) == 0: + return + nvars_only_cfuncs = [ + cfunc for cfunc in self.composite_functions if not (cfunc.vars_only) + ] + self.composite_functions = nvars_only_cfuncs + first_vars_only_cfunc = vars_only_cfuncs[0] + ref_count1 = sys.getrefcount(first_vars_only_cfunc) + del vars_only_cfuncs + gc.collect() # garbage collection.. + # # check reference count on one of the vars only cfunc + # # will only delete if no references available + ref_count2 = sys.getrefcount(first_vars_only_cfunc) + + print( + f"clear_vars_only_cfuncs: ref count1 {ref_count1} ref count2 {ref_count2}" + ) + print(f"need to have 2 for the object to be truly deleted") + return + def evaluate_composite_functions(self, compute_grad=True): """ compute the values and gradients of composite functions @@ -963,6 +1007,85 @@ def write_design_variables_file(self, comm, filename, root=0): return + def print_memory_size(self, comm, root: int, starting_message=""): + """print the memory of the funtofem model and its various constituents""" + + def print_list(m_list, name): + number = len(m_list) + size = sys.getsizeof(m_list) + sum([sys.getsizeof(item) for item in m_list]) + print(f"\t{number} {name} have size {size}") + return size + + if comm.rank == root: + print(f"\nFuntofem model size {starting_message}") + + # size of the whole funtofem model + model_size = sys.getsizeof(self) + print(f"\tmodel_size = {model_size}") + + # size of scenarios list + scenario_size = print_list( + m_list=self.scenarios, + name="scenarios", + ) + + # size of bodies list + body_size = print_list( + m_list=self.bodies, + name="bodies", + ) + + # size of regular functions + functions = [ + func for scenario in self.scenarios for func in scenario.functions + ] + function_size = print_list( + m_list=functions, + name="functions", + ) + + # size of variables + variables = [ + var for scenario in self.scenarios for var in scenario.variables + ] + [var for body in self.bodies for var in body.variables] + variable_size = print_list( + m_list=variables, + name="variables", + ) + + # size of composite funtions, vars only + vars_only_composite_functions = [ + cfunc for cfunc in self.composite_functions if cfunc.vars_only + ] + vars_only_cfunc_size = print_list( + m_list=vars_only_composite_functions, + name="vars only composite functions", + ) + + # size of composite funtions, not vars only + not_vars_only_composite_functions = [ + cfunc for cfunc in self.composite_functions if not (cfunc.vars_only) + ] + nvars_only_cfunc_size = print_list( + m_list=not_vars_only_composite_functions, + name="not vars only composite functions", + ) + + total_size = ( + model_size + + scenario_size + + body_size + + function_size + + variable_size + + vars_only_cfunc_size + + nvars_only_cfunc_size + ) + print(f"\ttotal size = {total_size}") + + print("\n", flush=True) + + return + def read_functions_file(self, comm, filename, root=0, **kwargs): """ Read the functions variables file funtofem.out diff --git a/funtofem/model/variable.py b/funtofem/model/variable.py index 0bb721ca..5da686e0 100644 --- a/funtofem/model/variable.py +++ b/funtofem/model/variable.py @@ -50,6 +50,7 @@ def __init__( active=True, coupled=False, id=0, + state=False, ): """ Variable object for FUNtoFEM @@ -70,6 +71,8 @@ def __init__( whether or not the design variable is coupled id: int id number of the design variable + state: bool + whether this variable is used as a state variable (for TACS only right now) Examples -------- @@ -82,6 +85,7 @@ def __init__( self.upper = upper self.active = active self.coupled = coupled + self.state = state self.id = id self.scale = scale self.analysis_type = analysis_type @@ -96,6 +100,7 @@ def set_bounds( scale=None, active=None, coupled=None, + state=None, ): """ Update the one or more of the attributes of the design variable @@ -114,6 +119,8 @@ def set_bounds( whether or not the design variable is active coupled: bool whether or not the design variable is coupled + state: bool + whether this variable is used as a state variable (for TACS only right now) """ if value is not None: @@ -128,6 +135,8 @@ def set_bounds( self.active = active if coupled is not None: self.coupled = coupled + if state is not None: + self.state = state # return the object for method cascading return self diff --git a/funtofem/optimization/optimization_manager.py b/funtofem/optimization/optimization_manager.py index 579eb05a..fe8f6512 100644 --- a/funtofem/optimization/optimization_manager.py +++ b/funtofem/optimization/optimization_manager.py @@ -46,8 +46,9 @@ def __init__( design_out_file=None, hot_start_file=None, debug: bool = False, - sparse: bool = False, + sparse: bool = True, write_checkpoints=False, + plot_hist=True, ): """ Constructs the optimization manager class using a funtofem model and driver @@ -77,6 +78,7 @@ def __init__( self.debug = debug self.sparse = sparse self.write_checkpoints = write_checkpoints + self.plot_hist = plot_hist # optimization meta data self._iteration = 0 @@ -133,7 +135,7 @@ def __init__( # add all variables (for off-scenario variables to derivatives dict for each function) to analysis functions for func in self.model.get_functions(): - for var in self.model.get_variables(): + for var in self.model.get_variables(optim=True): func.derivatives[var] = 0.0 # initialize funcs, sens in case of failure on first design iteration of hot start @@ -143,7 +145,7 @@ def __init__( self._sens = {} for func in self.model.get_functions(optim=True): self._sens[func.full_name] = {} - for var in self.model.get_variables(): + for var in self.model.get_variables(optim=True): if self.sparse: self._sens[func.full_name][self.SPARSE_VARS_GROUP] = None else: @@ -176,7 +178,7 @@ def _gatekeeper(self, x_dict): if self.sparse: # all vars in same group regular_dict = { var.name: float(x_dict[self.SPARSE_VARS_GROUP][ivar]) - for ivar, var in enumerate(self.model.get_variables()) + for ivar, var in enumerate(self.model.get_variables(optim=True)) } else: regular_dict = {key: float(x_dict[key]) for key in x_dict} @@ -240,13 +242,14 @@ def _gatekeeper(self, x_dict): for func in self.model.get_functions(optim=True) if func._plot } - forward_res = self.driver.solvers.flow.get_forward_residual() - forward_steps = self.driver.solvers.flow._last_forward_step - adjoint_res = self.driver.solvers.flow.get_adjoint_residual() - adjoint_steps = self.driver.solvers.flow._last_adjoint_step - self._design_hdl.write( - f"Forward resid {forward_res:2.5e} in {forward_steps} steps, Adjoint resid {adjoint_res:2.5e} in {adjoint_steps} steps and {_total_time:.4f} {time_units}\n" - ) + if self.driver.solvers.flow is not None: + forward_res = self.driver.solvers.flow.get_forward_residual() + forward_steps = self.driver.solvers.flow._last_forward_step + adjoint_res = self.driver.solvers.flow.get_adjoint_residual() + adjoint_steps = self.driver.solvers.flow._last_adjoint_step + self._design_hdl.write( + f"Forward resid {forward_res:2.5e} in {forward_steps} steps, Adjoint resid {adjoint_res:2.5e} in {adjoint_steps} steps and {_total_time:.4f} {time_units}\n" + ) self._design_hdl.write(f"Functions = {plot_funcs}\n") self._design_hdl.flush() @@ -284,7 +287,8 @@ def _gatekeeper(self, x_dict): if not func._plot: continue self._func_history[func.plot_name] += [func.value.real] - self._plot_history() + if self.plot_hist: + self._plot_history() return fail def _run_complete_analysis(self): @@ -294,12 +298,12 @@ def _run_complete_analysis(self): # update the model design variables if self.sparse: - variables = self.model.get_variables() + variables = self.model.get_variables(optim=True) for ivar, val in enumerate(self._x_dict[self.SPARSE_VARS_GROUP]): var = variables[ivar] var.value = float(val) else: - for var in self.model.get_variables(): + for var in self.model.get_variables(optim=True): for var_key in self._x_dict: if var.full_name == var_key: # assumes here that only pyoptsparse single variables (no var groups are made) @@ -317,6 +321,8 @@ def _run_complete_analysis(self): self._sens = {} # get only functions with optim=True, set with func.optimize() method (can method cascade it) for func in self.model.get_functions(optim=True): + if func.vars_only and self.sparse: + continue # for linear gradients just shown up front self._funcs[func.full_name] = func.value.real self._sens[func.full_name] = {} # if self.sparse and isinstance(func, CompositeFunction) and func.vars_only: @@ -325,22 +331,42 @@ def _run_complete_analysis(self): # ) # linear constraints define their jacobians up front if self.sparse: self._sens[func.full_name][self.SPARSE_VARS_GROUP] = np.array( - [func.derivatives[var].real for var in func.derivatives] + [ + func.derivatives[var].real + for var in self.model.get_variables(optim=True) + ] ) else: - for var in self.model.get_variables(): + for var in self.model.get_variables(optim=True): self._sens[func.full_name][var.full_name] = ( func.get_gradient_component(var).real ) return + def register_sparse_constraint(self, opt_problem, comp_func): + """ + register adjacency composite functions as sparse constraints + need composite function to be vars_only type + """ + assert comp_func.vars_only # the composite function + variables = [var for var in comp_func.derivatives] + opt_problem.addCon( + comp_func.full_name, + lower=comp_func.lower, + upper=comp_func.upper, + scale=comp_func.scale, + linear=True, + wrt=[self.SPARSE_VARS_GROUP], + jac={self.SPARSE_VARS_GROUP: comp_func.sparse_gradient}, + ) + return + def register_to_problem(self, opt_problem): """ add funtofem model variables and functions to a pyoptsparse optimization problem """ if self.sparse: - variables = self.model.get_variables() - values = np.array([var.value for var in variables]) + variables = self.model.get_variables(optim=True) opt_problem.addVarGroup( self.SPARSE_VARS_GROUP, len(variables), @@ -351,7 +377,7 @@ def register_to_problem(self, opt_problem): scale=np.array([var.scale for var in variables]), ) else: - for var in self.model.get_variables(): + for var in self.model.get_variables(optim=True): opt_problem.addVar( var.full_name, lower=var.lower, @@ -360,7 +386,7 @@ def register_to_problem(self, opt_problem): scale=var.scale, ) - for func in self.model.get_functions(optim=True): + for func in self.model.get_functions(optim=True, include_vars_only=True): if func._objective: opt_problem.addObj(func.full_name, scale=func.scale) else: @@ -375,6 +401,8 @@ def register_to_problem(self, opt_problem): wrt=[self.SPARSE_VARS_GROUP], jac={self.SPARSE_VARS_GROUP: func.sparse_gradient}, ) + # print(f"func sparse_gradient = {func.sparse_gradient}") + # exit() else: opt_problem.addCon( func.full_name, @@ -383,6 +411,10 @@ def register_to_problem(self, opt_problem): scale=func.scale, ) + # now remove all vars_only composite functions from the model + # if self.sparse: + # self.model.clear_vars_only_composite_functions() + return def eval_functions(self, x_dict):