Merge 750988b into c10d6d5

openmdao/core/group.py

@@ -1896,9 +1896,39 @@ def _solve_linear(self, vec_names, mode, rel_systems):
         rel_systems : set of str
             Set of names of relevant systems based on the current linear solve.
         """
-        vec_names = [v for v in vec_names if v in self._rel_vec_names]
+        if self._owns_approx_jac:
+            # No subsolves if we are approximating our jacobian. Instead, we behave like an
+            # ExplicitComponent and pass on the values in the derivatives vectors.
+            for vec_name in vec_names:
+                if vec_name in self._rel_vec_names:
+                    d_outputs = self._vectors['output'][vec_name]
+                    d_residuals = self._vectors['residual'][vec_name]
+
+                    if mode == 'fwd':
+                        if self._has_resid_scaling:
+                            with self._unscaled_context(outputs=[d_outputs],
+                                                        residuals=[d_residuals]):
+                                d_outputs.set_vec(d_residuals)
+                        else:
+                            d_outputs.set_vec(d_residuals)
+
+                        # ExplicitComponent jacobian defined with -1 on diagonal.
+                        d_outputs *= -1.0
 
-        self._linear_solver.solve(vec_names, mode, rel_systems)
+                    else:  # rev
+                        if self._has_resid_scaling:
+                            with self._unscaled_context(outputs=[d_outputs],
+                                                        residuals=[d_residuals]):
+                                d_residuals.set_vec(d_outputs)
+                        else:
+                            d_residuals.set_vec(d_outputs)
+
+                        # ExplicitComponent jacobian defined with -1 on diagonal.
+                        d_residuals *= -1.0
+
+        else:
+            vec_names = [v for v in vec_names if v in self._rel_vec_names]
+            self._linear_solver.solve(vec_names, mode, rel_systems)
 
     def _linearize(self, jac, sub_do_ln=True):
         """

openmdao/core/tests/test_approx_derivs.py

@@ -810,6 +810,86 @@ def test_opt_with_linear_constraint(self):
 
         assert_near_equal(p['circle.area'], np.pi, 1e-6)
 
+    def test_bug_subsolve(self):
+        # There was a bug where a group with an approximation was still performing a linear
+        # solve on its subsystems, which led to partials declared with 'val' corrupting the
+        # results.
+
+        class DistParab(om.ExplicitComponent):
+
+            def initialize(self):
+
+                self.options.declare('arr_size', types=int, default=10,
+                                     desc="Size of input and output vectors.")
+
+            def setup(self):
+                arr_size = self.options['arr_size']
+
+                self.add_input('x', val=np.ones(arr_size))
+                self.add_output('f_xy', val=np.ones(arr_size))
+
+                self.declare_partials('f_xy', 'x')
+
+            def compute(self, inputs, outputs):
+                x = inputs['x']
+                outputs['f_xy'] = x**2
+
+        class NonDistComp(om.ExplicitComponent):
+
+            def initialize(self):
+                self.options.declare('arr_size', types=int, default=10,
+                                     desc="Size of input and output vectors.")
+
+            def setup(self):
+                arr_size = self.options['arr_size']
+
+                self.add_input('f_xy', val=np.ones(arr_size))
+
+                self.add_output('g', val=np.ones(arr_size))
+
+                # Make this wrong to see if it shows up in the answer.
+                mat = np.array([7.0, 13, 27])
+
+                row_col = np.arange(arr_size)
+                self.declare_partials('g', ['f_xy'], rows=row_col, cols=row_col, val=mat)
+                #self.declare_partials('g', ['f_xy'])
+
+            def compute(self, inputs, outputs):
+                x = inputs['f_xy']
+                outputs['g'] = x * np.array([3.5, -1.0, 5.0])
+
+        size = 3
+
+        prob = om.Problem()
+        model = prob.model
+
+        ivc = om.IndepVarComp()
+        ivc.add_output('x', np.ones((size, )))
+
+        model.add_subsystem('p', ivc, promotes=['*'])
+        sub = model.add_subsystem('sub', om.Group(), promotes=['*'])
+
+        sub.add_subsystem("parab", DistParab(arr_size=size), promotes=['*'])
+        sub.add_subsystem("ndp", NonDistComp(arr_size=size), promotes=['*'])
+
+        model.add_design_var('x', lower=-50.0, upper=50.0)
+        model.add_constraint('g', lower=0.0)
+
+        sub.approx_totals(method='fd')
+
+        prob.setup()
+
+        prob.run_model()
+
+        of = ['sub.ndp.g']
+        totals = prob.driver._compute_totals(of=of, wrt=['p.x'], return_format='dict')
+        assert_near_equal(totals['sub.ndp.g']['p.x'], np.diag([7.0, -2.0, 10.0]), 1e-6)
+
+        totals = prob.check_totals()
+
+        for key, val in totals.items():
+            assert_near_equal(val['rel error'][0], 0.0, 1e-6)
+
 
 @unittest.skipUnless(MPI and PETScVector, "MPI and PETSc are required.")
 class TestGroupFiniteDifferenceMPI(unittest.TestCase):