diff --git a/benchmark/benchmark_balanced_field.py b/benchmark/benchmark_balanced_field.py
index cfc4b2f26..332103315 100644
--- a/benchmark/benchmark_balanced_field.py
+++ b/benchmark/benchmark_balanced_field.py
@@ -20,7 +20,7 @@ def _run_balanced_field_length_problem(tx=dm.GaussLobatto, timeseries=True, sim=
     p.driver.options['optimizer'] = optimizer
     p.driver.options['print_results'] = False
     if optimizer == 'IPOPT':
-        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['print_level'] = 0
 
     # First Phase: Brake release to V1 - both engines operable
     br_to_v1 = dm.Phase(ode_class=BalancedFieldODEComp,
diff --git a/benchmark/benchmark_racecar.py b/benchmark/benchmark_racecar.py
index d30fe1124..b1b5636c4 100644
--- a/benchmark/benchmark_racecar.py
+++ b/benchmark/benchmark_racecar.py
@@ -129,7 +129,7 @@ def _run_racecar_problem(transcription, timeseries=False):
     p.driver.opt_settings['compl_inf_tol'] = 1e-3
     p.driver.opt_settings['acceptable_iter'] = 0
     p.driver.opt_settings['tol'] = 1e-3
-    p.driver.opt_settings['print_level'] = 5
+    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'
diff --git a/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py b/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py
index 3eb1f6968..37d55b62a 100644
--- a/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py
+++ b/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py
@@ -29,7 +29,7 @@ def test_balanced_field_length_for_docs(self):
         p.driver.options['optimizer'] = optimizer
         p.driver.options['print_results'] = False
         if optimizer == 'IPOPT':
-            p.driver.opt_settings['print_level'] = 5
+            p.driver.opt_settings['print_level'] = 0
             p.driver.opt_settings['derivative_test'] = 'first-order'
 
         # First Phase: Brake release to V1 - both engines operable
diff --git a/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py b/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py
index 2811178ec..04aeb4a6f 100644
--- a/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py
+++ b/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py
@@ -158,7 +158,7 @@ def _run_problem(self, tx):
         p.driver.options['optimizer'] = 'IPOPT'
         p.driver.options['print_results'] = True
 
-        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['print_level'] = 0
         p.driver.opt_settings['mu_strategy'] = 'adaptive'
 
         p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
diff --git a/dymos/examples/balanced_field/test/test_balanced_field_length.py b/dymos/examples/balanced_field/test/test_balanced_field_length.py
index 6a23e3010..7e373c76a 100644
--- a/dymos/examples/balanced_field/test/test_balanced_field_length.py
+++ b/dymos/examples/balanced_field/test/test_balanced_field_length.py
@@ -18,7 +18,7 @@ def _make_problem(self):
 
         p.driver = om.pyOptSparseDriver()
         p.driver.options['optimizer'] = 'IPOPT'
-        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['print_level'] = 0
         p.driver.opt_settings['derivative_test'] = 'first-order'
 
         p.driver.declare_coloring()
diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py
index 6c13227ab..c3dbc8c14 100644
--- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py
+++ b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py
@@ -310,7 +310,7 @@ def test_brachistochrone_for_docs_coloring_demo_solve_segments(self):
         #
         p = om.Problem(model=om.Group())
         p.driver = om.pyOptSparseDriver(optimizer='IPOPT')
-        p.driver.opt_settings['print_level'] = 4
+        p.driver.opt_settings['print_level'] = 0
         # p.driver.declare_coloring()
 
         #
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
index 275423a57..c38057867 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
@@ -27,7 +27,7 @@ def _make_problem(transcription='gauss-lobatto', num_segments=8, transcription_o
     elif optimizer == 'IPOPT':
         p.driver.opt_settings['mu_init'] = 1e-3
         p.driver.opt_settings['max_iter'] = 500
-        p.driver.opt_settings['print_level'] = 5
+        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'
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py b/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
index a5eaf23c1..8b363b183 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
@@ -1,308 +1,312 @@
 """Unit Tests for the code that does automatic report generation"""
 import unittest
 import pathlib
-import sys
 import os
-from io import StringIO
+from packaging.version import Version
 
 import openmdao.api as om
-from openmdao.test_suite.components.paraboloid import Paraboloid
-from openmdao.test_suite.components.sellar_feature import SellarMDA
 import openmdao.core.problem
-from openmdao.core.constants import _UNDEFINED
-from openmdao.utils.assert_utils import assert_warning
-from openmdao.utils.general_utils import set_pyoptsparse_opt
-from openmdao.utils.reports_system import set_default_reports_dir, _reports_dir, register_report, \
-    list_reports, clear_reports, run_n2_report, setup_default_reports, report_function
+import openmdao.utils.reports_system as reports_system
 from openmdao.utils.testing_utils import use_tempdirs
-from openmdao.utils.mpi import MPI
 from openmdao.utils.tests.test_hooks import hooks_active
 from openmdao.visualization.n2_viewer.n2_viewer import _default_n2_filename
 from openmdao.visualization.scaling_viewer.scaling_report import _default_scaling_filename
+from openmdao import __version__ as openmdao_version
 
-try:
-    from openmdao.vectors.petsc_vector import PETScVector
-except ImportError:
-    PETScVector = None
-
-OPT, OPTIMIZER = set_pyoptsparse_opt('SLSQP')
-
-if OPTIMIZER:
-    from openmdao.drivers.pyoptsparse_driver import pyOptSparseDriver
 
 import dymos as dm
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
 
 
-@use_tempdirs
-class TestSubproblemReportToggle(unittest.TestCase):
+def setup_model_radau(do_reports):
+    p = om.Problem(model=om.Group())
 
-    def setUp(self):
-        self.n2_filename = _default_n2_filename
-        self.scaling_filename = _default_scaling_filename
+    p.driver = om.ScipyOptimizeDriver()
+    p.driver.declare_coloring(tol=1.0E-12)
 
-        # set things to a known initial state for all the test runs
-        openmdao.core.problem._problem_names = []  # need to reset these to simulate separate runs
-        os.environ.pop('OPENMDAO_REPORTS', None)
-        os.environ.pop('OPENMDAO_REPORTS_DIR', None)
-        # We need to remove the TESTFLO_RUNNING environment variable for these tests to run.
-        # The reports code checks to see if TESTFLO_RUNNING is set and will not do anything if set
-        # But we need to remember whether it was set so we can restore it
-        self.testflo_running = os.environ.pop('TESTFLO_RUNNING', None)
-        clear_reports()
-        set_default_reports_dir(_reports_dir)
+    t = dm.Radau(num_segments=10, order=3)
 
-        self.count = 0
+    traj = dm.Trajectory()
+    phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
+    p.model.add_subsystem('traj0', traj)
+    traj.add_phase('phase0', phase)
 
-    def tearDown(self):
-        # restore what was there before running the test
-        if self.testflo_running is not None:
-            os.environ['TESTFLO_RUNNING'] = self.testflo_running
+    # p.model.add_subsystem('traj0', traj)
 
-    @hooks_active
-    def test_no_sim_reports(self):
-        setup_default_reports()
+    phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
 
-        p = om.Problem(model=om.Group())
+    phase.add_state('x', fix_initial=True, fix_final=False)
+    phase.add_state('y', fix_initial=True, fix_final=False)
 
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring(tol=1.0E-12)
+    # Note that by omitting the targets here Dymos will automatically attempt to connect
+    # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found.
+    phase.add_state('v', fix_initial=True, fix_final=False)
 
-        t = dm.Radau(num_segments=10, order=3)
+    phase.add_control('theta', continuity=True, rate_continuity=True,
+                      units='deg', lower=0.01, upper=179.9)
 
-        traj = dm.Trajectory()
-        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
-        p.model.add_subsystem('traj0', traj)
-        traj.add_phase('phase0', phase)
+    phase.add_parameter('g', targets=['g'], units='m/s**2')
 
-        # p.model.add_subsystem('traj0', traj)
+    phase.add_boundary_constraint('x', loc='final', equals=10)
+    phase.add_boundary_constraint('y', loc='final', equals=5)
+    # Minimize time at the end of the phase
+    phase.add_objective('time_phase', loc='final', scaler=10)
 
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+    p.setup()
 
-        phase.add_state('x', fix_initial=True, fix_final=False)
-        phase.add_state('y', fix_initial=True, fix_final=False)
+    phase.set_simulate_options(method='RK23')
 
-        # Note that by omitting the targets here Dymos will automatically attempt to connect
-        # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found.
-        phase.add_state('v', fix_initial=True, fix_final=False)
+    p['traj0.phase0.t_initial'] = 0.0
+    p['traj0.phase0.t_duration'] = 2.0
 
-        phase.add_control('theta',
-                          continuity=True, rate_continuity=True,
-                          units='deg', lower=0.01, upper=179.9)
+    p['traj0.phase0.states:x'] = phase.interp('x', [0, 10])
+    p['traj0.phase0.states:y'] = phase.interp('y', [10, 5])
+    p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9])
+    p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100])
+    p['traj0.phase0.parameters:g'] = 9.80665
 
