Repository: OpenMDAO/dymos

Issue:
`phase.load_case` should use the new `phase.set_xxx_vals` API
### Description

`phase.load_case` was not using the set_xxx_vals API.
As a result, under the Birkhoff transcription, the initial and final state values were not being properly set.

### Example

N/A

### Dymos Version

1.10.1-dev

### Relevant environment information

_No response_

Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/68912280/badge)](https://coveralls.io/builds/68912280)

coverage: 92.532% (-0.2%) from 92.757%
when pulling **73a6363672317bd0530b9f5d5e2dd5267ef16572 on robfalck:load_case_set_state_vals**
into **d61cd9792489e4cbd428d9f5b9fbe4799a75797f on OpenMDAO:master**.


Patch:
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 58047247..26088aaf 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1,7 +1,6 @@
-from collections.abc import Iterable, Callable, Sequence
+from collections.abc import Iterable, Callable
 from copy import deepcopy
 import inspect
-from os import path
 import warnings
 
 import numpy as np
@@ -29,7 +28,7 @@ from ..utils.indexing import get_constraint_flat_idxs
 from ..utils.introspection import configure_time_introspection, _configure_constraint_introspection, \
     configure_controls_introspection, configure_parameters_introspection, \
     configure_timeseries_output_introspection, classify_var, configure_timeseries_expr_introspection
-from ..utils.misc import _unspecified
+from ..utils.misc import _unspecified, create_subprob
 from ..utils.lgl import lgl
 
 
@@ -58,6 +57,12 @@ class Phase(om.Group):
         A phase instance from which the initialized phase should copy its data.
     **kwargs : dict
         Dictionary of optional phase arguments.
+
+    Attributes
+    ----------
+    sim_prob : Problem or None
+        The OpenMDAO problem used for trajectory simulation.
+        This is None unless the simulate method has been called.
     """
     def __init__(self, from_phase=None, **kwargs):
         _kwargs = kwargs.copy()
@@ -69,10 +74,6 @@ class Phase(om.Group):
 
         # These are the options which will be set at setup time.
         # Prior to setup, the options are placed into the user_*_options dictionaries.
-        # self.time_options = TimeOptionsDictionary()
-        # self.state_options = {}
-        # self.control_options = {}
-        # self.parameter_options = {}
         self.refine_options = GridRefinementOptionsDictionary()
         self.simulate_options = SimulateOptionsDictionary()
         self.timeseries_ec_vars = {}
@@ -87,6 +88,7 @@ class Phase(om.Group):
                                            'subset': 'all',
                                            'outputs': {}}}
         self._objectives = {}
+        self.sim_prob = None
 
         super(Phase, self).__init__(**_kwargs)
 
@@ -2699,7 +2701,7 @@ class Phase(om.Group):
             prob.set_val(prob_path, val)
 
     def simulate(self, times_per_seg=None, method=_unspecified, atol=_unspecified, rtol=_unspecified,
-                 first_step=_unspecified, max_step=_unspecified, record_file=None):
+                 first_step=_unspecified, max_step=_unspecified, record_file=None, reports=False):
         """
         Simulate the Phase using scipy.integrate.solve_ivp.
 
@@ -2721,6 +2723,8 @@ class Phase(om.Group):
         record_file : str or None
             If a string, the file to which the result of the simulation will be saved.
             If None, no record of the simulation will be saved.
+        reports : bool or None or str or Sequence
+            Reports setting for the subproblems run under simualate.
 
         Returns
         -------
@@ -2729,7 +2733,9 @@ class Phase(om.Group):
             can be interrogated to obtain timeseries outputs in the same manner as other Phases
             to obtain results at the requested times.
         """
-        sim_prob = om.Problem(model=om.Group())
+        self.sim_prob = sim_prob = create_subprob(base_name=f'{self.name}_simulation',
+                                                  comm=self.comm,
+                                                  reports=reports)
 
         sim_phase = self.get_simulation_phase(times_per_seg=times_per_seg, method=method, atol=atol, rtol=rtol,
                                               first_step=first_step, max_step=max_step)
@@ -2740,7 +2746,10 @@ class Phase(om.Group):
             rec = om.SqliteRecorder(record_file)
             sim_prob.add_recorder(rec)
 
-        sim_prob.setup(check=True)
+        if om_version <= (3, 42, 2):
+            sim_prob.setup(check=True)
+        else:
+            sim_prob.setup(check=True, parent=self)
         sim_prob.final_setup()
 
         sim_phase.set_vals_from_phase(from_phase=self)
@@ -3064,8 +3073,7 @@ class Phase(om.Group):
             prev_time_val = np.flip(prev_time_val, axis=0)
             unique_idxs = np.flip(unique_idxs, axis=0)
 
-        self.set_val('t_initial', t_initial, units=prev_time_units)
-        self.set_val('t_duration', t_duration, units=prev_time_units)
+        self.set_time_val(initial=t_initial, duration=t_duration, units=prev_time_units)
 
         # Interpolate the timeseries state outputs from the previous solution onto the new grid.
         if not isinstance(self, dm.AnalyticPhase):
@@ -3081,19 +3089,12 @@ class Phase(om.Group):
 
                 prev_state_val = prev_vars[prev_state_path]['val']
                 prev_state_units = prev_vars[prev_state_path]['units']
-                interp_vals = self.interp(name=state_name,
-                                          xs=prev_time_val,
-                                          ys=prev_state_val[unique_idxs],
-                                          kind='slinear')
                 if options['lower'] is not None or options['upper'] is not None:
-                    interp_vals = interp_vals.clip(options['lower'], options['upper'])
-                self.set_val(f'states:{state_name}',
-                             interp_vals,
-                             units=prev_state_units)
-                try:
-                    self.set_val(f'initial_states:{state_name}', prev_state_val[0, ...], units=prev_state_units)
-                except KeyError:
-                    pass
+                    prev_state_val = prev_state_val.clip(options['lower'], options['upper'])
+                self.set_state_val(state_name,
+                                   vals=prev_state_val[unique_idxs],
+                                   time_vals=prev_time_val,
+                                   units=prev_state_units)
 
                 if options['fix_final']:
                     warning_message = f"{phase_name}.states:{state_name} specifies 'fix_final=True'. " \
@@ -3114,13 +3115,12 @@ class Phase(om.Group):
 
                 prev_control_val = prev_vars[prev_control_path]['val']
                 prev_control_units = prev_vars[prev_control_path]['units']
-                interp_vals = self.interp(name=control_name,
-                                          xs=prev_time_val,
-                                          ys=prev_control_val[unique_idxs],
-                                          kind='slinear')
                 if options['lower'] is not None or options['upper'] is not None:
-                    interp_vals = interp_vals.clip(options['lower'], options['upper'])
-                self.set_val(f'controls:{control_name}', interp_vals, units=prev_control_units)
+                    prev_control_val = prev_control_val.clip(options['lower'], options['upper'])
+                self.set_control_val(control_name,
+                                     vals=prev_control_val[unique_idxs],
+                                     time_vals=prev_time_val,
+                                     units=prev_control_units)
                 if options['fix_final']:
                     warning_message = f"{phase_name}.controls:{control_name} specifies 'fix_final=True'. " \
                                       f"If the given restart file has a" \
@@ -3132,7 +3132,7 @@ class Phase(om.Group):
             if f'{prev_phase_prom_path}.parameter_vals:{param_name}' in prev_vars:
                 prev_param_val = prev_vars[f'{prev_phase_prom_path}.parameter_vals:{param_name}']['val']
                 prev_param_units = prev_vars[f'{prev_phase_prom_path}.parameter_vals:{param_name}']['units']
-                self.set_val(f'parameters:{param_name}', prev_param_val[0, ...], units=prev_param_units)
+                self.set_parameter_val(param_name, prev_param_val[0, ...], units=prev_param_units)
             else:
                 issue_warning(f'Unable to find "{prev_phase_prom_path}.parameter_vals:{param_name}" '
                               f'in data from case being loaded.')
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 4060f5a4..f54c093c 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -11,6 +11,7 @@ from openmdao.utils.units import unit_conversion
 import numpy as np
 import networkx as nx
 
+import openmdao
 import openmdao.api as om
 from openmdao.utils.mpi import MPI
 
@@ -21,11 +22,14 @@ from .phase_linkage_comp import PhaseLinkageComp
 from ..phase.analytic_phase import AnalyticPhase
 from ..phase.options import TrajParameterOptionsDictionary
 from ..transcriptions.common import ParameterComp
-from ..utils.misc import get_rate_units, _unspecified, _none_or_unspecified
+from ..utils.misc import create_subprob, get_rate_units, _unspecified, _none_or_unspecified
 from ..utils.introspection import get_promoted_vars, get_source_metadata, _get_common_metadata
 from .._options import options as dymos_options
 
 
+om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
+
+
 class Trajectory(om.Group):
     """
     Class definition for a Trajectory group.
@@ -44,6 +48,9 @@ class Trajectory(om.Group):
         A dictionary of parameter names and their associated TrajectoryParameterOptionsDictionary
     phases : om.Group or om.ParallelGroup
         The Group which contains phases for this Trajectory.
+    sim_prob : Problem or None
+        The OpenMDAO problem used for trajectory simulation.
+        This is None unless the simulate method has been called.
 
     _linkages : OrderedDict
         A dictionary containing phase linkage information for the Trajectory.
@@ -59,6 +66,7 @@ class Trajectory(om.Group):
         self._phases = {}
         self._phase_graph = nx.DiGraph()
         self._has_connected_phases = False
+        self.sim_prob = None
 
         self.phases = om.ParallelGroup() if self.options['parallel_phases'] else om.Group()
 
@@ -1455,7 +1463,6 @@ class Trajectory(om.Group):
         for name, phs in self._phases.items():
             if phs.simulate_options is None:
                 continue
-
             sim_phs = phs.get_simulation_phase(times_per_seg=times_per_seg, method=method,
                                                atol=atol, rtol=rtol, first_step=first_step,
                                                max_step=max_step, reports=reports)
@@ -1466,7 +1473,9 @@ class Trajectory(om.Group):
 
         sim_traj.parameter_options.update(self.parameter_options)
 
-        sim_prob = om.Problem(model=om.Group(), reports=reports, comm=self.comm)
+        self.sim_prob = sim_prob = create_subprob(base_name=f'{self.name}_simulation',
+                                                  comm=self.comm,
+                                                  reports=reports)
 
         traj_name = self.name if self.name else 'sim_traj'
         sim_prob.model.add_subsystem(traj_name, sim_traj)
@@ -1483,7 +1492,10 @@ class Trajectory(om.Group):
             # fault of the user.
             warnings.filterwarnings(action='ignore', category=om.UnusedOptionWarning)
             warnings.filterwarnings(action='ignore', category=om.SetupWarning)
-            sim_prob.setup()
+            if om_version <= (3, 42, 2):
+                sim_prob.setup(check=True)
+            else:
+                sim_prob.setup(check=True, parent=self)
             sim_prob.final_setup()
 
         # Assign trajectory parameter values
diff --git a/dymos/transcriptions/explicit_shooting/ode_integration_comp.py b/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
index 7f3984d0..d53f605f 100644
--- a/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
+++ b/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
@@ -1,10 +1,16 @@
 import numpy as np
-import openmdao.api as om
 from scipy.integrate import solve_ivp
 
+import openmdao
+import openmdao.api as om
+
 from ..._options import options as dymos_options
 
 from .ode_evaluation_group import ODEEvaluationGroup
+from dymos.utils.misc import create_subprob
+
+
+om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
 
 
 class ODEIntegrationComp(om.ExplicitComponent):
@@ -93,7 +99,10 @@ class ODEIntegrationComp(om.ExplicitComponent):
                              'transcription where the number of nodes per segment can exceed 20 to 30.')
 
     def _setup_subprob(self):
