Repository: OpenMDAO/OpenMDAO

Issue:
Backtracking solvers don't respect bounds added during configure in some situations.
### Description

When a nonlinear solver is attached to a Group that is not root, and that Group adds bounded outputs to descendent components during configure, the `System._has_bounds` flag is not properly propagated up the system tree.

This results in the backtracking solver effectively doing nothing because the `_has_bounds` flag of its system is False.

OpenMDAO does propagate this information up to the root system in `Group._setup` but this is only called on the root system only, and intermediate groups don't acquire flags like `_has_bounds` from their decendents if they are not the root sysstem.

### Example

The following test fails because the `_has_bounds` flag is only propagated up to the root system in `Group._setup`, not to intermediate groups that may have their own bounds-enforcing solvers.

```python
    def test_boundsenforce_configure(self):

        class _TestGroup(om.Group):

            def setup(self):
                self.add_subsystem('exec_comp', om.ExecComp(), promotes=['*'])
                self.add_subsystem('implicit_comp', om.InputResidsComp(), promotes=['*'])

            def configure(self):
                ec = self._get_subsystem('exec_comp')
                ic = self._get_subsystem('implicit_comp')

                ec.add_expr('resid_y = y - x - 2.0 * z',
                            x = {'shape': (3, 1)},
                            y = {'shape': (3, 1)},
                            z = {'shape': (3, 1)},
                            resid_y = {'shape': (3, 1)})
                ec.add_expr('resid_z = x * z + z - 4.0',
                            resid_z = {'shape': (3, 1)})
                
                ic.add_input('resid_y', shape=(3, 1))
                ic.add_input('resid_z', shape=(3, 1))
                ic.add_output('z', shape=(3, 1), lower=1.5,
                              upper=np.array([2.6, 2.5, 2.65]).reshape((3,1)))
                ic.add_output('y', shape=(3, 1))
        
        p = om.Problem()
        G_outer = p.model.add_subsystem('G_outer', om.Group())
        G_outer.add_subsystem('G_inner', _TestGroup())
        G_outer.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
        G_outer.nonlinear_solver.options['maxiter'] = 10
        G_outer.linear_solver = om.ScipyKrylov()
        G_outer.nonlinear_solver.linesearch = om.BoundsEnforceLS()

        p.setup()

        # # Test lower bounds: should go to the lower bound and stall
        p['G_outer.G_inner.x'] = 2.0
        p['G_outer.G_inner.y'] = 0.
        p['G_outer.G_inner.z'] = 1.6
        p.run_model()

        for ind in range(3):
            assert_near_equal(p['G_outer.G_inner.z'][ind], [1.5], 1e-8)
```

### OpenMDAO Version

3.37.1-dev

### Relevant environment information

_No response_

Hints from issue comments:

Solution patch:
diff --git a/openmdao/core/component.py b/openmdao/core/component.py
index 63149484e..8fca7d382 100644
--- a/openmdao/core/component.py
+++ b/openmdao/core/component.py
@@ -157,8 +157,9 @@ class Component(System):
                                   'across multiple processes')
         self.options.declare('run_root_only', types=bool, default=False,
                              desc='If True, call compute, compute_partials, linearize, '
-                                  'apply_linear, apply_nonlinear, and compute_jacvec_product '
-                                  'only on rank 0 and broadcast the results to the other ranks.')
+                                  'apply_linear, apply_nonlinear, solve_linear, solve_nonlinear, '
+                                  'and compute_jacvec_product only on rank 0 and broadcast the '
+                                  'results to the other ranks.')
         self.options.declare('always_opt', types=bool, default=False,
                              desc='If True, force nonlinear operations on this component to be '
                                   'included in the optimization loop even if this component is not '
diff --git a/openmdao/core/group.py b/openmdao/core/group.py
index f2890caab..7eaa0393a 100644
--- a/openmdao/core/group.py
+++ b/openmdao/core/group.py
@@ -799,13 +799,17 @@ class Group(System):
         self._has_resid_scaling = False
         self._has_bounds = False
 
-        for subsys in self.system_iter(include_self=True, recurse=True):
-            subsys._apply_output_solver_options()
-
-            self._has_output_scaling |= subsys._has_output_scaling
-            self._has_output_adder |= subsys._has_output_adder
-            self._has_resid_scaling |= subsys._has_resid_scaling
-            self._has_bounds |= subsys._has_bounds
+        _has_applied_options = set()
+        for grp in self.system_iter(include_self=True, recurse=True, depth_first=True, typ=Group):
+            for subsys in grp.system_iter(include_self=True, recurse=False):
+                if subsys.pathname not in _has_applied_options:
+                    subsys._apply_output_solver_options()
+                    _has_applied_options.add(subsys.pathname)
+
+                grp._has_output_scaling |= subsys._has_output_scaling
+                grp._has_output_adder |= subsys._has_output_adder
+                grp._has_resid_scaling |= subsys._has_resid_scaling
+                grp._has_bounds |= subsys._has_bounds
 
         # promoted names must be known to determine implicit connections so this must be
         # called after _setup_var_data, and _setup_var_data will have to be partially redone
diff --git a/openmdao/core/implicitcomponent.py b/openmdao/core/implicitcomponent.py
index 0f5361aa7..02feac3db 100644
--- a/openmdao/core/implicitcomponent.py
+++ b/openmdao/core/implicitcomponent.py
@@ -284,6 +284,33 @@ class ImplicitComponent(Component):
                 finally:
                     d_inputs.read_only = d_outputs.read_only = d_residuals.read_only = False
 
+    def _solve_linear_wrapper(self, *args):
+        """
+        Call solve_linear based on the value of the "run_root_only" option.
+
+        Parameters
+        ----------
+        *args : list
+            List of positional arguments.
+        """
+        d_outputs, d_residuals, mode = args
+        if self._run_root_only():
+            if self.comm.rank == 0:
+                self.solve_linear(d_outputs, d_residuals, mode)
+                if mode == 'fwd':
+                    self.comm.bcast(d_outputs.asarray(), root=0)
+                else:  # rev
+                    self.comm.bcast((d_residuals.asarray()), root=0)
+            else:
+                if mode == 'fwd':
+                    new_outs = self.comm.bcast(None, root=0)
+                    d_outputs.set_val(new_outs)
+                else:  # rev
+                    new_res = self.comm.bcast(None, root=0)
+                    d_residuals.set_val(new_res)
+        else:
+            self.solve_linear(d_outputs, d_residuals, mode)
+
     def _solve_linear(self, mode, scope_out=_UNDEFINED, scope_in=_UNDEFINED):
         """
         Apply inverse jac product. The model is assumed to be in a scaled state.
@@ -314,7 +341,7 @@ class ImplicitComponent(Component):
 
                 try:
                     with self._call_user_function('solve_linear'):
-                        self.solve_linear(d_outputs, d_residuals, mode)
+                        self._solve_linear_wrapper(d_outputs, d_residuals, mode)
                 finally:
                     d_outputs.read_only = d_residuals.read_only = False
 
diff --git a/openmdao/core/system.py b/openmdao/core/system.py
index 8710736d1..639e13ca5 100644
--- a/openmdao/core/system.py
+++ b/openmdao/core/system.py
@@ -3170,7 +3170,7 @@ class System(object, metaclass=SystemMetaclass):
 
         return ()  # for dynamic coloring or no coloring
 
-    def system_iter(self, include_self=False, recurse=True, typ=None):
+    def system_iter(self, include_self=False, recurse=True, typ=None, depth_first=False):
         """
         Yield a generator of local subsystems of this system.
 
@@ -3183,20 +3183,32 @@ class System(object, metaclass=SystemMetaclass):
         typ : type
             If not None, only yield Systems that match that are instances of the
             given type.
+        depth_first : bool
+            If recurse is True, this specifies whether subsystems are returned
+            in depth-first order (if True) or bredth-first order (if False).
 
         Yields
         ------
         type or None
         """
-        if include_self and (typ is None or isinstance(self, typ)):
-            yield self
+        if not recurse or not depth_first:
+            if include_self and (typ is None or isinstance(self, typ)):
+                yield self
 
-        for s in self._subsystems_myproc:
-            if typ is None or isinstance(s, typ):
-                yield s
-            if recurse:
-                for sub in s.system_iter(recurse=True, typ=typ):
+            for s in self._subsystems_myproc:
+                if typ is None or isinstance(s, typ):
+                    yield s
+                if recurse:
+                    for sub in s.system_iter(recurse=True, typ=typ):
+                        yield sub
+        else:
+            for s in self._subsystems_myproc:
+                for sub in s.system_iter(recurse=True, typ=typ, depth_first=True):
                     yield sub
+                if typ is None or isinstance(s, typ):
+                    yield s
+            if include_self and (typ is None or isinstance(self, typ)):
+                yield self
 
     def _all_subsystem_iter(self):
         """


Test patch:
diff --git a/openmdao/core/tests/test_run_root_only.py b/openmdao/core/tests/test_run_root_only.py
index 91e0d62be..2e70e2f04 100644
--- a/openmdao/core/tests/test_run_root_only.py
+++ b/openmdao/core/tests/test_run_root_only.py
@@ -70,6 +70,7 @@ class ImplCompArrayMatVecCounted(TestImplCompArrayMatVec):
         self.napply_linears = 0
         self.nsolve_nonlinears = 0
         self.napply_nonlinears = 0
+        self.nsolve_linears = 0
 
     def apply_nonlinear(self, inputs, outputs, residuals):
         super().apply_nonlinear(inputs, outputs, residuals)
@@ -84,6 +85,10 @@ class ImplCompArrayMatVecCounted(TestImplCompArrayMatVec):
         super().apply_linear(inputs, outputs, d_inputs, d_outputs, d_residuals, mode)
         self.napply_linears += 1
 
+    def solve_linear(self, d_outputs, d_residuals, mode):
+        super().solve_linear(d_outputs, d_residuals, mode)
+        self.nsolve_linears += 1
+
 
 @unittest.skipUnless(MPI is not None and PETScVector is not None, "MPI and PETSc are required.")
 class TestRunRootOnly(unittest.TestCase):
@@ -295,6 +300,11 @@ class TestRunRootOnly(unittest.TestCase):
         else:
             self.assertEqual(comp.napply_linears, 0)
 
+        if prob.comm.rank == 0:
+            self.assertEqual(comp.nsolve_linears, 2)
+        else:
+            self.assertEqual(comp.nsolve_linears, 0)
+
         np.testing.assert_allclose(J['comp.x', 'indeps.x'], np.eye(2))
 
 
diff --git a/openmdao/solvers/linesearch/tests/test_backtracking.py b/openmdao/solvers/linesearch/tests/test_backtracking.py
index f426c8cd2..829af2c3e 100644
--- a/openmdao/solvers/linesearch/tests/test_backtracking.py
+++ b/openmdao/solvers/linesearch/tests/test_backtracking.py
@@ -842,6 +842,51 @@ class TestFeatureLineSearch(unittest.TestCase):
         for ind in range(3):
             assert_near_equal(top['comp.z'][ind], [1.5], 1e-8)
 
+    def test_boundsenforce_configure(self):
+
+        class _TestGroup(om.Group):
+
+            def setup(self):
+                self.add_subsystem('exec_comp', om.ExecComp(), promotes=['*'])
+                self.add_subsystem('implicit_comp', om.InputResidsComp(), promotes=['*'])
+
+            def configure(self):
+                ec = self._get_subsystem('exec_comp')
+                ic = self._get_subsystem('implicit_comp')
+
+                ec.add_expr('resid_y = y - x - 2.0 * z',
+                            x = {'shape': (3, 1)},
+                            y = {'shape': (3, 1)},
+                            z = {'shape': (3, 1)},
+                            resid_y = {'shape': (3, 1)})
+                ec.add_expr('resid_z = x * z + z - 4.0',
+                            resid_z = {'shape': (3, 1)})
+                
+                ic.add_input('resid_y', shape=(3, 1))
+                ic.add_input('resid_z', shape=(3, 1))
+                ic.add_output('z', shape=(3, 1), lower=1.5,
+                              upper=np.array([2.6, 2.5, 2.65]).reshape((3,1)))
+                ic.add_output('y', shape=(3, 1))
+        
+        p = om.Problem()
+        G_outer = p.model.add_subsystem('G_outer', om.Group())
+        G_outer.add_subsystem('G_inner', _TestGroup())
+        G_outer.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
+        G_outer.nonlinear_solver.options['maxiter'] = 10
+        G_outer.linear_solver = om.ScipyKrylov()
+        G_outer.nonlinear_solver.linesearch = om.BoundsEnforceLS()
+
+        p.setup()
+
+        # # Test lower bounds: should go to the lower bound and stall
+        p['G_outer.G_inner.x'] = 2.0
+        p['G_outer.G_inner.y'] = 0.
+        p['G_outer.G_inner.z'] = 1.6
+        p.run_model()
+
+        for ind in range(3):
+            assert_near_equal(p['G_outer.G_inner.z'][ind], [1.5], 1e-8)
+
     def test_feature_armijogoldsteinls_basic(self):
 
         top = om.Problem()