From 5ff3d3649eaa82b6e6e506c3b7b755cbcae83eb7 Mon Sep 17 00:00:00 2001 From: Andrew Thelen Date: Mon, 10 Jul 2023 17:38:50 -0400 Subject: [PATCH 1/5] modified mphys meld builder for multiple bodies --- funtofem/mphys/mphys_meld.py | 468 +++++++++++++++++++++-------------- 1 file changed, 276 insertions(+), 192 deletions(-) diff --git a/funtofem/mphys/mphys_meld.py b/funtofem/mphys/mphys_meld.py index 846cfa70..935f5d70 100644 --- a/funtofem/mphys/mphys_meld.py +++ b/funtofem/mphys/mphys_meld.py @@ -13,14 +13,11 @@ class MeldDispXfer(om.ExplicitComponent): """ def initialize(self): - self.options.declare("xfer_object", recordable=False) self.options.declare("struct_ndof") self.options.declare("struct_nnodes") self.options.declare("aero_nnodes") self.options.declare("check_partials") - - self.meld = None - self.initialized_meld = False + self.options.declare("bodies", recordable=False) self.struct_ndof = None self.struct_nnodes = None @@ -28,12 +25,11 @@ def initialize(self): self.check_partials = False def setup(self): - self.meld = self.options["xfer_object"] - self.struct_ndof = self.options["struct_ndof"] self.struct_nnodes = self.options["struct_nnodes"] self.aero_nnodes = self.options["aero_nnodes"] self.check_partials = self.options["check_partials"] + self.bodies = self.options["bodies"] # self.set_check_partial_options(wrt='*',method='cs',directional=True) @@ -74,24 +70,28 @@ def setup(self): # self.declare_partials('u_aero',['x_struct0','x_aero0','u_struct']) def compute(self, inputs, outputs): - x_s0 = np.array(inputs["x_struct0"], dtype=TransferScheme.dtype) - x_a0 = np.array(inputs["x_aero0"], dtype=TransferScheme.dtype) - u_a = np.array(outputs["u_aero"], dtype=TransferScheme.dtype) - u_s = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - for i in range(3): - u_s[i::3] = inputs["u_struct"][i :: self.struct_ndof] + for body in self.bodies: + + 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) + u_a = np.array(outputs["u_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) - self.meld.setStructNodes(x_s0) - self.meld.setAeroNodes(x_a0) + u_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + for i in range(3): + u_s[i::3] = inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] - if not self.initialized_meld: - self.meld.initialize() - self.initialized_meld = True + body.meld.setStructNodes(x_s0) + body.meld.setAeroNodes(x_a0) - self.meld.transferDisps(u_s, u_a) + if not body.initialized_meld: + body.meld.initialize() + body.initialized_meld = True + + body.meld.transferDisps(u_s, u_a) + + outputs["u_aero"][body.aero_coord_indices] = u_a - outputs["u_aero"] = u_a def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ @@ -101,67 +101,70 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): D = u_a - g(u_s,x_a0,x_s0) So explicit partials below for u_a are negative partials of D """ - if self.check_partials: - x_s0 = np.array(inputs["x_struct0"], dtype=TransferScheme.dtype) - x_a0 = np.array(inputs["x_aero0"], dtype=TransferScheme.dtype) - self.meld.setStructNodes(x_s0) - self.meld.setAeroNodes(x_a0) - u_s = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - for i in range(3): - u_s[i::3] = inputs["u_struct"][i :: self.struct_ndof] - u_a = np.zeros(self.aero_nnodes * 3, dtype=TransferScheme.dtype) - self.meld.transferDisps(u_s, u_a) - - if mode == "fwd": - if "u_aero" in d_outputs: - if "u_struct" in d_inputs: - d_in = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - for i in range(3): - d_in[i::3] = d_inputs["u_struct"][i :: self.struct_ndof] - prod = np.zeros(self.aero_nnodes * 3, dtype=TransferScheme.dtype) - self.meld.applydDduS(d_in, prod) - d_outputs["u_aero"] -= np.array(prod, dtype=float) - - if "x_aero0" in d_inputs: - if self.check_partials: - pass - else: - raise ValueError( - "MELD forward mode requested but not implemented" - ) - - if "x_struct0" in d_inputs: - if self.check_partials: - pass - else: - raise ValueError( - "MELD forward mode requested but not implemented" - ) - if mode == "rev": - if "u_aero" in d_outputs: - du_a = np.array(d_outputs["u_aero"], dtype=TransferScheme.dtype) - if "u_struct" in d_inputs: - # du_a/du_s^T * psi = - dD/du_s^T psi - prod = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - self.meld.applydDduSTrans(du_a, prod) - for i in range(3): - d_inputs["u_struct"][i :: self.struct_ndof] -= np.array( - prod[i::3], dtype=np.float64 + for body in self.bodies: + + if self.check_partials: + 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) + u_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + for i in range(3): + u_s[i::3] = inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + u_a = np.zeros(len(body.aero_coord_indices), dtype=TransferScheme.dtype) + body.meld.transferDisps(u_s, u_a) + + if mode == "fwd": + if "u_aero" in d_outputs: + if "u_struct" in d_inputs: + d_in = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + for i in range(3): + d_in[i::3] = d_inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + prod = np.zeros(len(body.aero_coord_indices), dtype=TransferScheme.dtype) + body.meld.applydDduS(d_in, prod) + d_outputs["u_aero"][body.aero_coord_indices] -= np.array(prod, dtype=float) + + if "x_aero0" in d_inputs: + if self.check_partials: + pass + else: + raise ValueError( + "MELD forward mode requested but not implemented" + ) + + if "x_struct0" in d_inputs: + if self.check_partials: + pass + else: + raise ValueError( + "MELD forward mode requested but not implemented" + ) + + if mode == "rev": + if "u_aero" in d_outputs: + du_a = np.array(d_outputs["u_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + if "u_struct" in d_inputs: + # du_a/du_s^T * psi = - dD/du_s^T psi + prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + body.meld.applydDduSTrans(du_a, prod) + for i in range(3): + d_inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] -= 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: + prod = np.zeros( + d_inputs["x_aero0"][body.aero_coord_indices].size, dtype=TransferScheme.dtype ) + body.meld.applydDdxA0(du_a, prod) + d_inputs["x_aero0"][body.aero_coord_indices] -= np.array(prod, dtype=float) - # du_a/dx_a0^T * psi = - psi^T * dD/dx_a0 in F2F terminology - if "x_aero0" in d_inputs: - prod = np.zeros( - d_inputs["x_aero0"].size, dtype=TransferScheme.dtype - ) - self.meld.applydDdxA0(du_a, prod) - d_inputs["x_aero0"] -= np.array(prod, dtype=float) - - if "x_struct0" in d_inputs: - prod = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - self.meld.applydDdxS0(du_a, prod) - d_inputs["x_struct0"] -= np.array(prod, dtype=float) + if "x_struct0" in d_inputs: + prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + body.meld.applydDdxS0(du_a, prod) + d_inputs["x_struct0"][body.struct_coord_indices] -= np.array(prod, dtype=float) class MeldLoadXfer(om.ExplicitComponent): @@ -170,14 +173,11 @@ class MeldLoadXfer(om.ExplicitComponent): """ def initialize(self): - self.options.declare("xfer_object", recordable=False) self.options.declare("struct_ndof") self.options.declare("struct_nnodes") self.options.declare("aero_nnodes") self.options.declare("check_partials") - - self.meld = None - self.initialized_meld = False + self.options.declare("bodies", recordable=False) self.struct_ndof = None self.struct_nnodes = None @@ -185,13 +185,11 @@ def initialize(self): self.check_partials = False def setup(self): - # get the transfer scheme object - self.meld = self.options["xfer_object"] - self.struct_ndof = self.options["struct_ndof"] self.struct_nnodes = self.options["struct_nnodes"] self.aero_nnodes = self.options["aero_nnodes"] self.check_partials = self.options["check_partials"] + self.bodies = self.options["bodies"] # self.set_check_partial_options(wrt='*',method='cs',directional=True) @@ -241,25 +239,28 @@ def setup(self): # self.declare_partials('f_struct',['x_struct0','x_aero0','u_struct','f_aero']) def compute(self, inputs, outputs): - if self.check_partials: - x_s0 = np.array(inputs["x_struct0"], dtype=TransferScheme.dtype) - x_a0 = np.array(inputs["x_aero0"], dtype=TransferScheme.dtype) - self.meld.setStructNodes(x_s0) - self.meld.setAeroNodes(x_a0) - f_a = np.array(inputs["f_aero"], dtype=TransferScheme.dtype) - f_s = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - - u_s = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - for i in range(3): - u_s[i::3] = inputs["u_struct"][i :: self.struct_ndof] - u_a = np.zeros(inputs["f_aero"].size, dtype=TransferScheme.dtype) - self.meld.transferDisps(u_s, u_a) - - self.meld.transferLoads(f_a, f_s) outputs["f_struct"][:] = 0.0 - for i in range(3): - outputs["f_struct"][i :: self.struct_ndof] = f_s[i::3] + for body in self.bodies: + + if self.check_partials: + 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) + f_a = np.array(inputs["f_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + f_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + + u_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + for i in range(3): + u_s[i::3] = inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + u_a = np.zeros(inputs["f_aero"][body.aero_coord_indices].size, dtype=TransferScheme.dtype) + body.meld.transferDisps(u_s, u_a) + + body.meld.transferLoads(f_a, f_s) + + for i in range(3): + outputs["f_struct"][body.struct_dof_indices[i :: self.struct_ndof]] = f_s[i::3] def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ @@ -269,90 +270,142 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): L = f_s - g(f_a,u_s,x_a0,x_s0) So explicit partials below for f_s are negative partials of L """ - if self.check_partials: - x_s0 = np.array(inputs["x_struct0"], dtype=TransferScheme.dtype) - x_a0 = np.array(inputs["x_aero0"], dtype=TransferScheme.dtype) - self.meld.setStructNodes(x_s0) - self.meld.setAeroNodes(x_a0) - f_a = np.array(inputs["f_aero"], dtype=TransferScheme.dtype) - f_s = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - - u_s = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - for i in range(3): - u_s[i::3] = inputs["u_struct"][i :: self.struct_ndof] - u_a = np.zeros(inputs["f_aero"].size, dtype=TransferScheme.dtype) - self.meld.transferDisps(u_s, u_a) - self.meld.transferLoads(f_a, f_s) - - if mode == "fwd": - if "f_struct" in d_outputs: - if "u_struct" in d_inputs: - d_in = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - for i in range(3): - d_in[i::3] = d_inputs["u_struct"][i :: self.struct_ndof] - prod = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - self.meld.applydLduS(d_in, prod) - for i in range(3): - d_outputs["f_struct"][i :: self.struct_ndof] -= np.array( - prod[i::3], dtype=float - ) - if "f_aero" in d_inputs: - # df_s/df_a psi = - dL/df_a * psi = -dD/du_s^T * psi - prod = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - df_a = np.array(d_inputs["f_aero"], dtype=TransferScheme.dtype) - self.meld.applydDduSTrans(df_a, prod) + for body in self.bodies: + + if self.check_partials: + 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) + f_a = np.array(inputs["f_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + f_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + + u_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + for i in range(3): + u_s[i::3] = inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + u_a = np.zeros(inputs["f_aero"][body.aero_coord_indices].size, dtype=TransferScheme.dtype) + body.meld.transferDisps(u_s, u_a) + body.meld.transferLoads(f_a, f_s) + + if mode == "fwd": + if "f_struct" in d_outputs: + if "u_struct" in d_inputs: + d_in = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + for i in range(3): + d_in[i::3] = d_inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + body.meld.applydLduS(d_in, prod) + for i in range(3): + d_outputs["f_struct"][body.struct_dof_indices[i :: self.struct_ndof]] -= np.array( + prod[i::3], dtype=float + ) + + if "f_aero" in d_inputs: + # df_s/df_a psi = - dL/df_a * psi = -dD/du_s^T * psi + prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + df_a = np.array(d_inputs["f_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + body.meld.applydDduSTrans(df_a, prod) + for i in range(3): + d_outputs["f_struct"][body.struct_dof_indices[i :: self.struct_ndof]] -= np.array( + prod[i::3], dtype=float + ) + + if "x_aero0" in d_inputs: + if self.check_partials: + pass + else: + raise ValueError("forward mode requested but not implemented") + + if "x_struct0" in d_inputs: + if self.check_partials: + pass + else: + raise ValueError("forward mode requested but not implemented") + + if mode == "rev": + if "f_struct" in d_outputs: + d_out = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) for i in range(3): - d_outputs["f_struct"][i :: self.struct_ndof] -= np.array( - prod[i::3], dtype=float - ) - - if "x_aero0" in d_inputs: - if self.check_partials: - pass - else: - raise ValueError("forward mode requested but not implemented") - - if "x_struct0" in d_inputs: - if self.check_partials: - pass - else: - raise ValueError("forward mode requested but not implemented") - - if mode == "rev": - if "f_struct" in d_outputs: - d_out = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - for i in range(3): - d_out[i::3] = d_outputs["f_struct"][i :: self.struct_ndof] - - if "u_struct" in d_inputs: - d_in = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - # df_s/du_s^T * psi = - dL/du_s^T * psi - self.meld.applydLduSTrans(d_out, d_in) - - for i in range(3): - d_inputs["u_struct"][i :: self.struct_ndof] -= np.array( - d_in[i::3], dtype=float - ) - - if "f_aero" in d_inputs: - # df_s/df_a^T psi = - dL/df_a^T * psi = -dD/du_s * psi - prod = np.zeros(self.aero_nnodes * 3, dtype=TransferScheme.dtype) - self.meld.applydDduS(d_out, prod) - d_inputs["f_aero"] -= np.array(prod, dtype=float) - - if "x_aero0" in d_inputs: - # df_s/dx_a0^T * psi = - psi^T * dL/dx_a0 in F2F terminology - prod = np.zeros(self.aero_nnodes * 3, dtype=TransferScheme.dtype) - self.meld.applydLdxA0(d_out, prod) - d_inputs["x_aero0"] -= np.array(prod, dtype=float) - - if "x_struct0" in d_inputs: - # df_s/dx_s0^T * psi = - psi^T * dL/dx_s0 in F2F terminology - prod = np.zeros(self.struct_nnodes * 3, dtype=TransferScheme.dtype) - self.meld.applydLdxS0(d_out, prod) - d_inputs["x_struct0"] -= np.array(prod, dtype=float) + d_out[i::3] = d_outputs["f_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + + if "u_struct" in d_inputs: + d_in = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + # df_s/du_s^T * psi = - dL/du_s^T * psi + body.meld.applydLduSTrans(d_out, d_in) + + for i in range(3): + d_inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] -= np.array( + d_in[i::3], dtype=float + ) + + if "f_aero" in d_inputs: + # df_s/df_a^T psi = - dL/df_a^T * psi = -dD/du_s * psi + prod = np.zeros(len(body.aero_coord_indices), dtype=TransferScheme.dtype) + body.meld.applydDduS(d_out, prod) + d_inputs["f_aero"][body.aero_coord_indices] -= np.array(prod, dtype=float) + + if "x_aero0" in d_inputs: + # df_s/dx_a0^T * psi = - psi^T * dL/dx_a0 in F2F terminology + prod = np.zeros(len(body.aero_coord_indices), dtype=TransferScheme.dtype) + body.meld.applydLdxA0(d_out, prod) + d_inputs["x_aero0"][body.aero_coord_indices] -= np.array(prod, dtype=float) + + if "x_struct0" in d_inputs: + # df_s/dx_s0^T * psi = - psi^T * dL/dx_s0 in F2F terminology + prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + body.meld.applydLdxS0(d_out, prod) + d_inputs["x_struct0"][body.struct_coord_indices] -= np.array(prod, dtype=float) + +class MeldBodyInstance: + def __init__( + self, + comm, + isym=-1, + n=200, + beta=0.5, + aero_gid=None, # list of aero grid IDs to be couple + aero_nnodes=None, + struct_gid=None, # list of struct grid IDs to be coupled + struct_nnodes=None, + struct_ndof=None, + linearized=False + ): + # determine input/output indices for the current body + if struct_gid is not None: # use grid ids from struct_builder's get_tagged_indices + self.struct_coord_indices = [] # for coordinates + self.struct_dof_indices = [] # for dof + for i in range(3): # accounts for xyz of each grid coordinate + self.struct_coord_indices = np.hstack([self.struct_coord_indices, struct_gid*3+i]) + for i in range(struct_ndof): # accounts for each structural dof of each grid coordinate + self.struct_dof_indices = np.hstack([self.struct_dof_indices, struct_gid*struct_ndof+i]) + self.struct_coord_indices = list(np.sort(self.struct_coord_indices).astype('int')) + self.struct_dof_indices = list(np.sort(self.struct_dof_indices).astype('int')) + + else: # use all indices + self.struct_coord_indices = list(np.arange(struct_nnodes*3)) + self.struct_dof_indices = list(np.arange(struct_nnodes*struct_ndof)) + + if aero_gid is not None: # use grid ids from aero_builder's get_tagged_indices + self.aero_coord_indices = [] + for i in range(3): + self.aero_coord_indices = np.hstack([self.aero_coord_indices, aero_gid*3+i]) + self.aero_coord_indices = list(np.sort(self.aero_coord_indices).astype('int')) + + else: # use all indices + self.aero_coord_indices = list(np.arange(aero_nnodes*3)) + + # new MELD instance + if linearized: + self.meld = TransferScheme.pyLinearizedMELD( + comm, comm, 0, comm, 0, isym, n, beta + ) + else: + self.meld = TransferScheme.pyMELD( + comm, comm, 0, comm, 0, isym, n, beta + ) + self.initialized_meld = False class MeldBuilder(Builder): def __init__( @@ -364,6 +417,7 @@ def __init__( beta=0.5, check_partials=False, linearized=False, + body_tags=[] ): self.aero_builder = aero_builder self.struct_builder = struct_builder @@ -372,36 +426,66 @@ def __init__( self.beta = beta self.check_partials = check_partials self.linearized = linearized + self.body_tags = body_tags + + if len(self.body_tags)>0: # make into lists, potentially for different bodies + if not hasattr(self.n,'__len__'): + self.n = [self.n]*len(self.body_tags) + if not hasattr(self.beta,'__len__'): + self.beta = [self.beta]*len(self.body_tags) + if not hasattr(self.linearized,'__len__'): + self.linearized = [self.linearized]*len(self.body_tags) def initialize(self, comm): self.nnodes_aero = self.aero_builder.get_number_of_nodes() self.nnodes_struct = self.struct_builder.get_number_of_nodes() self.ndof_struct = self.struct_builder.get_ndof() - if self.linearized: - self.meld = TransferScheme.pyLinearizedMELD( - comm, comm, 0, comm, 0, self.isym, self.n, self.beta - ) - else: - self.meld = TransferScheme.pyMELD( - comm, comm, 0, comm, 0, self.isym, self.n, self.beta - ) + self.bodies = [] + + if len(self.body_tags)>0: # body tags given + for i in range(len(self.body_tags)): + aero_gid = self.aero_builder.get_tagged_indices(self.body_tags[i]["aero"]) + struct_gid = self.struct_builder.get_tagged_indices(self.body_tags[i]["struct"]) + self.bodies += [MeldBodyInstance( + comm, + isym=self.isym, + n=self.n[i], + beta=self.beta[i], + aero_gid=aero_gid, + struct_gid=struct_gid, + struct_ndof=self.ndof_struct, + linearized=self.linearized[i] + )] + + else: # default: couple all nodes with a single MELD instance + self.bodies = [MeldBodyInstance( + comm, + isym=self.isym, + n=self.n, + beta=self.beta, + aero_nnodes=self.nnodes_aero, + struct_nnodes=self.nnodes_struct, + struct_ndof=self.ndof_struct, + linearized=self.linearized + )] + def get_coupling_group_subsystem(self, scenario_name=None): disp_xfer = MeldDispXfer( - xfer_object=self.meld, struct_ndof=self.ndof_struct, struct_nnodes=self.nnodes_struct, aero_nnodes=self.nnodes_aero, check_partials=self.check_partials, + bodies=self.bodies ) load_xfer = MeldLoadXfer( - xfer_object=self.meld, struct_ndof=self.ndof_struct, struct_nnodes=self.nnodes_struct, aero_nnodes=self.nnodes_aero, check_partials=self.check_partials, + bodies=self.bodies ) return disp_xfer, load_xfer From d5abb7ffe55ba8ca4434511b1163544ddab4c1d9 Mon Sep 17 00:00:00 2001 From: Andrew Thelen Date: Tue, 11 Jul 2023 12:31:38 -0400 Subject: [PATCH 2/5] added docstring to MeldBodyInstance; corrected misspelling --- funtofem/mphys/mphys_meld.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/funtofem/mphys/mphys_meld.py b/funtofem/mphys/mphys_meld.py index 935f5d70..897ad0c7 100644 --- a/funtofem/mphys/mphys_meld.py +++ b/funtofem/mphys/mphys_meld.py @@ -358,15 +358,19 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): d_inputs["x_struct0"][body.struct_coord_indices] -= np.array(prod, dtype=float) class MeldBodyInstance: + """ + Class that helps split OpenMDAO input/output indices for multiple bodies, + with a separate MELD instance for each. + """ def __init__( self, comm, isym=-1, n=200, beta=0.5, - aero_gid=None, # list of aero grid IDs to be couple + aero_gid=None, # list of aero grid IDs to couple aero_nnodes=None, - struct_gid=None, # list of struct grid IDs to be coupled + struct_gid=None, # list of struct grid IDs to couple struct_nnodes=None, struct_ndof=None, linearized=False @@ -470,7 +474,6 @@ def initialize(self, comm): linearized=self.linearized )] - def get_coupling_group_subsystem(self, scenario_name=None): disp_xfer = MeldDispXfer( struct_ndof=self.ndof_struct, From 153e92341ca54344dc19d23a58979a5823f825c7 Mon Sep 17 00:00:00 2001 From: Andrew Thelen Date: Mon, 17 Jul 2023 09:25:50 -0400 Subject: [PATCH 3/5] reformat MELD MPhys builder; added try loops around get_tagged_indices now that default now raises NotImplementedError --- funtofem/mphys/mphys_meld.py | 381 +++++++++++++++++++++++------------ 1 file changed, 256 insertions(+), 125 deletions(-) diff --git a/funtofem/mphys/mphys_meld.py b/funtofem/mphys/mphys_meld.py index 897ad0c7..5c1740c6 100644 --- a/funtofem/mphys/mphys_meld.py +++ b/funtofem/mphys/mphys_meld.py @@ -70,16 +70,23 @@ def setup(self): # self.declare_partials('u_aero',['x_struct0','x_aero0','u_struct']) def compute(self, inputs, outputs): - for body in self.bodies: - - 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) - u_a = np.array(outputs["u_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + 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 + ) + u_a = np.array( + outputs["u_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype + ) u_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) for i in range(3): - u_s[i::3] = inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + u_s[i::3] = inputs["u_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] body.meld.setStructNodes(x_s0) body.meld.setAeroNodes(x_a0) @@ -92,7 +99,6 @@ def compute(self, inputs, outputs): outputs["u_aero"][body.aero_coord_indices] = u_a - def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ The explicit component is defined as: @@ -103,27 +109,42 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ for body in self.bodies: - if self.check_partials: - 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) + 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) u_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) for i in range(3): - u_s[i::3] = inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + u_s[i::3] = inputs["u_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] u_a = np.zeros(len(body.aero_coord_indices), dtype=TransferScheme.dtype) body.meld.transferDisps(u_s, u_a) if mode == "fwd": if "u_aero" in d_outputs: if "u_struct" in d_inputs: - d_in = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + d_in = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) for i in range(3): - d_in[i::3] = d_inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] - prod = np.zeros(len(body.aero_coord_indices), dtype=TransferScheme.dtype) + d_in[i::3] = d_inputs["u_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] + prod = np.zeros( + len(body.aero_coord_indices), dtype=TransferScheme.dtype + ) body.meld.applydDduS(d_in, prod) - d_outputs["u_aero"][body.aero_coord_indices] -= np.array(prod, dtype=float) + d_outputs["u_aero"][body.aero_coord_indices] -= np.array( + prod, dtype=float + ) if "x_aero0" in d_inputs: if self.check_partials: @@ -143,28 +164,40 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): if mode == "rev": if "u_aero" in d_outputs: - du_a = np.array(d_outputs["u_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + du_a = np.array( + d_outputs["u_aero"][body.aero_coord_indices], + dtype=TransferScheme.dtype, + ) if "u_struct" in d_inputs: # du_a/du_s^T * psi = - dD/du_s^T psi - prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + prod = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) body.meld.applydDduSTrans(du_a, prod) for i in range(3): - d_inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] -= np.array( - prod[i::3], dtype=np.float64 - ) + d_inputs["u_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] -= 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: prod = np.zeros( - d_inputs["x_aero0"][body.aero_coord_indices].size, dtype=TransferScheme.dtype + d_inputs["x_aero0"][body.aero_coord_indices].size, + dtype=TransferScheme.dtype, ) body.meld.applydDdxA0(du_a, prod) - d_inputs["x_aero0"][body.aero_coord_indices] -= np.array(prod, dtype=float) + d_inputs["x_aero0"][body.aero_coord_indices] -= np.array( + prod, dtype=float + ) if "x_struct0" in d_inputs: - prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + prod = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) body.meld.applydDdxS0(du_a, prod) - d_inputs["x_struct0"][body.struct_coord_indices] -= np.array(prod, dtype=float) + d_inputs["x_struct0"][body.struct_coord_indices] -= np.array( + prod, dtype=float + ) class MeldLoadXfer(om.ExplicitComponent): @@ -239,28 +272,41 @@ def setup(self): # self.declare_partials('f_struct',['x_struct0','x_aero0','u_struct','f_aero']) def compute(self, inputs, outputs): - outputs["f_struct"][:] = 0.0 for body in self.bodies: - if self.check_partials: - 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) + 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) - f_a = np.array(inputs["f_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + f_a = np.array( + inputs["f_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype + ) f_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) u_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) for i in range(3): - u_s[i::3] = inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] - u_a = np.zeros(inputs["f_aero"][body.aero_coord_indices].size, dtype=TransferScheme.dtype) + u_s[i::3] = inputs["u_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] + u_a = np.zeros( + inputs["f_aero"][body.aero_coord_indices].size, + dtype=TransferScheme.dtype, + ) body.meld.transferDisps(u_s, u_a) body.meld.transferLoads(f_a, f_s) for i in range(3): - outputs["f_struct"][body.struct_dof_indices[i :: self.struct_ndof]] = f_s[i::3] + outputs["f_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] = f_s[i::3] def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ @@ -272,133 +318,195 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ for body in self.bodies: - if self.check_partials: - 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) + 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) - f_a = np.array(inputs["f_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + f_a = np.array( + inputs["f_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype + ) f_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) u_s = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) for i in range(3): - u_s[i::3] = inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] - u_a = np.zeros(inputs["f_aero"][body.aero_coord_indices].size, dtype=TransferScheme.dtype) + u_s[i::3] = inputs["u_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] + u_a = np.zeros( + inputs["f_aero"][body.aero_coord_indices].size, + dtype=TransferScheme.dtype, + ) body.meld.transferDisps(u_s, u_a) body.meld.transferLoads(f_a, f_s) if mode == "fwd": if "f_struct" in d_outputs: if "u_struct" in d_inputs: - d_in = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + d_in = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) for i in range(3): - d_in[i::3] = d_inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] - prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + d_in[i::3] = d_inputs["u_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] + prod = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) body.meld.applydLduS(d_in, prod) for i in range(3): - d_outputs["f_struct"][body.struct_dof_indices[i :: self.struct_ndof]] -= np.array( - prod[i::3], dtype=float - ) + d_outputs["f_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] -= np.array(prod[i::3], dtype=float) if "f_aero" in d_inputs: # df_s/df_a psi = - dL/df_a * psi = -dD/du_s^T * psi - prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) - df_a = np.array(d_inputs["f_aero"][body.aero_coord_indices], dtype=TransferScheme.dtype) + prod = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) + df_a = np.array( + d_inputs["f_aero"][body.aero_coord_indices], + dtype=TransferScheme.dtype, + ) body.meld.applydDduSTrans(df_a, prod) for i in range(3): - d_outputs["f_struct"][body.struct_dof_indices[i :: self.struct_ndof]] -= np.array( - prod[i::3], dtype=float - ) + d_outputs["f_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] -= np.array(prod[i::3], dtype=float) if "x_aero0" in d_inputs: if self.check_partials: pass else: - raise ValueError("forward mode requested but not implemented") + raise ValueError( + "forward mode requested but not implemented" + ) if "x_struct0" in d_inputs: if self.check_partials: pass else: - raise ValueError("forward mode requested but not implemented") + raise ValueError( + "forward mode requested but not implemented" + ) if mode == "rev": if "f_struct" in d_outputs: - d_out = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + d_out = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) for i in range(3): - d_out[i::3] = d_outputs["f_struct"][body.struct_dof_indices[i :: self.struct_ndof]] + d_out[i::3] = d_outputs["f_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] if "u_struct" in d_inputs: - d_in = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + d_in = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) # df_s/du_s^T * psi = - dL/du_s^T * psi body.meld.applydLduSTrans(d_out, d_in) for i in range(3): - d_inputs["u_struct"][body.struct_dof_indices[i :: self.struct_ndof]] -= np.array( - d_in[i::3], dtype=float - ) + d_inputs["u_struct"][ + body.struct_dof_indices[i :: self.struct_ndof] + ] -= np.array(d_in[i::3], dtype=float) if "f_aero" in d_inputs: # df_s/df_a^T psi = - dL/df_a^T * psi = -dD/du_s * psi - prod = np.zeros(len(body.aero_coord_indices), dtype=TransferScheme.dtype) + prod = np.zeros( + len(body.aero_coord_indices), dtype=TransferScheme.dtype + ) body.meld.applydDduS(d_out, prod) - d_inputs["f_aero"][body.aero_coord_indices] -= np.array(prod, dtype=float) + d_inputs["f_aero"][body.aero_coord_indices] -= np.array( + prod, dtype=float + ) if "x_aero0" in d_inputs: # df_s/dx_a0^T * psi = - psi^T * dL/dx_a0 in F2F terminology - prod = np.zeros(len(body.aero_coord_indices), dtype=TransferScheme.dtype) + prod = np.zeros( + len(body.aero_coord_indices), dtype=TransferScheme.dtype + ) body.meld.applydLdxA0(d_out, prod) - d_inputs["x_aero0"][body.aero_coord_indices] -= np.array(prod, dtype=float) + d_inputs["x_aero0"][body.aero_coord_indices] -= np.array( + prod, dtype=float + ) if "x_struct0" in d_inputs: # df_s/dx_s0^T * psi = - psi^T * dL/dx_s0 in F2F terminology - prod = np.zeros(len(body.struct_coord_indices), dtype=TransferScheme.dtype) + prod = np.zeros( + len(body.struct_coord_indices), dtype=TransferScheme.dtype + ) body.meld.applydLdxS0(d_out, prod) - d_inputs["x_struct0"][body.struct_coord_indices] -= np.array(prod, dtype=float) + d_inputs["x_struct0"][body.struct_coord_indices] -= np.array( + prod, dtype=float + ) + class MeldBodyInstance: """ Class that helps split OpenMDAO input/output indices for multiple bodies, - with a separate MELD instance for each. + with a separate MELD instance for each """ + def __init__( self, comm, isym=-1, n=200, beta=0.5, - aero_gid=None, # list of aero grid IDs to couple + aero_gid=None, # list of aero grid IDs to couple aero_nnodes=None, - struct_gid=None, # list of struct grid IDs to couple + struct_gid=None, # list of struct grid IDs to couple struct_nnodes=None, struct_ndof=None, - linearized=False + linearized=False, ): - # determine input/output indices for the current body - if struct_gid is not None: # use grid ids from struct_builder's get_tagged_indices - self.struct_coord_indices = [] # for coordinates - self.struct_dof_indices = [] # for dof - for i in range(3): # accounts for xyz of each grid coordinate - self.struct_coord_indices = np.hstack([self.struct_coord_indices, struct_gid*3+i]) - for i in range(struct_ndof): # accounts for each structural dof of each grid coordinate - self.struct_dof_indices = np.hstack([self.struct_dof_indices, struct_gid*struct_ndof+i]) - self.struct_coord_indices = list(np.sort(self.struct_coord_indices).astype('int')) - self.struct_dof_indices = list(np.sort(self.struct_dof_indices).astype('int')) - - else: # use all indices - self.struct_coord_indices = list(np.arange(struct_nnodes*3)) - self.struct_dof_indices = list(np.arange(struct_nnodes*struct_ndof)) - - if aero_gid is not None: # use grid ids from aero_builder's get_tagged_indices + if ( + struct_gid is not None + ): # use grid ids from struct_builder's get_tagged_indices + self.struct_coord_indices = [] # for coordinates + self.struct_dof_indices = [] # for dof + for i in range(3): # accounts for xyz of each grid coordinate + self.struct_coord_indices = np.hstack( + [self.struct_coord_indices, struct_gid * 3 + i] + ) + for i in range( + struct_ndof + ): # accounts for each structural dof of each grid coordinate + self.struct_dof_indices = np.hstack( + [self.struct_dof_indices, struct_gid * struct_ndof + i] + ) + self.struct_coord_indices = list( + np.sort(self.struct_coord_indices).astype("int") + ) + self.struct_dof_indices = list( + np.sort(self.struct_dof_indices).astype("int") + ) + + else: # use all indices + self.struct_coord_indices = list(np.arange(struct_nnodes * 3)) + self.struct_dof_indices = list(np.arange(struct_nnodes * struct_ndof)) + + if aero_gid is not None: # use grid ids from aero_builder's get_tagged_indices self.aero_coord_indices = [] for i in range(3): - self.aero_coord_indices = np.hstack([self.aero_coord_indices, aero_gid*3+i]) - self.aero_coord_indices = list(np.sort(self.aero_coord_indices).astype('int')) + self.aero_coord_indices = np.hstack( + [self.aero_coord_indices, aero_gid * 3 + i] + ) + self.aero_coord_indices = list( + np.sort(self.aero_coord_indices).astype("int") + ) - else: # use all indices - self.aero_coord_indices = list(np.arange(aero_nnodes*3)) + else: # use all indices + self.aero_coord_indices = list(np.arange(aero_nnodes * 3)) # new MELD instance if linearized: @@ -406,11 +514,10 @@ def __init__( comm, comm, 0, comm, 0, isym, n, beta ) else: - self.meld = TransferScheme.pyMELD( - comm, comm, 0, comm, 0, isym, n, beta - ) + self.meld = TransferScheme.pyMELD(comm, comm, 0, comm, 0, isym, n, beta) self.initialized_meld = False + class MeldBuilder(Builder): def __init__( self, @@ -421,7 +528,7 @@ def __init__( beta=0.5, check_partials=False, linearized=False, - body_tags=[] + body_tags=[], ): self.aero_builder = aero_builder self.struct_builder = struct_builder @@ -432,13 +539,13 @@ def __init__( self.linearized = linearized self.body_tags = body_tags - if len(self.body_tags)>0: # make into lists, potentially for different bodies - if not hasattr(self.n,'__len__'): - self.n = [self.n]*len(self.body_tags) - if not hasattr(self.beta,'__len__'): - self.beta = [self.beta]*len(self.body_tags) - if not hasattr(self.linearized,'__len__'): - self.linearized = [self.linearized]*len(self.body_tags) + if len(self.body_tags) > 0: # make into lists, potentially for different bodies + if not hasattr(self.n, "__len__"): + self.n = [self.n] * len(self.body_tags) + if not hasattr(self.beta, "__len__"): + self.beta = [self.beta] * len(self.body_tags) + if not hasattr(self.linearized, "__len__"): + self.linearized = [self.linearized] * len(self.body_tags) def initialize(self, comm): self.nnodes_aero = self.aero_builder.get_number_of_nodes() @@ -447,32 +554,56 @@ def initialize(self, comm): self.bodies = [] - if len(self.body_tags)>0: # body tags given + if len(self.body_tags) > 0: # body tags given for i in range(len(self.body_tags)): - aero_gid = self.aero_builder.get_tagged_indices(self.body_tags[i]["aero"]) - struct_gid = self.struct_builder.get_tagged_indices(self.body_tags[i]["struct"]) - self.bodies += [MeldBodyInstance( - comm, - isym=self.isym, - n=self.n[i], - beta=self.beta[i], - aero_gid=aero_gid, - struct_gid=struct_gid, - struct_ndof=self.ndof_struct, - linearized=self.linearized[i] - )] - - else: # default: couple all nodes with a single MELD instance - self.bodies = [MeldBodyInstance( - comm, - isym=self.isym, - n=self.n, - beta=self.beta, - aero_nnodes=self.nnodes_aero, - struct_nnodes=self.nnodes_struct, - struct_ndof=self.ndof_struct, - linearized=self.linearized - )] + try: + aero_gid = self.aero_builder.get_tagged_indices( + self.body_tags[i]["aero"] + ) + except: + if comm.rank == 0: + print( + "get_tagged_indices has not been implemented in the aero builder; all nodes will be used" + ) + aero_gid = None + try: + struct_gid = self.struct_builder.get_tagged_indices( + self.body_tags[i]["struct"] + ) + except: + if comm.rank == 0: + print( + "get_tagged_indices has not been implemented in the struct builder; all nodes will be used" + ) + struct_gid = None + self.bodies += [ + MeldBodyInstance( + comm, + isym=self.isym, + n=self.n[i], + beta=self.beta[i], + aero_gid=aero_gid, + aero_nnodes=self.nnodes_aero, + struct_gid=struct_gid, + struct_nnodes=self.nnodes_struct, + struct_ndof=self.ndof_struct, + linearized=self.linearized[i], + ) + ] + + else: # default: couple all nodes with a single MELD instance + self.bodies = [ + MeldBodyInstance( + comm, + isym=self.isym, + n=self.n, + beta=self.beta, + aero_nnodes=self.nnodes_aero, + struct_nnodes=self.nnodes_struct, + struct_ndof=self.ndof_struct, + linearized=self.linearized, + ) + ] def get_coupling_group_subsystem(self, scenario_name=None): disp_xfer = MeldDispXfer( @@ -480,7 +611,7 @@ def get_coupling_group_subsystem(self, scenario_name=None): struct_nnodes=self.nnodes_struct, aero_nnodes=self.nnodes_aero, check_partials=self.check_partials, - bodies=self.bodies + bodies=self.bodies, ) load_xfer = MeldLoadXfer( @@ -488,7 +619,7 @@ def get_coupling_group_subsystem(self, scenario_name=None): struct_nnodes=self.nnodes_struct, aero_nnodes=self.nnodes_aero, check_partials=self.check_partials, - bodies=self.bodies + bodies=self.bodies, ) return disp_xfer, load_xfer From 417a0199cfab96016c846b1c6903aa9bae4ce8d4 Mon Sep 17 00:00:00 2001 From: Andrew Thelen Date: Mon, 17 Jul 2023 15:28:37 -0400 Subject: [PATCH 4/5] add atleast_1d to get_tagged_indices call; parameter info to MeldBodyInstance docstring --- funtofem/mphys/mphys_meld.py | 47 ++++++++++++++++++++++++++++-------- 1 file changed, 37 insertions(+), 10 deletions(-) diff --git a/funtofem/mphys/mphys_meld.py b/funtofem/mphys/mphys_meld.py index 5c1740c6..85b6e792 100644 --- a/funtofem/mphys/mphys_meld.py +++ b/funtofem/mphys/mphys_meld.py @@ -453,6 +453,29 @@ class MeldBodyInstance: """ Class that helps split OpenMDAO input/output indices for multiple bodies, with a separate MELD instance for each + + Parameters + ---------- + comm : :class:`~mpi4py.MPI.Comm` + The communicator object created for this xfer object instance. + isym : int + Whether to search for symmetries in the geometry for transfer + n : int + The number of nearest neighbors included in the transfers + beta : float + The exponential decay factor used to average loads, displacements, etc. + aero_gid : list + List of aerodynamic node IDs on the current body and rank + aero_nnodes : int + Number of aerodynamic nodes on the current body and rank + struct_gid : list + List of structural node IDs on the current body and rank + struct_nnodes : int + Number of structural nodes on the current body and rank + struct_ndof : int + Number of degrees of freedom at each structural node + linearized : bool + Whether to use linearized MELD """ def __init__( @@ -461,9 +484,9 @@ def __init__( isym=-1, n=200, beta=0.5, - aero_gid=None, # list of aero grid IDs to couple + aero_gid=None, aero_nnodes=None, - struct_gid=None, # list of struct grid IDs to couple + struct_gid=None, struct_nnodes=None, struct_ndof=None, linearized=False, @@ -474,13 +497,15 @@ def __init__( ): # use grid ids from struct_builder's get_tagged_indices self.struct_coord_indices = [] # for coordinates self.struct_dof_indices = [] # for dof - for i in range(3): # accounts for xyz of each grid coordinate + for i in range( + 3 + ): # account for xyz of each node (i.e., [x_0, y_0, z_0, ..., x_N, y_N, z_N]) self.struct_coord_indices = np.hstack( [self.struct_coord_indices, struct_gid * 3 + i] ) for i in range( struct_ndof - ): # accounts for each structural dof of each grid coordinate + ): # account for each structural DOF of each node (i.e., [u_0, v_0, w_0, ..., u_N, v_N, w_N]) self.struct_dof_indices = np.hstack( [self.struct_dof_indices, struct_gid * struct_ndof + i] ) @@ -557,20 +582,22 @@ def initialize(self, comm): if len(self.body_tags) > 0: # body tags given for i in range(len(self.body_tags)): try: - aero_gid = self.aero_builder.get_tagged_indices( - self.body_tags[i]["aero"] + aero_gid = np.atleast_1d( + self.aero_builder.get_tagged_indices(self.body_tags[i]["aero"]) ) - except: + except NotImplementedError: if comm.rank == 0: print( "get_tagged_indices has not been implemented in the aero builder; all nodes will be used" ) aero_gid = None try: - struct_gid = self.struct_builder.get_tagged_indices( - self.body_tags[i]["struct"] + struct_gid = np.atleast_1d( + self.struct_builder.get_tagged_indices( + self.body_tags[i]["struct"] + ) ) - except: + except NotImplementedError: if comm.rank == 0: print( "get_tagged_indices has not been implemented in the struct builder; all nodes will be used" From 0f5555c489ee97b4ecd1abb8d8dfaf937acd9073 Mon Sep 17 00:00:00 2001 From: Andrew Thelen Date: Mon, 17 Jul 2023 16:58:07 -0400 Subject: [PATCH 5/5] cleaned up body index calculation; changed *_gid to *_node_ids --- funtofem/mphys/mphys_meld.py | 78 +++++++++++++++++------------------- 1 file changed, 36 insertions(+), 42 deletions(-) diff --git a/funtofem/mphys/mphys_meld.py b/funtofem/mphys/mphys_meld.py index 85b6e792..20199ff1 100644 --- a/funtofem/mphys/mphys_meld.py +++ b/funtofem/mphys/mphys_meld.py @@ -457,21 +457,21 @@ class MeldBodyInstance: Parameters ---------- comm : :class:`~mpi4py.MPI.Comm` - The communicator object created for this xfer object instance. + The communicator object created for this xfer object instance isym : int Whether to search for symmetries in the geometry for transfer n : int The number of nearest neighbors included in the transfers beta : float The exponential decay factor used to average loads, displacements, etc. - aero_gid : list + aero_node_ids : list List of aerodynamic node IDs on the current body and rank aero_nnodes : int - Number of aerodynamic nodes on the current body and rank - struct_gid : list + Total number of aerodynamic nodes on the current rank + struct_node_ids : list List of structural node IDs on the current body and rank struct_nnodes : int - Number of structural nodes on the current body and rank + Total number of structural nodes on the current rank struct_ndof : int Number of degrees of freedom at each structural node linearized : bool @@ -484,51 +484,45 @@ def __init__( isym=-1, n=200, beta=0.5, - aero_gid=None, + aero_node_ids=None, aero_nnodes=None, - struct_gid=None, + struct_node_ids=None, struct_nnodes=None, struct_ndof=None, linearized=False, ): # determine input/output indices for the current body if ( - struct_gid is not None + struct_node_ids is not None ): # use grid ids from struct_builder's get_tagged_indices - self.struct_coord_indices = [] # for coordinates - self.struct_dof_indices = [] # for dof - for i in range( - 3 - ): # account for xyz of each node (i.e., [x_0, y_0, z_0, ..., x_N, y_N, z_N]) - self.struct_coord_indices = np.hstack( - [self.struct_coord_indices, struct_gid * 3 + i] - ) - for i in range( - struct_ndof - ): # account for each structural DOF of each node (i.e., [u_0, v_0, w_0, ..., u_N, v_N, w_N]) - self.struct_dof_indices = np.hstack( - [self.struct_dof_indices, struct_gid * struct_ndof + i] - ) - self.struct_coord_indices = list( - np.sort(self.struct_coord_indices).astype("int") - ) - self.struct_dof_indices = list( - np.sort(self.struct_dof_indices).astype("int") + # account for xyz of each node + self.struct_coord_indices = np.zeros(len(struct_node_ids) * 3, dtype=int) + for ii in range(3): + self.struct_coord_indices[ii::3] = 3 * struct_node_ids + ii + self.struct_coord_indices = list(self.struct_coord_indices) + + # account for structural DOFs of each node + self.struct_dof_indices = np.zeros( + len(struct_node_ids) * struct_ndof, dtype=int ) + for ii in range(struct_ndof): + self.struct_dof_indices[ii::struct_ndof] = ( + struct_node_ids * struct_ndof + ii + ) + self.struct_dof_indices = list(self.struct_dof_indices) else: # use all indices self.struct_coord_indices = list(np.arange(struct_nnodes * 3)) self.struct_dof_indices = list(np.arange(struct_nnodes * struct_ndof)) - if aero_gid is not None: # use grid ids from aero_builder's get_tagged_indices - self.aero_coord_indices = [] - for i in range(3): - self.aero_coord_indices = np.hstack( - [self.aero_coord_indices, aero_gid * 3 + i] - ) - self.aero_coord_indices = list( - np.sort(self.aero_coord_indices).astype("int") - ) + if ( + aero_node_ids is not None + ): # use grid ids from aero_builder's get_tagged_indices + # account for xyz of each node + self.aero_coord_indices = np.zeros(len(aero_node_ids) * 3, dtype=int) + for ii in range(3): + self.aero_coord_indices[ii::3] = 3 * aero_node_ids + ii + self.aero_coord_indices = list(self.aero_coord_indices) else: # use all indices self.aero_coord_indices = list(np.arange(aero_nnodes * 3)) @@ -582,7 +576,7 @@ def initialize(self, comm): if len(self.body_tags) > 0: # body tags given for i in range(len(self.body_tags)): try: - aero_gid = np.atleast_1d( + aero_node_ids = np.atleast_1d( self.aero_builder.get_tagged_indices(self.body_tags[i]["aero"]) ) except NotImplementedError: @@ -590,9 +584,9 @@ def initialize(self, comm): print( "get_tagged_indices has not been implemented in the aero builder; all nodes will be used" ) - aero_gid = None + aero_node_ids = None try: - struct_gid = np.atleast_1d( + struct_node_ids = np.atleast_1d( self.struct_builder.get_tagged_indices( self.body_tags[i]["struct"] ) @@ -602,16 +596,16 @@ def initialize(self, comm): print( "get_tagged_indices has not been implemented in the struct builder; all nodes will be used" ) - struct_gid = None + struct_node_ids = None self.bodies += [ MeldBodyInstance( comm, isym=self.isym, n=self.n[i], beta=self.beta[i], - aero_gid=aero_gid, + aero_node_ids=aero_node_ids, aero_nnodes=self.nnodes_aero, - struct_gid=struct_gid, + struct_node_ids=struct_node_ids, struct_nnodes=self.nnodes_struct, struct_ndof=self.ndof_struct, linearized=self.linearized[i],