-        self._eval_subprob = p = om.Problem(comm=self.comm, reports=self._reports)
+        self._eval_subprob = p = create_subprob(base_name=f'{self.pathname}_subprob',
+                                                comm=self.comm,
+                                                reports=self._reports)
+
         p.model.add_subsystem('ode_eval',
                               ODEEvaluationGroup(ode_class=self.options['ode_class'],
                                                  time_options=self.time_options,
@@ -107,7 +116,10 @@ class ODEIntegrationComp(om.ExplicitComponent):
                               promotes_inputs=['*'],
                               promotes_outputs=['*'])
 
-        p.setup()
+        if om_version <= (3, 34, 2):
+            p.setup()
+        else:
+            p.setup(parent=self)
         p.final_setup()
 
     def _set_segment_index(self, idx):
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index 750874a6..e360fa5e 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -686,7 +686,7 @@ class PseudospectralBase(TranscriptionBase):
         if np.isscalar(vals):
             interp_vals = vals
         else:
-            interp_vals = phase.interp(name, vals, time_vals,
+            interp_vals = phase.interp(name, ys=vals, xs=time_vals,
                                        nodes='state_input',
                                        kind=interpolation_kind)
         input_data[f'states:{name}'] = interp_vals
diff --git a/dymos/utils/misc.py b/dymos/utils/misc.py
index fd55f736..60240ba5 100644
--- a/dymos/utils/misc.py
+++ b/dymos/utils/misc.py
@@ -2,9 +2,11 @@ from collections.abc import Iterable
 
 import numpy as np
 
+import openmdao.api as om
+from openmdao.core.constants import _ReprClass
+
 from .constants import INF_BOUND
 from .indexing import get_desvar_indices
-from openmdao.core.constants import _ReprClass
 
 
 # unique object to check if default is given (when None is an allowed value)
@@ -212,3 +214,44 @@ def GroupWrapperConfig(comp_class, config_io_args=None):
             self.configure_io(*args)
 
     return WrappedClass
+
+
+def create_subprob(base_name, comm, reports=False):
+    """
+    Create a new problem using basename possibly appended with unique identifiers if name collisions occur.
+
+    Parameters
+    ----------
+    base_name : str
+        The base name of the problem. This may be appended by `_{int}` to obtain a unique problem name.
+        In the event of running under MPI, an 8-character hash may further append the name to ensure
+        it is unique.
+    comm : comm
+        The MPI comm to be used by the subproblem.
+    reports : bool or None or str or Sequence
+        Reports setting for the subproblems run under simualate.
+
+    Returns
+    -------
+    Problem
+        The instantiated OpenMDAO problem instance.
+    """
+    from openmdao.core.problem import _problem_names
+
+    # Find a unique sim problem name. This mostly causes problems
+    # when many simulations are being run in a single process, as in testing.
+    i = 0
+    sim_prob_name = f'{base_name}_{i}'
+    while sim_prob_name in _problem_names:
+        i += 1
+        sim_prob_name = f'{base_name}_{i}'
+
+    try:
+        p = om.Problem(comm=comm, reports=reports, name=sim_prob_name)
+    except ValueError:
+        # Testing under MPI, we still might have name collisions. In that case, add a random hash
+        # to the end of the problem name.
+        import hashlib
+        str_hash = hashlib.sha256(used_for_security=False)[:8]
+        p = om.Problem(comm=comm, reports=reports, name=f'{sim_prob_name}_{str_hash}')
+    return p


Test patch:
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py b/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
index a01a8ad0..e4899487 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
@@ -4,21 +4,24 @@ import pathlib
 import os
 from packaging.version import Version
 
+import openmdao
 import openmdao.api as om
 import openmdao.core.problem
 from openmdao.utils.testing_utils import use_tempdirs
 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
 
 
 import dymos as dm
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
 
 
-def setup_model_radau(do_reports):
-    p = om.Problem(model=om.Group())
+om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
+
+
+def setup_model_radau(do_reports, probname):
+    p = om.Problem(model=om.Group(), name=probname)
 
     p.driver = om.ScipyOptimizeDriver()
     p.driver.declare_coloring(tol=1.0E-12)
@@ -64,16 +67,13 @@ def setup_model_radau(do_reports):
     phase.set_control_val('theta', [5, 100.5])
     phase.set_parameter_val('g', 9.80665)
 
-    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)
+    dm.run_problem(p, run_driver=True, simulate=True, simulate_kwargs={'reports': do_reports})
 
     return p
 
 
-def setup_model_shooting(do_reports):
-    prob = om.Problem()
+def setup_model_shooting(do_reports, probname):
+    prob = om.Problem(name=probname)
 
     prob.driver = om.ScipyOptimizeDriver()
     prob.driver.declare_coloring(tol=1.0E-12)
@@ -112,198 +112,97 @@ def setup_model_shooting(do_reports):
     phase.set_control_val('theta', [0.01, 90], units='deg')
     phase.set_parameter_val('g', 1.0)
 
-    dm.run_problem(prob, run_driver=True, simulate=False)
+    dm.run_problem(prob, run_driver=True, simulate=True)
 
     return prob
 
 
-# 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
-
-    @use_tempdirs
-    class TestSubproblemReportToggle(unittest.TestCase):
-
-        def setUp(self):
-            self.n2_filename = _default_n2_filename
-            self.scaling_filename = _default_scaling_filename
-
-            # 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()
-
-            self.count = 0
-
-        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_no_sim_reports(self):
-            p = setup_model_radau(do_reports=False)
-
-            report_subdirs = sorted([e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()])
-
-            # Test that a report subdir was made
-            self.assertEqual(len(report_subdirs), 1)
-
-            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')
-
-        @hooks_active
-        def test_make_sim_reports(self):
-            p = setup_model_radau(do_reports=True)
-
-            report_subdirs = sorted([e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()])
-
-            # Test that a report subdir was made
-            # There is the nominal problem, the simulation problem, and a subproblem for the simulation.
-            self.assertEqual(len(report_subdirs), 3)
-
-            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')
-
-        @hooks_active
-        def test_explicitshooting_no_subprob_reports(self):
-            p = setup_model_shooting(do_reports=False)
-
-            report_subdirs = sorted([e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()])
-
-            # Test that a report subdir was made
-            self.assertEqual(len(report_subdirs), 1)
-
-            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')
-
-        @hooks_active
-        def test_explicitshooting_make_subprob_reports(self):
-            p = setup_model_shooting(do_reports=True)
-
-            report_subdirs = sorted([e for e in pathlib.Path(get_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), 2)
-
-            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')
-
-            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')
-
-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
-
-    @use_tempdirs
-    class TestSubproblemReportToggle(unittest.TestCase):
+@use_tempdirs
+class TestSubproblemReportToggle(unittest.TestCase):
 
-        def setUp(self):
-            self.n2_filename = _default_n2_filename
-            self.scaling_filename = _default_scaling_filename
+    def setUp(self):
+        self.n2_filename = _default_n2_filename
+        self.scaling_filename = _default_scaling_filename
 
-            # 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)
+        # 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)
 
-            self.count = 0
+        self.count = 0
 
-        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
+    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_no_sim_reports(self):
-            setup_default_reports()
+    @unittest.skipIf(om_version <= (3, 34, 2), 'Requires OpenMDAO version later than 3.34.2')
+    @hooks_active
+    def test_no_sim_reports(self):
+        p = setup_model_radau(do_reports=False, probname='test_no_sim_reports')
 
-            p = setup_model_radau(do_reports=False)
+        main_outputs_dir = p.get_outputs_dir()
+        main_reports_dir = p.get_reports_dir()
 
-            problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name)
-            report_subdirs = sorted([e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()])
+        sim_outputs_dir = main_outputs_dir / 'traj0_simulation_out'
+        sim_reports_dir = sim_outputs_dir / 'reports'
 
-            # Test that a report subdir was made
-            self.assertEqual(len(report_subdirs), 1)
+        self.assertFalse(sim_reports_dir.exists())
 
-            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')
+    @unittest.skipIf(om_version <= (3, 34, 2), 'Requires OpenMDAO version later than 3.34.2')
+    @hooks_active
+    def test_make_sim_reports(self):
+        p = setup_model_radau(do_reports=True, probname='test_make_sim_reports')
 
-        @hooks_active
-        def test_make_sim_reports(self):
-            setup_default_reports()
+        main_reports_dir = p.get_reports_dir()
 
-            p = setup_model_radau(do_reports=True)
+        traj = p.model._get_subsystem('traj0')
+        sim_reports_dir = traj.sim_prob.get_reports_dir()
 
-            report_subdirs = sorted([e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()])
+        self.assertTrue((main_reports_dir / self.n2_filename).exists())
+        self.assertTrue(sim_reports_dir.exists())
+        self.assertTrue((sim_reports_dir / self.n2_filename).exists())
 
-            # 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)
+    @unittest.skipIf(om_version <= (3, 34, 2), 'Requires OpenMDAO version later than 3.34.2')
+    @hooks_active
+    def test_explicitshooting_no_subprob_reports(self):
+        p = setup_model_shooting(do_reports=False,
+                                 probname='test_explicitshooting_no_subprob_reports')
 
-            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')
+        main_reports_dir = p.get_reports_dir()
+        subprob_reports_dir = p.model.phase0.integrator._eval_subprob.get_reports_dir()
 
-        @hooks_active
-        def test_explicitshooting_no_subprob_reports(self):
-            setup_default_reports()
+        main_reports = os.listdir(main_reports_dir)
 
-            p = setup_model_shooting(do_reports=False)
+        self.assertFalse(subprob_reports_dir.exists())
 
-            problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name)
-            report_subdirs = sorted([e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()])
+        self.assertIn(self.n2_filename, main_reports)
+        self.assertIn(self.scaling_filename, main_reports)
 
