Repository: OpenMDAO/dymos (python)

Issue:
Solver-based shooting with a Krylov solver fails due to the DirectSolver attached to `indep_states`
### Description

When using Newton solver-based shooting with a Krylov linear solver, the `LinearRunOnce` preconditioner fails with the following error.
It seems the error comes from the `DirectSolver` attached to `indep_states` component.

Dymos 1.9.0 didn't cause this error because `indep_states` didn't have a `DirectSolver` for the shooting settings. I believe the behavior changed since PR #991 .

### Error message
```
Traceback (most recent call last):
  File "/Users/shugo/dymos_bug/brachistochrone_shooting.py", line 77, in <module>
    dm.run_problem(p, run_driver=False)
  File "/Users/shugo/packages/dymos/dymos/run_problem.py", line 100, in run_problem
    failed = problem.run_model()
             ^^^^^^^^^^^^^^^^^^^
  File "/Users/shugo/packages/openmdao-dev/openmdao/core/problem.py", line 654, in run_model
    self.model.run_solve_nonlinear()
  File "/Users/shugo/packages/openmdao-dev/openmdao/core/system.py", line 4552, in run_solve_nonlinear
    self._solve_nonlinear()
  File "/Users/shugo/packages/openmdao-dev/openmdao/core/group.py", line 3593, in _solve_nonlinear
    self._nonlinear_solver._solve_with_cache_check()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/nonlinear/nonlinear_runonce.py", line 26, in _solve_with_cache_check
    self.solve()  # don't use caching
    ^^^^^^^^^^^^
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/nonlinear/nonlinear_runonce.py", line 45, in solve
    self._gs_iter()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/solver.py", line 803, in _gs_iter
    subsys._solve_nonlinear()
  File "/Users/shugo/packages/openmdao-dev/openmdao/core/group.py", line 3593, in _solve_nonlinear
    self._nonlinear_solver._solve_with_cache_check()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/nonlinear/nonlinear_runonce.py", line 26, in _solve_with_cache_check
    self.solve()  # don't use caching
    ^^^^^^^^^^^^
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/nonlinear/nonlinear_runonce.py", line 45, in solve
    self._gs_iter()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/solver.py", line 803, in _gs_iter
    subsys._solve_nonlinear()
  File "/Users/shugo/packages/openmdao-dev/openmdao/core/group.py", line 3593, in _solve_nonlinear
    self._nonlinear_solver._solve_with_cache_check()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/nonlinear/nonlinear_runonce.py", line 26, in _solve_with_cache_check
    self.solve()  # don't use caching
    ^^^^^^^^^^^^
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/nonlinear/nonlinear_runonce.py", line 45, in solve
    self._gs_iter()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/solver.py", line 803, in _gs_iter
    subsys._solve_nonlinear()
  File "/Users/shugo/packages/openmdao-dev/openmdao/core/group.py", line 3593, in _solve_nonlinear
    self._nonlinear_solver._solve_with_cache_check()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/solver.py", line 845, in _solve_with_cache_check
    self.solve()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/solver.py", line 592, in solve
    raise err
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/solver.py", line 588, in solve
    self._solve()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/solver.py", line 674, in _solve
    self._single_iteration()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/nonlinear/newton.py", line 227, in _single_iteration
    system._linearize(my_asm_jac, sub_do_ln=do_sub_ln)
  File "/Users/shugo/packages/openmdao-dev/openmdao/core/group.py", line 3822, in _linearize
    subsys._linear_solver._linearize()
  File "/Users/shugo/packages/openmdao-dev/openmdao/solvers/linear/direct.py", line 277, in _linearize
    matrix = self._assembled_jac._int_mtx._matrix
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: 'NoneType' object has no attribute '_matrix'
```

### Workarounds
- Switch back to Dymos 1.9.0
- Comment out the DirectSolver line added in PR #991 .
- Remove the `DirectSolver` in runscript. For example, after the setup, `traj.phases.phase0.indep_states.linear_solver = None`.

### Potential fix
I believe attaching a`DirectSolver` for `indep_states` is not necessary for shooting cases, because Dymos has a `DirectSolver` at the phase level anyway.


### Example

Here is a Brachistochrone example. I made a few changes to the collocation example in the Dymos docs:
**1. Use solver-based shooting**
```python
phase0 = dm.Phase(ode_class=BrachistochroneODE, transcription=dm.GaussLobatto(num_segments=10, compressed=True, solve_segments='forward'))
phase = traj.add_phase('phase0', phase0)
```

**2. Change the phase-level linear solver to PETScKrylov**
```python
phase0.linear_solver = om.PETScKrylov(assemble_jac=True, iprint=2)
phase0.linear_solver.precon = om.LinearRunOnce(iprint=-1)
```

