Merge pull request #66 from anilyil/distributed_io

Updated variable declarations for the new distributed API in openmdao 3.9.0

.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
 }

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')

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'])

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

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'])

mphys/geo_disp.py

@@ -6,21 +6,23 @@ 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):
         nnodes = self.options['number_of_nodes']
         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'])
 

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)