-        phase.add_parameter('g', targets=['g'], units='m/s**2')
+    if do_reports:
+        dm.run_problem(p, run_driver=True, simulate=True, simulate_kwargs={'reports': True})
+    else:
+        dm.run_problem(p, run_driver=True, simulate=True)
 
-        phase.add_boundary_constraint('x', loc='final', equals=10)
-        phase.add_boundary_constraint('y', loc='final', equals=5)
-        # Minimize time at the end of the phase
-        phase.add_objective('time_phase', loc='final', scaler=10)
+    return p
 
-        p.setup()
 
-        phase.set_simulate_options(method='RK23')
+def setup_model_shooting(do_reports):
+    prob = om.Problem()
 
-        p['traj0.phase0.t_initial'] = 0.0
-        p['traj0.phase0.t_duration'] = 2.0
+    prob.driver = om.ScipyOptimizeDriver()
+    prob.driver.declare_coloring(tol=1.0E-12)
 
-        p['traj0.phase0.states:x'] = phase.interp('x', [0, 10])
-        p['traj0.phase0.states:y'] = phase.interp('y', [10, 5])
-        p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9])
-        p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100])
-        p['traj0.phase0.parameters:g'] = 9.80665
+    tx = dm.ExplicitShooting(num_segments=3, grid='gauss-lobatto',
+                             method='rk4', order=5,
+                             num_steps_per_segment=5,
+                             compressed=False,
+                             subprob_reports=do_reports)
 