### Complete script
```python
import openmdao.api as om
import dymos as dm
from dymos.examples.plotting import plot_results
from dymos.examples.brachistochrone import BrachistochroneODE
import matplotlib.pyplot as plt

#
# Initialize the Problem and the optimization driver
#
p = om.Problem(model=om.Group())
# p.driver = om.ScipyOptimizeDriver()
# p.driver.declare_coloring()

#
# Create a trajectory and add a phase to it
#
traj = p.model.add_subsystem('traj', dm.Trajectory())

# -------------------------------
# forward solver-based shooting
# -------------------------------
phase0 = dm.Phase(ode_class=BrachistochroneODE, transcription=dm.GaussLobatto(num_segments=10, compressed=True, solve_segments='forward'))
phase = traj.add_phase('phase0', phase0)

#
# Set the variables
#
phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))

phase.add_state('x', fix_initial=True, fix_final=False)

phase.add_state('y', fix_initial=True, fix_final=False)

phase.add_state('v', fix_initial=True, fix_final=False)

phase.add_control('theta', continuity=True, rate_continuity=True,
                  units='deg', lower=0.01, upper=179.9)

phase.add_parameter('g', units='m/s**2', val=9.80665)

#
# Minimize time at the end of the phase
#
phase.add_objective('time', loc='final', scaler=10)

# ---------------------------------------------
# change linear solver for shooting
# ---------------------------------------------
# linear solver
phase0.linear_solver = om.PETScKrylov(assemble_jac=True, iprint=2)
phase0.linear_solver.precon = om.LinearRunOnce(iprint=-1)

#
# Setup the Problem
#
p.setup()

# ---------------------------------------------
# Remove DirectSolver that causes an error
# ---------------------------------------------
# traj.phases.phase0.indep_states.linear_solver = None

#
# Set the initial values
#
p['traj.phase0.t_initial'] = 0.0
p['traj.phase0.t_duration'] = 2.0

p.set_val('traj.phase0.states:x', phase.interp('x', ys=[0, 10]))
p.set_val('traj.phase0.states:y', phase.interp('y', ys=[10, 5]))
p.set_val('traj.phase0.states:v', phase.interp('v', ys=[0, 9.9]))
p.set_val('traj.phase0.controls:theta', phase.interp('theta', ys=[5, 100.5]))

#
# Solve for the optimal trajectory
#
dm.run_problem(p, run_driver=False)

# plot results
sol = om.CaseReader('dymos_solution.db').get_case('final')
plot_results([('traj.phase0.timeseries.x', 'traj.phase0.timeseries.y',
               'x', 'y'),
              ('traj.phase0.timeseries.x', 'traj.phase0.timeseries.v',
               'x', 'v')],
             title='plot',
             p_sol=sol, p_sim=None)
plt.show()
```


### Dymos Version

1.9.1

### Relevant environment information

OpenMDAO 3.28.0

Pull request:
Fixed an issue where applying DirectSolver to StateIndependentsComp was breaking when used with other linear solvers under MPI.
### Summary

Previously dymos automatically assigned a DirectSolver to StateIndependentsComp when necessary, but this was problematic depending on the layout of the systems on MPI (mostly due to the presence of distributed systems within the ODE.

Dymos needs solvers on phases in the following conditions:
- `solve_segments` is True for one or more states in pseudospectral and Birkhoff transcriptions.
- pseudospectral phases with `input_initial` set to `True` for one or more states
- `duration_balance` is used and therefore a residual is associated with the time duration of a phase.
**By default**, Dymos will automatically assign `NewtonSolver` and `DirectSolver` to the phase if one of these conditions is met.

Furthermore, if a Trajectory has its phases under a ParallelGroup and there are explicit connections between these phases, some form of solver is required to ensure that the inputs/ouputs between the phases are consistent.
**By default**, Dymos will add NonlinearBlockJac and PETScKrylov as the default solvers in `trajectory.phases` when this is the case.

The user can bypass this behavior by assigning some other solver to the phase or `Trajectory.phases`. If the user really truly wants to use a RunOnce solver in these situations, then the new Trajectory and Phase option `auto_solvers` should be set to `False`.

Adds new options to Trajectory:
- `parallel_phases` dictates whether the top level phase group is a ParallelGroup or a standard Group.  Setting this to `False` allows the user to force phases to run in series when operating in parallel.
- `auto_solvers` is `True` by default which causes dymos to attempt to place solvers where it thinks they are necessary. Setting this to `False` will disable this behavior.

New options for Phase
- `auto_solvers` is `True` by default which causes dymos to attempt to place solvers where it thinks they are necessary. Setting this to `False` will disable this behavior.

### Related Issues

- Resolves #1009 
- Supercedes #1011 

### Backwards incompatibilities

None

### New Dependencies

None


Patch:
diff --git a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
index 98da7d17..22d3b3bd 100644
--- a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
+++ b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
@@ -4,6 +4,7 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "tags": [
      "active-ipynb",
      "remove-input",
@@ -76,6 +77,7 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "tags": [
      "output_scroll"
     ]
@@ -148,7 +150,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "class CannonballODE(om.ExplicitComponent):\n",
@@ -232,7 +236,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "def initial_guess(t_dur, gam0=np.pi/3):\n",
@@ -265,10 +271,12 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": false
+   },
    "outputs": [],
    "source": [
-    "p = om.Problem(model=om.Group())\n",
+    "p = om.Problem(name='cannonball_implicit_duration', model=om.Group())\n",
     "\n",
     "p.driver = om.pyOptSparseDriver()\n",
     "p.driver.options['optimizer'] = 'IPOPT'\n",
@@ -321,6 +329,14 @@
     "# A nonlinear solver is used to find the duration of required to satisfy the condition.\n",
     "phase.set_duration_balance('h', val=0.0)\n",
     "\n",
+    "# Prior to setup, this phase has the default NonlinearRunOnce and LinearRunOnce solvers.\n",
+    "# Because of the implicit behavior introduced by set_duration_balance, it will automatically\n",
+    "# use a NewtonSolver with an Armijo Goldstein linesearch.\n",
+    "# In this case that linesearch causes extrapolation of the atmospheric model into invalid regions.\n",
+    "# To account for this, we explicitly assign a NewtonSolver with a different linesearch.\n",
+    "phase.nonlinear_solver = om.NewtonSolver(iprint=0, solve_subsystems=True)\n",
+    "phase.nonlinear_solver.linesearch = None\n",
+    "\n",
     "phase.add_objective('r', loc='final', scaler=-1.0)\n",
     "\n",
     "p.model.connect('size_comp.mass', 'traj.phase.parameters:m')\n",
@@ -350,7 +366,7 @@
     "#####################################################\n",
     "# Run the optimization and final explicit simulation\n",
     "#####################################################\n",
-    "dm.run_problem(p)"
+    "dm.run_problem(p, simulate=True, make_plots=True)"
    ]
   },
   {
@@ -363,60 +379,21 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": false
+   },
    "outputs": [],
    "source": [
-    "exp_out = traj.simulate()\n",
-    "\n",
-    "rad = p.get_val('radius', units='m')[0]\n",
-    "print(f'optimal radius: {rad} m ')\n",
-    "mass = p.get_val('size_comp.mass', units='kg')[0]\n",
-    "print(f'cannonball mass: {mass} kg ')\n",
-    "angle = p.get_val('traj.phase.timeseries.gam', units='deg')[0, 0]\n",
-    "print(f'launch angle: {angle} deg')\n",
-    "max_range = p.get_val('traj.phase.timeseries.r')[-1, 0]\n",
-    "print(f'maximum range: {max_range} m')\n",
-    "\n",
-    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))\n",
-    "\n",
-    "time_imp = p.get_val('traj.phase.timeseries.time')\n",
-    "\n",
-    "time_exp = exp_out.get_val('traj.phase.timeseries.time')\n",
-    "\n",
-    "r_imp = p.get_val('traj.phase.timeseries.r')\n",
-    "\n",
-    "r_exp = exp_out.get_val('traj.phase.timeseries.r')\n",
-    "\n",
-    "h_imp = p.get_val('traj.phase.timeseries.h')\n",
-    "\n",
-    "h_exp = exp_out.get_val('traj.phase.timeseries.h')\n",
-    "\n",
-    "axes.plot(r_imp, h_imp, 'bo')\n",
-    "\n",
-    "axes.plot(r_exp, h_exp, 'r--')\n",
-    "\n",
-    "axes.set_xlabel('range (m)')\n",
-    "axes.set_ylabel('altitude (m)')\n",
-    "\n",
-    "fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))\n",
-    "states = ['r', 'h', 'v', 'gam']\n",
-    "for i, state in enumerate(states):\n",
-    "    x_imp = p.get_val(f'traj.phase.timeseries.{state}')\n",
-    "\n",
-    "    x_exp = exp_out.get_val(f'traj.phase.timeseries.{state}')\n",
-    "\n",
-    "    axes[i].set_ylabel(state)\n",
-    "\n",
-    "    axes[i].plot(time_imp, x_imp, 'bo')\n",
-    "    axes[i].plot(time_exp, x_exp, 'r--')\n",
+    "from IPython.display import IFrame\n",
     "\n",