-            # Test that a report subdir was made
-            self.assertEqual(len(report_subdirs), 1)
+    @unittest.skipIf(om_version <= (3, 34, 2), 'Requires OpenMDAO version later than 3.34.2')
+    @hooks_active
+    def test_explicitshooting_make_subprob_reports(self):
+        p = setup_model_shooting(do_reports=True,
+                                 probname='test_explicitshooting_make_subprob_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')
+        main_reports_dir = p.get_reports_dir()
+        subprob_reports_dir = p.model.phase0.integrator._eval_subprob.get_reports_dir()
 
-        @hooks_active
-        def test_explicitshooting_make_subprob_reports(self):
-            setup_default_reports()
+        main_reports = os.listdir(main_reports_dir)
+        subprob_reports = os.listdir(subprob_reports_dir)
 
-            p = setup_model_shooting(do_reports=True)
+        self.assertIn(self.n2_filename, main_reports)
+        self.assertIn(self.n2_filename, subprob_reports)
 
-            problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name)
-            report_subdirs = sorted([e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()])
+        self.assertIn(self.scaling_filename, main_reports)
 
-            # Test that a report subdir was made
-            # There is the nominal problem and a subproblem for integration
-            self.assertEqual(len(report_subdirs), 2)
+        # The subprob has no optimization, so should not have a scaling report
+        self.assertNotIn(self.scaling_filename, subprob_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')
 
-            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')
+if __name__ == '__main__':
+    unittest.main()
diff --git a/dymos/test/test_load_case.py b/dymos/test/test_load_case.py
index 7cbab5fe..1f07393a 100644
--- a/dymos/test/test_load_case.py
+++ b/dymos/test/test_load_case.py
@@ -184,6 +184,76 @@ class TestLoadCase(unittest.TestCase):
                           q.model.phase0.interp(xs=time_q, ys=theta_q, nodes='all'),
                           tolerance=1.0E-2)
 
