From f67e476ec14429bbeb00cf1560f89f506f9ec449 Mon Sep 17 00:00:00 2001
From: Anil Yildirim <anily@umich.edu>
Date: Mon, 3 May 2021 12:07:45 +0300
Subject: [PATCH 01/18] updated adflow wrapper to latest distributed api

---
 mphys/mphys_adflow.py | 73 +++++++++++++++++++------------------------
 1 file changed, 33 insertions(+), 40 deletions(-)

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

From 6d8d47d3bf23b013fcdbfe70c9311656c5492cdc Mon Sep 17 00:00:00 2001
From: Anil Yildirim <anily@umich.edu>
Date: Mon, 3 May 2021 14:42:56 +0300
Subject: [PATCH 02/18] updated dvgeo

---
 mphys/mphys_dvgeo.py | 33 ++++++++++++++++-----------------
 1 file changed, 16 insertions(+), 17 deletions(-)

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

From 4029835332dbcaa154f1619d5714ab84f209461c Mon Sep 17 00:00:00 2001
From: Anil Yildirim <anily@umich.edu>
Date: Mon, 3 May 2021 15:16:43 +0300
Subject: [PATCH 03/18] updated geo_disp, meld, rlt, tacs, and meld thermal

---
 mphys/geo_disp.py          |  4 +++-
 mphys/mphys_meld.py        | 13 +++++++----
 mphys/mphys_meldthermal.py | 19 +++++++--------
 mphys/mphys_rlt.py         | 47 ++++++++++++++++++++++++++++++--------
 mphys/mphys_tacs.py        | 45 ++++++++++++++----------------------
 5 files changed, 74 insertions(+), 54 deletions(-)

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/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_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..44819e58 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'])
@@ -578,8 +571,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 +595,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 +647,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 +665,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')
 

From bd07684533342f9e3c86207afc4ed89f257cd884 Mon Sep 17 00:00:00 2001
From: Anil Yildirim <anily@umich.edu>
Date: Mon, 3 May 2021 15:26:23 +0300
Subject: [PATCH 04/18] updated distributed declarations in tests

---
 tests/integration_tests/test_adflow_derivs.py | 5 +++--
 tests/integration_tests/test_tacs_derivs.py   | 8 ++------
 2 files changed, 5 insertions(+), 8 deletions(-)

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

From 8fab08d81f38657d93efcdf49cbc4391951153d6 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Tue, 22 Jun 2021 18:32:46 -0400
Subject: [PATCH 05/18] Add distributed converter with test. Test fails in
 parallel due to OM bug currently

---
 mphys/distributed_converter.py                | 72 +++++++++++++++++++
 .../unit_tests/test_distributed_converter.py  | 70 ++++++++++++++++++
 2 files changed, 142 insertions(+)
 create mode 100644 mphys/distributed_converter.py
 create mode 100644 tests/unit_tests/test_distributed_converter.py

