diff --git a/funtofem/mphys/mphys_meld.py b/funtofem/mphys/mphys_meld.py index b5b9e375..c58bc4f8 100644 --- a/funtofem/mphys/mphys_meld.py +++ b/funtofem/mphys/mphys_meld.py @@ -26,6 +26,11 @@ def initialize(self): self.options.declare("aero_nnodes") self.options.declare("check_partials") self.options.declare("bodies", recordable=False) + self.options.declare( + "use_reference_coordinates", + types=bool, + desc="Use separate aero and struct reference coordinates for transfer scheme initialization (same variable name with '_ref' appended)", + ) self.struct_ndof = None self.struct_nnodes = None @@ -63,6 +68,21 @@ def setup(self): desc="structural node displacements", tags=["mphys_coupling"], ) + if self.options["use_reference_coordinates"]: + self.add_input( + X_STRUCT0 + "_ref", + shape_by_conn=True, + distributed=True, + desc="baseline structural node coordinates with which to initialize MELD", + tags=["mphys_coordinates"], + ) + self.add_input( + X_AERO0 + "_ref", + shape_by_conn=True, + distributed=True, + desc="baseline aero surface node coordinates with which to initialize MELD", + tags=["mphys_coordinates"], + ) # outputs self.add_output( @@ -77,8 +97,31 @@ def setup(self): # partials # self.declare_partials(U_AERO,[X_STRUCT0,X_AERO0,U_STRUCT]) + def _initialize_xfer(self, inputs, body): + if self.options["use_reference_coordinates"]: + x_s0 = np.array( + inputs[X_STRUCT0 + "_ref"][body.struct_coord_indices], + dtype=TransferScheme.dtype, + ) + x_a0 = np.array( + inputs[X_AERO0 + "_ref"][body.aero_coord_indices], + dtype=TransferScheme.dtype, + ) + else: + x_s0 = np.array( + inputs[X_STRUCT0][body.struct_coord_indices], + dtype=TransferScheme.dtype, + ) + x_a0 = np.array( + inputs[X_AERO0][body.aero_coord_indices], dtype=TransferScheme.dtype + ) + body.meld.setStructNodes(x_s0) + body.meld.setAeroNodes(x_a0) + body.meld.initialize() + def compute(self, inputs, outputs): for body in self.bodies: + x_s0 = np.array( inputs[X_STRUCT0][body.struct_coord_indices], dtype=TransferScheme.dtype, @@ -100,7 +143,7 @@ def compute(self, inputs, outputs): body.meld.setAeroNodes(x_a0) if not body.initialized_meld: - body.meld.initialize() + self._initialize_xfer(inputs, body) body.initialized_meld = True body.meld.transferDisps(u_s, u_a) @@ -556,6 +599,7 @@ def __init__( check_partials=False, linearized=False, body_tags=None, + use_reference_coordinates=False, ): self.aero_builder = aero_builder self.struct_builder = struct_builder @@ -565,6 +609,7 @@ def __init__( self.under_check_partials = check_partials self.linearized = linearized self.body_tags = body_tags if body_tags is not None else [] + self.use_reference_coordinates = use_reference_coordinates if len(self.body_tags) > 0: # make into lists, potentially for different bodies if not hasattr(self.n, "__len__"): @@ -641,6 +686,7 @@ def get_coupling_group_subsystem(self, scenario_name=None): aero_nnodes=self.nnodes_aero, check_partials=self.under_check_partials, bodies=self.bodies, + use_reference_coordinates=self.use_reference_coordinates, ) load_xfer = MeldLoadXfer( diff --git a/funtofem/mphys/mphys_meld_lfd.py b/funtofem/mphys/mphys_meld_lfd.py index db6a03cd..fec50df6 100644 --- a/funtofem/mphys/mphys_meld_lfd.py +++ b/funtofem/mphys/mphys_meld_lfd.py @@ -4,10 +4,16 @@ import numpy as np import openmdao.api as om -from .mphys_meld import MeldBuilder +from mphys import MPhysVariables from funtofem import TransferScheme +from .mphys_meld import MeldBuilder + +# Set MPhys variable names +X_STRUCT0 = MPhysVariables.Structures.COORDINATES +X_AERO0 = MPhysVariables.Aerodynamics.Surface.COORDINATES_INITIAL + class ModeTransfer(om.ExplicitComponent): def initialize(self): @@ -16,14 +22,29 @@ def initialize(self): self.options.declare("ndof_struct") self.options.declare("nnodes_aero") self.options.declare("meld") + self.options.declare( + "use_reference_coordinates", + types=bool, + desc="Use separate aero and struct reference coordinates for transfer scheme initialization (same variable name with '_ref' appended)", + ) + self.options.declare( + "ase_mode", + types=bool, + desc="Add _motion onto the end of mode_shapes_aero, so that they can later be combined with control surface modes", + ) self.first_pass = True def setup(self): # self.set_check_partial_options(wrt='*',method='cs',directional=True) + if self.options["ase_mode"]: + self.mode_shapes_aero_name = "mode_shapes_aero_motion" + else: + self.mode_shapes_aero_name = "mode_shapes_aero" + self.add_input( - "x_struct0", + X_STRUCT0, shape_by_conn=True, distributed=True, tags=["mphys_coordinates"], @@ -35,8 +56,23 @@ def setup(self): tags="mphys_coupling", ) self.add_input( - "x_aero0", shape_by_conn=True, distributed=True, tags=["mphys_coordinates"] + X_AERO0, shape_by_conn=True, distributed=True, tags=["mphys_coordinates"] ) + if self.options["use_reference_coordinates"]: + self.add_input( + X_STRUCT0 + "_ref", + shape_by_conn=True, + distributed=True, + desc="baseline structural node coordinates with which to initialize MELD", + tags=["mphys_coordinates"], + ) + self.add_input( + X_AERO0 + "_ref", + shape_by_conn=True, + distributed=True, + desc="baseline aero surface node coordinates with which to initialize MELD", + tags=["mphys_coordinates"], + ) nmodes = self.options["nmodes"] self.nnodes_aero = self.options["nnodes_aero"] @@ -45,17 +81,28 @@ def setup(self): aero_mode_size = (self.nnodes_aero * 3, nmodes) self.add_output( - "mode_shapes_aero", + self.mode_shapes_aero_name, shape=aero_mode_size, distributed=True, tags=["mphys_coupling"], ) + def _initialize_xfer(self, inputs, meld): + if self.options["use_reference_coordinates"]: + aero_X = np.array(inputs[X_AERO0 + "_ref"], dtype=TransferScheme.dtype) + struct_X = np.array(inputs[X_STRUCT0 + "_ref"], dtype=TransferScheme.dtype) + else: + aero_X = np.array(inputs[X_AERO0], dtype=TransferScheme.dtype) + struct_X = np.array(inputs[X_STRUCT0], dtype=TransferScheme.dtype) + meld.setStructNodes(struct_X) + meld.setAeroNodes(aero_X) + meld.initialize() + def compute(self, inputs, outputs): meld = self.options["meld"] nmodes = self.options["nmodes"] - aero_X = np.array(inputs["x_aero0"], dtype=TransferScheme.dtype) - struct_X = np.array(inputs["x_struct0"], dtype=TransferScheme.dtype) + aero_X = np.array(inputs[X_AERO0], dtype=TransferScheme.dtype) + struct_X = np.array(inputs[X_STRUCT0], dtype=TransferScheme.dtype) aero_modes = np.zeros((aero_X.size, nmodes), dtype=TransferScheme.dtype) struct_modes = inputs["mode_shapes_struct"].reshape((-1, nmodes)) @@ -64,7 +111,7 @@ def compute(self, inputs, outputs): meld.setStructNodes(struct_X) if self.first_pass: - meld.initialize() + self._initialize_xfer(inputs, meld) self.first_pass = False struct_mode = np.zeros(self.nnodes_struct * 3, dtype=TransferScheme.dtype) @@ -76,13 +123,13 @@ def compute(self, inputs, outputs): meld.transferDisps(struct_mode, aero_mode) aero_modes[:, mode] = aero_mode - outputs["mode_shapes_aero"] = aero_modes.copy() + outputs[self.mode_shapes_aero_name] = aero_modes.copy() def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): meld = self.options["meld"] nmodes = self.options["nmodes"] - x_s0 = np.array(inputs["x_struct0"], dtype=TransferScheme.dtype) - x_a0 = np.array(inputs["x_aero0"], dtype=TransferScheme.dtype) + x_s0 = np.array(inputs[X_STRUCT0], dtype=TransferScheme.dtype) + x_a0 = np.array(inputs[X_AERO0], dtype=TransferScheme.dtype) meld.setStructNodes(x_s0) meld.setAeroNodes(x_a0) @@ -94,7 +141,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): u_a = np.zeros(self.nnodes_aero * 3, dtype=TransferScheme.dtype) meld.transferDisps(u_s, u_a) if mode == "fwd": - if "mode_shapes_aero" in d_outputs: + if self.mode_shapes_aero_name in d_outputs: if "mode_shapes_struct" in d_inputs: d_in = np.zeros( self.nnodes_struct * 3, dtype=TransferScheme.dtype @@ -107,13 +154,13 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): self.nnodes_aero * 3, dtype=TransferScheme.dtype ) meld.applydDduS(d_in, prod) - d_outputs["mode_shapes_aero"][:, imode] -= np.array( + d_outputs[self.mode_shapes_aero_name][:, imode] -= np.array( prod, dtype=float ) if mode == "rev": - if "mode_shapes_aero" in d_outputs: + if self.mode_shapes_aero_name in d_outputs: du_a = np.array( - d_outputs["mode_shapes_aero"][:, imode], + d_outputs[self.mode_shapes_aero_name][:, imode], dtype=TransferScheme.dtype, ) if "mode_shapes_struct" in d_inputs: @@ -128,19 +175,19 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): ] -= np.array(prod[i::3], dtype=np.float64) # du_a/dx_a0^T * psi = - psi^T * dD/dx_a0 in F2F terminology - if "x_aero0" in d_inputs: + if X_AERO0 in d_inputs: prod = np.zeros( - d_inputs["x_aero0"].size, dtype=TransferScheme.dtype + d_inputs[X_AERO0].size, dtype=TransferScheme.dtype ) meld.applydDdxA0(du_a, prod) - d_inputs["x_aero0"] -= np.array(prod, dtype=float) + d_inputs[X_AERO0] -= np.array(prod, dtype=float) - if "x_struct0" in d_inputs: + if X_STRUCT0 in d_inputs: prod = np.zeros( self.nnodes_struct * 3, dtype=TransferScheme.dtype ) meld.applydDdxS0(du_a, prod) - d_inputs["x_struct0"] -= np.array(prod, dtype=float) + d_inputs[X_STRUCT0] -= np.array(prod, dtype=float) class MeldLfdBuilder(MeldBuilder): @@ -153,15 +200,30 @@ def __init__( n=200, beta=0.5, check_partials=False, + use_reference_coordinates=False, + ase_mode=False, ): self.nmodes = nmodes - super().__init__(aero_builder, struct_builder, isym, n, beta, check_partials) + self.ase_mode = ase_mode + super().__init__( + aero_builder, + struct_builder, + isym, + n, + beta, + check_partials, + False, + None, + use_reference_coordinates, + ) - def get_post_coupling_subsystem(self): + def get_post_coupling_subsystem(self, scenario_name=None): return ModeTransfer( nmodes=self.nmodes, nnodes_struct=self.nnodes_struct, ndof_struct=self.ndof_struct, nnodes_aero=self.nnodes_aero, - meld=self.meld, + meld=self.bodies[0].meld, # TODO: implement multi-body mode transfer + use_reference_coordinates=self.use_reference_coordinates, + ase_mode=self.ase_mode, )