+    def test_load_case_radau_to_birkhoff(self):
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+        import dymos as dm
+
+        p = setup_problem(dm.Radau(num_segments=20))
+
+        # Solve for the optimal trajectory
+        dm.run_problem(p)
+
+        # Load the solution
+        case = om.CaseReader('dymos_solution.db').get_case('final')
+
+        # create a problem with a different transcription with a different number of variables
+        q = setup_problem(dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=50)))
+
+        # Fill q with junk so that we can be sure load_case worked
+        q['phase0.t_initial'] = -88
+        q['phase0.t_duration'] = 88
+
+        q['phase0.states:x'] = -88
+        q['phase0.states:y'] = -88
+        q['phase0.states:v'] = -88
+
+        q['phase0.initial_states:x'] = -88
+        q['phase0.initial_states:y'] = -88
+        q['phase0.initial_states:v'] = -88
+
+        q['phase0.final_states:x'] = -88
+        q['phase0.final_states:y'] = -88
+        q['phase0.final_states:v'] = -88
+
+        # Load the values from the previous solution
+        q.load_case(case)
+
+        # Run the model to ensure we find the same output values as those that we recorded
+        q.run_model()
+
+        time_p = case.get_val('phase0.timeseries.time')
+        theta_p = case.get_val('phase0.timeseries.theta')
+
+        time_q = q.get_val('phase0.timeseries.time')
+        theta_q = q.get_val('phase0.timeseries.theta')
+
+        x_p = case.get_val('phase0.timeseries.x')
+        y_p = case.get_val('phase0.timeseries.y')
+        v_p = case.get_val('phase0.timeseries.v')
+
+        x0_q = q.get_val('phase0.initial_states:x')
+        xf_q = q.get_val('phase0.final_states:x')
+
+        y0_q = q.get_val('phase0.initial_states:y')
+        yf_q = q.get_val('phase0.final_states:y')
+
+        v0_q = q.get_val('phase0.initial_states:v')
+        vf_q = q.get_val('phase0.final_states:v')
+
+        assert_near_equal(q.model.phase0.interp(xs=time_p, ys=theta_p, nodes='all'),
+                          q.model.phase0.interp(xs=time_q, ys=theta_q, nodes='all'),
+                          tolerance=1.0E-2)
+
+        assert_near_equal(x_p[0, ...], x0_q, tolerance=1.0E-5)
+        assert_near_equal(x_p[-1, ...], xf_q, tolerance=1.0E-5)
+
+        assert_near_equal(y_p[0, ...], y0_q, tolerance=1.0E-5)
+        assert_near_equal(y_p[-1, ...], yf_q, tolerance=1.0E-5)
+
+        assert_near_equal(v_p[0, ...], v0_q, tolerance=1.0E-5)
+        assert_near_equal(v_p[-1, ...], vf_q, tolerance=1.0E-5)
+
     def test_load_case_warn_fix_final_states(self):
         import openmdao.api as om
         from openmdao.utils.assert_utils import assert_warnings