-        dm.run_problem(p, run_driver=True, simulate=True)
+    phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
 
-        problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name)
-        report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()]
+    phase.set_time_options(units='s', fix_initial=True, duration_bounds=(1.0, 10.0))
 
-        # Test that a report subdir was made
-        self.assertEqual(len(report_subdirs), 1)
+    # automatically discover states
+    phase.set_state_options('x', fix_initial=True)
+    phase.set_state_options('y', fix_initial=True)
+    phase.set_state_options('v', fix_initial=True)
 
-        path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename)
-        self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
-        path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename)
-        self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
+    phase.add_parameter('g', val=1.0, units='m/s**2', opt=True, lower=1, upper=9.80665)
+    phase.add_control('theta', val=45.0, units='deg', opt=True, lower=1.0E-6, upper=179.9,
+                      ref=90., rate2_continuity=True)
 
-    @hooks_active
-    def test_make_sim_reports(self):
-        setup_default_reports()
+    phase.add_boundary_constraint('x', loc='final', equals=10.0)
+    phase.add_boundary_constraint('y', loc='final', equals=5.0)
 
-        p = om.Problem(model=om.Group())
+    prob.model.add_subsystem('phase0', phase)
 
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring(tol=1.0E-12)
+    phase.add_objective('time', loc='final')
 
-        t = dm.Radau(num_segments=10, order=3)
+    prob.setup(force_alloc_complex=True)
 
-        traj = dm.Trajectory()
-        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
-        p.model.add_subsystem('traj0', traj)
-        traj.add_phase('phase0', phase)
+    prob.set_val('phase0.t_initial', 0.0)
+    prob.set_val('phase0.t_duration', 2)
+    prob.set_val('phase0.states:x', 0.0)
+    prob.set_val('phase0.states:y', 10.0)
+    prob.set_val('phase0.states:v', 1.0E-6)
+    prob.set_val('phase0.parameters:g', 1.0, units='m/s**2')
+    prob.set_val('phase0.controls:theta', phase.interp('theta', ys=[0.01, 90]), units='deg')
 
-        # p.model.add_subsystem('traj0', traj)
+    dm.run_problem(prob, run_driver=True, simulate=False)
 
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+    return prob
 
-        phase.add_state('x', fix_initial=True, fix_final=False)
-        phase.add_state('y', fix_initial=True, fix_final=False)
 
-        # Note that by omitting the targets here Dymos will automatically attempt to connect
-        # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found.
-        phase.add_state('v', fix_initial=True, fix_final=False)
+# reports API between 3.18 and 3.19, so handle it here in order to be able to test against older
+# versions of openmdao
+if Version(openmdao_version) > Version("3.18"):
+    from openmdao.utils.reports_system import get_reports_dir, clear_reports
 
-        phase.add_control('theta',
-                          continuity=True, rate_continuity=True,
-                          units='deg', lower=0.01, upper=179.9)
+    @use_tempdirs
+    class TestSubproblemReportToggle(unittest.TestCase):
 
