Fixed an issue where applying DirectSolver to StateIndependentsComp was breaking when used with other linear solvers under MPI.

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"

dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py

@@ -9,14 +9,20 @@
 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')
 
 
 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 @@ def test_brachistochrone_bounded_solve_segments(self):
 
             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)

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 @@ def compute(self, inputs, outputs):
         # 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 @@ def compute(self, inputs, outputs):
         #####################################################
         # 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()