-    "plt.show()"
+    "IFrame(src='reports/cannonball_implicit_duration/traj_results_report.html', width=700, height=1000)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "tags": [
      "remove-input",
      "remove-output"
diff --git a/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py b/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py
index c6026d70..75621f3b 100644
--- a/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py
+++ b/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py
@@ -10,7 +10,7 @@ from dymos.examples.finite_burn_orbit_raise.finite_burn_eom import FiniteBurnODE
 def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=False,
               connected=False, default_nonlinear_solver=None, default_linear_solver=None):
     """
-    Build a traejctory for the finite burn orbit raise problem.
+    Build a trajectory for the finite burn orbit raise problem.
 
     Parameters
     ----------
@@ -31,8 +31,13 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F
     t = {'gauss-lobatto': dm.GaussLobatto(num_segments=5, order=transcription_order, compressed=compressed),
          'radau': dm.Radau(num_segments=5, order=transcription_order, compressed=compressed)}
 
-    traj = dm.Trajectory(default_nonlinear_solver=default_nonlinear_solver,
-                         default_linear_solver=default_linear_solver)
+    traj = dm.Trajectory()
+
+    if default_nonlinear_solver is not None:
+        traj.phases.nonlinear_solver = default_nonlinear_solver
+
+    if default_linear_solver is not None:
+        traj.phases.linear_solver = default_linear_solver
 
     traj.add_parameter('c', opt=False, val=1.5, units='DU/TU',
                        targets={'burn1': ['c'], 'burn2': ['c']})
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 73a5f8b1..d036fafc 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -114,6 +114,8 @@ class Phase(om.Group):
                              desc='Keyword arguments provided when initializing the ODE System')
         self.options.declare('transcription', types=TranscriptionBase,
                              desc='Transcription technique of the optimal control problem.')
+        self.options.declare('auto_solvers', types=bool, default=True,
+                             desc='If True, attempt to automatically assign solvers if necessary.')
 
     def add_state(self, name, units=_unspecified, shape=_unspecified,
                   rate_source=_unspecified, targets=_unspecified,
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 31cce9b7..55886fe7 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -61,7 +61,7 @@ class Trajectory(om.Group):
         self._has_connected_phases = False
 
         self.parameter_options = {}
-        self.phases = om.ParallelGroup()
+        self.phases = om.ParallelGroup() if self.options['parallel_phases'] else om.Group()
 
     def initialize(self):
         """
@@ -69,17 +69,12 @@ class Trajectory(om.Group):
         """
         self.options.declare('sim_mode', types=bool, default=False,
                              desc='Used internally by Dymos when invoking simulate on a trajectory')
-        self.options.declare('default_nonlinear_solver',
-                             types=(om.NonlinearBlockJac, om.NewtonSolver, om.BroydenSolver),
-                             default=None, allow_none=True, recordable=False,
-                             desc='A nonlinear solver to be used when Phases are connected but being '
-                                  'run in parallel. If not specified, Dymos will automatically use '
-                                  'NonlinearBlockJac in this situation.')
-        self.options.declare('default_linear_solver', default=None, allow_none=True, recordable=False,
-                             types=(om.PETScKrylov, om.LinearBlockJac, om.NonlinearBlockGS),
-                             desc='A linear solver to be used when Phases are connected but being '
-                                  'run in parallel. If not specified, Dymos will automatically use '
-                                  'PETScKrylov in this situation.')
+        self.options.declare('parallel_phases',
+                             default=True, types=bool,
+                             desc='If True, the top-level container of all phases will be a ParallelGroup,'
+                                  'otherwise it will be a standard OpenMDAO Group.')
+        self.options.declare('auto_solvers', types=bool, default=True,
+                             desc='If True, attempt to automatically assign solvers if necessary.')
 
     def add_phase(self, name, phase, **kwargs):
         """