-        phase.add_parameter('g', targets=['g'], units='m/s**2')
+        def setUp(self):
+            self.n2_filename = _default_n2_filename
+            self.scaling_filename = _default_scaling_filename
 
-        phase.add_boundary_constraint('x', loc='final', equals=10)
-        phase.add_boundary_constraint('y', loc='final', equals=5)
-        # Minimize time at the end of the phase
-        phase.add_objective('time_phase', loc='final', scaler=10)
+            # set things to a known initial state for all the test runs
+            openmdao.core.problem._problem_names = []  # need to reset these to simulate separate runs
+            os.environ.pop('OPENMDAO_REPORTS', None)
+            os.environ.pop('OPENMDAO_REPORTS_DIR', None)
+            # We need to remove the TESTFLO_RUNNING environment variable for these tests to run.
+            # The reports code checks to see if TESTFLO_RUNNING is set and will not do anything if set
+            # But we need to remember whether it was set so we can restore it
+            self.testflo_running = os.environ.pop('TESTFLO_RUNNING', None)
+            clear_reports()
 
-        p.setup()
+            self.count = 0
 
-        phase.set_simulate_options(method='RK23')
+        def tearDown(self):
+            # restore what was there before running the test
+            if self.testflo_running is not None:
+                os.environ['TESTFLO_RUNNING'] = self.testflo_running
 
-        p['traj0.phase0.t_initial'] = 0.0
-        p['traj0.phase0.t_duration'] = 2.0
+        @hooks_active
+        def test_no_sim_reports(self):
+            p = setup_model_radau(do_reports=False)
 
-        p['traj0.phase0.states:x'] = phase.interp('x', [0, 10])
-        p['traj0.phase0.states:y'] = phase.interp('y', [10, 5])
-        p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9])
-        p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100])
-        p['traj0.phase0.parameters:g'] = 9.80665
+            report_subdirs = [e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()]
 
-        dm.run_problem(p, run_driver=True, simulate=True, simulate_kwargs={'reports': True})
+            # Test that a report subdir was made
+            self.assertEqual(len(report_subdirs), 1)
 
-        problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name)
-        report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()]
+            path = pathlib.Path(report_subdirs[0]).joinpath(self.n2_filename)
+            self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
+            path = pathlib.Path(report_subdirs[0]).joinpath(self.scaling_filename)
+            self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
 
-        # Test that a report subdir was made
-        # # There is the nominal problem, the simulation problem, and a subproblem for each segment in the simulation.
-        self.assertEqual(len(report_subdirs), 12)
+        @hooks_active
+        def test_make_sim_reports(self):
+            p = setup_model_radau(do_reports=True)
 
-        for subdir in report_subdirs:
-            path = pathlib.Path(subdir).joinpath(self.n2_filename)
-            self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
+            report_subdirs = [e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()]
 
-    @hooks_active
-    def test_explicitshooting_no_subprob_reports(self):
-        setup_default_reports()
+            # Test that a report subdir was made
+            # # There is the nominal problem, the simulation problem, and a subproblem for each segment in the simulation.
+            self.assertEqual(len(report_subdirs), 12)
 
-        prob = om.Problem()
+            for subdir in report_subdirs:
+                path = pathlib.Path(subdir).joinpath(self.n2_filename)
+                self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
 
-        prob.driver = om.ScipyOptimizeDriver()
-        prob.driver.declare_coloring(tol=1.0E-12)
+        @hooks_active
+        def test_explicitshooting_no_subprob_reports(self):
+            p = setup_model_shooting(do_reports=False)
 
-        tx = dm.ExplicitShooting(num_segments=3, grid='gauss-lobatto',
-                                 method='rk4', order=5,
-                                 num_steps_per_segment=5,
-                                 compressed=False)
+            report_subdirs = [e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()]
 
-        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+            # Test that a report subdir was made
+            self.assertEqual(len(report_subdirs), 1)
 
-        phase.set_time_options(units='s', fix_initial=True, duration_bounds=(1.0, 10.0))
+            path = pathlib.Path(report_subdirs[0]).joinpath(self.n2_filename)
+            self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
+            path = pathlib.Path(report_subdirs[0]).joinpath(self.scaling_filename)
+            self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
 
-        # automatically discover states
-        phase.set_state_options('x', fix_initial=True)
-        phase.set_state_options('y', fix_initial=True)
-        phase.set_state_options('v', fix_initial=True)
+        @hooks_active
+        def test_explicitshooting_make_subprob_reports(self):
+            p = setup_model_shooting(do_reports=True)
 
