diff --git a/.vscode/settings.json b/.vscode/settings.json index f5d797e4..66978c95 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -15,5 +15,6 @@ ], "python.testing.pytestEnabled": false, "python.testing.nosetestsEnabled": false, - "python.testing.unittestEnabled": true + "python.testing.unittestEnabled": true, + "python.analysis.useImportHeuristic": true } diff --git a/examples/aero_only/agard/run_fun3d.py b/examples/aero_only/agard/run_fun3d.py index 6e215ef3..a501015d 100644 --- a/examples/aero_only/agard/run_fun3d.py +++ b/examples/aero_only/agard/run_fun3d.py @@ -10,14 +10,13 @@ class Top(Multipoint): def setup(self): # FUN3D options - project_rootname = 'agard_debug' boundary_tags = [3] mach = 0.5, - aoa = 0.0 + aoa = 1.0 yaw = 0.0 q_inf = 3000. - reynolds = 1.0 + reynolds = 0.0 dvs = self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*']) dvs.add_output('aoa', val=aoa, units='deg') @@ -26,7 +25,7 @@ def setup(self): dvs.add_output('q_inf', q_inf) dvs.add_output('reynolds', reynolds) - aero_builder = Fun3dSfeBuilder(project_rootname, boundary_tags) + aero_builder = Fun3dSfeBuilder(boundary_tags) aero_builder.initialize(self.comm) self.add_subsystem('mesh',aero_builder.get_mesh_coordinate_subsystem()) @@ -38,10 +37,12 @@ def setup(self): prob = om.Problem() prob.model = Top() -prob.setup() +prob.setup(mode='rev') om.n2(prob, show_browser=False, outfile='mphys_fun3d.html') prob.run_model() if MPI.COMM_WORLD.rank == 0: - print('C_L, C_D =',prob['cruise.C_L'], prob['cruise.C_D']) + print('C_L, C_D =',prob['cruise.C_L'], prob['cruise.C_D'], prob['cruise.C_X'],prob['cruise.C_Z']) + +prob.check_totals(of='cruise.C_L',wrt='aoa') diff --git a/examples/aerostructural/mach_tutorial_wing/vlm_meld_tacs/run_analysis_modal.py b/examples/aerostructural/mach_tutorial_wing/vlm_meld_tacs/run_analysis_modal.py new file mode 100644 index 00000000..bafe633a --- /dev/null +++ b/examples/aerostructural/mach_tutorial_wing/vlm_meld_tacs/run_analysis_modal.py @@ -0,0 +1,90 @@ +#rst Imports +from __future__ import print_function, division +import numpy as np +from mpi4py import MPI + +import openmdao.api as om + +from mphys import Multipoint +from mphys.scenario_aerostructural import ScenarioAeroStructural +from mphys.mphys_vlm import VlmBuilder +from mphys.mphys_modal_solver import ModalBuilder +from mphys.mphys_meld import MeldBuilder + +import tacs_setup + +class Top(Multipoint): + def setup(self): + # VLM + mesh_file = 'wing_VLM.dat' + mach = 0.85 + aoa0 = 2.0 + aoa1 = 5.0 + q_inf = 3000. + vel = 178. + nu = 3.5E-5 + + aero_builder = VlmBuilder(mesh_file) + aero_builder.initialize(self.comm) + + dvs = self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*']) + dvs.add_output('aoa', val=[aoa0,aoa1], units='deg') + dvs.add_output('mach', mach) + dvs.add_output('q_inf', q_inf) + dvs.add_output('vel', vel) + dvs.add_output('nu', nu) + + self.add_subsystem('mesh_aero',aero_builder.get_mesh_coordinate_subsystem()) + + # TACS + tacs_options = {'add_elements': tacs_setup.add_elements, + 'get_funcs' : tacs_setup.get_funcs, + 'mesh_file' : 'wingbox_Y_Z_flip.bdf', + 'f5_writer' : tacs_setup.f5_writer } + + nmodes = 50 + struct_builder = ModalBuilder(tacs_options, nmodes) + struct_builder.initialize(self.comm) + ndv_struct = struct_builder.get_ndv() + + self.add_subsystem('mesh_struct',struct_builder.get_mesh_coordinate_subsystem()) + + dvs.add_output('dv_struct', np.array(ndv_struct*[0.002])) + self.connect('dv_struct','mesh_struct.dv_struct') + + # MELD setup + isym = 1 + ldxfer_builder = MeldBuilder(aero_builder, struct_builder, isym=isym) + ldxfer_builder.initialize(self.comm) + + for iscen, scenario in enumerate(['cruise','maneuver']): + nonlinear_solver = om.NonlinearBlockGS(maxiter=25, iprint=2, use_aitken=True, rtol = 1E-14, atol=1E-14) + linear_solver = om.LinearBlockGS(maxiter=25, iprint=2, use_aitken=True, rtol = 1e-14, atol=1e-14) + self.mphys_add_scenario(scenario,ScenarioAeroStructural(aero_builder=aero_builder, + struct_builder=struct_builder, + ldxfer_builder=ldxfer_builder), + nonlinear_solver, linear_solver) + + for discipline in ['aero','struct']: + self.mphys_connect_scenario_coordinate_source('mesh_%s' % discipline, scenario, discipline) + + for modal_var in ['mode_shape', 'modal_stiffness']: + self.connect(f'mesh_struct.{modal_var}',f'{scenario}.{modal_var}') + for dv in ['q_inf','vel','nu','mach', 'dv_struct']: + self.connect(dv, f'{scenario}.{dv}') + self.connect('aoa', f'{scenario}.aoa', src_indices=[iscen]) + +################################################################################ +# OpenMDAO setup +################################################################################ +prob = om.Problem() +prob.model = Top() + +prob.setup() +om.n2(prob, show_browser=False, outfile='mphys_as_vlm_modal.html') + +prob.run_model() + +if MPI.COMM_WORLD.rank == 0: + print('C_L =',prob['cruise.C_L']) + print('func_struct =',prob['maneuver.func_struct']) diff --git a/mphys/__init__.py b/mphys/__init__.py index 64e35410..78470ab4 100644 --- a/mphys/__init__.py +++ b/mphys/__init__.py @@ -1,3 +1,4 @@ #!/usr/bin/env python from .builder import Builder from .multipoint import Multipoint, MultipointParallel +from .distributed_converter import DistributedConverter, DistributedVariableDescription diff --git a/mphys/distributed_converter.py b/mphys/distributed_converter.py new file mode 100644 index 00000000..304b6085 --- /dev/null +++ b/mphys/distributed_converter.py @@ -0,0 +1,83 @@ +import openmdao.api as om + + +class DistributedVariableDescription: + """ + Attributes of a variable for conversion from serial to distributed + or distributed to serial + """ + + def __init__(self, name: str, shape: tuple, tags=[]): + self.name = name + self.shape = shape + self.tags = tags + + +class DistributedConverter(om.ExplicitComponent): + """ + An ExplicitComponent to convert from distributed to serial and serial to distributed variables. + Mphys requires the coupling inputs and outputs to be distributed variables, so this + class is provided to help with those conversions. + For each mphys variable, a {variable}_serial version is created for the nonparallel solver to connect to and the + distributed version will have the full vector on the root processor and zero length on the other processors. + Given a list of distributed inputs in the options, the component will add variables to the inputs as distributed and + produce {variable}_serial as outputs. + Given a list of distributed outputs in the options, the component will add variables to the outputs as distributed and + add {variable}_serial as inputs. + + """ + + def initialize(self): + self.options.declare( + 'distributed_inputs', default=[], + desc='List of DistributedVariableDescription objects that will be converted from distributed to serial') + self.options.declare( + 'distributed_outputs', default=[], + desc='List of DistributedVariableDescription objects that will be converted from serial to distributed') + + def setup(self): + for input in self.options['distributed_inputs']: + self.add_input(input.name, shape_by_conn=True, tags=input.tags, distributed=True) + self.add_output(f'{input.name}_serial', shape=input.shape, distributed=False) + + for output in self.options['distributed_outputs']: + shape = output.shape if self.comm.Get_rank() == 0 else 0 + self.add_input(f'{output.name}_serial', shape_by_conn=True, distributed=False) + self.add_output(output.name, shape=shape, tags=output.tags, distributed=True) + + def compute(self, inputs, outputs): + for input in self.options['distributed_inputs']: + if self.comm.Get_rank() == 0: + outputs[f'{input.name}_serial'] = inputs[input.name] + outputs[f'{input.name}_serial'] = self.comm.bcast(outputs[f'{input.name}_serial']) + + for output in self.options['distributed_outputs']: + if self.comm.Get_rank() == 0: + outputs[output.name] = inputs[f'{output.name}_serial'] + + def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): + if mode == 'fwd': + for input in self.options['distributed_inputs']: + if input.name in d_inputs and f'{input.name}_serial' in d_outputs: + if self.comm.Get_rank() == 0: + d_outputs[f'{input.name}_serial'] += d_inputs[input.name] + d_outputs[f'{input.name}_serial'] = self.comm.bcast( + d_outputs[f'{input.name}_serial']) + + for output in self.options['distributed_outputs']: + if output.name in d_outputs and f'{output.name}_serial' in d_inputs: + if self.comm.Get_rank() == 0: + d_outputs[output.name] += d_inputs[f'{output.name}_serial'] + + if mode == 'rev': + for input in self.options['distributed_inputs']: + if input.name in d_inputs and f'{input.name}_serial' in d_outputs: + if self.comm.Get_rank() == 0: + d_inputs[input.name] += d_outputs[f'{input.name}_serial'] + + for output in self.options['distributed_outputs']: + if output.name in d_outputs and f'{output.name}_serial' in d_inputs: + if self.comm.Get_rank() == 0: + d_inputs[f'{output.name}_serial'] += d_outputs[output.name] + d_inputs[f'{output.name}_serial'] = self.comm.bcast( + d_inputs[f'{output.name}_serial']) diff --git a/mphys/geo_disp.py b/mphys/geo_disp.py index 8feb25f6..c12874bd 100644 --- a/mphys/geo_disp.py +++ b/mphys/geo_disp.py @@ -6,7 +6,6 @@ class GeoDisp(om.ExplicitComponent): displacements to the geometry-changed aerodynamic surface """ def initialize(self): - self.options['distributed'] = True self.options.declare('number_of_nodes') def setup(self): @@ -14,13 +13,16 @@ def setup(self): local_size = nnodes * 3 self.add_input('x_aero0', shape_by_conn=True, + distributed=True, desc='aerodynamic surface with geom changes', tags=['mphys_coordinates']) self.add_input('u_aero', shape_by_conn=True, + distributed=True, desc='aerodynamic surface displacements', tags=['mphys_coupling']) self.add_output('x_aero', shape=local_size, + distributed=True, desc='deformed aerodynamic surface', tags=['mphys_coupling']) diff --git a/mphys/integrated_forces.py b/mphys/integrated_forces.py index ec7792bd..17151cbc 100644 --- a/mphys/integrated_forces.py +++ b/mphys/integrated_forces.py @@ -110,7 +110,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): if mode == 'fwd': if 'aoa' in d_inputs: - daoa_rad = np.pi / 180.0 * d_inputs['aoa'] + daoa_rad = d_inputs['aoa'] if 'Lift' in d_outputs or 'C_L' in d_outputs: d_lift_d_aoa = ( - fx_total * np.cos(aoa) * daoa_rad - fz_total * np.sin(aoa) * daoa_rad ) @@ -128,7 +128,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): d_outputs['C_D'] += d_drag_d_aoa / (q_inf * area) if 'yaw' in d_inputs: - dyaw_rad = np.pi / 180.0 * d_inputs['yaw'] + dyaw_rad = d_inputs['yaw'] if 'Drag' in d_outputs or 'C_D' in d_outputs: d_drag_d_yaw = ( fx_total * np.cos(aoa) * (-np.sin(yaw) * dyaw_rad) - fy_total * np.cos(yaw) * dyaw_rad @@ -271,27 +271,24 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): if 'Lift' in d_outputs or 'C_L' in d_outputs: d_lift = d_outputs['Lift'] if 'Lift' in d_outputs else 0.0 d_cl = d_outputs['C_L'] if 'C_L' in d_outputs else 0.0 - d_lift_d_aoa_rad = ( - fx_total * np.cos(aoa) + d_inputs['aoa'] += ( - fx_total * np.cos(aoa) - fz_total * np.sin(aoa) ) * (d_lift + d_cl / (q_inf * area)) - d_inputs['aoa'] += d_lift_d_aoa_rad * np.pi / 180.0 if 'Drag' in d_outputs or 'C_D' in d_outputs: d_drag = d_outputs['Drag'] if 'Drag' in d_outputs else 0.0 d_cd = d_outputs['C_D'] if 'C_D' in d_outputs else 0.0 - d_drag_d_aoa_rad = ( fx_total * (-np.sin(aoa)) * np.cos(yaw) + d_inputs['aoa'] += ( fx_total * (-np.sin(aoa)) * np.cos(yaw) + fz_total * ( np.cos(aoa)) * np.cos(yaw) ) * (d_drag + d_cd / (q_inf * area)) - d_inputs['aoa'] += d_drag_d_aoa_rad * np.pi / 180.0 if 'yaw' in d_inputs: if 'Drag' in d_outputs or 'C_D' in d_outputs: d_drag = d_outputs['Drag'] if 'Drag' in d_outputs else 0.0 d_cd = d_outputs['C_D'] if 'C_D' in d_outputs else 0.0 - d_drag_d_yaw_rad = ( fx_total * np.cos(aoa) * (-np.sin(yaw)) + d_inputs['yaw'] += ( fx_total * np.cos(aoa) * (-np.sin(yaw)) - fy_total * np.cos(yaw) + fz_total * np.sin(aoa) * (-np.sin(yaw)) ) * (d_drag + d_cd / (q_inf * area)) - d_inputs['yaw'] += d_drag_d_yaw_rad * np.pi / 180.0 if 'ref_area' in d_inputs: d_nondim = - 1.0 / (q_inf * area**2.0) @@ -423,8 +420,8 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): nnodes = 3 prob = Problem() ivc = IndepVarComp() - ivc.add_output('aoa',val=45.0) - ivc.add_output('yaw',val=135.0) + ivc.add_output('aoa',val=45.0, units='deg') + ivc.add_output('yaw',val=135.0, units='deg') ivc.add_output('ref_area',val=0.2) ivc.add_output('moment_center',shape=3,val=np.zeros(3)) ivc.add_output('ref_length', val = 3.0) diff --git a/mphys/mphys_adflow.py b/mphys/mphys_adflow.py index ff04fe9d..8f25793e 100644 --- a/mphys/mphys_adflow.py +++ b/mphys/mphys_adflow.py @@ -21,7 +21,6 @@ class ADflowMesh(ExplicitComponent): def initialize(self): self.options.declare("aero_solver", recordable=False) - self.options["distributed"] = True def setup(self): @@ -31,11 +30,17 @@ def setup(self): coord_size = self.x_a0.size self.add_output( - "x_aero0", shape=coord_size, desc="initial aerodynamic surface node coordinates", tags=["mphys_coordinates"] + "x_aero0", + distributed=True, + shape=coord_size, + desc="initial aerodynamic surface node coordinates", + tags=["mphys_coordinates"], ) def mphys_add_coordinate_input(self): - self.add_input("x_aero0_points", shape_by_conn=True, desc="aerodynamic surface with geom changes") + self.add_input( + "x_aero0_points", distributed=True, shape_by_conn=True, desc="aerodynamic surface with geom changes" + ) # return the promoted name and coordinates return "x_aero0_points", self.x_a0 @@ -139,8 +144,6 @@ def initialize(self): # self.options.declare('use_OM_KSP', default=False, types=bool, # desc="uses OpenMDAO's PestcKSP linear solver with ADflow's preconditioner to solve the adjoint.") - self.options["distributed"] = True - def setup(self): # self.set_check_partial_options(wrt='*',directional=True) @@ -152,8 +155,8 @@ def setup(self): # state inputs and outputs local_volume_coord_size = solver.mesh.getSolverGrid().size - self.add_input("x_aero", shape_by_conn=True, tags=["mphys_coupling"]) - self.add_output("adflow_vol_coords", shape=local_volume_coord_size, tags=["mphys_coupling"]) + self.add_input("x_aero", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_output("adflow_vol_coords", distributed=True, shape=local_volume_coord_size, tags=["mphys_coupling"]) # self.declare_partials(of='adflow_vol_coords', wrt='x_aero') @@ -202,8 +205,6 @@ def initialize(self): self.options.declare("restart_failed_analysis", default=False) self.options.declare("err_on_convergence_fail", default=False) - self.options["distributed"] = True - # testing flag used for unit-testing to prevent the call to actually solve # NOT INTENDED FOR USERS!!! FOR TESTING ONLY self._do_solve = True @@ -228,8 +229,8 @@ def setup(self): # state inputs and outputs local_state_size = solver.getStateSize() - self.add_input("adflow_vol_coords", shape_by_conn=True, tags=["mphys_coupling"]) - self.add_output("adflow_states", shape=local_state_size, tags=["mphys_coupling"]) + self.add_input("adflow_vol_coords", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_output("adflow_states", distributed=True, shape=local_state_size, tags=["mphys_coupling"]) # self.declare_partials(of='adflow_states', wrt='*') @@ -260,7 +261,9 @@ def set_ap(self, ap): for (args, kwargs) in self.ap_vars: name = args[0] size = args[1] - self.add_input(name, shape=size, val=kwargs["value"], units=kwargs["units"], tags=["mphys_input"]) + self.add_input( + name, distributed=False, shape=size, val=kwargs["value"], units=kwargs["units"], tags=["mphys_input"] + ) if self.comm.rank == 0: print("%s (%s)" % (name, kwargs["units"])) @@ -456,19 +459,17 @@ class ADflowForces(ExplicitComponent): def initialize(self): self.options.declare("aero_solver", recordable=False) - self.options["distributed"] = True - def setup(self): # self.set_check_partial_options(wrt='*',directional=True) self.solver = self.options["aero_solver"] solver = self.solver - self.add_input("adflow_vol_coords", shape_by_conn=True, tags=["mphys_coupling"]) - self.add_input("adflow_states", shape_by_conn=True, tags=["mphys_coupling"]) + self.add_input("adflow_vol_coords", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_input("adflow_states", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) local_surface_coord_size = solver.mesh.getSurfaceCoordinates().size - self.add_output("f_aero", shape=local_surface_coord_size, tags=["mphys_coupling"]) + self.add_output("f_aero", distributed=True, shape=local_surface_coord_size, tags=["mphys_coupling"]) # self.declare_partials(of='f_aero', wrt='*') @@ -494,7 +495,7 @@ def set_ap(self, ap): for (args, kwargs) in self.ap_vars: name = args[0] size = args[1] - self.add_input(name, shape=size, units=kwargs["units"], tags=["mphys_input"]) + self.add_input(name, distributed=False, shape=size, units=kwargs["units"], tags=["mphys_input"]) # if self.comm.rank == 0: # print('%s (%s)'%(name, kwargs['units'])) @@ -565,8 +566,6 @@ class AdflowHeatTransfer(ExplicitComponent): def initialize(self): self.options.declare("aero_solver") - self.options["distributed"] = True - def setup(self): # self.set_check_partial_options(wrt='*',directional=True) @@ -575,11 +574,16 @@ def setup(self): local_nodes, _ = solver._getSurfaceSize(solver.allIsothermalWallsGroup) - self.add_input("adflow_vol_coords", shape_by_conn=True, tags=["mphys_coupling"]) - self.add_input("adflow_states", shape_by_conn=True, tags=["mphys_coupling"]) + self.add_input("adflow_vol_coords", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_input("adflow_states", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) self.add_output( - "q_convect", val=np.ones(local_nodes) * -499, shape=local_nodes, units="W/m**2", tags=["mphys_coupling"] + "q_convect", + distributed=True, + val=np.ones(local_nodes) * -499, + shape=local_nodes, + units="W/m**2", + tags=["mphys_coupling"], ) # self.declare_partials(of='f_aero', wrt='*') @@ -604,7 +608,7 @@ def set_ap(self, ap): for (args, kwargs) in self.ap_vars: name = args[0] size = args[1] - self.add_input(name, shape=size, units=kwargs["units"], tags=["mphys_input"]) + self.add_input(name, distributed=False, shape=size, units=kwargs["units"], tags=["mphys_input"]) if self.comm.rank == 0: print(name) @@ -747,19 +751,8 @@ def setup(self): n1 = np.sum(n_list[:irank]) n2 = np.sum(n_list[: irank + 1]) - self.add_input( - "adflow_vol_coords", - src_indices=np.arange(n1, n2, dtype=int), - shape=local_coord_size, - tags=["mphys_coupling"], - ) - self.add_input( - "adflow_states", src_indices=np.arange(s1, s2, dtype=int), shape=local_state_size, tags=["mphys_coupling"] - ) - - # TODO shape by conn does not work for these - # self.add_input('adflow_vol_coords', shape_by_conn=True, tags=['mphys_coupling']) - # self.add_input('adflow_states', shape_by_conn=True, tags=['mphys_coupling']) + self.add_input("adflow_vol_coords", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_input("adflow_states", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) # self.declare_partials(of=f_name, wrt='*') @@ -785,7 +778,7 @@ def mphys_set_ap(self, ap): for (args, kwargs) in self.ap_vars: name = args[0] size = args[1] - self.add_input(name, shape=size, units=kwargs["units"], tags=["mphys_input"]) + self.add_input(name, distributed=False, shape=size, units=kwargs["units"], tags=["mphys_input"]) # if self.comm.rank == 0: # print('%s with units %s'%(name, kwargs['units'])) @@ -807,7 +800,7 @@ def mphys_set_ap(self, ap): if self.comm.rank == 0: print("%s (%s)" % (f_name, units)) - self.add_output(f_name, shape=1, units=units, tags=["mphys_result"]) + self.add_output(f_name, distributed=False, shape=1, units=units, tags=["mphys_result"]) # def mphys_add_prop_funcs(self, prop_funcs): # save this list @@ -838,7 +831,7 @@ def mphys_add_funcs(self, funcs): # if self.comm.rank == 0: # print("%s (%s)" % (f_name, units)) - self.add_output(f_name, shape=1, units=units, tags=["mphys_result"]) + self.add_output(f_name, distributed=False, shape=1, units=units, tags=["mphys_result"]) def _set_states(self, inputs): self.solver.setStates(inputs["adflow_states"]) diff --git a/mphys/mphys_dvgeo.py b/mphys/mphys_dvgeo.py index 21c4cd0c..532173c6 100644 --- a/mphys/mphys_dvgeo.py +++ b/mphys/mphys_dvgeo.py @@ -18,7 +18,6 @@ def initialize(self): self.options.declare("ffd_file", default=None) self.options.declare("vsp_file", default=None) self.options.declare("vsp_options", default=None) - self.options["distributed"] = True def setup(self): @@ -60,14 +59,14 @@ def nom_add_discipline_coords(self, discipline, points=None): if points is None: # no pointset info is provided, just do a generic i/o. We will add these points during the first compute - self.add_input("x_%s_in" % discipline, shape_by_conn=True) - self.add_output("x_%s0" % discipline, copy_shape="x_%s_in" % discipline) + self.add_input("x_%s_in" % discipline, distributed=True, shape_by_conn=True) + self.add_output("x_%s0" % discipline, distributed=True, copy_shape="x_%s_in" % discipline) else: # we are provided with points. we can do the full initialization now self.nom_addPointSet(points, "x_%s0" % discipline, add_output=False) - self.add_input("x_%s_in" % discipline, val=points.flatten()) - self.add_output("x_%s0" % discipline, val=points.flatten()) + self.add_input("x_%s_in" % discipline, distributed=True, val=points.flatten()) + self.add_output("x_%s0" % discipline, distributed=True, val=points.flatten()) def nom_addPointSet(self, points, ptName, add_output=True, **kwargs): # add the points to the dvgeo object @@ -76,7 +75,7 @@ def nom_addPointSet(self, points, ptName, add_output=True, **kwargs): if add_output: # add an output to the om component - self.add_output(ptName, val=points.flatten()) + self.add_output(ptName, distributed=True, val=points.flatten()) def nom_add_point_dict(self, point_dict): # add every pointset in the dict, and set the ptset name as the key @@ -85,14 +84,14 @@ def nom_add_point_dict(self, point_dict): def nom_addGeoDVGlobal(self, dvName, value, func): # define the input - self.add_input(dvName, shape=len(value)) + self.add_input(dvName, distributed=False, shape=len(value)) # call the dvgeo object and add this dv self.DVGeo.addGeoDVGlobal(dvName, value, func) def nom_addGeoDVLocal(self, dvName, axis="y"): nVal = self.DVGeo.addGeoDVLocal(dvName, axis=axis) - self.add_input(dvName, shape=nVal) + self.add_input(dvName, distributed=False, shape=nVal) return nVal def nom_addVSPVariable(self, component, group, parm, **kwargs): @@ -107,31 +106,31 @@ def nom_addVSPVariable(self, component, group, parm, **kwargs): val = self.DVGeo.DVs[dvName].value.copy() # add the input with the correct value, VSP DVs always have a size of 1 - self.add_input(dvName, shape=1, val=val) + self.add_input(dvName, distributed=False, shape=1, val=val) def nom_addThicknessConstraints2D(self, name, leList, teList, nSpan=10, nChord=10): self.DVCon.addThicknessConstraints2D(leList, teList, nSpan, nChord, lower=1.0, name=name) comm = self.comm if comm.rank == 0: - self.add_output(name, val=np.ones((nSpan * nChord,)), shape=nSpan * nChord) + self.add_output(name, distributed=True, val=np.ones((nSpan * nChord,)), shape=nSpan * nChord) else: - self.add_output(name, shape=(0,)) + self.add_output(name, distributed=True, shape=(0,)) def nom_addThicknessConstraints1D(self, name, ptList, nCon, axis): self.DVCon.addThicknessConstraints1D(ptList, nCon, axis, name=name) comm = self.comm if comm.rank == 0: - self.add_output(name, val=np.ones(nCon), shape=nCon) + self.add_output(name, distributed=True, val=np.ones(nCon), shape=nCon) else: - self.add_output(name, shape=(0)) + self.add_output(name, distributed=True, shape=(0)) def nom_addVolumeConstraint(self, name, leList, teList, nSpan=10, nChord=10): self.DVCon.addVolumeConstraint(leList, teList, nSpan=nSpan, nChord=nChord, name=name) comm = self.comm if comm.rank == 0: - self.add_output(name, val=1.0) + self.add_output(name, distributed=True, val=1.0) else: - self.add_output(name, shape=0) + self.add_output(name, distributed=True, shape=0) def nom_add_LETEConstraint(self, name, volID, faceID): self.DVCon.addLeTeConstraints(volID, faceID, name=name) @@ -140,9 +139,9 @@ def nom_add_LETEConstraint(self, name, volID, faceID): nCon = len(conobj.indSetA) comm = self.comm if comm.rank == 0: - self.add_output(name, val=np.zeros((nCon,)), shape=nCon) + self.add_output(name, distributed=True, val=np.zeros((nCon,)), shape=nCon) else: - self.add_output(name, shape=0) + self.add_output(name, distributed=True, shape=0) return nCon def nom_addRefAxis(self, **kwargs): diff --git a/mphys/mphys_fun3d.py b/mphys/mphys_fun3d.py index 19a42905..e764bb49 100644 --- a/mphys/mphys_fun3d.py +++ b/mphys/mphys_fun3d.py @@ -16,13 +16,12 @@ class Fun3dMesh(om.IndepVarComp): def initialize(self): self.options.declare('meshdef_solver') self.options.declare('boundary_tag_list') - self.options['distributed'] = True def setup(self): boundary_tag_list = self.options['boundary_tag_list'] meshdef = self.options['meshdef_solver'] x_aero0 = meshdef.get_boundary_node_coordinates(boundary_tag_list, owned_only = True, concatenate=True) - self.add_output('x_aero0', val=x_aero0.flatten(), desc='initial aerodynamic surface node coordinates', + self.add_output('x_aero0', distributed=True, val=x_aero0.flatten(), desc='initial aerodynamic surface node coordinates', tags=['mphys_coordinates']) class Fun3dFsiSolverGroup(om.Group): diff --git a/mphys/mphys_meld.py b/mphys/mphys_meld.py index ea2306da..d8c0e7fb 100644 --- a/mphys/mphys_meld.py +++ b/mphys/mphys_meld.py @@ -15,8 +15,6 @@ def initialize(self): self.options.declare('aero_nnodes') self.options.declare('check_partials') - self.options['distributed'] = True - self.meld = None self.initialized_meld = False @@ -37,17 +35,21 @@ def setup(self): # inputs self.add_input('x_struct0', shape_by_conn=True, + distributed=True, desc='initial structural node coordinates', tags=['mphys_coordinates']) self.add_input('x_aero0', shape_by_conn=True, + distributed=True, desc='initial aero surface node coordinates', tags=['mphys_coordinates']) self.add_input('u_struct', shape_by_conn=True, + distributed=True, desc='structural node displacements', tags=['mphys_coupling']) # outputs self.add_output('u_aero', shape = self.aero_nnodes*3, + distributed=True, val=np.zeros(self.aero_nnodes*3), desc='aerodynamic surface displacements', tags=['mphys_coupling']) @@ -148,8 +150,6 @@ def initialize(self): self.options.declare('aero_nnodes') self.options.declare('check_partials') - self.options['distributed'] = True - self.meld = None self.initialized_meld = False @@ -174,20 +174,25 @@ def setup(self): # inputs self.add_input('x_struct0', shape_by_conn=True, + distributed=True, desc='initial structural node coordinates', tags=['mphys_coordinates']) self.add_input('x_aero0', shape_by_conn=True, + distributed=True, desc='initial aero surface node coordinates', tags=['mphys_coordinates']) self.add_input('u_struct', shape_by_conn=True, + distributed=True, desc='structural node displacements', tags=['mphys_coupling']) self.add_input('f_aero', shape_by_conn=True, + distributed=True, desc='aerodynamic force vector', tags=['mphys_coupling']) # outputs self.add_output('f_struct', shape = struct_nnodes*struct_ndof, + distributed=True, desc='structural force vector', tags=['mphys_coupling']) diff --git a/mphys/mphys_meldthermal.py b/mphys/mphys_meldthermal.py index 440d9314..e1da0d7c 100644 --- a/mphys/mphys_meldthermal.py +++ b/mphys/mphys_meldthermal.py @@ -20,8 +20,6 @@ def initialize(self): self.options.declare('check_partials') self.options.declare('mapping') - self.options['distributed'] = True - self.meldThermal = None self.initialized_meld = False @@ -40,14 +38,15 @@ def setup(self): conv_nnodes = self.conv_nnodes # inputs - self.add_input('x_struct0', shape_by_conn=True, desc='initial structural node coordinates') - self.add_input('x_aero0', shape_by_conn=True, desc='initial aerodynamic surface node coordinates') - self.add_input('T_conduct', shape_by_conn=True, desc='conductive node displacements') + self.add_input('x_struct0', distributed=True, shape_by_conn=True, desc='initial structural node coordinates') + self.add_input('x_aero0', distributed=True, shape_by_conn=True, desc='initial aerodynamic surface node coordinates') + self.add_input('T_conduct', distributed=True, shape_by_conn=True, desc='conductive node displacements') # outputs print('T_convect', conv_nnodes) self.add_output('T_convect', shape = conv_nnodes, + distributed=True, val=np.ones(conv_nnodes)*301, desc='conv surface temperatures') @@ -98,8 +97,6 @@ def initialize(self): self.options.declare('check_partials') self.options.declare('mapping') - self.options['distributed'] = True - self.meldThermal = None self.initialized_meld = False @@ -118,14 +115,14 @@ def setup(self): self.check_partials= self.options['check_partials'] # inputs - self.add_input('x_struct0', shape_by_conn=True, desc='initial structural node coordinates') - self.add_input('x_aero0', shape_by_conn=True, desc='initial aerodynamic surface node coordinates') - self.add_input('q_convect', shape_by_conn=True, desc='initial conv heat transfer rate') + self.add_input('x_struct0', distributed=True, shape_by_conn=True, desc='initial structural node coordinates') + self.add_input('x_aero0', distributed=True, shape_by_conn=True, desc='initial aerodynamic surface node coordinates') + self.add_input('q_convect', distributed=True, shape_by_conn=True, desc='initial conv heat transfer rate') print('q_conduct', self.cond_nnodes) # outputs - self.add_output('q_conduct', shape = self.cond_nnodes, desc='heat transfer rate on the conduction mesh at the interface') + self.add_output('q_conduct', distributed=True, shape = self.cond_nnodes, desc='heat transfer rate on the conduction mesh at the interface') def compute(self, inputs, outputs): diff --git a/mphys/mphys_modal_solver.py b/mphys/mphys_modal_solver.py index 534f21ef..d7b7173a 100644 --- a/mphys/mphys_modal_solver.py +++ b/mphys/mphys_modal_solver.py @@ -1,30 +1,27 @@ #!/usr/bin/env python -from __future__ import print_function import numpy as np from mpi4py import MPI -from tacs import TACS, elements, functions +from tacs import TACS import openmdao.api as om +from mphys import Builder +from mphys.mphys_tacs import TACSFuncsGroup class ModalDecomp(om.ExplicitComponent): def initialize(self): - self.options.declare('struct_solver', default = None, desc='tacs object') + self.options.declare('tacs_assembler', default = None, desc='tacs object') self.options.declare('ndv', default = 1, desc='number of design variables in tacs') self.options.declare('nmodes', default = 15, desc = 'number of modes to kept') - self.options['distributed'] = True def setup(self): - - # TACS assembler setup - self.tacs = self.options['struct_solver'] + self.tacs_assembler = self.options['tacs_assembler'] self.ndv = self.options['ndv'] self.nmodes = self.options['nmodes'] if self.comm.rank == 0: - # create some TACS bvecs that will be needed later - self.xpts = self.tacs.createNodeVec() - self.tacs.getNodes(self.xpts) + self.xpts = self.tacs_assembler.createNodeVec() + self.tacs_assembler.getNodes(self.xpts) - self.vec = self.tacs.createVec() + self.vec = self.tacs_assembler.createVec() # OpenMDAO setup node_size = self.xpts.getArray().size @@ -34,29 +31,30 @@ def setup(self): node_size = 0 self.state_size = 0 - self.add_input('dv_struct',shape=self.ndv, desc='structural design variables') + self.add_input('dv_struct', shape_by_conn=True, desc='structural design variables', tags=['mphys_input']) # instead of using 2d arrays for the mode shapes, we use a flattened array with modes back to back. # this is because in OpenMDAO version 2.10.1, the src_indices option for 2D arrays with empty i/o on # some procs is broken. Flattened works at least for the analysis. - self.add_output('mode_shape', shape=self.nmodes*self.state_size, desc='structural mode shapes') + self.add_output('mode_shape', shape=self.nmodes*self.state_size, + distributed=True, + desc='structural mode shapes', + tags=['mphys_coupling']) - if self.comm.rank == 0: - output_shape = self.nmodes - else: - output_shape = 0 + output_shape = self.nmodes + self.add_output('modal_mass', shape=output_shape, desc='modal mass', tags=['mphys_coupling']) + self.add_output('modal_stiffness', shape=output_shape, desc='modal stiffness', tags=['mphys_coupling']) - self.add_output('modal_mass', shape=output_shape, desc='modal mass') - self.add_output('modal_stiffness', shape=output_shape, desc='modal stiffness') - - self.add_output('x_struct0', shape = node_size, desc = 'undeformed nodal coordinates') + self.add_output('x_struct0', shape = node_size, + distributed = True, + desc = 'undeformed nodal coordinates', tags=['mphys_coodinates']) def compute(self,inputs,outputs): if self.comm.rank == 0: - self.tacs.setDesignVars(np.array(inputs['dv_struct'],dtype=TACS.dtype)) + self.tacs_assembler.setDesignVars(np.array(inputs['dv_struct'],dtype=TACS.dtype)) - kmat = self.tacs.createFEMat() - self.tacs.assembleMatType(TACS.PY_STIFFNESS_MATRIX,kmat,TACS.PY_NORMAL) + kmat = self.tacs_assembler.createFEMat() + self.tacs_assembler.assembleMatType(TACS.PY_STIFFNESS_MATRIX,kmat,TACS.PY_NORMAL) pc = TACS.Pc(kmat) subspace = 100 restarts = 2 @@ -66,10 +64,10 @@ def compute(self,inputs,outputs): sigma_hz = 1.0 sigma = 2.0*np.pi*sigma_hz - mmat = self.tacs.createFEMat() - self.tacs.assembleMatType(TACS.PY_MASS_MATRIX,mmat,TACS.PY_NORMAL) + mmat = self.tacs_assembler.createFEMat() + self.tacs_assembler.assembleMatType(TACS.PY_MASS_MATRIX,mmat,TACS.PY_NORMAL) - self.freq = TACS.FrequencyAnalysis(self.tacs, sigma, mmat, kmat, self.gmres, + self.freq = TACS.FrequencyAnalysis(self.tacs_assembler, sigma, mmat, kmat, self.gmres, num_eigs=self.nmodes, eig_tol=1e-12) self.freq.solve() @@ -78,62 +76,47 @@ def compute(self,inputs,outputs): eig, err = self.freq.extractEigenvector(imode,self.vec) outputs['modal_mass'][imode] = 1.0 outputs['modal_stiffness'][imode] = eig - for idof in range(3): - # put every mode back to back instead of using a 2d array bec. the pseudo parallelism breaks that way... This should be fixed in OM and we can go back to using 2D arrays - outputs['mode_shape'][imode*self.state_size:(imode+1)*self.state_size] = self.vec.getArray() + outputs['mode_shape'][imode*self.state_size:(imode+1)*self.state_size] = self.vec.getArray() + + outputs['modal_mass'] = self.comm.bcast(outputs['modal_mass']) + outputs['modal_stiffness'] = self.comm.bcast(outputs['modal_stiffness']) class ModalSolver(om.ExplicitComponent): """ - Steady Modal structural solver + Steady modal structural solver K z - mf = 0 """ def initialize(self): self.options.declare('nmodes',default=1) - self.options['distributed'] = True def setup(self): - # adjust the indices for procs - if self.comm.rank == 0: - nmodes = self.options['nmodes'] - src_indices = np.arange(0, nmodes, dtype=int) - input_shape = nmodes - else: - nmodes = 0 - src_indices = np.zeros(0) - input_shape=0 - - self.add_input('modal_stiffness', shape=input_shape, - src_indices=src_indices, - val=np.ones(nmodes), desc = 'modal stiffness') - self.add_input('mf', shape=input_shape, - src_indices=src_indices, - val=np.ones(nmodes), - desc = 'modal force') + nmodes = self.options['nmodes'] + + self.add_input('modal_stiffness', shape_by_conn=True, desc = 'modal stiffness', tags=['mphys_coupling']) + self.add_input('mf', shape_by_conn=True, desc = 'modal force') self.add_output('z', shape=nmodes, val=np.ones(nmodes), desc = 'modal displacement') def compute(self,inputs,outputs): - if self.comm.rank == 0: - k = inputs['modal_stiffness'] - outputs['z'] = inputs['mf'] / k + k = inputs['modal_stiffness'] + outputs['z'] = inputs['mf'] / k def compute_jacvec_product(self,inputs,d_inputs,d_outputs,mode): - if self.comm.rank == 0: - k = inputs['modal_stiffness'] - if mode == 'fwd': - if 'z' in d_outputs: - if 'mf' in d_inputs: - d_outputs['z'] += d_inputs['mf'] / k - if 'modal_stiffness' in d_inputs: - d_outputs['z'] += - inputs['mf'] / (k**2.0) * d_inputs['modal_stiffness'] - if mode == 'rev': - if 'z' in d_outputs: - if 'mf' in d_inputs: - d_inputs['mf'] += d_outputs['z'] / k - if 'modal_stiffness' in d_inputs: - d_inputs['modal_stiffness'] += - inputs['mf'] / (k**2.0) * d_outputs['z'] + k = inputs['modal_stiffness'] + if mode == 'fwd': + if 'z' in d_outputs: + if 'mf' in d_inputs: + d_outputs['z'] += d_inputs['mf'] / k + if 'modal_stiffness' in d_inputs: + d_outputs['z'] += - inputs['mf'] / (k**2.0) * d_inputs['modal_stiffness'] + if mode == 'rev': + if 'z' in d_outputs: + if 'mf' in d_inputs: + d_inputs['mf'] += d_outputs['z'] / k + if 'modal_stiffness' in d_inputs: + d_inputs['modal_stiffness'] += - inputs['mf'] / (k**2.0) * d_outputs['z'] class ModalForces(om.ExplicitComponent): def initialize(self): @@ -144,29 +127,14 @@ def setup(self): self.nmodes = self.options['nmodes'] self.mode_size = self.options['mode_size'] - # adjust the indices for procs - if self.comm.rank == 0: - src_indices = np.arange(0, self.nmodes*self.mode_size, dtype=int) - input_shape = self.nmodes*self.mode_size - else: - src_indices = np.zeros(0) - input_shape=0 - self.add_input( 'mode_shape', - shape=input_shape, - src_indices=src_indices, + shape_by_conn=True, + distributed=True, desc='structural mode shapes' ) - if self.comm.rank == 0: - src_indices = np.arange(0, self.mode_size, dtype=int) - input_shape = self.mode_size - else: - src_indices = np.zeros((0)) - input_shape=0 - - self.add_input('f_struct',shape=input_shape, src_indices=src_indices,desc = 'nodal force') + self.add_input('f_struct',shape_by_conn=True, distributed=True, desc = 'nodal force', tags=['mphys_coupling']) self.add_output('mf',shape=self.nmodes, desc = 'modal force') @@ -175,17 +143,20 @@ def compute(self,inputs,outputs): outputs['mf'][:] = 0.0 for imode in range(self.nmodes): outputs['mf'][imode] = np.sum(inputs['mode_shape'][imode*self.mode_size:(imode+1)*self.mode_size] * inputs['f_struct'][:]) + outputs['mf'] = self.comm.bcast(outputs['mf']) def compute_jacvec_product(self,inputs,d_inputs,d_outputs,mode): - if self.comm.rank == 0: - if mode=='fwd': - if 'mf' in d_outputs: - if 'f_struct' in d_inputs: + if mode=='fwd': + if 'mf' in d_outputs: + if 'f_struct' in d_inputs: + if self.comm.rank == 0: for imode in range(self.options['nmodes']): d_outputs['mf'][imode] += np.sum(inputs['mode_shape'][imode*self.mode_size:(imode+1)*self.mode_size] * d_inputs['f_struct'][:]) - if mode=='rev': - if 'mf' in d_outputs: - if 'f_struct' in d_inputs: + d_outputs['mf'] = self.comm.bcast(d_outputs['mf']) + if mode=='rev': + if 'mf' in d_outputs: + if 'f_struct' in d_inputs: + if self.comm.rank == 0: for imode in range(self.options['nmodes']): d_inputs['f_struct'][:] += inputs['mode_shape'][imode*self.mode_size:(imode+1)*self.mode_size] * d_outputs['mf'][imode] @@ -198,32 +169,10 @@ def setup(self): self.nmodes = self.options['nmodes'] self.mode_size = self.options['mode_size'] - # adjust the indices for procs - if self.comm.rank == 0: - src_indices = np.arange(0, self.nmodes*self.mode_size, dtype=int) - input_shape = self.nmodes*self.mode_size - else: - src_indices = np.zeros(0) - input_shape=0 - - self.add_input( - 'mode_shape', - shape=input_shape, - src_indices=src_indices, - desc='structural mode shapes' - ) - - if self.comm.rank == 0: - src_indices = np.arange(0, self.nmodes, dtype=int) - input_shape = self.nmodes - else: - src_indices = np.zeros(0) - input_shape = 0 - - self.add_input('z',shape=input_shape, src_indices=src_indices, desc='modal displacement') + self.add_input('mode_shape', shape_by_conn=True, distributed=True, desc='structural mode shapes', tags=['mphys_coupling']) + self.add_input('z', shape_by_conn=True, desc='modal displacement') - # its important that we set this to zero since this displacement value is used for the first iteration of the aero - self.add_output('u_struct', shape=self.mode_size, val = np.zeros(self.mode_size), desc = 'nodal displacement') + self.add_output('u_struct', shape=self.mode_size, val = np.zeros(self.mode_size), desc = 'nodal displacement', tags=['mphys_coupling']) def compute(self,inputs,outputs): if self.comm.rank == 0: @@ -232,27 +181,34 @@ def compute(self,inputs,outputs): outputs['u_struct'][:] += inputs['mode_shape'][imode*self.mode_size:(imode+1)*self.mode_size] * inputs['z'][imode] def compute_jacvec_product(self,inputs,d_inputs,d_outputs,mode): - if self.comm.rank == 0: - if mode=='fwd': - if 'u_struct' in d_outputs: - if 'z' in d_inputs: + if mode=='fwd': + if 'u_struct' in d_outputs: + if 'z' in d_inputs: + if self.comm.rank == 0: for imode in range(self.options['nmodes']): d_outputs['u_struct'][:] += inputs['mode_shape'][imode*self.mode_size:(imode+1)*self.mode_size] * d_inputs['z'][imode] - if mode=='rev': - if 'u_struct' in d_outputs: - if 'z' in d_inputs: + if 'mode_shape' in d_inputs: + if self.comm.rank == 0: + for imode in range(self.options['nmodes']): + d_outputs['u_struct'][:] += d_inputs['mode_shape'][imode*self.mode_size:(imode+1)*self.mode_size] * inputs['z'][imode] + if mode=='rev': + if 'u_struct' in d_outputs: + if 'z' in d_inputs: + if self.comm.rank == 0: for imode in range(self.options['nmodes']): d_inputs['z'][imode] += np.sum(inputs['mode_shape'][imode*self.mode_size:(imode+1)*self.mode_size] * d_outputs['u_struct'][:]) + d_inputs['z'] = self.comm.bcast(d_inputs['z']) + if 'mode_shape' in d_inputs: + if self.comm.rank == 0: + for imode in range(self.options['nmodes']): + d_inputs['mode_shape'][imode*self.mode_size:(imode+1)*self.mode_size] += inputs['z'][imode] * d_outputs['u_struct'][:] class ModalGroup(om.Group): def initialize(self): - self.options.declare('solver') - self.options.declare('solver_objects') self.options.declare('nmodes') self.options.declare('nnodes') self.options.declare('ndof') self.options.declare('check_partials') - self.options.declare('as_coupling') def setup(self): nmodes = self.options['nmodes'] @@ -273,118 +229,79 @@ def setup(self): promotes_outputs=['u_struct'] ) -# self.add_subsystem('funcs', TacsFunctions( -# struct_solver=self.struct_solver, -# struct_objects=self.struct_objects, -# check_partials=self.check_partials), -# promotes_inputs=['x_struct0', 'dv_struct'] -# ) -# -# self.add_subsystem('mass', TacsMass( -# struct_solver=self.struct_solver, -# struct_objects=self.struct_objects, -# check_partials=self.check_partials), -# promotes_inputs=['x_struct0', 'dv_struct'] -# ) - - def configure(self): - if self.comm.rank == 0: - self.connect('modal_forces.mf', 'modal_solver.mf') - self.connect('modal_solver.z', 'modal_disps.z') - -class ModalBuilder(object): + self.connect('modal_forces.mf', 'modal_solver.mf') + self.connect('modal_solver.z', 'modal_disps.z') - def __init__(self, options,nmodes=15,check_partials=False): +class ModalBuilder(Builder): + def __init__(self, options, nmodes, include_tacs_functions = True, check_partials=False): self.options = options self.nmodes = nmodes self.check_partials = check_partials + self.include_tacs_functions = include_tacs_functions - # api level method for all builders - def init_solver(self, comm): - - # save the comm bec. we will use this to control parallel output sizes + def initialize(self, comm): self.comm = comm solver_dict={} - # only root proc will do useful stuff if comm.rank == 0: mesh = TACS.MeshLoader(MPI.COMM_SELF) mesh.scanBDFFile(self.options['mesh_file']) ndof, ndv = self.options['add_elements'](mesh) - tacs = mesh.createTACS(ndof) + assembler = mesh.createTACS(ndof) - nnodes = int(tacs.createNodeVec().getArray().size / 3) + self.number_of_nodes = int(assembler.createNodeVec().getArray().size / 3) - mat = tacs.createFEMat() + mat = assembler.createFEMat() pc = TACS.Pc(mat) subspace = 100 restarts = 2 gmres = TACS.KSM(mat, pc, subspace, restarts) - solver_dict['nnodes'] = nnodes - solver_dict['get_funcs'] = self.options['get_funcs'] - - # check if the user provided a load function - if 'load_function' in self.options: - solver_dict['load_function'] = self.options['load_function'] - - # put the rest of the stuff in a tuple solver_objects = [mat, pc, gmres, solver_dict] - # on other procs, we just place dummy data else: - solver_dict['nnodes'] = 0 - solver_objects = np.zeros((0)) - tacs = None + self.number_of_nodes = 0 + solver_objects = [] + assembler = None ndv = 0 ndof = 0 - # ndv and ndof should be same on all procs - ndv = comm.bcast(ndv, root=0) - ndof = comm.bcast(ndof, root=0) - solver_dict['ndv'] = ndv - solver_dict['ndof'] = ndof + self.ndv = comm.bcast(ndv, root=0) + self.ndof = comm.bcast(ndof, root=0) + if 'f5_writer' in self.options.keys(): + solver_objects[3]['f5_writer'] = self.options['f5_writer'] + solver_objects[3]['ndv'] = self.ndv + solver_objects[3]['get_funcs'] = self.options['get_funcs'] - # set the attributes on all procs - self.solver_dict=solver_dict - self.solver = tacs + self.tacs_assembler = assembler self.solver_objects = solver_objects - # api level method for all builders - def get_solver(self): - return self.solver - - # api level method for all builders - def get_element(self, **kwargs): - return ModalGroup(solver=self.solver, solver_objects=self.solver_objects, - nmodes=self.nmodes, - nnodes=self.solver_dict['nnodes'], - ndof=self.solver_dict['ndof'], - check_partials=self.check_partials, **kwargs) - - def get_mesh_element(self): - return ModalDecomp(struct_solver=self.solver, - ndv=self.solver_dict['ndv'], + def get_coupling_group_subsystem(self): + return ModalGroup(nmodes=self.nmodes, + nnodes=self.number_of_nodes, + ndof=self.ndof, + check_partials=self.check_partials) + + def get_mesh_coordinate_subsystem(self): + return ModalDecomp(tacs_assembler=self.tacs_assembler, + ndv=self.ndv, nmodes=self.nmodes) - def get_mesh_connections(self): - return { - # because we dont have a solver or funcs key, - # mphys just assume that these will be connected - # to the solver. - 'modal_stiffness': 'modal_stiffness', - 'mode_shape': 'mode_shape', - } + def get_post_coupling_subsystem(self): + return TACSFuncsGroup( + tacs_assembler=self.tacs_assembler, + solver_objects=self.solver_objects, + check_partials=True) # always want the function component to load the state in def get_ndof(self): - return self.solver_dict['ndof'] + return self.ndof - def get_nnodes(self): - return self.solver_dict['nnodes'] + def get_number_of_nodes(self): + return self.number_of_nodes def get_ndv(self): - return self.solver_dict['ndv'] + return self.ndv diff --git a/mphys/mphys_rlt.py b/mphys/mphys_rlt.py index f02d7b2c..34cb1f2c 100644 --- a/mphys/mphys_rlt.py +++ b/mphys/mphys_rlt.py @@ -25,8 +25,6 @@ def initialize(self): # Flag used to prevent warning for fwd derivative d(u_a)/d(x_a0) self.options.declare("check_partials", default=False) - self.options["distributed"] = True - def setup(self): # get the inputs @@ -62,19 +60,31 @@ def setup(self): # Inputs self.add_input( - "x_struct0", shape_by_conn=True, desc="initial structural node coordinates", tags=["mphys_coordinates"] + "x_struct0", + distributed=True, + shape_by_conn=True, + desc="initial structural node coordinates", + tags=["mphys_coordinates"], ) self.add_input( "x_aero0", + distributed=True, shape_by_conn=True, desc="Initial aerodynamic surface node coordinates", tags=["mphys_coordinates"], ) - self.add_input("u_struct", shape_by_conn=True, desc="Structural node displacements", tags=["mphys_coupling"]) + self.add_input( + "u_struct", + distributed=True, + shape_by_conn=True, + desc="Structural node displacements", + tags=["mphys_coupling"], + ) # Outputs self.add_output( "u_aero", + distributed=True, shape=total_dof_aero, val=np.zeros(total_dof_aero), desc="Aerodynamic surface displacements", @@ -162,8 +172,6 @@ def initialize(self): # Flag used to prevent warning for fwd derivative d(u_a)/d(x_a0) self.options.declare("check_partials", default=True) - self.options["distributed"] = True - # Set everything we need to None before setup self.transfer = None self.tacs = None @@ -207,19 +215,38 @@ def setup(self): # Inputs self.add_input( - "x_struct0", shape_by_conn=True, desc="initial structural node coordinates", tags=["mphys_coordinates"] + "x_struct0", + distributed=True, + shape_by_conn=True, + desc="initial structural node coordinates", + tags=["mphys_coordinates"], ) self.add_input( "x_aero0", + distributed=True, shape_by_conn=True, desc="Initial aerodynamic surface node coordinates", tags=["mphys_coordinates"], ) - self.add_input("u_struct", shape_by_conn=True, desc="Structural node displacements", tags=["mphys_coupling"]) - self.add_input("f_aero", shape_by_conn=True, desc="Aerodynamic force vector", tags=["mphys_coupling"]) + self.add_input( + "u_struct", + distributed=True, + shape_by_conn=True, + desc="Structural node displacements", + tags=["mphys_coupling"], + ) + self.add_input( + "f_aero", distributed=True, shape_by_conn=True, desc="Aerodynamic force vector", tags=["mphys_coupling"] + ) # Outputs - self.add_output("f_struct", shape=total_dof_struct, desc="structural force vector", tags=["mphys_coupling"]) + self.add_output( + "f_struct", + distributed=True, + shape=total_dof_struct, + desc="structural force vector", + tags=["mphys_coupling"], + ) # TODO disable for now for the modal solver stuff. # Partials diff --git a/mphys/mphys_tacs.py b/mphys/mphys_tacs.py index 092dcf81..366f94ca 100644 --- a/mphys/mphys_tacs.py +++ b/mphys/mphys_tacs.py @@ -10,14 +10,13 @@ class TacsMesh(om.IndepVarComp): """ def initialize(self): self.options.declare('tacs_assembler', default = None, desc='the tacs object itself', recordable=False) - self.options['distributed'] = True def setup(self): tacs_assembler = self.options['tacs_assembler'] xpts = tacs_assembler.createNodeVec() x = xpts.getArray() tacs_assembler.getNodes(xpts) - self.add_output('x_struct0', val=x, shape=x.size, desc='structural node coordinates', tags=['mphys_coordinates']) + self.add_output('x_struct0', distributed=True, val=x, shape=x.size, desc='structural node coordinates', tags=['mphys_coordinates']) class TacsSolver(om.ImplicitComponent): """ @@ -33,8 +32,6 @@ def initialize(self): self.options.declare('struct_objects', recordable=False) self.options.declare('check_partials') - self.options['distributed'] = True - self.tacs_assembler = None self.pc = None @@ -84,13 +81,13 @@ def setup(self): self.ndof = int(state_size/(node_size/3)) # inputs - self.add_input('dv_struct', shape=ndv, desc='tacs design variables', tags=['mphys_input']) - self.add_input('x_struct0', shape_by_conn=True, desc='structural node coordinates',tags=['mphys_coordinates']) - self.add_input('f_struct', shape_by_conn=True, desc='structural load vector', tags=['mphys_coupling']) + self.add_input('dv_struct', distributed=False, shape=ndv, desc='tacs design variables', tags=['mphys_input']) + self.add_input('x_struct0', distributed=True, shape_by_conn=True, desc='structural node coordinates',tags=['mphys_coordinates']) + self.add_input('f_struct', distributed=True, shape_by_conn=True, desc='structural load vector', tags=['mphys_coupling']) # outputs # its important that we set this to zero since this displacement value is used for the first iteration of the aero - self.add_output('u_struct', shape=state_size, val = np.zeros(state_size),desc='structural state vector', tags=['mphys_coupling']) + self.add_output('u_struct', distributed=True, shape=state_size, val = np.zeros(state_size),desc='structural state vector', tags=['mphys_coupling']) # partials #self.declare_partials('u_struct',['dv_struct','x_struct0','f_struct']) @@ -318,8 +315,6 @@ def initialize(self): self.options.declare('struct_objects') self.options.declare('check_partials') - self.options['distributed'] = True - self.tacs_assembler = None self.pc = None @@ -363,11 +358,11 @@ def setup(self): # inputs # self.add_input('dv_struct', shape=ndv , desc='tacs design variables') - self.add_input('x_struct0', shape_by_conn=True, desc='structural node coordinates', tags=['mphys_coordinates']) - self.add_input('q_conduct', shape_by_conn=True, desc='structural load vector', tags=['mphys_coupling']) + self.add_input('x_struct0', distributed=True, shape_by_conn=True, desc='structural node coordinates', tags=['mphys_coordinates']) + self.add_input('q_conduct', distributed=True, shape_by_conn=True, desc='structural load vector', tags=['mphys_coupling']) # outputs - self.add_output('T_conduct', shape=surface_nodes.size//3, val = np.ones(surface_nodes.size//3)*300,desc='temperature vector', tags=['mphys_coupling']) + self.add_output('T_conduct', distributed=True, shape=surface_nodes.size//3, val = np.ones(surface_nodes.size//3)*300,desc='temperature vector', tags=['mphys_coupling']) # partials #self.declare_partials('u_struct',['dv_struct','x_struct0','f_struct']) @@ -440,8 +435,6 @@ def initialize(self): self.options.declare('struct_objects', recordable=False) self.options.declare('check_partials') - self.options['distributed'] = True - self.ans = None self.tacs_assembler = None @@ -475,9 +468,9 @@ def setup(self): # OpenMDAO part of setup # TODO move the dv_struct to an external call where we add the DVs - self.add_input('dv_struct', shape = ndv, desc='tacs design variables', tags=['mphys_input']) - self.add_input('x_struct0', shape_by_conn=True, desc='structural node coordinates',tags=['mphys_coordinates']) - self.add_input('u_struct', shape_by_conn=True, desc='structural state vector', tags=['mphys_coupling']) + self.add_input('dv_struct', distributed=False, shape = ndv, desc='tacs design variables', tags=['mphys_input']) + self.add_input('x_struct0', distributed=True, shape_by_conn=True, desc='structural node coordinates',tags=['mphys_coordinates']) + self.add_input('u_struct', distributed=True, shape_by_conn=True, desc='structural state vector', tags=['mphys_coupling']) # Remove the mass function from the func list if it is there # since it is not dependent on the structural state @@ -488,7 +481,7 @@ def setup(self): self.func_list = func_no_mass if len(self.func_list) > 0: - self.add_output('func_struct', shape=len(self.func_list), desc='structural function values', tags=['mphys_result']) + self.add_output('func_struct', distributed=False, shape=len(self.func_list), desc='structural function values', tags=['mphys_result']) # declare the partials #self.declare_partials('f_struct',['dv_struct','x_struct0','u_struct']) @@ -537,37 +530,36 @@ def compute_jacvec_product(self,inputs, d_inputs, d_outputs, mode): if not self.check_partials: raise ValueError('TACS forward mode requested but not implemented') if mode == 'rev': + # always update internal because same tacs object could be used by multiple scenarios + # and we need to load this scenario's state back into TACS before doing derivatives self._update_internal(inputs) - #if 'func_struct' in d_outputs: - # Next few lines temporary until func_struct can be declared a serial output - if True: - if 'func_struct' in d_outputs: - proc_contribution = d_outputs['func_struct'][:] - else: - proc_contribution = np.zeros(len(self.func_list)) - d_func = self.comm.allreduce(proc_contribution) / self.comm.size + if 'func_struct' in d_outputs: + proc_contribution = d_outputs['func_struct'][:] + else: # not sure why OM would call this method without func_struct, but here for safety + proc_contribution = np.zeros(len(self.func_list)) + d_func = self.comm.allreduce(proc_contribution) / self.comm.size - for ifunc, func in enumerate(self.func_list): - self.tacs_assembler.evalFunctions([func]) - if 'dv_struct' in d_inputs: - dvsens = np.zeros(d_inputs['dv_struct'].size,dtype=TACS.dtype) - self.tacs_assembler.evalDVSens(func, dvsens) - d_inputs['dv_struct'][:] += np.array(dvsens,dtype=float) * d_func[ifunc] + for ifunc, func in enumerate(self.func_list): + self.tacs_assembler.evalFunctions([func]) + if 'dv_struct' in d_inputs: + dvsens = np.zeros(d_inputs['dv_struct'].size,dtype=TACS.dtype) + self.tacs_assembler.evalDVSens(func, dvsens) + d_inputs['dv_struct'][:] += np.array(dvsens,dtype=float) * d_func[ifunc] - if 'x_struct0' in d_inputs: - xpt_sens = self.xpt_sens - xpt_sens_array = xpt_sens.getArray() - self.tacs_assembler.evalXptSens(func, xpt_sens) + if 'x_struct0' in d_inputs: + xpt_sens = self.xpt_sens + xpt_sens_array = xpt_sens.getArray() + self.tacs_assembler.evalXptSens(func, xpt_sens) - d_inputs['x_struct0'][:] += np.array(xpt_sens_array,dtype=float) * d_func[ifunc] + d_inputs['x_struct0'][:] += np.array(xpt_sens_array,dtype=float) * d_func[ifunc] - if 'u_struct' in d_inputs: - prod = self.tacs_assembler.createVec() - self.tacs_assembler.evalSVSens(func,prod) - prod_array = prod.getArray() + if 'u_struct' in d_inputs: + prod = self.tacs_assembler.createVec() + self.tacs_assembler.evalSVSens(func,prod) + prod_array = prod.getArray() - d_inputs['u_struct'][:] += np.array(prod_array,dtype=float) * d_func[ifunc] + d_inputs['u_struct'][:] += np.array(prod_array,dtype=float) * d_func[ifunc] class TacsMass(om.ExplicitComponent): """ @@ -578,8 +570,6 @@ def initialize(self): self.options.declare('struct_objects', recordable=False) self.options.declare('check_partials') - self.options['distributed'] = True - self.ans = None self.tacs_assembler = None @@ -604,10 +594,10 @@ def setup(self): self.xpt_sens = tacs_assembler.createNodeVec() # OpenMDAO part of setup - self.add_input('dv_struct', shape=ndv, desc='tacs design variables', tags=['mphys_input']) - self.add_input('x_struct0', shape_by_conn=True, desc='structural node coordinates', tags=['mphys_coordinates']) + self.add_input('dv_struct', distributed=False, shape=ndv, desc='tacs design variables', tags=['mphys_input']) + self.add_input('x_struct0', distributed=True, shape_by_conn=True, desc='structural node coordinates', tags=['mphys_coordinates']) - self.add_output('mass', 0.0, desc = 'structural mass', tags=['mphys_result']) + self.add_output('mass', val=0.0, distributed=False, desc = 'structural mass', tags=['mphys_result']) #self.declare_partials('mass',['dv_struct','x_struct0']) def _update_internal(self,inputs): @@ -656,8 +646,6 @@ def initialize(self): self.options.declare('load_function', default = None, desc='function that prescribes the loads', recordable=False) self.options.declare('tacs_assembler', recordable=False) - self.options['distributed'] = True - self.ndof = 0 def setup(self): @@ -676,8 +664,8 @@ def setup(self): self.ndof = int(state_size / ( node_size / 3 )) # OpenMDAO setup - self.add_input('x_struct0', shape_by_conn=True, desc='structural node coordinates', tags=['mphys_coordinates']) - self.add_output('f_struct', shape=state_size, desc='structural load', tags=['mphys_coupling']) + self.add_input('x_struct0', distributed=True, shape_by_conn=True, desc='structural node coordinates', tags=['mphys_coordinates']) + self.add_output('f_struct', distributed=True, shape=state_size, desc='structural load', tags=['mphys_coupling']) #self.declare_partials('f_struct','x_struct0') @@ -803,11 +791,10 @@ def initialize(self, comm): self.tacs_assembler = tacs_assembler self.solver_objects = [mat, pc, gmres, solver_dict] - def get_coupling_group_subsystem(self, **kwargs): + def get_coupling_group_subsystem(self): return TacsGroup(tacs_assembler=self.tacs_assembler, solver_objects=self.solver_objects, - check_partials=self.check_partials, - **kwargs) + check_partials=self.check_partials) def get_mesh_coordinate_subsystem(self): return TacsMesh(tacs_assembler=self.tacs_assembler) diff --git a/mphys/mphys_vlm.py b/mphys/mphys_vlm.py index 3c58d409..da3d95b5 100644 --- a/mphys/mphys_vlm.py +++ b/mphys/mphys_vlm.py @@ -1,7 +1,7 @@ import numpy as np import openmdao.api as om -from mphys import Builder +from mphys import Builder, DistributedConverter, DistributedVariableDescription from vlm_solver import Vlm @@ -9,7 +9,7 @@ class VlmMesh(om.IndepVarComp): def initialize(self): self.options.declare('x_aero0') def setup(self): - self.add_output('x_aero0', val=self.options['x_aero0'], tags=['mphys_coordinates']) + self.add_output('x_aero0', val=self.options['x_aero0'])#, tags=['mphys_coordinates']) class VlmSolver(om.ImplicitComponent): def initialize(self): @@ -208,33 +208,73 @@ def compute_partials(self,inputs,partials): partials['C_D','Cp'] = self.vlm.CD_Cp -class VlmGroup(om.Group): +class VlmMeshGroup(om.Group): + def initialize(self): + self.options.declare('x_aero0') + + def setup(self): + x_aero0 = self.options['x_aero0'] + self.add_subsystem('vlm_mesh', VlmMesh(x_aero0 = x_aero0)) + + vars = [DistributedVariableDescription(name='x_aero0', + shape=(x_aero0.size), + tags =['mphys_coordinates'])] + + self.add_subsystem('distributor',DistributedConverter(distributed_outputs=vars), + promotes_outputs=[var.name for var in vars]) + for var in vars: + self.connect(f'vlm_mesh.{var.name}', f'distributor.{var.name}_serial') + + +class VlmSolverGroup(om.Group): def initialize(self): self.options.declare('solver') + self.options.declare('n_aero') self.options.declare('complex_step', default=False) def setup(self): self.solver = self.options['solver'] complex_step = self.options['complex_step'] + n_aero = self.options['n_aero'] + + in_vars = [DistributedVariableDescription(name='x_aero', + shape=(n_aero*3), + tags =['mphys_coupling'])] + out_vars = [DistributedVariableDescription(name='f_aero', + shape=(n_aero*3), + tags =['mphys_coupling'])] + + self.add_subsystem('collector',DistributedConverter(distributed_inputs=in_vars), + promotes_inputs=[var.name for var in in_vars]) self.add_subsystem('aero_solver', VlmSolver( solver=self.solver, complex_step = complex_step), - promotes_inputs=['aoa','mach','x_aero'], + promotes_inputs=['aoa','mach'], promotes_outputs=['Cp'] ) self.add_subsystem('aero_forces', VlmForces( solver=self.solver, complex_step=complex_step), - promotes_inputs=['q_inf','Cp','x_aero'], - promotes_outputs=['f_aero'] + promotes_inputs=['q_inf','Cp'], ) self.add_subsystem('aero_coefficients', VlmCoefficients( solver=self.solver, complex_step=complex_step), - promotes_inputs=['mach','vel','nu','Cp','x_aero'], + promotes_inputs=['mach','vel','nu','Cp'], promotes_outputs=['C_L','C_D'] ) + self.add_subsystem('distributor',DistributedConverter(distributed_outputs=out_vars), + promotes_outputs=[var.name for var in out_vars]) + + connection_dest = ['aero_solver', 'aero_forces', 'aero_coefficients'] + for var in in_vars: + for dest in connection_dest: + self.connect(f'collector.{var.name}_serial', f'{dest}.{var.name}') + + for var in out_vars: + self.connect(f'aero_forces.{var.name}', f'distributor.{var.name}_serial') + class VlmBuilder(Builder): def __init__(self, meshfile, compute_traction=False, complex_step=False): @@ -246,15 +286,17 @@ def initialize(self, comm): self.solver = Vlm(compute_traction=self.compute_traction) self.solver.read_mesh(self.meshfile) self.x_aero0 = self.solver.xa + self.n_aero = self.x_aero0.size // 3 + self.comm = comm def get_mesh_coordinate_subsystem(self): - return VlmMesh(x_aero0=self.x_aero0) + return VlmMeshGroup(x_aero0=self.x_aero0) def get_coupling_group_subsystem(self): - return VlmGroup(solver=self.solver, complex_step=self.complex_step) + return VlmSolverGroup(solver=self.solver, n_aero= self.n_aero, complex_step=self.complex_step) def get_number_of_nodes(self): - return self.x_aero0.size // 3 + return self.n_aero if self.comm.Get_rank() == 0 else 0 def get_ndof(self): - return self.x_aero0.size // 3 + return 3 diff --git a/tests/integration_tests/test_adflow_derivs.py b/tests/integration_tests/test_adflow_derivs.py index 7c058660..4942ca6b 100644 --- a/tests/integration_tests/test_adflow_derivs.py +++ b/tests/integration_tests/test_adflow_derivs.py @@ -84,10 +84,11 @@ def setup(self): # ivc for aero DVs self.add_subsystem("dv_aero", om.IndepVarComp(), promotes=["*"]) # ivc for mesh coordinate dvs. we need a separate distributed component for the shape by conn stuff to work - dv_coords = self.add_subsystem("dv_coords", om.IndepVarComp(distributed=False), promotes=["*"]) + dv_coords = self.add_subsystem("dv_coords", om.IndepVarComp(), promotes=["*"]) + # TODO this only works for a serial run. must be updated for parallel x_a = adflow_builder.solver.getSurfaceCoordinates(groupName="allWalls").flatten() - dv_coords.add_output("x_aero", val=x_a) + dv_coords.add_output("x_aero", val=x_a, distributed=False) # normally we would have a mesh comp but we just do the parallel ivc for the test. self.mphys_add_scenario("cruise", ScenarioAerodynamic(aero_builder=adflow_builder)) diff --git a/tests/integration_tests/test_tacs_derivs.py b/tests/integration_tests/test_tacs_derivs.py index 5ce2ab2d..42882f9e 100644 --- a/tests/integration_tests/test_tacs_derivs.py +++ b/tests/integration_tests/test_tacs_derivs.py @@ -67,10 +67,6 @@ def get_funcs(tacs): ms = functions.StructuralMass(tacs) return [ks, ms] -class Forces(om.IndepVarComp): - def initialize(self): - self.options['distributed'] = True - class Top(Multipoint): def setup(self): @@ -88,8 +84,8 @@ def setup(self): dvs = self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*']) dvs.add_output('dv_struct', np.array(ndv_struct * [0.01])) - forces = self.add_subsystem('forces', Forces(), promotes=['*']) - forces.add_output('f_struct', np.ones(f_size)) + forces = self.add_subsystem('forces', om.IndepVarComp(), promotes=['*']) + forces.add_output('f_struct', val=np.ones(f_size), distributed=True) self.add_subsystem('mesh', struct_builder.get_mesh_coordinate_subsystem()) self.mphys_add_scenario('analysis', ScenarioStructural(struct_builder=struct_builder)) diff --git a/tests/unit_tests/test_distributed_converter.py b/tests/unit_tests/test_distributed_converter.py new file mode 100644 index 00000000..4053ad3c --- /dev/null +++ b/tests/unit_tests/test_distributed_converter.py @@ -0,0 +1,70 @@ +import unittest +from mpi4py import MPI +import numpy as np +import openmdao.api as om + +from mphys import DistributedConverter, DistributedVariableDescription + +from common_methods import CommonMethods +from openmdao.utils.assert_utils import assert_near_equal + + +class TestDistributedConverter(unittest.TestCase): + N_PROCS = 1 #TODO should be 2 or more but there is a bug in OM currently + + def setUp(self): + self.common = CommonMethods() + self.prob = om.Problem() + + vars_in = [ + DistributedVariableDescription('in1', shape=(4, 4), tags=['mphys_coupling']), + DistributedVariableDescription('in2', shape=(10), tags=['mphys_coordinates']) + ] + + vars_out = [ + DistributedVariableDescription('out1', shape=(5), tags=['mphys_coupling']), + DistributedVariableDescription('out2', shape=(15), tags=['mphys_result']) + ] + + inputs = self.prob.model.add_subsystem('inputs', om.IndepVarComp()) + + in1_shape = (4, 4) if MPI.COMM_WORLD.Get_rank() == 0 else 0 + in2_shape = 10 if MPI.COMM_WORLD.Get_rank() == 0 else 0 + inputs.add_output('in1', val=np.ones(in1_shape, dtype=float), distributed=True) + inputs.add_output('in2', val=np.arange(in2_shape, dtype=float), distributed=True) + + out1_shape = 5 + out2_shape = 15 + inputs.add_output('out1', val=np.ones(out1_shape, dtype=float), distributed=False) + inputs.add_output('out2', val=np.arange(out2_shape, dtype=float), distributed=False) + + self.prob.model.add_subsystem('converter', DistributedConverter( + distributed_inputs=vars_in, distributed_outputs=vars_out)) + + for var in ['in1', 'in2']: + self.prob.model.connect(f'inputs.{var}', f'converter.{var}') + for var in ['out1', 'out2']: + self.prob.model.connect(f'inputs.{var}', f'converter.{var}_serial') + + self.prob.setup(force_alloc_complex=True) + + def test_run_model(self): + self.common.test_run_model(self, write_n2=False) + + def test_check_partials(self): + partials = self.prob.check_partials(compact_print=True, method='cs') + tol = 1e-9 + for in_var in ['in1', 'in2']: + rel_error = partials['converter'][(f'{in_var}_serial', in_var)]['rel error'] + assert_near_equal(rel_error.reverse, 0.0, tolerance=tol) + assert_near_equal(rel_error.forward, 0.0, tolerance=tol) + assert_near_equal(rel_error.forward_reverse, 0.0, tolerance=tol) + for out_var in ['out1', 'out2']: + rel_error = partials['converter'][(out_var, f'{out_var}_serial')]['rel error'] + assert_near_equal(rel_error.reverse, 0.0, tolerance=tol) + assert_near_equal(rel_error.forward, 0.0, tolerance=tol) + assert_near_equal(rel_error.forward_reverse, 0.0, tolerance=tol) + + +if __name__ == '__main__': + unittest.main()