@@ -1018,27 +1013,30 @@ class Trajectory(om.Group):
         These solvers can be changed through the
         'default_nonlinear_solver' and 'default_linear_solver' options.
         """
-        if self._has_connected_phases and MPI and self.comm.size > 1:
+        if not self.options['auto_solvers']:
+            return
 
+        warn = False
+        if self._has_connected_phases and isinstance(self.phases, om.ParallelGroup) and self.comm.size > 1:
+            msg = (f'{self.msginfo}\n'
+                   f'  Non-default solvers are required\n'
+                   f'    connected phases in parallel: True\n')
             if isinstance(self.phases.nonlinear_solver, om.NonlinearRunOnce):
-                if self.options['default_nonlinear_solver'] is None:
-                    msg = f'{self.msginfo}: Setting phases.nonlinear_solver to `om.NonlinearBlockJac(iprint=0)`.\n' \
-                          f'Connected phases in parallel require a non-default nonlinear solver.\n' \
-                          f'Use {self.pathname}.options[\'default_nonlinear_solver\'] to explicitly set the solver.'
-                    om.issue_warning(msg)
-                    self.phases.nonlinear_solver = om.NonlinearBlockJac(iprint=0)
-                else:
-                    self.phases.nonlinear_solver = self.options['default_nonlinear_solver']
+                warn = True
+                msg += (f'  Setting \n'
+                        f'    {self.pathname}.phases.nonlinear_solver = om.NonlinearBlockJac(iprint=0)\n'
+                        f'  Explicitly set {self.pathname}.phases.nonlinear_solver to override.\n')
+                self.phases.nonlinear_solver = om.NonlinearBlockJac(iprint=0)
 
             if isinstance(self.phases.linear_solver, om.LinearRunOnce):
-                if self.options['default_linear_solver'] is None:
-                    msg = f'{self.msginfo}: Setting phases.linear_solver to `om.PETScKrylov()`.\n' \
-                          f'Connected phases in parallel require a non-default linear solver.\n' \
-                          f'Use {self.pathname}.options[\'default_linear_solver\'] to explicitly set the solver.'
-                    om.issue_warning(msg)
-                    self.phases.linear_solver = om.PETScKrylov()
-                else:
-                    self.phases.linear_solver = self.options['default_linear_solver']
+                warn = True
+                msg += (f'  Setting\n'
+                        f'    {self.pathname}.phases.linear_solver = om.PETScKrylov(iprint=0)\n'
+                        f'  Explicitly set {self.pathname}.phases.linear_solver to override.\n')
+                self.phases.linear_solver = om.PETScKrylov(iprint=0)
+
+            if warn:
+                om.issue_warning(msg)
 
     def configure(self):
         """
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index 1c7d2824..e51342eb 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -10,7 +10,6 @@ from .explicit_shooting_continuity_comp import ExplicitShootingContinuityComp
 from ..transcription_base import TranscriptionBase
 from ..grid_data import BirkhoffGrid, GaussLobattoGrid, RadauGrid, UniformGrid
 from .ode_integration_comp import ODEIntegrationComp
-from ..._options import options as dymos_options
 from ...utils.misc import get_rate_units, CoerceDesvar
 from ...utils.indexing import get_src_indices_by_row
 from ...utils.introspection import get_promoted_vars, get_source_metadata, get_targets, _get_targets_metadata
@@ -675,7 +674,7 @@ class ExplicitShooting(TranscriptionBase):
         """
         pass
 
-    def configure_solvers(self, phase):
+    def configure_solvers(self, phase, requires_solvers=None):
         """
         Setup the solvers for this transcription.
 
@@ -683,8 +682,10 @@ class ExplicitShooting(TranscriptionBase):
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
+        requires_solvers : None
+            Required to extend TranscriptionBase.configure_solvers but not used.
         """
-        pass
+        super().configure_solvers(phase, requires_solvers=None)
 
     def get_parameter_connections(self, name, phase):
         """
diff --git a/dymos/transcriptions/pseudospectral/birkhoff.py b/dymos/transcriptions/pseudospectral/birkhoff.py
index df914ead..9840ab15 100644
--- a/dymos/transcriptions/pseudospectral/birkhoff.py
+++ b/dymos/transcriptions/pseudospectral/birkhoff.py
@@ -343,7 +343,7 @@ class Birkhoff(TranscriptionBase):
         """
         pass
 
-    def configure_solvers(self, phase):
+    def configure_solvers(self, phase, requires_solvers=None):
         """
         Configure the solvers.
 