-        phase.add_parameter('g', val=1.0, units='m/s**2', opt=True, lower=1, upper=9.80665)
-        phase.add_control('theta', val=45.0, units='deg', opt=True, lower=1.0E-6, upper=179.9,
-                          ref=90., rate2_continuity=True)
+            report_subdirs = [e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()]
 
-        phase.add_boundary_constraint('x', loc='final', equals=10.0)
-        phase.add_boundary_constraint('y', loc='final', equals=5.0)
+            # Test that a report subdir was made
+            # There is the nominal problem, a subproblem for integration, and a subproblem for the derivatives.
+            self.assertEqual(len(report_subdirs), 3)
 
-        prob.model.add_subsystem('phase0', phase)
+            path = pathlib.Path(report_subdirs[0]).joinpath(self.n2_filename)
+            self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
+            path = pathlib.Path(report_subdirs[0]).joinpath(self.scaling_filename)
+            self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
 
-        phase.add_objective('time', loc='final')
+            for subdir in report_subdirs:
+                path = pathlib.Path(subdir).joinpath(self.n2_filename)
+                self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
 
-        prob.setup(force_alloc_complex=True)
+else:  # old OM versions before reports API changed...
+    from openmdao.utils.reports_system import set_default_reports_dir, _reports_dir, clear_reports, \
+        setup_default_reports
 
-        prob.set_val('phase0.t_initial', 0.0)
-        prob.set_val('phase0.t_duration', 2)
-        prob.set_val('phase0.states:x', 0.0)
-        prob.set_val('phase0.states:y', 10.0)
-        prob.set_val('phase0.states:v', 1.0E-6)
-        prob.set_val('phase0.parameters:g', 1.0, units='m/s**2')
-        prob.set_val('phase0.controls:theta', phase.interp('theta', ys=[0.01, 90]), units='deg')
+    @use_tempdirs
+    class TestSubproblemReportToggle(unittest.TestCase):
 
-        dm.run_problem(prob, run_driver=True, simulate=False)
+        def setUp(self):
+            self.n2_filename = _default_n2_filename
+            self.scaling_filename = _default_scaling_filename
 
-        problem_reports_dir = pathlib.Path(_reports_dir).joinpath(prob._name)
-        report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()]
+            # set things to a known initial state for all the test runs
+            openmdao.core.problem._problem_names = []  # need to reset these to simulate separate runs
+            os.environ.pop('OPENMDAO_REPORTS', None)
+            os.environ.pop('OPENMDAO_REPORTS_DIR', None)
+            # We need to remove the TESTFLO_RUNNING environment variable for these tests to run.
+            # The reports code checks to see if TESTFLO_RUNNING is set and will not do anything if set
+            # But we need to remember whether it was set so we can restore it
+            self.testflo_running = os.environ.pop('TESTFLO_RUNNING', None)
+            clear_reports()
+            set_default_reports_dir(_reports_dir)
 
-        # Test that a report subdir was made
-        self.assertEqual(len(report_subdirs), 1)
+            self.count = 0
 
-        path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename)
-        self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
-        path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename)
-        self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
+        def tearDown(self):
+            # restore what was there before running the test
+            if self.testflo_running is not None:
+                os.environ['TESTFLO_RUNNING'] = self.testflo_running
 
-    @hooks_active
-    def test_explicitshooting_make_subprob_reports(self):
-        setup_default_reports()
+        @hooks_active
+        def test_no_sim_reports(self):
+            setup_default_reports()
 
-        prob = om.Problem()
+            p = setup_model_radau(do_reports=False)
 
-        prob.driver = om.ScipyOptimizeDriver()
-        prob.driver.declare_coloring(tol=1.0E-12)
+            problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name)
+            report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()]
 
-        tx = dm.ExplicitShooting(num_segments=3, grid='gauss-lobatto',
-                                 method='rk4', order=5,
-                                 num_steps_per_segment=5,
-                                 compressed=False,
-                                 subprob_reports=True)
+            # Test that a report subdir was made
+            self.assertEqual(len(report_subdirs), 1)
 
-        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+            path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename)
+            self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
+            path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename)
+            self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
 
-        phase.set_time_options(units='s', fix_initial=True, duration_bounds=(1.0, 10.0))
+        @hooks_active
+        def test_make_sim_reports(self):
+            setup_default_reports()
 