diff --git a/mphys/distributed_converter.py b/mphys/distributed_converter.py
new file mode 100644
index 00000000..2ca422dc
--- /dev/null
+++ b/mphys/distributed_converter.py
@@ -0,0 +1,72 @@
+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.
+    """
+
+    def initialize(self):
+        self.options.declare('distributed_inputs', default=[])
+        self.options.declare('distributed_outputs', default=[])
+
+    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/tests/unit_tests/test_distributed_converter.py b/tests/unit_tests/test_distributed_converter.py
new file mode 100644
index 00000000..409985ae
--- /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 = 2
+
+    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()

From 0fa898b1bb8eceaf8c279ddc63bddda199f144d2 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Tue, 22 Jun 2021 18:34:32 -0400
Subject: [PATCH 06/18] Turn on heuristic imports for linting in examples

---
 .vscode/settings.json | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

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
 }

From b1b1f654bf5ffb2bbf0090a4efc8ca8116b9739c Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Tue, 22 Jun 2021 18:54:46 -0400
Subject: [PATCH 07/18] Add distributed converters to VLM solver subsystems

---
 mphys/mphys_vlm.py | 63 ++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 52 insertions(+), 11 deletions(-)

diff --git a/mphys/mphys_vlm.py b/mphys/mphys_vlm.py
index 3c58d409..f18e6fec 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,16 @@ 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
 
     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
 
     def get_ndof(self):
-        return self.x_aero0.size // 3
+        return 3

From 0eb22507d7d12fad688be94ebc1269925700029b Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 08:30:10 -0400
Subject: [PATCH 08/18] Add dimensions to angle inputs of integrated force
 calculator

---
 examples/aero_only/agard/run_fun3d.py | 13 +++++++------
 mphys/integrated_forces.py            | 17 +++++++----------
 2 files changed, 14 insertions(+), 16 deletions(-)

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

From 69c60bd922ed1c7d0439093d68a9032bd482bd08 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 08:45:54 -0400
Subject: [PATCH 09/18] Add distributed converter to __init__.py

---
 mphys/__init__.py | 1 +
 1 file changed, 1 insertion(+)

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

From e4e0e7284ca4a9a51850428abf5abc712bf2a279 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 08:46:40 -0400
Subject: [PATCH 10/18] VLM return the correct number of nodes for parallel
 runs

---
 mphys/mphys_vlm.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/mphys/mphys_vlm.py b/mphys/mphys_vlm.py
index f18e6fec..da3d95b5 100644
--- a/mphys/mphys_vlm.py
+++ b/mphys/mphys_vlm.py
@@ -287,6 +287,7 @@ def initialize(self, comm):
         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 VlmMeshGroup(x_aero0=self.x_aero0)
@@ -295,7 +296,7 @@ def get_coupling_group_subsystem(self):
         return VlmSolverGroup(solver=self.solver, n_aero= self.n_aero, complex_step=self.complex_step)
 
     def get_number_of_nodes(self):
-        return self.n_aero
+        return self.n_aero if self.comm.Get_rank() == 0 else 0
 
     def get_ndof(self):
         return 3

From fef5fb0be501014c9dfa466c7a62cf56194c895f Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 08:52:16 -0400
Subject: [PATCH 11/18] Add some description of the distributed converter

---
 mphys/distributed_converter.py | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/mphys/distributed_converter.py b/mphys/distributed_converter.py
index 2ca422dc..304b6085 100644
--- a/mphys/distributed_converter.py
+++ b/mphys/distributed_converter.py
@@ -18,11 +18,22 @@ 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=[])
-        self.options.declare('distributed_outputs', default=[])
+        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']:

From 4a81b6238df9ad7c79d80d591dba47cb097b4c10 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 08:58:54 -0400
Subject: [PATCH 12/18] Make the converter test 1 proc until OM bug is fixed

---
 tests/unit_tests/test_distributed_converter.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/unit_tests/test_distributed_converter.py b/tests/unit_tests/test_distributed_converter.py
index 409985ae..28f486d9 100644
--- a/tests/unit_tests/test_distributed_converter.py
+++ b/tests/unit_tests/test_distributed_converter.py
@@ -9,7 +9,7 @@
 from openmdao.utils.assert_utils import assert_near_equal
 
 class TestDistributedConverter(unittest.TestCase):
-    N_PROCS = 2
+    N_PROCS = 1 #TODO should be 2 or more but there is a bug in OM currently
 
     def setUp(self):
         self.common = CommonMethods()

From dfaa487401f6717c2097eec7623df29c078e2176 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 09:00:21 -0400
Subject: [PATCH 13/18] formatting

---
 .../unit_tests/test_distributed_converter.py  | 26 +++++++++----------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/tests/unit_tests/test_distributed_converter.py b/tests/unit_tests/test_distributed_converter.py
index 28f486d9..4053ad3c 100644
--- a/tests/unit_tests/test_distributed_converter.py
+++ b/tests/unit_tests/test_distributed_converter.py
@@ -8,8 +8,9 @@
 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
+    N_PROCS = 1  #TODO should be 2 or more but there is a bug in OM currently
 
     def setUp(self):
         self.common = CommonMethods()
@@ -52,18 +53,17 @@ def test_run_model(self):
 
     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)
-
+        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__':

From bd3fbbecac6d7827f9787c71c72eb6dbf57d764e Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 12:50:27 -0400
Subject: [PATCH 14/18] Some updates to the modal solver

---
 .../vlm_meld_tacs/run_analysis_modal.py       |  90 +++++
 mphys/mphys_modal_solver.py                   | 310 +++++++-----------
 2 files changed, 207 insertions(+), 193 deletions(-)
 create mode 100644 examples/aerostructural/mach_tutorial_wing/vlm_meld_tacs/run_analysis_modal.py

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..534eaccf
--- /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']:
+                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.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/mphys_modal_solver.py b/mphys/mphys_modal_solver.py
index 534f21ef..51b32695 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')
-
-        if self.comm.rank == 0:
-            output_shape = self.nmodes
-        else:
-            output_shape = 0
+        self.add_output('mode_shape', shape=self.nmodes*self.state_size,
+                        distributed=True,
+                        desc='structural mode shapes',
+                        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')
+        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('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
+        nmodes = self.options['nmodes']
 
-        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')
+        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]
 
@@ -206,24 +177,10 @@ def setup(self):
             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('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')
 
-        self.add_input('z',shape=input_shape, src_indices=src_indices, 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 +189,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 +237,78 @@ 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)
+        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=self.check_partials
+        )
 
     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

From cd5a3bee105a2288f7bf9952b0ece285647233f7 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 14:06:16 -0400
Subject: [PATCH 15/18] TACS functions always update state for modal

---
 mphys/mphys_modal_solver.py | 13 +++----------
 1 file changed, 3 insertions(+), 10 deletions(-)

diff --git a/mphys/mphys_modal_solver.py b/mphys/mphys_modal_solver.py
index 51b32695..d7b7173a 100644
--- a/mphys/mphys_modal_solver.py
+++ b/mphys/mphys_modal_solver.py
@@ -169,14 +169,6 @@ 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_by_conn=True, distributed=True, desc='structural mode shapes', tags=['mphys_coupling'])
         self.add_input('z', shape_by_conn=True, desc='modal displacement')
 
@@ -280,6 +272,8 @@ def initialize(self, comm):
 
         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']
 
@@ -301,8 +295,7 @@ def get_post_coupling_subsystem(self):
         return TACSFuncsGroup(
             tacs_assembler=self.tacs_assembler,
             solver_objects=self.solver_objects,
-            check_partials=self.check_partials
-        )
+            check_partials=True) # always want the function component to load the state in
 
     def get_ndof(self):
         return self.ndof

From 6b6eab03cfd934dc2b35e565a4745c588ff1152f Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 14:06:35 -0400
Subject: [PATCH 16/18] Connect dv_struct for modal functions

---
 .../mach_tutorial_wing/vlm_meld_tacs/run_analysis_modal.py    | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

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
index 534eaccf..bafe633a 100644
--- 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
@@ -70,7 +70,7 @@ def setup(self):
 
             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']:
+            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])
 
@@ -81,7 +81,7 @@ def setup(self):
 prob.model = Top()
 
 prob.setup()
-om.n2(prob, show_browser=False, outfile='mphys_as_vlm.html')
+om.n2(prob, show_browser=False, outfile='mphys_as_vlm_modal.html')
 
 prob.run_model()
 

From da863f8b1741de6094b404dd4893b24a0bd9ee72 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 14:07:16 -0400
Subject: [PATCH 17/18] Remove unnecessary "if True:"

---
 mphys/mphys_tacs.py | 52 ++++++++++++++++++++++-----------------------
 1 file changed, 25 insertions(+), 27 deletions(-)

diff --git a/mphys/mphys_tacs.py b/mphys/mphys_tacs.py
index 44819e58..366f94ca 100644
--- a/mphys/mphys_tacs.py
+++ b/mphys/mphys_tacs.py
@@ -530,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):
     """
@@ -792,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)

From 672e8980cb071d41b11e4b8b014f9614b7dd1359 Mon Sep 17 00:00:00 2001
From: Kevin Jacobson <kevin.e.jacobson@nasa.gov>
Date: Wed, 23 Jun 2021 14:09:44 -0400
Subject: [PATCH 18/18] Update fun3d mesh for distributed variables

---
 mphys/mphys_fun3d.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

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