@@ -351,35 +351,17 @@ class Birkhoff(TranscriptionBase):
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
+        requires_solvers : dict[str: bool]
+            A dictionary mapping a string descriptor of a reason why a solver is required,
+            and whether a solver is required.
         """
         ode_iter_group = phase._get_subsystem('ode_iter_group')
+        req_solvers = {'implicit outputs': ode_iter_group._implicit_outputs}
 
-        if ode_iter_group._implicit_outputs or self._implicit_duration:
-            # Only override the solvers if the user hasn't set them to something else.
-            if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce):
-                msg = f'{phase.msginfo}: Phase requires non-default nonlinear solver due to the use of\n' \
-                      f'solve_segments or the use of set_duration_balance.\n' \
-                      f'Setting nonlinear solver to ' \
-                      f'om.NewtonSolver(solve_subsystems=True, maxiter=100, iprint=2, stall_limit=3)\n' \
-                      f'by default.' \
-                      f'To avoid this warning explicitly assign a nonlinear solver to this phase.'
-
-                om.issue_warning(msg)
-                newton = phase.nonlinear_solver = om.NewtonSolver()
-                newton.options['solve_subsystems'] = True
-                newton.options['maxiter'] = 100
-                newton.options['iprint'] = 2
-                newton.options['stall_limit'] = 3
-                newton.linesearch = om.BoundsEnforceLS()
-
-            if isinstance(phase.linear_solver, om.LinearRunOnce):
-                msg = f'{phase.msginfo}: Phase requires non-default linear solver due to the use of\n' \
-                      f'solve_segments or the use of set_duration_balance.\n' \
-                      f'Setting linear solver to om.DirectSolver() by default.' \
-                      f'To avoid this warning explicitly assign a linear solver to this phase.'
-
-                om.issue_warning(msg)
-                phase.linear_solver = om.DirectSolver()
+        if requires_solvers is not None:
+            req_solvers.update(requires_solvers)
+
+        super().configure_solvers(phase, requires_solvers=req_solvers)
 
     def configure_timeseries_outputs(self, phase):
         """
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index c9b9f962..41144006 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -20,6 +20,12 @@ class PseudospectralBase(TranscriptionBase):
     **kwargs : dict
         Dictionary of optional arguments.
     """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self.any_solved_segs = False
+        self.any_connected_opt_segs = False
+
     def initialize(self):
         """
         Declare transcription options.
@@ -80,8 +86,6 @@ class PseudospectralBase(TranscriptionBase):
         """
         grid_data = self.grid_data
 
-        self.any_solved_segs = False
-        self.any_connected_opt_segs = False
         for options in phase.state_options.values():
             # Transcription solve_segments overrides state solve_segments if its not set
             if options['solve_segments'] is None:
@@ -95,7 +99,6 @@ class PseudospectralBase(TranscriptionBase):
         if self.any_solved_segs or self.any_connected_opt_segs:
             indep = StateIndependentsComp(grid_data=grid_data,
                                           state_options=phase.state_options)
-            indep.linear_solver = om.DirectSolver()
         else:
             indep = om.IndepVarComp()
 
@@ -495,7 +498,7 @@ class PseudospectralBase(TranscriptionBase):
         """
         pass
 
-    def configure_solvers(self, phase):
+    def configure_solvers(self, phase, requires_solvers=None):
         """
         Configure the solvers.
 
@@ -503,19 +506,15 @@ class PseudospectralBase(TranscriptionBase):
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
-        """
-        if self.any_solved_segs or self._implicit_duration:
-            # Only override the solvers if the user hasn't set them to something else.
-            if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce):
-                newton = phase.nonlinear_solver = om.NewtonSolver()
-                newton.options['solve_subsystems'] = True
-                newton.options['maxiter'] = 100
-                newton.options['iprint'] = 2
-                newton.options['stall_limit'] = 3
-                newton.linesearch = om.BoundsEnforceLS()
-
-            if isinstance(phase.linear_solver, om.LinearRunOnce):
-                phase.linear_solver = om.DirectSolver()
+        requires_solvers : dict[str: bool]
+            A dictionary mapping a string descriptor of a reason why a solver is required,
+            and whether a solver is required.
+        """
+        req_solvers = {'solved segments': self.any_solved_segs,
+                       'input initial': self.any_connected_opt_segs}
+        if requires_solvers is not None:
+            req_solvers.update(requires_solvers)
+        super().configure_solvers(phase, requires_solvers=req_solvers)
 
     def setup_timeseries_outputs(self, phase):
         """
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index db64a0a0..077fb190 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -399,17 +399,53 @@ class TranscriptionBase(object):
         """
         raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement method setup_solvers.')
 
-    def configure_solvers(self, phase):
+    def configure_solvers(self, phase, requires_solvers=None):
         """