-        # automatically discover states
-        phase.set_state_options('x', fix_initial=True)
-        phase.set_state_options('y', fix_initial=True)
-        phase.set_state_options('v', fix_initial=True)
+            p = setup_model_radau(do_reports=True)
 
-        phase.add_parameter('g', val=1.0, units='m/s**2', opt=True, lower=1, upper=9.80665)
-        phase.add_control('theta', val=45.0, units='deg', opt=True, lower=1.0E-6, upper=179.9,
-                          ref=90., rate2_continuity=True)
+            report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()]
 
-        phase.add_boundary_constraint('x', loc='final', equals=10.0)
-        phase.add_boundary_constraint('y', loc='final', equals=5.0)
+            # Test that a report subdir was made
+            # # There is the nominal problem, the simulation problem, and a subproblem for each segment in the simulation.
+            self.assertEqual(len(report_subdirs), 12)
 
-        prob.model.add_subsystem('phase0', phase)
+            for subdir in report_subdirs:
+                path = pathlib.Path(subdir).joinpath(self.n2_filename)
+                self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
 
-        phase.add_objective('time', loc='final')
+        @hooks_active
+        def test_explicitshooting_no_subprob_reports(self):
+            setup_default_reports()
 
-        prob.setup(force_alloc_complex=True)
+            p = setup_model_shooting(do_reports=False)
 
-        prob.set_val('phase0.t_initial', 0.0)
-        prob.set_val('phase0.t_duration', 2)
-        prob.set_val('phase0.states:x', 0.0)
-        prob.set_val('phase0.states:y', 10.0)
-        prob.set_val('phase0.states:v', 1.0E-6)
-        prob.set_val('phase0.parameters:g', 1.0, units='m/s**2')
-        prob.set_val('phase0.controls:theta', phase.interp('theta', ys=[0.01, 90]), units='deg')
+            problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name)
+            report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()]
 
-        dm.run_problem(prob, run_driver=True, simulate=False)
+            # Test that a report subdir was made
+            self.assertEqual(len(report_subdirs), 1)
 
-        problem_reports_dir = pathlib.Path(_reports_dir).joinpath(prob._name)
-        report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()]
+            path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename)
+            self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
+            path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename)
+            self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
 
-        # Test that a report subdir was made
-        # There is the nominal problem, a subproblem for integration, and a subproblem for the derivatives.
-        self.assertEqual(len(report_subdirs), 3)
+        @hooks_active
+        def test_explicitshooting_make_subprob_reports(self):
+            setup_default_reports()
 
-        path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename)
-        self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
-        path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename)
-        self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
+            p = setup_model_shooting(do_reports=True)
 
-        for subdir in report_subdirs:
-            path = pathlib.Path(subdir).joinpath(self.n2_filename)
+            problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name)
+            report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()]
+
+            # Test that a report subdir was made
+            # There is the nominal problem, a subproblem for integration, and a subproblem for the derivatives.
+            self.assertEqual(len(report_subdirs), 3)
+
+            path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename)
             self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
+            path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename)
+            self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found')
+
+            for subdir in report_subdirs:
+                path = pathlib.Path(subdir).joinpath(self.n2_filename)
+                self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found')
diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py
index 3e9c6b0fa..516ec16c9 100644
--- a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py
+++ b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py
@@ -129,7 +129,7 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP
     elif optimizer == 'IPOPT':
         p.driver.opt_settings['hessian_approximation'] = 'limited-memory'
         p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
-        p.driver.opt_settings['print_level'] = 4
+        p.driver.opt_settings['print_level'] = 0
         p.driver.opt_settings['linear_solver'] = 'mumps'
         p.driver.opt_settings['mu_strategy'] = 'adaptive'
         # p.driver.opt_settings['derivative_test'] = 'first-order'
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 356a5239f..29815fb6c 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
@@ -38,8 +38,8 @@ def test_ex_two_burn_orbit_raise_connected(self):
         sim_case2 = om.CaseReader('dymos_simulation.db').get_case('final')
 
         # Verify that the second case has the same inputs and outputs
-        assert_cases_equal(case1, p, tol=1.0E-9)
-        assert_cases_equal(sim_case1, sim_case2, tol=1.0E-9)
+        assert_cases_equal(case1, p, tol=1.0E-8)
+        assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8)
 
     def test_restart_from_solution_radau(self):
         optimizer = 'IPOPT'
@@ -94,8 +94,8 @@ def test_ex_two_burn_orbit_raise_connected(self):
         sim_case2 = om.CaseReader('dymos_simulation.db').get_case('final')
 
         # Verify that the second case has the same inputs and outputs