-        Configure the solvers for this transcription.
+        Configure the solvers.
 
         Parameters
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
-        """
-        raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement method '
-                                  f'configure_solvers.')
+        requires_solvers : dict[str: bool]
+            A dictionary mapping a string descriptor of a reason why a solver is required,
+            and whether a solver is required.
+        """
+        if not phase.options['auto_solvers']:
+            return
+
+        req_solvers = {'implicit duration': self._implicit_duration}
+
+        if requires_solvers is not None:
+            req_solvers.update(requires_solvers)
+
+        reasons_txt = '\n'.join(f'    {key}: {val}' for key, val in req_solvers.items())
+        warn = False
+
+        if any(req_solvers.values()):
+            msg = (f'{phase.msginfo}\n'
+                   f'  Non-default solvers are required\n'
+                   f'{reasons_txt}\n')
+            if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce):
+                msg += (f'  Setting `{phase.pathname}.nonlinear_solver = om.NewtonSolver(iprint=0, '
+                        f'solve_subsystems=True, maxiter=1000, stall_limit=3)`\n'
+                        f'  Explicitly set {phase.pathname}.nonlinear_solver to override.\n')
+                warn = True
+                phase.nonlinear_solver = om.NewtonSolver(iprint=0)
+                phase.nonlinear_solver.options['solve_subsystems'] = True
+                phase.nonlinear_solver.options['maxiter'] = 1000
+                phase.nonlinear_solver.options['stall_limit'] = 3
+                phase.nonlinear_solver.linesearch = om.ArmijoGoldsteinLS()
+
+            if isinstance(phase.linear_solver, om.LinearRunOnce):
+                warn = True
+                msg += (f'  Setting `{phase.pathname}.linear_solver = om.DirectSolver(iprint=2)`\n'
+                        f'  Explicitly set {phase.pathname}.linear_solver to override.\n')
+                phase.linear_solver = om.DirectSolver(iprint=0)
+
+            if warn:
+                msg += f'  Set `{phase.pathname}.options["auto_solvers"] = False` to disable this behavior.'
+                om.issue_warning(msg)
 
     def setup_timeseries_outputs(self, phase):
         """


Test patch:
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
index f0909489..aa9110de 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
@@ -9,6 +9,11 @@ import dymos.examples.brachistochrone.test.ex_brachistochrone_vector_states as e
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
+try:
+    from petsc4py import PETSc
+except ImportError:
+    PETSc = None
+
 from openmdao.utils.general_utils import set_pyoptsparse_opt
 
 OPT, OPTIMIZER = set_pyoptsparse_opt('SLSQP')
@@ -16,7 +21,8 @@ OPT, OPTIMIZER = set_pyoptsparse_opt('SLSQP')
 
 def _make_problem(transcription='gauss-lobatto', num_segments=8, transcription_order=3,
                   compressed=True, optimizer='SLSQP', force_alloc_complex=False,
-                  solve_segments=False, y_bounds=None):
+                  solve_segments=False, y_bounds=None, phase_nonlinear_solver=None,
+                  phase_linear_solver=None):
     p = om.Problem(model=om.Group())
 
     p.driver = om.pyOptSparseDriver()
@@ -44,6 +50,12 @@ def _make_problem(transcription='gauss-lobatto', num_segments=8, transcription_o
 
     traj = dm.Trajectory()
     phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
+
+    if phase_nonlinear_solver is not None:
+        phase.nonlinear_solver = phase_nonlinear_solver
+    if phase_linear_solver is not None:
+        phase.linear_solver = phase_linear_solver
+
     p.model.add_subsystem('traj0', traj)
     traj.add_phase('phase0', phase)
 
@@ -461,3 +473,90 @@ class TestBrachistochroneSolveSegments(unittest.TestCase):
 
             warnings.simplefilter('always')
             self.assertIn(expected_warning, [str(w.message) for w in ctx])
+
+    @unittest.skipIf(PETSc is None, 'PETSc is not available.')
+    def test_brachistochrone_solve_segments_radau_krylov_solver(self, optimizer='SLSQP'):
+        #
+        # Initialize the Problem and the optimization driver
+        #
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = optimizer
+        if optimizer == 'SNOPT':
+            p.driver.opt_settings['iSumm'] = 6
+            p.driver.opt_settings['Verify level'] = 3
+        elif optimizer == 'IPOPT':
+            p.driver.opt_settings['mu_init'] = 1e-3
+            p.driver.opt_settings['max_iter'] = 500
+            p.driver.opt_settings['print_level'] = 0
+            p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'  # for faster convergence
+            p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
+            p.driver.opt_settings['mu_strategy'] = 'monotone'
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        #
+        # Create a trajectory and add a phase to it
+        #
+        traj = p.model.add_subsystem('traj0', dm.Trajectory())
+
+        # -------------------------------
+        # forward solver-based shooting
+        # -------------------------------
+        phase0 = dm.Phase(ode_class=BrachistochroneODE,
+                          transcription=dm.GaussLobatto(num_segments=10, compressed=True,
+                                                        solve_segments='forward'))
+        phase = traj.add_phase('phase0', phase0)
+
+        #
+        # Set the variables
+        #
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=False)
+
+        phase.add_state('y', fix_initial=True, fix_final=False)
+
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_control('theta', continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665)
+
+        #
+        # Minimize time at the end of the phase
+        #
+        phase.add_objective('time', loc='final', scaler=10)
+        phase.add_boundary_constraint('x', loc='final', equals=10)
+        phase.add_boundary_constraint('y', loc='final', equals=5)
+
+        # ---------------------------------------------
+        # change linear solver for shooting
+        # ---------------------------------------------
+        # linear solver
+        phase0.linear_solver = om.PETScKrylov(assemble_jac=False, iprint=-1)
+        phase0.linear_solver.precon = om.LinearRunOnce(iprint=-1)
+
+        #
+        # Setup the Problem
+        #
+        p.setup()
+
+        #
+        # Set the initial values
+        #
+        p['traj0.phase0.t_initial'] = 0.0
+        p['traj0.phase0.t_duration'] = 2.0
+
+        p.set_val('traj0.phase0.states:x', phase.interp('x', ys=[0, 10]))
+        p.set_val('traj0.phase0.states:y', phase.interp('y', ys=[10, 5]))
+        p.set_val('traj0.phase0.states:v', phase.interp('v', ys=[0, 9.9]))
+        p.set_val('traj0.phase0.controls:theta', phase.interp('theta', ys=[5, 100.5]))
+
+        #
+        # Solve for the optimal trajectory
+        #
+        dm.run_problem(p, run_driver=True)
+
+        self.assert_results(p)
diff --git a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
index bc607b7a..fe61c251 100644
--- a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
+++ b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
@@ -1,14 +1,10 @@
 import unittest
 
-import matplotlib.pyplot as plt
 import numpy as np
 
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 
-plt.switch_backend('Agg')
-
-
 def initial_guess(t_dur, gam0=np.pi/3):
     t = np.linspace(0, t_dur, num=100)
     g = -9.80665
@@ -206,6 +202,12 @@ class TestTwoPhaseCannonballForDocs(unittest.TestCase):
         # converge to the initial point.
         phase.set_duration_balance('h', val=0.0)
 
+        # In this problem, the default ArmijoGoldsteinLS has issues with extrapolating
+        # the states and causes the optimization to fail.
+        # Using the default linesearch or BoundsEnforceLS work better here.
+        phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, iprint=2)
+        phase.nonlinear_solver.linesearch = None
+
         phase.add_objective('r', loc='final', scaler=-1.0)
 
         p.model.connect('size_comp.mass', 'traj.phase.parameters:m')
@@ -235,60 +237,11 @@ class TestTwoPhaseCannonballForDocs(unittest.TestCase):
         #####################################################
         # Run the optimization and final explicit simulation
         #####################################################
-        dm.run_problem(p)
+        dm.run_problem(p, simulate=True)
 
         assert_near_equal(p.get_val('traj.phase.states:r')[-1],
                           3183.25, tolerance=1.0)
 
-        exp_out = traj.simulate()
-
-        #############################################
-        # Plot the results
-        #############################################
-        rad = p.get_val('radius', units='m')[0]
-        print(f'optimal radius: {rad} m ')
-        mass = p.get_val('size_comp.mass', units='kg')[0]
-        print(f'cannonball mass: {mass} kg ')
-        angle = p.get_val('traj.phase.timeseries.gam', units='deg')[0, 0]
-        print(f'launch angle: {angle} deg')
-        max_range = p.get_val('traj.phase.timeseries.r')[-1, 0]
-        print(f'maximum range: {max_range} m')
-
-        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))
-
-        time_imp = p.get_val('traj.phase.timeseries.time')
-
-        time_exp = exp_out.get_val('traj.phase.timeseries.time')
-
-        r_imp = p.get_val('traj.phase.timeseries.r')
-
-        r_exp = exp_out.get_val('traj.phase.timeseries.r')
-
-        h_imp = p.get_val('traj.phase.timeseries.h')
-
-        h_exp = exp_out.get_val('traj.phase.timeseries.h')
-
-        axes.plot(r_imp, h_imp, 'bo')
-
-        axes.plot(r_exp, h_exp, 'b--')
-
-        axes.set_xlabel('range (m)')
-        axes.set_ylabel('altitude (m)')
-
-        fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))
-        states = ['r', 'h', 'v', 'gam']
-        for i, state in enumerate(states):
-            x_imp = p.get_val(f'traj.phase.timeseries.{state}')
-
-            x_exp = exp_out.get_val(f'traj.phase.timeseries.{state}')
-
-            axes[i].set_ylabel(state)
-
-            axes[i].plot(time_imp, x_imp, 'bo')
-            axes[i].plot(time_exp, x_exp, 'b--')
-
-        plt.show()
-
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()
diff --git a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py
index 78c3a2ef..80fa1139 100644
--- a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py
+++ b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py
@@ -1,13 +1,8 @@
 import unittest
 
-import matplotlib.pyplot as plt
-
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 
-plt.switch_backend('Agg')
-
-
 @use_tempdirs
 class TestTwoPhaseCannonballForDocs(unittest.TestCase):
 
@@ -257,75 +252,11 @@ class TestTwoPhaseCannonballForDocs(unittest.TestCase):
         #####################################################
         # Run the optimization and final explicit simulation
         #####################################################
-        dm.run_problem(p)
+        dm.run_problem(p, simulate=True)
 
         assert_near_equal(p.get_val('traj.descent.states:r')[-1],
                           3183.25, tolerance=1.0E-2)
 
-        # use the explicit simulation to check the final collocation solution accuracy
-        exp_out = traj.simulate()
-
-        #############################################
-        # Plot the results
-        #############################################
-        rad = p.get_val('radius', units='m')[0]
-        print(f'optimal radius: {rad} m ')
-        mass = p.get_val('size_comp.mass', units='kg')[0]
-        print(f'cannonball mass: {mass} kg ')
-        angle = p.get_val('traj.ascent.timeseries.gam', units='deg')[0, 0]
-        print(f'launch angle: {angle} deg')
-        max_range = p.get_val('traj.descent.timeseries.r')[-1, 0]
-        print(f'maximum range: {max_range} m')
-
-        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))
-
-        time_imp = {'ascent': p.get_val('traj.ascent.timeseries.time'),
-                    'descent': p.get_val('traj.descent.timeseries.time')}
-
-        time_exp = {'ascent': exp_out.get_val('traj.ascent.timeseries.time'),
-                    'descent': exp_out.get_val('traj.descent.timeseries.time')}
-
-        r_imp = {'ascent': p.get_val('traj.ascent.timeseries.r'),
-                 'descent': p.get_val('traj.descent.timeseries.r')}
-
-        r_exp = {'ascent': exp_out.get_val('traj.ascent.timeseries.r'),
-                 'descent': exp_out.get_val('traj.descent.timeseries.r')}
-
-        h_imp = {'ascent': p.get_val('traj.ascent.timeseries.h'),
-                 'descent': p.get_val('traj.descent.timeseries.h')}
-
-        h_exp = {'ascent': exp_out.get_val('traj.ascent.timeseries.h'),
-                 'descent': exp_out.get_val('traj.descent.timeseries.h')}
-
-        axes.plot(r_imp['ascent'], h_imp['ascent'], 'bo')
-
-        axes.plot(r_imp['descent'], h_imp['descent'], 'ro')
-
-        axes.plot(r_exp['ascent'], h_exp['ascent'], 'b--')
-
-        axes.plot(r_exp['descent'], h_exp['descent'], 'r--')
-
-        axes.set_xlabel('range (m)')
-        axes.set_ylabel('altitude (m)')
-
-        fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))
-        states = ['r', 'h', 'v', 'gam']
-        for i, state in enumerate(states):
-            x_imp = {'ascent': p.get_val(f'traj.ascent.timeseries.{state}'),
-                     'descent': p.get_val(f'traj.descent.timeseries.{state}')}
-
-            x_exp = {'ascent': exp_out.get_val(f'traj.ascent.timeseries.{state}'),
-                     'descent': exp_out.get_val(f'traj.descent.timeseries.{state}')}
-
-            axes[i].set_ylabel(state)
-
-            axes[i].plot(time_imp['ascent'], x_imp['ascent'], 'bo')
-            axes[i].plot(time_imp['descent'], x_imp['descent'], 'ro')
-            axes[i].plot(time_exp['ascent'], x_exp['ascent'], 'b--')
-            axes[i].plot(time_exp['descent'], x_exp['descent'], 'r--')
-
-        plt.show()
-
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()
diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py b/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py
index c3d11c8d..ad55821b 100644
--- a/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py
+++ b/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py
@@ -114,7 +114,6 @@ class TestExampleTwoBurnOrbitRaiseConnected(unittest.TestCase):
                                          show_output=False, restart='dymos_solution.db',
                                          connected=True, solution_record_file='dymos_solution2.db',
                                          simulation_record_file='dymos_simulation2.db')
-#
 
         case2 = om.CaseReader('dymos_solution2.db').get_case('final')
         sim_case2 = om.CaseReader('dymos_simulation2.db').get_case('final')
@@ -127,41 +126,30 @@ class TestExampleTwoBurnOrbitRaiseConnected(unittest.TestCase):
     def test_restart_from_solution_radau_to_connected(self):
         optimizer = 'IPOPT'
 
-        expected_warnings = \
-            [(om.OpenMDAOWarning,
-              "'traj' <class Trajectory>: Setting phases.nonlinear_solver to `om.NonlinearBlockJac(iprint=0)`.\n"
-              "Connected phases in parallel require a non-default nonlinear solver.\n"
-              "Use traj.options[\'default_nonlinear_solver\'] to explicitly set the solver."),
-             (om.OpenMDAOWarning,
-              "'traj' <class Trajectory>: Setting phases.linear_solver to `om.PETScKrylov()`.\n"
-              "Connected phases in parallel require a non-default linear solver.\n"
-              "Use traj.options[\'default_linear_solver\'] to explicitly set the solver.")]
-
-        with assert_warnings(expected_warnings):
-            p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
-                                             compressed=False, optimizer=optimizer,
-                                             show_output=False, connected=True)
+        p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
+                                         compressed=False, optimizer=optimizer,
+                                         show_output=False, connected=True)
 
-            if p.model.traj.phases.burn2 in p.model.traj.phases._subsystems_myproc:
-                assert_near_equal(p.get_val('traj.burn2.states:deltav')[0], 0.3995,
-                                  tolerance=4.0E-3)
+        if p.model.traj.phases.burn2 in p.model.traj.phases._subsystems_myproc:
+            assert_near_equal(p.get_val('traj.burn2.states:deltav')[0], 0.3995,
+                              tolerance=4.0E-3)
 
-            case1 = om.CaseReader('dymos_solution.db').get_case('final')
-            sim_case1 = om.CaseReader('dymos_simulation.db').get_case('final')
+        case1 = om.CaseReader('dymos_solution.db').get_case('final')
+        sim_case1 = om.CaseReader('dymos_simulation.db').get_case('final')
 
-            # Run again without an actual optimizer
-            p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
-                                             compressed=False, optimizer=optimizer, run_driver=False,
-                                             show_output=False, restart='dymos_solution.db',
-                                             connected=True, solution_record_file='dymos_solution2.db',
-                                             simulation_record_file='dymos_simulation2.db')
+        # Run again without an actual optimizer
+        p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
+                                         compressed=False, optimizer=optimizer, run_driver=False,
+                                         show_output=False, restart='dymos_solution.db',
+                                         connected=True, solution_record_file='dymos_solution2.db',
+                                         simulation_record_file='dymos_simulation2.db')
 
-            case2 = om.CaseReader('dymos_solution2.db').get_case('final')
-            sim_case2 = om.CaseReader('dymos_simulation2.db').get_case('final')
+        case2 = om.CaseReader('dymos_solution2.db').get_case('final')
+        sim_case2 = om.CaseReader('dymos_simulation2.db').get_case('final')
 
-            # Verify that the second case has the same inputs and outputs
-            assert_cases_equal(case1, case2, tol=1.0E-8)
-            assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8)
+        # Verify that the second case has the same inputs and outputs
+        assert_cases_equal(case1, case2, tol=1.0E-8)
+        assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8)
 
 
 if __name__ == '__main__':  # pragma: no cover
diff --git a/dymos/trajectory/test/test_trajectory.py b/dymos/trajectory/test/test_trajectory.py
index 45ef246d..c65d9080 100644
--- a/dymos/trajectory/test/test_trajectory.py
+++ b/dymos/trajectory/test/test_trajectory.py
@@ -1779,7 +1779,7 @@ class TestInvalidLinkages(unittest.TestCase):
         p.model.add_design_var('radius', lower=0.01, upper=0.10,
                                ref0=0.01, ref=0.10, units='m')
 
-        traj = p.model.add_subsystem('traj', dm.Trajectory())
+        traj = p.model.add_subsystem('traj', dm.Trajectory(parallel_phases=True))
 
         transcription = dm.Radau(num_segments=5, order=3, compressed=True)
 
@@ -1856,9 +1856,6 @@ class TestInvalidLinkages(unittest.TestCase):
         p.model.connect('size_comp.mass', 'traj.parameters:m')
         p.model.connect('size_comp.S', 'traj.parameters:S')
 
-        # A linear solver at the top level can improve performance.
-        p.model.linear_solver = om.DirectSolver()
-
         # Finish Problem Setup
         p.setup()