-        assert_cases_equal(case1, p, tol=1.0E-9)
-        assert_cases_equal(sim_case1, sim_case2, tol=1.0E-9)
+        assert_cases_equal(case1, p, tol=1.0E-8)
+        assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8)
 
     def test_restart_from_solution_radau_to_connected(self):
         optimizer = 'IPOPT'
diff --git a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py
index 9a69a655c..6c4298bdd 100644
--- a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py
+++ b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py
@@ -42,7 +42,7 @@ def make_problem(self, transcription=GaussLobatto, optimizer='SLSQP', numseg=30)
         elif optimizer == 'IPOPT':
             p.driver.opt_settings['hessian_approximation'] = 'limited-memory'
             # p.driver.opt_settings['nlp_scaling_method'] = 'user-scaling'
-            p.driver.opt_settings['print_level'] = 4
+            p.driver.opt_settings['print_level'] = 0
             p.driver.opt_settings['linear_solver'] = 'mumps'
             p.driver.declare_coloring()
         else:
diff --git a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py
index 29d2eab31..52c8b2e70 100644
--- a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py
+++ b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py
@@ -27,7 +27,7 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto',
         p.driver.opt_settings['Major step limit'] = 0.5
     elif optimizer == 'IPOPT':
         p.driver.opt_settings['tol'] = 1.0E-5
-        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['print_level'] = 0
         p.driver.opt_settings['mu_strategy'] = 'monotone'
         p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
         p.driver.opt_settings['mu_init'] = 0.01
diff --git a/dymos/examples/racecar/doc/test_doc_racecar.py b/dymos/examples/racecar/doc/test_doc_racecar.py
index 32ad1af3f..a7261c24b 100644
--- a/dymos/examples/racecar/doc/test_doc_racecar.py
+++ b/dymos/examples/racecar/doc/test_doc_racecar.py
@@ -137,7 +137,7 @@ def test_racecar_for_docs(self):
         p.driver.opt_settings['compl_inf_tol'] = 1e-3
         p.driver.opt_settings['acceptable_iter'] = 0
         p.driver.opt_settings['tol'] = 1e-3
-        p.driver.opt_settings['print_level'] = 5
+        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'
diff --git a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py
index 2a79813a3..992f34740 100644
--- a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py
+++ b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py
@@ -34,7 +34,7 @@ def make_problem(self, transcription=Radau, optimizer='SLSQP', numseg=30):
             p.driver.opt_settings['Verify level'] = 3
         elif optimizer == 'IPOPT':
             p.driver.opt_settings['max_iter'] = 500
-            p.driver.opt_settings['print_level'] = 4
+            p.driver.opt_settings['print_level'] = 0
             p.driver.opt_settings['tol'] = 1.0E-6
             p.driver.opt_settings['acceptable_tol'] = 1.0E-5
         traj = p.model.add_subsystem('traj', Trajectory())
diff --git a/dymos/examples/robot_arm/test/test_robot_arm.py b/dymos/examples/robot_arm/test/test_robot_arm.py
index 84438e6dd..69a561189 100644
--- a/dymos/examples/robot_arm/test/test_robot_arm.py
+++ b/dymos/examples/robot_arm/test/test_robot_arm.py
@@ -38,7 +38,7 @@ def make_problem(self, transcription=Radau, optimizer='SLSQP', numseg=30):
         elif optimizer == 'IPOPT':
             p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
             p.driver.opt_settings['max_iter'] = 500
-            p.driver.opt_settings['print_level'] = 4
+            p.driver.opt_settings['print_level'] = 0
         traj = p.model.add_subsystem('traj', Trajectory())
 
         phase = traj.add_phase('phase', Phase(ode_class=RobotArmODE,
diff --git a/dymos/examples/shuttle_reentry/test/test_reentry.py b/dymos/examples/shuttle_reentry/test/test_reentry.py
index 2b0a3b7bc..76d9a0e6f 100644
--- a/dymos/examples/shuttle_reentry/test/test_reentry.py
+++ b/dymos/examples/shuttle_reentry/test/test_reentry.py
@@ -26,7 +26,7 @@ def make_problem(self, constrained=True, transcription=GaussLobatto, optimizer='
         if optimizer == 'IPOPT':
             p.driver.opt_settings['max_iter'] = 500
             p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
-            p.driver.opt_settings['print_level'] = 5
+            p.driver.opt_settings['print_level'] = 0
             p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
             p.driver.opt_settings['tol'] = 1.0E-7
             p.driver.opt_settings['mu_strategy'] = 'monotone'
@@ -152,7 +152,7 @@ def test_reentry_mixed_controls(self):
         p.driver.declare_coloring(tol=1.0E-12)
         p.driver.options['optimizer'] = 'IPOPT'
         p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
-        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['print_level'] = 0
         p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
         p.driver.opt_settings['mu_strategy'] = 'monotone'
 
diff --git a/dymos/examples/vanderpol/vanderpol_dymos.py b/dymos/examples/vanderpol/vanderpol_dymos.py
index c9acfa75a..34a5e2572 100644
--- a/dymos/examples/vanderpol/vanderpol_dymos.py
+++ b/dymos/examples/vanderpol/vanderpol_dymos.py
@@ -20,7 +20,7 @@ def vanderpol(transcription='gauss-lobatto', num_segments=40, transcription_orde
         if optimizer == 'SNOPT':
             p.driver.opt_settings['iSumm'] = 6  # show detailed SNOPT output
         elif optimizer == 'IPOPT':
-            p.driver.opt_settings['print_level'] = 4
+            p.driver.opt_settings['print_level'] = 0
     p.driver.declare_coloring()
 
     # define a Trajectory object and add to model
diff --git a/dymos/examples/water_rocket/doc/test_doc_water_rocket.py b/dymos/examples/water_rocket/doc/test_doc_water_rocket.py
index db405f857..41fefd39e 100644
--- a/dymos/examples/water_rocket/doc/test_doc_water_rocket.py
+++ b/dymos/examples/water_rocket/doc/test_doc_water_rocket.py
@@ -25,7 +25,7 @@ def test_water_rocket_height_for_docs(self):
         traj = p.model.add_subsystem('traj', traj)
 
         p.driver = om.pyOptSparseDriver(optimizer='IPOPT', print_results=False)
-        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['print_level'] = 0
         p.driver.opt_settings['max_iter'] = 1000
         p.driver.opt_settings['mu_strategy'] = 'monotone'
         p.driver.declare_coloring(tol=1.0E-12)
@@ -64,7 +64,7 @@ def test_water_rocket_range_for_docs(self):
         traj = p.model.add_subsystem('traj', traj)
 
         p.driver = om.pyOptSparseDriver(optimizer='IPOPT')
-        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['print_level'] = 0
         p.driver.opt_settings['max_iter'] = 1000
         p.driver.opt_settings['mu_strategy'] = 'monotone'
         p.driver.declare_coloring(tol=1.0E-12)
diff --git a/dymos/test/test_run_problem.py b/dymos/test/test_run_problem.py
index 779747211..5f4067bc8 100755
--- a/dymos/test/test_run_problem.py
+++ b/dymos/test/test_run_problem.py
@@ -35,7 +35,7 @@ def test_run_HS_problem_radau(self):
         elif optimizer == 'IPOPT':
             p.driver.opt_settings['hessian_approximation'] = 'limited-memory'
             # p.driver.opt_settings['nlp_scaling_method'] = 'user-scaling'
-            p.driver.opt_settings['print_level'] = 4
+            p.driver.opt_settings['print_level'] = 0
             p.driver.opt_settings['max_iter'] = 200
             p.driver.opt_settings['linear_solver'] = 'mumps'
 
@@ -102,7 +102,7 @@ def test_run_HS_problem_gl(self):
         elif optimizer == 'IPOPT':
             p.driver.opt_settings['hessian_approximation'] = 'limited-memory'
             # p.driver.opt_settings['nlp_scaling_method'] = 'user-scaling'
-            p.driver.opt_settings['print_level'] = 4
+            p.driver.opt_settings['print_level'] = 0
             p.driver.opt_settings['linear_solver'] = 'mumps'
 
         traj = p.model.add_subsystem('traj', dm.Trajectory())
diff --git a/dymos/trajectory/test/test_trajectory_parameters.py b/dymos/trajectory/test/test_trajectory_parameters.py
index 89446dc9b..76487855d 100644
--- a/dymos/trajectory/test/test_trajectory_parameters.py
+++ b/dymos/trajectory/test/test_trajectory_parameters.py
@@ -175,7 +175,7 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP
     elif optimizer == 'IPOPT':
         p.driver.opt_settings['hessian_approximation'] = 'limited-memory'
         p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
-        p.driver.opt_settings['print_level'] = 4
+        p.driver.opt_settings['print_level'] = 0
         p.driver.opt_settings['linear_solver'] = 'mumps'
         p.driver.opt_settings['mu_strategy'] = 'adaptive'
         # p.driver.opt_settings['derivative_test'] = 'first-order'