Repository: OpenMDAO/dymos
Pull request #1106 (2024-08-20 19:21:02)

bokeh timeseries plots now comply with OpenMDAO's placement of recorder files.
### Summary

The latest OpenMDAO change puts recorder files in the problem outputs directory by default.  This change makes dymos' timeseries plots handle that correctly.

### Related Issues

- Resolves #

### Backwards incompatibilities

None

### New Dependencies

None


review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/OpenMDAO/dymos/pull/1106"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>

Patch:
diff --git a/.github/workflows/dymos_docs_workflow.yml b/.github/workflows/dymos_docs_workflow.yml
index a924e72f..f1aabc43 100644
--- a/.github/workflows/dymos_docs_workflow.yml
+++ b/.github/workflows/dymos_docs_workflow.yml
@@ -31,7 +31,7 @@ jobs:
             SCIPY: '1.13'
             PETSc: '3.19'
             PYOPTSPARSE: 'v2.11.0'
-            OPENMDAO: 'latest'
+            OPENMDAO: 'dev'
             OPTIONAL: '[docs]'
             JAX: '0.4.28'
             PUBLISH_DOCS: 1
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 3bf8b319..33924e04 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -6,3 +6,11 @@ repos:
     # Run the linter.
     - id: ruff
       types_or: [ python, pyi, jupyter ]
+-   repo: local
+    hooks:
+    - id: check-notebook-outputs
+      name: Check that notebooks have no output cells.
+      description: Verifies that documentation notebooks contain no output cells.
+      language: python
+      pass_filenames: false
+      entry: python ./docs/clean_notebooks.py --dryrun ./docs
\ No newline at end of file
diff --git a/benchmark/benchmark_balanced_field.py b/benchmark/benchmark_balanced_field.py
index 4b4ff108..6b0f25a2 100644
--- a/benchmark/benchmark_balanced_field.py
+++ b/benchmark/benchmark_balanced_field.py
@@ -262,8 +262,8 @@ class BenchmarkBalancedFieldLength(unittest.TestCase):
                                            solvesegs=False)
 
     def benchmark_radau_timeseries_nosim_nosolveseg(self):
-            _run_balanced_field_length_problem(tx=dm.Radau, timeseries=True, sim=False,
-                                               solvesegs=False)
+        _run_balanced_field_length_problem(tx=dm.Radau, timeseries=True, sim=False,
+                                           solvesegs=False)
 
     def benchmark_radau_notimeseries_nosim_solveseg(self):
         _run_balanced_field_length_problem(tx=dm.Radau, timeseries=False, sim=False,
diff --git a/benchmark/benchmark_brachistochrone.py b/benchmark/benchmark_brachistochrone.py
index 692680a8..87c9c510 100644
--- a/benchmark/benchmark_brachistochrone.py
+++ b/benchmark/benchmark_brachistochrone.py
@@ -14,7 +14,7 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     p = om.Problem(model=om.Group())
 
     p.driver = om.pyOptSparseDriver()
-    OPT, OPTIMIZER = set_pyoptsparse_opt(optimizer, fallback=True)
+    _, OPTIMIZER = set_pyoptsparse_opt(optimizer, fallback=True)
     p.driver.options['optimizer'] = OPTIMIZER
 
     if simul_derivs:
diff --git a/benchmark/benchmark_racecar.py b/benchmark/benchmark_racecar.py
index 4c230e73..413c7c72 100644
--- a/benchmark/benchmark_racecar.py
+++ b/benchmark/benchmark_racecar.py
@@ -7,7 +7,6 @@ import openmdao.api as om
 import dymos as dm
 
 from dymos.examples.racecar.combinedODE import CombinedODE
-from dymos.examples.racecar.spline import get_spline, get_track_points
 from dymos.examples.racecar.tracks import ovaltrack  # track curvature imports
 
 
@@ -17,10 +16,7 @@ def _run_racecar_problem(transcription, timeseries=False, make_plots=False):
 
     # generate nodes along the centerline for curvature calculation (different
     # than collocation nodes)
-    points = get_track_points(track)
 
-    # fit the centerline spline.
-    finespline, gates, gatesd, curv, slope = get_spline(points, s=0.0)
     # by default 10000 points
     s_final = track.get_total_length()
 
diff --git a/docs/build_source_docs.py b/docs/build_source_docs.py
index f7033ab6..189e41cb 100644
--- a/docs/build_source_docs.py
+++ b/docs/build_source_docs.py
@@ -15,7 +15,6 @@ packages = [
     'transcriptions.common',
     'transcriptions.explicit_shooting',
     'transcriptions.pseudospectral',
-    'transcriptions.solve_ivp',
     'utils',
 ]
 
diff --git a/docs/clean_notebooks.py b/docs/clean_notebooks.py
new file mode 100644
index 00000000..5c034399
--- /dev/null
+++ b/docs/clean_notebooks.py
@@ -0,0 +1,78 @@
+#!/usr/bin/env python
+
+import argparse
+import json
+import os
+import pathlib
+import subprocess
+
+
+EXCLUDE_DIRS = ('_build', '.ipynb_checkpoints')
+
+
+class bcolors:
+    FILE = '\033[96m'  # blue
+    CYAN = '\033[96m'
+    GREEN = '\033[92m'
+    WARNING = '\033[93m'
+    FAIL = '\033[91m'
+    ENDC = '\033[0m'
+    BOLD = '\033[1m'
+    UNDERLINE = '\033[4m'
+
+
+def _clear_notebook_outputs(path=os.getcwd(), dry_run=True):
+
+    # Find all the notebooks
+    notebooks = []
+    if pathlib.Path(path).is_dir():
+        for dirpath, dirs, files in os.walk(path, topdown=True):
+            dirs[:] = [d for d in dirs if d not in EXCLUDE_DIRS]
+            notebooks.extend(pathlib.Path(dirpath).glob('*.ipynb'))
+        if not notebooks:
+            print(f'{bcolors.FAIL}No notebooks found.{bcolors.ENDC}')
+            return
+    elif path.endswith('.ipynb'):
+        notebooks.append(path)
+    else:
+        print(f'{bcolors.FILE}{path}{bcolors.ENDC}{bcolors.FAIL} is not a notebook file.{bcolors.ENDC}')
+        return
+
+    num_cleaned = 0
+    for file in notebooks:
+        with open(file) as f:
+            json_data = json.load(f)
+            for i in json_data['cells']:
+                if 'execution_count' in i and i['execution_count'] is not None:
+                    break
+            else:
+                # Notebook clean, process next notebook.
+                continue
+        num_cleaned += 1
+        if dry_run:
+            print(f'Would clear outputs from {bcolors.FILE}{file}{bcolors.ENDC}. '
+                  f'{bcolors.WARNING}(dryrun = True){bcolors.ENDC}')
+        else:
+            print(f'Clearing {bcolors.FILE}{file}{bcolors.ENDC}...', end='')
+            subprocess.Popen(f"jupyter nbconvert --clear-output --inplace {file}",
+                             shell=True,
+                             stdout=subprocess.DEVNULL,
+                             stderr=subprocess.STDOUT).wait()
+            print('done')
+
+    if num_cleaned == 0:
+        print('No unclean notebooks found.')
+
+    return num_cleaned
+
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Clear notebook outputs in place.')
+    parser.add_argument('path', nargs='?', default='.', help='The or directory path to clean.')
+    parser.add_argument('-d', '--dryrun', action='store_true',
+                        help='Print notebooks with outputs but do not clean them.')
+    args = parser.parse_args()
+    num_cleaned = _clear_notebook_outputs(path=args.path, dry_run=args.dryrun)
+    if num_cleaned >0 and  args.dryrun:
+        exit(1)
diff --git a/docs/dymos_book/examples/balanced_field/balanced_field.ipynb b/docs/dymos_book/examples/balanced_field/balanced_field.ipynb
index 8c2842fb..878cfb79 100644
--- a/docs/dymos_book/examples/balanced_field/balanced_field.ipynb
+++ b/docs/dymos_book/examples/balanced_field/balanced_field.ipynb
@@ -720,7 +720,8 @@
    },
    "outputs": [],
    "source": [
-    "sim_case = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sim_prob = traj.sim_prob\n",
+    "sim_case = om.CaseReader(sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "fig, axes = plt.subplots(2, 1, sharex=True, gridspec_kw={'top': 0.92}, figsize=(12,6))\n",
     "for phase in ['br_to_v1', 'rto', 'v1_to_vr', 'rotate', 'climb']:\n",
@@ -813,7 +814,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb
index 42d06fa6..829c43ec 100644
--- a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb
+++ b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb
@@ -458,13 +458,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim_case = om.CaseReader('dymos_solution.db').get_case('final')\n",
+    "sol_case = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
     "\n",
     "fig, axes = plt.subplots(2, 1, sharex=True, gridspec_kw={'top': 0.92}, figsize=(12, 6))\n",
     "for phase in ['br_to_v1', 'rto', 'v1_to_vr', 'rotate', 'climb']:\n",
-    "    r = sim_case.get_val(f'traj.{phase}.timeseries.r', units='ft')\n",
-    "    v = sim_case.get_val(f'traj.{phase}.timeseries.v', units='kn')\n",
-    "    t = sim_case.get_val(f'traj.{phase}.timeseries.time', units='s')\n",
+    "    r = sol_case.get_val(f'traj.{phase}.timeseries.r', units='ft')\n",
+    "    v = sol_case.get_val(f'traj.{phase}.timeseries.v', units='kn')\n",
+    "    t = sol_case.get_val(f'traj.{phase}.timeseries.time', units='s')\n",
     "    axes[0].plot(t, r, '-', label=phase)\n",
     "    axes[1].plot(t, v, '-', label=phase)\n",
     "fig.suptitle('Balanced Field Length')\n",
@@ -474,11 +474,11 @@
     "axes[0].grid(True)\n",
     "axes[1].grid(True)\n",
     "\n",
-    "tv1 = sim_case.get_val('traj.br_to_v1.timeseries.time', units='s')[-1, 0]\n",
-    "v1 = sim_case.get_val('traj.br_to_v1.timeseries.v', units='kn')[-1, 0]\n",
+    "tv1 = sol_case.get_val('traj.br_to_v1.timeseries.time', units='s')[-1, 0]\n",
+    "v1 = sol_case.get_val('traj.br_to_v1.timeseries.v', units='kn')[-1, 0]\n",
     "\n",
-    "tf_rto = sim_case.get_val('traj.rto.timeseries.time', units='s')[-1, 0]\n",
-    "rf_rto = sim_case.get_val('traj.rto.timeseries.r', units='ft')[-1, 0]\n",
+    "tf_rto = sol_case.get_val('traj.rto.timeseries.time', units='s')[-1, 0]\n",
+    "rf_rto = sol_case.get_val('traj.rto.timeseries.r', units='ft')[-1, 0]\n",
     "\n",
     "axes[0].annotate(f'field length = {r[-1, 0]:5.1f} ft', xy=(t[-1, 0], r[-1, 0]),\n",
     "                 xycoords='data', xytext=(0.7, 0.5),\n",
@@ -551,7 +551,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/docs/dymos_book/examples/brachistochrone/brachistochrone_tandem_phases.ipynb b/docs/dymos_book/examples/brachistochrone/brachistochrone_tandem_phases.ipynb
index 37266ae9..7b3bc58f 100644
--- a/docs/dymos_book/examples/brachistochrone/brachistochrone_tandem_phases.ipynb
+++ b/docs/dymos_book/examples/brachistochrone/brachistochrone_tandem_phases.ipynb
@@ -320,7 +320,11 @@
     "i = np.array(0)\n",
     "legend_contents = []\n",
     "\n",
-    "sol_case = om.CaseReader('dymos_solution.db').get_case('final')\n",
+    "print(p.get_outputs_dir())\n",
+    "import os\n",
+    "print(os.listdir(os.getcwd()))\n",
+    "\n",
+    "sol_case = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
     "\n",
     "sol_x = sol_case.get_val('phase0.timeseries.x')\n",
     "sol_y = sol_case.get_val('phase0.timeseries.y')\n",
@@ -363,13 +367,6 @@
     "\n",
     "HTML(html_with_custom_height)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_controls.ipynb b/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_controls.ipynb
index 82d24ce0..a1e6870e 100644
--- a/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_controls.ipynb
+++ b/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_controls.ipynb
@@ -77,19 +77,19 @@
    ]
   },
   {
-    "cell_type": "code",
-    "execution_count": null,
-    "metadata": {
-     "tags": [
-      "active-ipynb",
-      "remove-input",
-      "remove-output"
-     ]
-    },
-    "outputs": [],
-    "source": [
-     "%matplotlib inline"
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
     ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
   },
   {
    "cell_type": "code",
@@ -193,8 +193,7 @@
     "# Test the results\n",
     "print(p.get_val('traj.phase0.timeseries.time')[-1])\n",
     "\n",
-    "# Check the validity of our results by using scipy.integrate.solve_ivp to\n",
-    "# integrate the solution.\n",
+    "# Check the validity of our results.\n",
     "sim_out = traj.simulate()\n",
     "\n",
     "# Plot the results\n",
@@ -273,7 +272,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   },
   "orphan": true
  },
diff --git a/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_init_duration_states.ipynb b/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_init_duration_states.ipynb
index 488d759a..5e719dab 100644
--- a/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_init_duration_states.ipynb
+++ b/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_init_duration_states.ipynb
@@ -61,19 +61,19 @@
    ]
   },
   {
-    "cell_type": "code",
-    "execution_count": null,
-    "metadata": {
-     "tags": [
-      "active-ipynb",
-      "remove-input",
-      "remove-output"
-     ]
-    },
-    "outputs": [],
-    "source": [
-     "%matplotlib inline"
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
     ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
   },
   {
    "cell_type": "code",
@@ -195,10 +195,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Check the validity of our results by using scipy.integrate.solve_ivp to\n",
-    "# integrate the solution.\n",
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "# Check the validity of our results.\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "# Plot the results\n",
     "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 4.5))\n",
@@ -259,7 +258,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/docs/dymos_book/examples/bryson_denham/bryson_denham.ipynb b/docs/dymos_book/examples/bryson_denham/bryson_denham.ipynb
index 9d50450e..0356c9f5 100644
--- a/docs/dymos_book/examples/bryson_denham/bryson_denham.ipynb
+++ b/docs/dymos_book/examples/bryson_denham/bryson_denham.ipynb
@@ -461,8 +461,8 @@
    },
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "t = sol.get_val('traj.phase0.timeseries.time')\n",
     "x = sol.get_val('traj.phase0.timeseries.x')\n",
diff --git a/docs/dymos_book/examples/cart_pole/cart_pole.ipynb b/docs/dymos_book/examples/cart_pole/cart_pole.ipynb
index 169e1295..c74a0ccb 100644
--- a/docs/dymos_book/examples/cart_pole/cart_pole.ipynb
+++ b/docs/dymos_book/examples/cart_pole/cart_pole.ipynb
@@ -329,8 +329,8 @@
     "print('objective value:', obj)\n",
     "\n",
     "# get optimization solution and simulation (time integration) results\n",
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "# plot time histories of x, x_dot, theta, theta_dot\n",
     "plot_results([('traj.phase.timeseries.time', 'traj.phase.timeseries.x', 'time (s)', 'x (m)'),\n",
diff --git a/docs/dymos_book/examples/double_integrator/double_integrator.ipynb b/docs/dymos_book/examples/double_integrator/double_integrator.ipynb
index 52432c1f..16013f2a 100644
--- a/docs/dymos_book/examples/double_integrator/double_integrator.ipynb
+++ b/docs/dymos_book/examples/double_integrator/double_integrator.ipynb
@@ -137,19 +137,19 @@
    ]
   },
   {
-    "cell_type": "code",
-    "execution_count": null,
-    "metadata": {
-     "tags": [
-      "active-ipynb",
-      "remove-input",
-      "remove-output"
-     ]
-    },
-    "outputs": [],
-    "source": [
-     "%matplotlib inline"
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
     ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
   },
   {
    "cell_type": "code",
@@ -221,8 +221,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.x',\n",
     "               'time (s)', 'x $(m)$'),\n",
diff --git a/docs/dymos_book/examples/finite_burn_orbit_raise/finite_burn_orbit_raise.ipynb b/docs/dymos_book/examples/finite_burn_orbit_raise/finite_burn_orbit_raise.ipynb
index a5391b4c..b87ff234 100644
--- a/docs/dymos_book/examples/finite_burn_orbit_raise/finite_burn_orbit_raise.ipynb
+++ b/docs/dymos_book/examples/finite_burn_orbit_raise/finite_burn_orbit_raise.ipynb
@@ -121,19 +121,19 @@
    ]
   },
   {
-    "cell_type": "code",
-    "execution_count": null,
-    "metadata": {
-     "tags": [
-      "active-ipynb",
-      "remove-input",
-      "remove-output"
-     ]
-    },
-    "outputs": [],
-    "source": [
-     "%matplotlib inline"
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
     ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
   },
   {
    "cell_type": "code",
@@ -307,8 +307,8 @@
    },
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "fig = plt.figure(figsize=(12, 6))\n",
     "fig.suptitle('Two Burn Orbit Raise Solution')\n",
diff --git a/docs/dymos_book/examples/hull/hull_problem.ipynb b/docs/dymos_book/examples/hull/hull_problem.ipynb
index 2f034de6..2dd17f45 100644
--- a/docs/dymos_book/examples/hull/hull_problem.ipynb
+++ b/docs/dymos_book/examples/hull/hull_problem.ipynb
@@ -447,8 +447,8 @@
    },
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "t = sol.get_val('traj.phase0.timeseries.time')\n",
     "x = sol.get_val('traj.phase0.timeseries.x')\n",
diff --git a/docs/dymos_book/examples/hypersensitive/hypersensitive.ipynb b/docs/dymos_book/examples/hypersensitive/hypersensitive.ipynb
index 1f8a9277..2365c203 100644
--- a/docs/dymos_book/examples/hypersensitive/hypersensitive.ipynb
+++ b/docs/dymos_book/examples/hypersensitive/hypersensitive.ipynb
@@ -264,8 +264,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.x',\n",
     "               'time (s)', 'x $(m)$'),\n",
diff --git a/docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb b/docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb
index f10dc27b..8b1bb5af 100644
--- a/docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb
+++ b/docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb
@@ -166,19 +166,19 @@
    ]
   },
   {
-    "cell_type": "code",
-    "execution_count": null,
-    "metadata": {
-     "tags": [
-      "active-ipynb",
-      "remove-input",
-      "remove-output"
-     ]
-    },
-    "outputs": [],
-    "source": [
-     "%matplotlib inline"
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
     ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
   },
   {
    "cell_type": "code",
@@ -295,8 +295,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.h',\n",
     "               'time (s)', 'altitude (m)'),\n",
diff --git a/docs/dymos_book/examples/mountain_car/mountain_car.ipynb b/docs/dymos_book/examples/mountain_car/mountain_car.ipynb
index 5ad084fa..429f5b3b 100644
--- a/docs/dymos_book/examples/mountain_car/mountain_car.ipynb
+++ b/docs/dymos_book/examples/mountain_car/mountain_car.ipynb
@@ -515,8 +515,8 @@
    },
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "t = sol.get_val('traj.phase0.timeseries.time')\n",
     "x = sol.get_val('traj.phase0.timeseries.x')\n",
diff --git a/docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb b/docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb
index adec8568..a9857816 100644
--- a/docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb
+++ b/docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb
@@ -142,19 +142,19 @@
    ]
   },
   {
-    "cell_type": "code",
-    "execution_count": null,
-    "metadata": {
-     "tags": [
-      "active-ipynb",
-      "remove-input",
-      "remove-output"
-     ]
-    },
-    "outputs": [],
-    "source": [
-     "%matplotlib inline"
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
     ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
   },
   {
    "cell_type": "code",
@@ -448,8 +448,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "#############################################\n",
     "# Plot the results\n",
@@ -560,7 +560,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/docs/dymos_book/examples/reentry/reentry.ipynb b/docs/dymos_book/examples/reentry/reentry.ipynb
index d2989ecf..bace185c 100644
--- a/docs/dymos_book/examples/reentry/reentry.ipynb
+++ b/docs/dymos_book/examples/reentry/reentry.ipynb
@@ -685,8 +685,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.alpha',\n",
     "               'time (s)', 'alpha (rad)'),\n",
diff --git a/docs/dymos_book/examples/ssto_earth/ssto_earth.ipynb b/docs/dymos_book/examples/ssto_earth/ssto_earth.ipynb
index f5d78dd7..df0b7eda 100644
--- a/docs/dymos_book/examples/ssto_earth/ssto_earth.ipynb
+++ b/docs/dymos_book/examples/ssto_earth/ssto_earth.ipynb
@@ -204,7 +204,7 @@
     "\n",
     "        # Setup partials\n",
     "        # Complex-step derivatives\n",
-    "        self.declare_coloring(wrt='*', method='cs', show_sparsity=True)\n",
+    "        self.declare_coloring(wrt='*', method='cs')\n",
     "\n",
     "    def compute(self, inputs, outputs):\n",
     "\n",
@@ -240,19 +240,19 @@
    ]
   },
   {
-    "cell_type": "code",
-    "execution_count": null,
-    "metadata": {
-     "tags": [
-      "active-ipynb",
-      "remove-input",
-      "remove-output"
-     ]
-    },
-    "outputs": [],
-    "source": [
-     "%matplotlib inline"
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
     ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
   },
   {
    "cell_type": "code",
@@ -345,8 +345,8 @@
    },
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "fig, [traj_ax, control_ax] = plt.subplots(nrows=2, ncols=1, figsize=(10, 8))\n",
     "\n",
diff --git a/docs/dymos_book/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb b/docs/dymos_book/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb
index ea8881ad..d7dcd3d8 100644
--- a/docs/dymos_book/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb
+++ b/docs/dymos_book/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb
@@ -369,8 +369,8 @@
    },
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "fig, [traj_ax, control_ax, param_ax] = plt.subplots(nrows=3, ncols=1, figsize=(10, 8))\n",
     "\n",
diff --git a/docs/dymos_book/examples/ssto_moon_polynomial_controls/ssto_moon_polynomial_controls.ipynb b/docs/dymos_book/examples/ssto_moon_polynomial_controls/ssto_moon_polynomial_controls.ipynb
index 830d0838..54ce7aaf 100644
--- a/docs/dymos_book/examples/ssto_moon_polynomial_controls/ssto_moon_polynomial_controls.ipynb
+++ b/docs/dymos_book/examples/ssto_moon_polynomial_controls/ssto_moon_polynomial_controls.ipynb
@@ -353,8 +353,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "fig, [traj_ax, control_ax, param_ax] = plt.subplots(nrows=3, ncols=1, figsize=(10, 8))\n",
     "\n",
@@ -459,7 +459,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/docs/dymos_book/examples/vanderpol/vanderpol.ipynb b/docs/dymos_book/examples/vanderpol/vanderpol.ipynb
index f1b8db32..06aabf37 100644
--- a/docs/dymos_book/examples/vanderpol/vanderpol.ipynb
+++ b/docs/dymos_book/examples/vanderpol/vanderpol.ipynb
@@ -394,35 +394,34 @@
     "from dymos.examples.plotting import plot_results\n",
     "\n",
     "# Display the results\n",
-    "def vanderpol_plots():\n",
-    "    sol = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "    sim = om.CaseReader('dymos_simulation.db').get_case('final')\n",
-    "    \n",
-    "    plot_results([('traj.phase0.timeseries.time',\n",
-    "                   'traj.phase0.timeseries.x1',\n",
-    "                   'time (s)',\n",
-    "                   '$x_1$ (V)'),\n",
-    "                 ('traj.phase0.timeseries.time',\n",
-    "                  'traj.phase0.timeseries.x0',\n",
-    "                  'time (s)',\n",
-    "                  '$x_0$ (V/s)'),\n",
-    "                  ('traj.phase0.timeseries.time',\n",
-    "                   'traj.phase0.timeseries.J',\n",
-    "                   'time (s)',\n",
-    "                   'J'),\n",
-    "                  ('traj.phase0.timeseries.x0',\n",
-    "                   'traj.phase0.timeseries.x1',\n",
-    "                   '$x_0$ (V/s)',\n",
-    "                   '$x_1$ (V)'),\n",
-    "                 ('traj.phase0.timeseries.time',\n",
-    "                  'traj.phase0.timeseries.u',\n",
-    "                  'time (s)',\n",
-    "                  'control u'),\n",
-    "                  ],\n",
-    "                 title='Van Der Pol Simulation',\n",
-    "                 p_sol=sol, p_sim=sim)\n",
-    "\n",
-    "vanderpol_plots()"
+    "\n",
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim_prob = p.model.traj.sim_prob\n",
+    "sim = om.CaseReader(sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
+    "\n",
+    "plot_results([('traj.phase0.timeseries.time',\n",
+    "                'traj.phase0.timeseries.x1',\n",
+    "                'time (s)',\n",
+    "                '$x_1$ (V)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "              'traj.phase0.timeseries.x0',\n",
+    "              'time (s)',\n",
+    "              '$x_0$ (V/s)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "                'traj.phase0.timeseries.J',\n",
+    "                'time (s)',\n",
+    "                'J'),\n",
+    "              ('traj.phase0.timeseries.x0',\n",
+    "                'traj.phase0.timeseries.x1',\n",
+    "                '$x_0$ (V/s)',\n",
+    "                '$x_1$ (V)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "              'traj.phase0.timeseries.u',\n",
+    "              'time (s)',\n",
+    "              'control u'),\n",
+    "              ],\n",
+    "              title='Van Der Pol Simulation',\n",
+    "              p_sol=sol, p_sim=sim)\n"
    ]
   },
   {
@@ -470,7 +469,33 @@
    },
    "outputs": [],
    "source": [
-    "vanderpol_plots()"
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim_prob = p.model.traj.sim_prob\n",
+    "sim = om.CaseReader(sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
+    "\n",
+    "plot_results([('traj.phase0.timeseries.time',\n",
+    "                'traj.phase0.timeseries.x1',\n",
+    "                'time (s)',\n",
+    "                '$x_1$ (V)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "              'traj.phase0.timeseries.x0',\n",
+    "              'time (s)',\n",
+    "              '$x_0$ (V/s)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "                'traj.phase0.timeseries.J',\n",
+    "                'time (s)',\n",
+    "                'J'),\n",
+    "              ('traj.phase0.timeseries.x0',\n",
+    "                'traj.phase0.timeseries.x1',\n",
+    "                '$x_0$ (V/s)',\n",
+    "                '$x_1$ (V)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "              'traj.phase0.timeseries.u',\n",
+    "              'time (s)',\n",
+    "              'control u'),\n",
+    "              ],\n",
+    "              title='Van Der Pol Simulation',\n",
+    "              p_sol=sol, p_sim=sim)"
    ]
   },
   {
@@ -516,7 +541,33 @@
    },
    "outputs": [],
    "source": [
-    "vanderpol_plots()"
+    "sol = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim_prob = p.model.traj.sim_prob\n",
+    "sim = om.CaseReader(sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')\n",
+    "\n",
+    "plot_results([('traj.phase0.timeseries.time',\n",
+    "                'traj.phase0.timeseries.x1',\n",
+    "                'time (s)',\n",
+    "                '$x_1$ (V)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "              'traj.phase0.timeseries.x0',\n",
+    "              'time (s)',\n",
+    "              '$x_0$ (V/s)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "                'traj.phase0.timeseries.J',\n",
+    "                'time (s)',\n",
+    "                'J'),\n",
+    "              ('traj.phase0.timeseries.x0',\n",
+    "                'traj.phase0.timeseries.x1',\n",
+    "                '$x_0$ (V/s)',\n",
+    "                '$x_1$ (V)'),\n",
+    "              ('traj.phase0.timeseries.time',\n",
+    "              'traj.phase0.timeseries.u',\n",
+    "              'time (s)',\n",
+    "              'control u'),\n",
+    "              ],\n",
+    "              title='Van Der Pol Simulation',\n",
+    "              p_sol=sol, p_sim=sim)"
    ]
   },
   {
@@ -556,7 +607,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/docs/dymos_book/examples/water_rocket/water_rocket.ipynb b/docs/dymos_book/examples/water_rocket/water_rocket.ipynb
index 8339d289..80756dc4 100644
--- a/docs/dymos_book/examples/water_rocket/water_rocket.ipynb
+++ b/docs/dymos_book/examples/water_rocket/water_rocket.ipynb
@@ -618,8 +618,8 @@
     "for key, entry in summary.items():\n",
     "    print(f'{key}: {entry.value:6.4f} {entry.unit}')\n",
     "    \n",
-    "sol_out = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim_out = om.CaseReader('dymos_simulation.db').get_case('final')"
+    "sol_out = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim_out = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')"
    ]
   },
   {
@@ -735,8 +735,8 @@
     "for key, entry in summary.items():\n",
     "    print(f'{key}: {entry.value:6.4f} {entry.unit}')\n",
     "\n",
-    "sol_out = om.CaseReader('dymos_solution.db').get_case('final')\n",
-    "sim_out = om.CaseReader('dymos_simulation.db').get_case('final')"
+    "sol_out = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
+    "sim_out = om.CaseReader(traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db').get_case('final')"
    ]
   },
   {
diff --git a/docs/dymos_book/faq/connect_scalar_parameters_to_ode.ipynb b/docs/dymos_book/faq/connect_scalar_parameters_to_ode.ipynb
index 1d1895e6..ffcf2328 100644
--- a/docs/dymos_book/faq/connect_scalar_parameters_to_ode.ipynb
+++ b/docs/dymos_book/faq/connect_scalar_parameters_to_ode.ipynb
@@ -91,19 +91,19 @@
    ]
   },
   {
-    "cell_type": "code",
-    "execution_count": null,
-    "metadata": {
-     "tags": [
-      "active-ipynb",
-      "remove-input",
-      "remove-output"
-     ]
-    },
-    "outputs": [],
-    "source": [
-     "%matplotlib inline"
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
     ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
   },
   {
    "cell_type": "code",
@@ -198,7 +198,8 @@
     "dm.run_problem(p, simulate=True)\n",
     "\n",
     "# Generate the explicitly simulated trajectory\n",
-    "exp_out = om.CaseReader('dymos_simulation.db').get_case('final')\n",
+    "sim_prob_dir = traj.sim_prob.get_outputs_dir()\n",
+    "exp_out = om.CaseReader(sim_prob_dir / 'dymos_simulation.db').get_case('final')\n",
     "\n",
     "# Extract the timeseries from the implicit solution and the explicit simulation\n",
     "x = p.get_val('traj.phase0.timeseries.x')\n",
@@ -270,7 +271,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.3"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
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 e296b300..a94e4e79 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
@@ -4,6 +4,7 @@ import openmdao.api as om
 
 import dymos as dm
 from dymos.examples.finite_burn_orbit_raise.finite_burn_eom import FiniteBurnODE
+from dymos.utils.misc import om_version
 
 
 def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=False,
@@ -307,8 +308,14 @@ if __name__ == '__main__':  # pragma: no cover
                                      compressed=False, optimizer=optimizer, simulate=True,
                                      connected=CONNECTED, show_output=False)
 
-    sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-    sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+    sol_db = 'dymos_solution.db'
+    sim_db = 'dymos_simulation.db'
+    if om_version()[0] > (3, 34, 2):
+        sol_db = p.get_outputs_dir() / sol_db
+        sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+    sol_case = om.CaseReader(sol_db).get_case('final')
+    sim_case = om.CaseReader(sim_db).get_case('final')
 
     # The last phase in this case is run in reverse time if CONNECTED=True,
     # so grab the correct index to test the resulting delta-V.
diff --git a/dymos/grid_refinement/refinement.py b/dymos/grid_refinement/refinement.py
index eb0b5ad4..348496dd 100644
--- a/dymos/grid_refinement/refinement.py
+++ b/dymos/grid_refinement/refinement.py
@@ -7,12 +7,11 @@ from .write_iteration import write_error, write_refine_iter
 
 from dymos.grid_refinement.error_estimation import check_error
 from dymos.load_case import load_case, find_phases
+from dymos.utils.misc import om_version
 
 import numpy as np
 import sys
 
-import openmdao
-
 
 def _refine_iter(problem, refine_iteration_limit=0, refine_method='hp', case_prefix=None, reset_iter_counts=True):
     """
@@ -62,8 +61,7 @@ def _refine_iter(problem, refine_iteration_limit=0, refine_method='hp', case_pre
                 for stream in f, sys.stdout:
                     write_refine_iter(stream, i, phases, refine_results)
 
-                om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
-                if om_version < (3, 27, 1):
+                if om_version()[0] < (3, 27, 1):
                     prev_soln = {'inputs': problem.model.list_inputs(out_stream=None, units=True, prom_name=True),
                                  'outputs': problem.model.list_outputs(out_stream=None, units=True, prom_name=True)}
 
@@ -85,7 +83,7 @@ def _refine_iter(problem, refine_iteration_limit=0, refine_method='hp', case_pre
                 failed = problem.run_driver(case_prefix=f'{_case_prefix}{refine_method}_{i}_')
 
             for stream in [f, sys.stdout]:
-                if i == refine_iteration_limit-1:
+                if i == refine_iteration_limit - 1:
                     print('Iteration limit exceeded. Unable to satisfy specified tolerance', file=stream)
                     failed = True
                 else:
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 5b24f21f..acf2964d 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -7,7 +7,6 @@ import numpy as np
 
 from scipy import interpolate
 
-import openmdao
 import openmdao.api as om
 from openmdao.utils.mpi import MPI
 from openmdao.utils.om_warnings import issue_warning
@@ -28,14 +27,10 @@ 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, create_subprob
+from ..utils.misc import _unspecified, create_subprob, om_version
 from ..utils.lgl import lgl
 
 
-om_dev_version = openmdao.__version__.endswith('dev')
-om_version = tuple(int(s) for s in openmdao.__version__.split('-')[0].split('.'))
-
-
 class Phase(om.Group):
     """
     The Phase object in Dymos.
@@ -2743,7 +2738,7 @@ class Phase(om.Group):
             rec = om.SqliteRecorder(record_file)
             sim_prob.add_recorder(rec)
 
-        if om_version <= (3, 42, 2):
+        if om_version()[0] <= (3, 42, 2):
             sim_prob.setup(check=True)
         else:
             sim_prob.setup(check=True, parent=self)
diff --git a/dymos/run_problem.py b/dymos/run_problem.py
index 04e3701c..8a3fbfd2 100755
--- a/dymos/run_problem.py
+++ b/dymos/run_problem.py
@@ -1,10 +1,12 @@
+import os
+import pathlib
 import warnings
 
-import openmdao
 import openmdao.api as om
 from openmdao.recorders.case import Case
 from ._options import options as dymos_options
 from dymos.trajectory.trajectory import Trajectory
+from dymos.utils.misc import om_version
 from dymos.visualization.timeseries_plots import timeseries_plots
 
 from .grid_refinement.refinement import _refine_iter
@@ -64,7 +66,7 @@ def run_problem(problem, refine_method='hp', refine_iteration_limit=0, run_drive
         If True and model has been run previously, reset all iteration counters.
     """
     if restart is not None:
-        if isinstance(restart, str):
+        if isinstance(restart, (str, pathlib.Path)):
             case = om.CaseReader(restart).get_case('final')
         elif isinstance(restart, Case):
             case = restart
@@ -82,8 +84,7 @@ def run_problem(problem, refine_method='hp', refine_iteration_limit=0, run_drive
     problem.final_setup()
 
     if restart is not None:
-        om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
-        if om_version < (3, 27, 1):
+        if om_version()[0] < (3, 27, 1):
             from dymos.load_case import load_case
             load_case(problem, case, deprecation_warning=False)
         else:
@@ -104,6 +105,7 @@ def run_problem(problem, refine_method='hp', refine_iteration_limit=0, run_drive
     problem.record(f'{_case_prefix}final')  # save case for potential restart
     problem.cleanup()
 
+    sims = {}
     if simulate:
         _simulate_kwargs = simulate_kwargs if simulate_kwargs is not None else {}
         if 'record_file' in _simulate_kwargs:
@@ -114,17 +116,42 @@ def run_problem(problem, refine_method='hp', refine_iteration_limit=0, run_drive
                              'argument "case_prefix", not part of the simulate_kwargs dictionary.')
         for subsys in problem.model.system_iter(include_self=True, recurse=True):
             if isinstance(subsys, Trajectory):
-                subsys.simulate(record_file=simulation_record_file, case_prefix=case_prefix, **_simulate_kwargs)
+                sim_prob = subsys.simulate(record_file=simulation_record_file,
+                                           case_prefix=case_prefix,
+                                           **_simulate_kwargs)
+                sims[subsys.pathname] = sim_prob
 
     if make_plots:
+
+        if om_version()[0] > (3, 34, 2):
+            outputs_dir = problem.get_outputs_dir()
+            if os.sep in str(solution_record_file):
+                _sol_record_file = solution_record_file
+            else:
+                _sol_record_file = outputs_dir / solution_record_file
+
+            if simulate:
+                sim_outputs_dir = list(sims.values())[0].get_outputs_dir()
+                if os.sep in str(simulation_record_file):
+                    _sim_record_file = simulation_record_file
+                else:
+                    _sim_record_file = sim_outputs_dir / simulation_record_file
+            else:
+                _sim_record_file = None
+        else:
+            _sol_record_file = solution_record_file
+            _sim_record_file = None if not simulate else simulation_record_file
+
         if dymos_options['plots'] == 'bokeh':
             from dymos.visualization.timeseries.bokeh_timeseries_report import make_timeseries_report
-            make_timeseries_report(prob=problem, solution_record_file=solution_record_file,
-                                   simulation_record_file=simulation_record_file)
+            make_timeseries_report(prob=problem,
+                                   solution_record_file=_sol_record_file,
+                                   simulation_record_file=_sim_record_file)
         else:
-            _sim_record_file = None if not simulate else simulation_record_file
             _plot_kwargs = plot_kwargs if plot_kwargs is not None else {}
-            timeseries_plots(solution_record_file, simulation_record_file=_sim_record_file,
-                             plot_dir=plot_dir, problem=problem, **_plot_kwargs)
+            plots_dir = problem.get_reports_dir() / 'plots'
+            timeseries_plots(_sol_record_file,
+                             simulation_record_file=_sim_record_file,
+                             plot_dir=plots_dir, **_plot_kwargs)
 
     return failed
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index a17e09eb..6e3ce590 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -11,7 +11,6 @@ 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
 
@@ -22,13 +21,11 @@ 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 create_subprob, get_rate_units, _unspecified, _none_or_unspecified
+from ..utils.misc import create_subprob, get_rate_units, om_version, \
+    _unspecified, _none_or_unspecified
 from ..utils.introspection import get_promoted_vars, get_source_metadata, _get_common_metadata
 
 
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
-
-
 class Trajectory(om.Group):
     """
     Class definition for a Trajectory group.
@@ -1372,7 +1369,7 @@ class Trajectory(om.Group):
                         name = options['constraint_name']
                     else:
                         name = options['name']
-                    _, shape, units, linear = tx._get_objective_src(name, loc, phase, ode_outputs=ode_outputs)
+                    _, _, units, _ = tx._get_objective_src(name, loc, phase, ode_outputs=ode_outputs)
 
                     equals = options['equals']
                     lower = options['lower']
@@ -1491,10 +1488,10 @@ class Trajectory(om.Group):
             # fault of the user.
             warnings.filterwarnings(action='ignore', category=om.UnusedOptionWarning)
             warnings.filterwarnings(action='ignore', category=om.SetupWarning)
-            if om_version <= (3, 42, 2):
-                sim_prob.setup(check=True)
+            if om_version()[0] <= (3, 34, 2):
+                sim_prob.setup(check=False)
             else:
-                sim_prob.setup(check=True, parent=self)
+                sim_prob.setup(check=False, parent=self)
             sim_prob.final_setup()
 
         # Assign trajectory parameter values
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index b0d52d41..06d36c11 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -27,7 +27,8 @@ class ExplicitShooting(TranscriptionBase):
     If we view integration as a function
 
     .. math::
-        \bar{x}_{f} = \mathcal{I}(\bar{x}_0, t_0, t_d, \bar{\theta}) = \bar{x}_{0} + \int_{t_0}^{t_0+t_d} \left( f_{ode}(\bar{x}, t, \bar{\theta}) \right) dt
+        \bar{x}_{f} = \mathcal{I}(\bar{x}_0, t_0, t_d, \bar{\theta}) \\
+        \bar{x}_{f} = \bar{x}_{0} + \int_{t_0}^{t_0+t_d} \left( f_{ode}(\bar{x}, t, \bar{\theta}) \right) dt
 
     then the inputs are the initial states ($\bar{x}$), the initial time and duration ($t_0$ and $t_d$), and some set
     of parameters that impact the ODE ($\theta$). For Dymos, $\theta$ may include the phase parameters, or the node values
diff --git a/dymos/transcriptions/explicit_shooting/ode_integration_comp.py b/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
index d53f605f..882094bf 100644
--- a/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
+++ b/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
@@ -1,16 +1,12 @@
 import numpy as np
 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('.')])
+from dymos.utils.misc import create_subprob, om_version
 
 
 class ODEIntegrationComp(om.ExplicitComponent):
@@ -116,10 +112,10 @@ class ODEIntegrationComp(om.ExplicitComponent):
                               promotes_inputs=['*'],
                               promotes_outputs=['*'])
 
-        if om_version <= (3, 34, 2):
-            p.setup()
+        if om_version()[0] <= (3, 34, 2):
+            p.setup(check=None)
         else:
-            p.setup(parent=self)
+            p.setup(check=None, parent=self)
         p.final_setup()
 
     def _set_segment_index(self, idx):
@@ -250,7 +246,7 @@ class ODEIntegrationComp(om.ExplicitComponent):
                            desc=f'value for parameter {param_name}')
 
             param_size = np.prod(options['shape'], dtype=int)
-            self.parameter_idxs[param_name] = np.s_[self.p_size:self.p_size+param_size]
+            self.parameter_idxs[param_name] = np.s_[self.p_size:self.p_size + param_size]
             self.p_size += param_size
 
     def _setup_controls(self):
@@ -365,7 +361,7 @@ class ODEIntegrationComp(om.ExplicitComponent):
         for param_name, options in self.parameter_options.items():
             param_size = np.prod(options['shape'], dtype=int)
             self._parameter_idxs_in_z[param_name] = np.s_[start_z: start_z + param_size]
-            self._parameter_idxs_in_theta[param_name] = np.s_[start_theta: start_theta+param_size]
+            self._parameter_idxs_in_theta[param_name] = np.s_[start_theta: start_theta + param_size]
             z_input_names.extend([self._param_input_names[param_name]] * param_size)
             start_z += param_size
             start_theta += param_size
@@ -376,7 +372,7 @@ class ODEIntegrationComp(om.ExplicitComponent):
                 control_param_shape = (num_input_nodes,) + options['shape']
                 control_param_size = np.prod(control_param_shape, dtype=int)
                 self._control_idxs_in_z[control_name] = np.s_[start_z:start_z + control_param_size]
-                self._control_idxs_in_theta[control_name] = np.s_[start_theta:start_theta+control_param_size]
+                self._control_idxs_in_theta[control_name] = np.s_[start_theta:start_theta + control_param_size]
                 z_input_names.extend([self._control_input_names[control_name]] * control_param_size)
                 start_z += control_param_size
                 start_theta += control_param_size
@@ -384,7 +380,7 @@ class ODEIntegrationComp(om.ExplicitComponent):
                 control_param_shape = (len(control_input_node_ptau),) + options['shape']
                 control_param_size = np.prod(control_param_shape, dtype=int)
                 self._control_idxs_in_z[control_name] = np.s_[start_z:start_z + control_param_size]
-                self._control_idxs_in_theta[control_name] = np.s_[start_theta:start_theta+control_param_size]
+                self._control_idxs_in_theta[control_name] = np.s_[start_theta:start_theta + control_param_size]
                 z_input_names.extend([self._control_input_names[control_name]] * control_param_size)
                 start_z += control_param_size
                 start_theta += control_param_size
@@ -616,7 +612,7 @@ class ODEIntegrationComp(om.ExplicitComponent):
         x_dot, f_x, f_t, f_theta = self.eval_ode(x, t, theta, eval_solution=True, eval_derivs=True)
 
         dh_dz = np.zeros((1, n_z))
-        dh_dz[0, n_x+1] = 1./td
+        dh_dz[0, n_x + 1] = 1. / td
         dt_dz_dot = dh_dz
 
         dx_dz_dot = f_x @ dx_dz + f_t @ dt_dz + f_theta @ dtheta_dz + x_dot @ dt_dz_dot
@@ -759,8 +755,8 @@ class ODEIntegrationComp(om.ExplicitComponent):
                 sol = solve_ivp(self._f_augmented, t_span=t_span_seg, t_eval=t_eval_seg, y0=y0,
                                 args=(theta, self._dtheta_dz), method=method, first_step=first_step, max_step=max_step,
                                 atol=atol, rtol=rtol)
-                dx_dz_out[row_seg_i:row_seg_i+nnps[i], :] = sol.y.T[:, n_x:n_x+n_x*n_z]
-                dt_dz_out[row_seg_i:row_seg_i+nnps[i], :] = sol.y.T[:, -n_z:]
+                dx_dz_out[row_seg_i:row_seg_i + nnps[i], :] = sol.y.T[:, n_x:n_x + n_x * n_z]
+                dt_dz_out[row_seg_i:row_seg_i + nnps[i], :] = sol.y.T[:, -n_z:]
             else:
                 sol = solve_ivp(self._f_primal, t_span=t_span_seg, t_eval=t_eval_seg, y0=y0, args=(theta,),
                                 method=method, first_step=first_step, max_step=max_step, atol=atol, rtol=rtol)
@@ -768,7 +764,7 @@ class ODEIntegrationComp(om.ExplicitComponent):
             if not sol.success:
                 raise om.AnalysisError(f'solve_ivp failed: {sol.message}')
 
-            x_out[row_seg_i:row_seg_i+nnps[i], :] = sol.y.T[:, :n_x]  # Save solution to the output nodes
+            x_out[row_seg_i:row_seg_i + nnps[i], :] = sol.y.T[:, :n_x]  # Save solution to the output nodes
             t_out[row_seg_i:row_seg_i + nnps[i], 0] = sol.t
             y0 = sol.y.T[-1, :]  # Set initial y for the next segment
             row_seg_i += nnps[i]  # Increment node associated with the start of the next segment
@@ -831,8 +827,8 @@ class ODEIntegrationComp(om.ExplicitComponent):
         dt_dz = self._dt_dz_out
         dx_dz = self._dx_dz_out
 
-        partials[t_name, 't_duration'] = dt_dz[:, self.x_size+1]
-        partials[f'{t_name}_phase', 't_duration'] = dt_dz[:, self.x_size+1]
+        partials[t_name, 't_duration'] = dt_dz[:, self.x_size + 1]
+        partials[f'{t_name}_phase', 't_duration'] = dt_dz[:, self.x_size + 1]
 
         for state_name in self.state_options:
             of = self._state_output_names[state_name]
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py b/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py
index 0c6a5da7..09cde3b3 100644
--- a/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py
@@ -1,14 +1,11 @@
 import numpy as np
 
-import openmdao
 import openmdao.api as om
 
 from openmdao.utils.general_utils import ensure_compatible
 
-from ...._options import options as dymos_options
-
-
-_om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
+from dymos._options import options as dymos_options
+from dymos.utils.misc import om_version
 
 
 class BirkhoffStateResidComp(om.ImplicitComponent):
@@ -52,7 +49,7 @@ class BirkhoffStateResidComp(om.ImplicitComponent):
         self.add_input(name, **kwargs)
         self.add_residual(resid_name, **kwargs)
 
-        if _om_version > (3, 31, 1):
+        if om_version()[0] > (3, 31, 1):
             self.declare_partials(of=resid_name, wrt=name, rows=ar, cols=ar, val=1.0)
         else:
             self.declare_partials(of='*', wrt='*', method='fd')
diff --git a/dymos/transcriptions/pseudospectral/gauss_lobatto.py b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
index 52bbd7ca..6389a50f 100644
--- a/dymos/transcriptions/pseudospectral/gauss_lobatto.py
+++ b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
@@ -538,7 +538,8 @@ class GaussLobatto(PseudospectralBase):
             src_units = meta['units']
             src_tags = meta['tags']
             if 'dymos.static_output' in src_tags:
-                raise RuntimeError(f'ODE output {var} is tagged with "dymos.static_output" and cannot be a timeseries output.')
+                raise RuntimeError(f'ODE output {var} is tagged with "dymos.static_output" '
+                                   'and cannot be a timeseries output.')
 
         src_idxs = om.slicer[node_idxs, ...]
 
diff --git a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
index d5005df7..4457440e 100644
--- a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
+++ b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
@@ -397,7 +397,8 @@ class Radau(PseudospectralBase):
             src_units = meta['units']
             src_tags = meta['tags']
             if 'dymos.static_output' in src_tags:
-                raise RuntimeError(f'ODE output {var} is tagged with "dymos.static_output" and cannot be a timeseries output.')
+                raise RuntimeError(f'ODE output {var} is tagged with "dymos.static_output" and '
+                                   'cannot be a timeseries output.')
 
         src_idxs = om.slicer[node_idxs, ...]
 
diff --git a/dymos/utils/misc.py b/dymos/utils/misc.py
index e8c37a15..c0e81a1f 100644
--- a/dymos/utils/misc.py
+++ b/dymos/utils/misc.py
@@ -1,5 +1,7 @@
+
 import numpy as np
 
+import openmdao
 import openmdao.api as om
 from openmdao.core.constants import _ReprClass
 
@@ -253,3 +255,26 @@ def create_subprob(base_name, comm, reports=False):
         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
+
+
+def om_version():
+    """
+    Return version infromation for OpenMDAO.
+
+    This information is useful for executing code that requires a specific
+    version of OpenMDAO. The tuple format returned by this function can
+    easily be compared using a statement like `if om_version()[0] < (3, 3, 0)`.
+
+    Returns
+    -------
+    tuple
+        The semantic version of OpenMDAO in a comparable tuple.
+    str
+        One of "dev" or "release".
+    """
+    try:
+        numeric, rel = openmdao.__version__.split('-')
+    except ValueError:
+        numeric = openmdao.__version__
+        rel = 'release'
+    return tuple([int(s) for s in numeric.split('.')]), rel
diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py
index efe12e46..7e43d72a 100644
--- a/dymos/visualization/timeseries/bokeh_timeseries_report.py
+++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py
@@ -245,13 +245,13 @@ def _load_data_sources(traj_and_phase_meta=None, solution_record_file=None, simu
     """
     data_dict = {}
 
-    if Path(solution_record_file).is_file():
+    if solution_record_file is not None and Path(solution_record_file).is_file():
         sol_cr = om.CaseReader(solution_record_file)
         sol_case = sol_cr.get_case('final')
     else:
         sol_case = None
 
-    if Path(simulation_record_file).is_file():
+    if simulation_record_file is not None and Path(simulation_record_file).is_file():
         sim_cr = om.CaseReader(simulation_record_file)
         sim_case = sim_cr.get_case('final')
     else:
diff --git a/dymos/visualization/timeseries_plots.py b/dymos/visualization/timeseries_plots.py
index bdc403da..37a02b30 100644
--- a/dymos/visualization/timeseries_plots.py
+++ b/dymos/visualization/timeseries_plots.py
@@ -136,7 +136,7 @@ def _mpl_timeseries_plots(time_units, var_units, phase_names, phases_node_path,
         plt.subplots_adjust(bottom=0.23, top=0.9, left=0.2)
 
         # save to file
-        plot_file_path = plot_dir_path.joinpath(f'{var_name.replace(":", "_")}.png')
+        plot_file_path = pathlib.Path(plot_dir_path).joinpath(f'{var_name.replace(":", "_")}.png')
         plt.savefig(plot_file_path, dpi=dpi)
         plt.close(fig)
         plotfiles.append(plot_file_path)
@@ -358,21 +358,9 @@ def timeseries_plots(solution_recorder_filename, simulation_record_file=None, pl
         for models with only static parameters that are uninteresting to plot.
     """
     # get ready to generate plot files
-
-    if not pathlib.Path(plot_dir).is_absolute():
-        if problem is None:
-            plot_dir_path = pathlib.Path.cwd().joinpath(plot_dir)
-        else:
-            try:
-                repdir = problem.get_reports_dir()
-            except AttributeError:
-                from openmdao.utils.reports_system import get_reports_dir
-                repdir = get_reports_dir(problem)
-            plot_dir_path = pathlib.Path(repdir).joinpath(plot_dir)
-    else:
-        plot_dir_path = plot_dir
-
-    plot_dir_path.mkdir(parents=True, exist_ok=True)
+    if problem is not None:
+        plot_dir = problem.get_reports_dir() / plot_dir
+    pathlib.Path(plot_dir).mkdir(parents=True, exist_ok=True)
 
     cr = om.CaseReader(solution_recorder_filename)
 
@@ -430,10 +418,10 @@ def timeseries_plots(solution_recorder_filename, simulation_record_file=None, pl
 
     if dymos_options['plots'] == 'bokeh':
         _bokeh_timeseries_plots(time_units, var_units, phase_names, phases_node_path,
-                                last_solution_case, last_simulation_case, plot_dir_path)
+                                last_solution_case, last_simulation_case, plot_dir)
     elif dymos_options['plots'] == 'matplotlib':
         fnames = _mpl_timeseries_plots(time_units, var_units, phase_names, phases_node_path,
-                                       last_solution_case, last_simulation_case, plot_dir_path,
+                                       last_solution_case, last_simulation_case, plot_dir,
                                        dpi, include_parameters)
         if (problem is not None) and make_html:
             for name in fnames:
diff --git a/pyproject.toml b/pyproject.toml
index 13246923..a7c048d9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -87,8 +87,12 @@ indent-width = 4
 indent-style = "space"
 line-ending = "lf"
 
+[tool.ruff.lint]
+select = ["E", "F"]
+preview = false
+
 [tool.ruff.lint.pycodestyle]
-max-line-length = 100
+max-line-length = 130
 
 [tool.ruff.lint.pydocstyle]
 convention = "numpy"
@@ -97,4 +101,6 @@ convention = "numpy"
 # Ignore `F401` (unused import) in api files.
 "__init__.py" = ["F401"]
 # Ignore import related conventions in docs
-"**/dymos_book/**" = ["E402"]
+"**/dymos_book/**" = ["E402", "E501"]
+# Ignore Test classifiers
+"test_*.py" = ["D", "E", "F"]


Test patch:
diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml
index 13556d19..acb26f65 100644
--- a/.github/workflows/dymos_tests_workflow.yml
+++ b/.github/workflows/dymos_tests_workflow.yml
@@ -305,6 +305,7 @@ jobs:
           pipdeptree
 
       - name: 'Upload environment artifact'
+        if: env.RUN_BUILD
         uses: actions/upload-artifact@v3
         with:
           name: ${{ matrix.NAME }}_environment
@@ -312,6 +313,7 @@ jobs:
           retention-days: 5
 
       - name: Check NumPy 2.0 Compatibility
+        if: env.RUN_BUILD
         run: |
           echo "============================================================="
           echo "Check Dymos code for NumPy 2.0 compatibility"
@@ -321,7 +323,7 @@ jobs:
           ruff check . --select NPY201
 
       - name: Perform linting with Ruff
-        if: ${{ matrix.NAME == 'baseline' }}
+        if: env.RUN_BUILD && matrix.NAME == 'baseline'
         run: |
           echo "============================================================="
           echo "Lint Dymos code per settings in pyproject.toml"
diff --git a/docs/dymos_book/test/test_jupyter_output_linting.py b/docs/dymos_book/test/test_jupyter_output_linting.py
index 8ac029b2..0aa55ea8 100644
--- a/docs/dymos_book/test/test_jupyter_output_linting.py
+++ b/docs/dymos_book/test/test_jupyter_output_linting.py
@@ -80,7 +80,7 @@ class LintJupyterOutputsTestCase(unittest.TestCase):
                   "        missing_packages.append('pyoptsparse')\n",
                   "if missing_packages:\n",
                   "    raise EnvironmentError('This notebook requires the following packages '\n",
-                  "                           'please install them and restart this notebook\\'s runtime: {\",\".join(missing_packages)}')"]
+                  "                           'please install them and restart this notebook\\'s runtime: {\",\".join(missing_packages)}')"]  # noqa: E501
 
         mpi_header = ['%pylab inline\n',
                       'from ipyparallel import Client, error\n',
diff --git a/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py b/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
index 41539630..a5337817 100644
--- a/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
+++ b/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
@@ -1,4 +1,5 @@
 import os
+import pathlib
 import unittest
 
 import openmdao.api as om
@@ -6,6 +7,7 @@ from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.testing_utils import use_tempdirs
 
 import dymos as dm
+from dymos.utils.misc import om_version
 from dymos.examples.aircraft_steady_flight.aircraft_ode import AircraftODE
 
 optimizer = os.environ.get('DYMOS_DEFAULT_OPT', 'SLSQP')
@@ -292,9 +294,7 @@ class TestAircraftCruise(unittest.TestCase):
         p['assumptions.mass_empty'] = 0.15E6
         p['assumptions.mass_payload'] = 84.02869 * 400
 
-        dm.run_problem(p, simulate=False, make_plots=True)
-
-        # TODO: Simulation is inaccurate with Birkhoff and high orders. due to Vandermonde interpolation.
+        dm.run_problem(p, simulate=True, make_plots=True)
 
         time = p.get_val('traj.phase0.timeseries.time')
         tas = p.get_val('traj.phase0.timeseries.TAS', units='km/s')
@@ -304,15 +304,19 @@ class TestAircraftCruise(unittest.TestCase):
         assert_near_equal(range, tas*time, tolerance=1.0E-4)
         assert_near_equal(mass_fuel[-1, ...], 0.0, tolerance=1.0E-4)
 
-        # case = om.CaseReader('dymos_simulation.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sim_out_dir = traj.sim_prob.get_outputs_dir()
+        else:
+            sim_out_dir = pathlib.Path.cwd()
+        case = om.CaseReader(sim_out_dir / 'dymos_simulation.db').get_case('final')
 
-        # time = case.get_val('traj.phase0.timeseries.time')
-        # tas = case.get_val('traj.phase0.timeseries.TAS', units='km/s')
-        # range = case.get_val('traj.phase0.timeseries.range')
-        # mass_fuel = case.get_val('traj.phase0.timeseries.mass_fuel')
+        time = case.get_val('traj.phase0.timeseries.time')
+        tas = case.get_val('traj.phase0.timeseries.TAS', units='km/s')
+        range = case.get_val('traj.phase0.timeseries.range')
+        mass_fuel = case.get_val('traj.phase0.timeseries.mass_fuel')
 
-        # assert_near_equal(range, tas*time, tolerance=1.0E-4)
-        # assert_near_equal(mass_fuel[-1, ...], 0.0, tolerance=1.0E-4)
+        assert_near_equal(range, tas*time, tolerance=1.0E-4)
+        assert_near_equal(mass_fuel[-1, ...], 0.0, tolerance=1.0E-4)
 
 
 if __name__ == '__main__':  # pragma: no cover
diff --git a/dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py b/dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py
index 7f98f2a6..4b7c8023 100644
--- a/dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py
+++ b/dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py
@@ -1,16 +1,14 @@
 import unittest
 
-import openmdao
 import openmdao.api as om
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 import dymos as dm
 from dymos.utils.lgl import lgl
+from dymos.utils.misc import om_version
 from dymos.examples.aircraft_steady_flight.aircraft_ode import AircraftODE
 
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
-
 
 def ex_aircraft_steady_flight(transcription, optimizer='SLSQP', use_boundary_constraints=False,
                               show_output=True):
@@ -132,7 +130,7 @@ class TestExSteadyAircraftFlight(unittest.TestCase):
                           726.85, tolerance=1.0E-2)
 
     @require_pyoptsparse(optimizer='IPOPT')
-    @unittest.skipIf(om_version < (3, 32, 2), 'Test requires OpenMDAO 3.32.2 or later')
+    @unittest.skipIf(om_version()[0] < (3, 32, 2), 'Test requires OpenMDAO 3.32.2 or later')
     def test_ex_aircraft_steady_flight_opt_birkhoff(self):
 
         tx = dm.Birkhoff(num_nodes=50, grid_type='cgl',
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 595e57a2..e4ccd25c 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
@@ -3,6 +3,8 @@ import unittest
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 from openmdao.utils.assert_utils import assert_near_equal
 
+from dymos.utils.misc import om_version
+
 
 SHOW_PLOTS = True
 
@@ -237,8 +239,14 @@ class TestBalancedFieldLengthForDocs(unittest.TestCase):
 
         dm.run_problem(p, run_driver=True, simulate=True, make_plots=True)
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol = om.CaseReader(sol_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         sol_r_f_climb = sol.get_val('traj.climb.timeseries.r')[-1, ...]
         sol_r_f_rto = sol.get_val('traj.rto.timeseries.r')[-1, ...]
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 9fd8ca32..9fc42825 100644
--- a/dymos/examples/balanced_field/test/test_balanced_field_length.py
+++ b/dymos/examples/balanced_field/test/test_balanced_field_length.py
@@ -10,6 +10,7 @@ from openmdao.utils.assert_utils import assert_near_equal
 
 import dymos as dm
 from dymos.examples.balanced_field.balanced_field_ode import BalancedFieldODEComp
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -230,14 +231,24 @@ class TestBalancedFieldLengthRestart(unittest.TestCase):
         p = self._make_problem()
         dm.run_problem(p, run_driver=True, simulate=False)
 
-        sol_results = om.CaseReader('dymos_solution.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+
+        sol_results = om.CaseReader(sol_db).get_case('final')
 
         assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01)
 
-        dm.run_problem(p, run_driver=True, simulate=True, restart='dymos_solution.db')
+        dm.run_problem(p, run_driver=True, simulate=True, restart=sol_db)
+
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
 
-        sol_results = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_results = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_results = om.CaseReader(sol_db).get_case('final')
+        sim_results = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01)
         assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01)
@@ -252,14 +263,20 @@ class TestBalancedFieldLengthRestart(unittest.TestCase):
         p = self._make_problem()
         dm.run_problem(p, run_driver=True, simulate=True)
 
-        sol_results = om.CaseReader('dymos_solution.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol_results = om.CaseReader(sol_db).get_case('final')
 
         assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01)
 
-        dm.run_problem(p, run_driver=True, simulate=True, restart='dymos_simulation.db')
+        dm.run_problem(p, run_driver=True, simulate=True, restart=sim_db)
 
-        sol_results = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_results = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_results = om.CaseReader(sol_db).get_case('final')
+        sim_results = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01)
         assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01)
diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py
index e1783e3b..6424cbfa 100644
--- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py
+++ b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py
@@ -12,6 +12,7 @@ except ImportError:
 
 
 from openmdao.utils.testing_utils import use_tempdirs
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -143,8 +144,11 @@ class TestBrachistochroneStaticGravity(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016,
                           tolerance=1.0E-3)
 
-        # Generate the explicitly simulated trajectory
-        exp_out = om.CaseReader('dymos_simulation.db').get_case('final')
+        # Load the explicitly simulated trajectory
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+        exp_out = om.CaseReader(sim_db).get_case('final')
 
         # Extract the timeseries from the implicit solution and the explicit simulation
         x = p.get_val('traj.phase0.timeseries.x')
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py b/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py
index 5e7f906d..0dd0ed9e 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py
@@ -5,6 +5,7 @@ import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 from openmdao.utils.assert_utils import assert_near_equal
 import dymos as dm
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -90,26 +91,31 @@ class TestBrachExecCompODE(unittest.TestCase):
 
         return p
 
-    def run_asserts(self):
+    def run_asserts(self, c):
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = c.get_outputs_dir() / sol_db
+            sim_db = c.model.traj0.sim_prob.get_outputs_dir() / sim_db
 
-        for db in ['dymos_solution.db', 'dymos_simulation.db']:
-            p = om.CaseReader(db).get_case('final')
+        for db in [sol_db, sim_db]:
+            c = om.CaseReader(db).get_case('final')
 
-            t_initial = p.get_val('traj0.phase0.timeseries.time')[0]
-            tf = p.get_val('traj0.phase0.timeseries.time')[-1]
+            t_initial = c.get_val('traj0.phase0.timeseries.time')[0]
+            tf = c.get_val('traj0.phase0.timeseries.time')[-1]
 
-            x0 = p.get_val('traj0.phase0.timeseries.x')[0]
-            xf = p.get_val('traj0.phase0.timeseries.x')[-1]
+            x0 = c.get_val('traj0.phase0.timeseries.x')[0]
+            xf = c.get_val('traj0.phase0.timeseries.x')[-1]
 
-            y0 = p.get_val('traj0.phase0.timeseries.y')[0]
-            yf = p.get_val('traj0.phase0.timeseries.y')[-1]
+            y0 = c.get_val('traj0.phase0.timeseries.y')[0]
+            yf = c.get_val('traj0.phase0.timeseries.y')[-1]
 
-            v0 = p.get_val('traj0.phase0.timeseries.v')[0]
-            vf = p.get_val('traj0.phase0.timeseries.v')[-1]
+            v0 = c.get_val('traj0.phase0.timeseries.v')[0]
+            vf = c.get_val('traj0.phase0.timeseries.v')[-1]
 
-            g = p.get_val('traj0.phase0.parameter_vals:g')[0]
+            g = c.get_val('traj0.phase0.parameter_vals:g')[0]
 
-            thetaf = p.get_val('traj0.phase0.timeseries.theta')[-1]
+            thetaf = c.get_val('traj0.phase0.timeseries.theta')[-1]
 
             assert_near_equal(t_initial, 0.0)
             assert_near_equal(x0, 0.0)
@@ -125,16 +131,16 @@ class TestBrachExecCompODE(unittest.TestCase):
             assert_near_equal(thetaf, 100.12, tolerance=0.01)
 
     def test_ex_brachistochrone_radau_uncompressed(self):
-        self._make_problem(transcription='radau-ps', compressed=False)
-        self.run_asserts()
+        p = self._make_problem(transcription='radau-ps', compressed=False)
+        self.run_asserts(p)
 
     def test_ex_brachistochrone_gl_uncompressed(self):
-        self._make_problem(transcription='gauss-lobatto', compressed=False)
-        self.run_asserts()
+        p = self._make_problem(transcription='gauss-lobatto', compressed=False)
+        self.run_asserts(p)
 
     def test_ex_brachistochrone_birkhoff_uncompressed(self):
-        self._make_problem(transcription='birkhoff', transcription_order=10)
-        self.run_asserts()
+        p = self._make_problem(transcription='birkhoff', transcription_order=10)
+        self.run_asserts(p)
 
 
 @use_tempdirs
@@ -305,26 +311,31 @@ class TestBrachCallableODE(unittest.TestCase):
 
         return p
 
-    def run_asserts(self):
+    def run_asserts(self, c):
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = c.get_outputs_dir() / sol_db
+            sim_db = c.model.traj0.sim_prob.get_outputs_dir() / sim_db
 
-        for db in ['dymos_solution.db', 'dymos_simulation.db']:
-            p = om.CaseReader(db).get_case('final')
+        for db in (sol_db, sim_db):
+            c = om.CaseReader(db).get_case('final')
 
-            t_initial = p.get_val('traj0.phase0.timeseries.time')[0]
-            tf = p.get_val('traj0.phase0.timeseries.time')[-1]
+            t_initial = c.get_val('traj0.phase0.timeseries.time')[0]
+            tf = c.get_val('traj0.phase0.timeseries.time')[-1]
 
-            x0 = p.get_val('traj0.phase0.timeseries.x')[0]
-            xf = p.get_val('traj0.phase0.timeseries.x')[-1]
+            x0 = c.get_val('traj0.phase0.timeseries.x')[0]
+            xf = c.get_val('traj0.phase0.timeseries.x')[-1]
 
-            y0 = p.get_val('traj0.phase0.timeseries.y')[0]
-            yf = p.get_val('traj0.phase0.timeseries.y')[-1]
+            y0 = c.get_val('traj0.phase0.timeseries.y')[0]
+            yf = c.get_val('traj0.phase0.timeseries.y')[-1]
 
-            v0 = p.get_val('traj0.phase0.timeseries.v')[0]
-            vf = p.get_val('traj0.phase0.timeseries.v')[-1]
+            v0 = c.get_val('traj0.phase0.timeseries.v')[0]
+            vf = c.get_val('traj0.phase0.timeseries.v')[-1]
 
-            g = p.get_val('traj0.phase0.parameter_vals:g')[0]
+            g = c.get_val('traj0.phase0.parameter_vals:g')[0]
 
-            thetaf = p.get_val('traj0.phase0.timeseries.theta')[-1]
+            thetaf = c.get_val('traj0.phase0.timeseries.theta')[-1]
 
             assert_near_equal(t_initial, 0.0)
             assert_near_equal(x0, 0.0)
@@ -340,10 +351,12 @@ class TestBrachCallableODE(unittest.TestCase):
             assert_near_equal(thetaf, 100.12, tolerance=0.01)
 
     def test_ex_brachistochrone_radau_uncompressed(self):
-        self._make_problem(transcription='radau-ps', compressed=False)
-        self.run_asserts()
+        p = self._make_problem(transcription='radau-ps', compressed=False)
+        self.run_asserts(p)
 
     def test_in_series(self):
-        self._make_problem(transcription='gauss-lobatto', compressed=False)
-        self._make_problem(transcription='radau-ps', compressed=False)
-        self.run_asserts()
+        p = self._make_problem(transcription='gauss-lobatto', compressed=False)
+        self.run_asserts(p)
+
+        p = self._make_problem(transcription='radau-ps', compressed=False)
+        self.run_asserts(p)
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_implicit_recorder.py b/dymos/examples/brachistochrone/test/test_brachistochrone_implicit_recorder.py
index 5068b2fe..e2322797 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_implicit_recorder.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_implicit_recorder.py
@@ -7,6 +7,7 @@ except ImportError:
     matplotlib = None
 
 from openmdao.utils.testing_utils import use_tempdirs
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -66,7 +67,10 @@ class TestBrachistochroneRecordingExample(unittest.TestCase):
         # Test the results
         assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
 
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+        case = om.CaseReader(sol_db).get_case('final')
 
         assert_near_equal(p['phase0.control_values:theta'],
                           case.get_val('phase0.timeseries.theta'))
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_integrated_control.py b/dymos/examples/brachistochrone/test/test_brachistochrone_integrated_control.py
index 7e72926f..d72fff9f 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_integrated_control.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_integrated_control.py
@@ -8,6 +8,7 @@ import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs
 
 import dymos as dm
+from dymos.utils.misc import om_version
 
 
 class BrachistochroneODE(om.ExplicitComponent):
@@ -206,8 +207,14 @@ class TestBrachistochroneIntegratedControl(unittest.TestCase):
         # Solve for the optimal trajectory
         dm.run_problem(p, simulate=True, make_plots=True)
 
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
 
         # Test the results
         assert_near_equal(sol_case.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_simulate_options.py b/dymos/examples/brachistochrone/test/test_brachistochrone_simulate_options.py
index 3237e351..9622f525 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_simulate_options.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_simulate_options.py
@@ -5,6 +5,8 @@ import numpy as np
 import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs
 
+from dymos.utils.misc import om_version
+
 
 class BrachistochroneODE(om.ExplicitComponent):
 
@@ -144,7 +146,11 @@ class TestBrachistochroneSimulate(unittest.TestCase):
                 # Run the driver to simulate the problem
                 dm.run_problem(p, run_driver=False, simulate=True)
 
-                sim = om.CaseReader('dymos_simulation.db').get_case('final')
+                sim_db = 'dymos_simulation.db'
+                if om_version()[0] > (3, 34, 2):
+                    sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+                sim = om.CaseReader(sim_db).get_case('final')
                 t = sim.get_val('traj.phase0.timeseries.time')
 
                 opt_times_per_seg = phase.simulate_options['times_per_seg']
@@ -157,6 +163,10 @@ class TestBrachistochroneSimulate(unittest.TestCase):
 
                 dm.run_problem(p, run_driver=False, simulate=True, simulate_kwargs={'times_per_seg': 7})
 
-                sim = om.CaseReader('dymos_simulation.db').get_case('final')
+                sim_db = 'dymos_simulation.db'
+                if om_version()[0] > (3, 34, 2):
+                    sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+                sim = om.CaseReader(sim_db).get_case('final')
                 t = sim.get_val('traj.phase0.timeseries.time')
                 self.assertEqual(7 * num_segments, t.size)
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py b/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
index d5fa4926..3e8973de 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py
@@ -13,9 +13,7 @@ from openmdao.visualization.scaling_viewer.scaling_report import _default_scalin
 
 import dymos as dm
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
-
-
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
+from dymos.utils.misc import om_version
 
 
 def setup_model_radau(do_reports, probname):
@@ -138,7 +136,7 @@ class TestSubproblemReportToggle(unittest.TestCase):
         if self.testflo_running is not None:
             os.environ['TESTFLO_RUNNING'] = self.testflo_running
 
-    @unittest.skipIf(om_version <= (3, 34, 2), 'Requires OpenMDAO version later than 3.34.2')
+    @unittest.skipIf(om_version()[0] <= (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')
@@ -150,7 +148,7 @@ class TestSubproblemReportToggle(unittest.TestCase):
 
         self.assertFalse(sim_reports_dir.exists())
 
-    @unittest.skipIf(om_version <= (3, 34, 2), 'Requires OpenMDAO version later than 3.34.2')
+    @unittest.skipIf(om_version()[0] <= (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')
@@ -164,7 +162,7 @@ class TestSubproblemReportToggle(unittest.TestCase):
         self.assertTrue(sim_reports_dir.exists())
         self.assertTrue((sim_reports_dir / self.n2_filename).exists())
 
-    @unittest.skipIf(om_version <= (3, 34, 2), 'Requires OpenMDAO version later than 3.34.2')
+    @unittest.skipIf(om_version()[0] <= (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,
@@ -180,7 +178,7 @@ class TestSubproblemReportToggle(unittest.TestCase):
         self.assertIn(self.n2_filename, main_reports)
         self.assertIn(self.scaling_filename, main_reports)
 
-    @unittest.skipIf(om_version <= (3, 34, 2), 'Requires OpenMDAO version later than 3.34.2')
+    @unittest.skipIf(om_version()[0] <= (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,
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_timeseries_feedback.py b/dymos/examples/brachistochrone/test/test_brachistochrone_timeseries_feedback.py
index 52ab3eef..9024da39 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_timeseries_feedback.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_timeseries_feedback.py
@@ -6,6 +6,8 @@ import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs
 from openmdao.utils.assert_utils import assert_near_equal
 
+from dymos.utils.misc import om_version
+
 
 class BrachistochroneODE(om.ExplicitComponent):
 
@@ -148,8 +150,14 @@ class TestBrachistochroneTimeseriesFeedback(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, run_driver=True, simulate=True)
 
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
 
         final_time_timeseries_sol = sol_case.get_val('traj.phase0.timeseries.time')[-1, ...]
         final_time_traj_param_sol = sol_case.get_val('traj.parameter_vals:final_time')[0, 0]
diff --git a/dymos/examples/brachistochrone/test/test_simulate_units.py b/dymos/examples/brachistochrone/test/test_simulate_units.py
index a9cf2106..cdb8fd33 100644
--- a/dymos/examples/brachistochrone/test/test_simulate_units.py
+++ b/dymos/examples/brachistochrone/test/test_simulate_units.py
@@ -8,6 +8,7 @@ from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneOD
 
 from openmdao.utils.testing_utils import use_tempdirs
 from openmdao.utils.assert_utils import assert_near_equal
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -85,8 +86,14 @@ class TestBrachistochroneSimulate_units(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, simulate=True)
 
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(sim_case.get_val('traj.phase0.parameter_vals:g', units='m/s**2')[0],
                           sol_case.get_val('traj.phase0.parameter_vals:g', units='m/s**2')[0])
diff --git a/dymos/examples/brachistochrone/test/test_state_rate_introspection.py b/dymos/examples/brachistochrone/test/test_state_rate_introspection.py
index c6e79dd5..77863324 100644
--- a/dymos/examples/brachistochrone/test/test_state_rate_introspection.py
+++ b/dymos/examples/brachistochrone/test/test_state_rate_introspection.py
@@ -7,6 +7,7 @@ from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 import dymos as dm
+from dymos.utils.misc import om_version
 from dymos.utils.testing_utils import assert_timeseries_near_equal
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
 
@@ -86,8 +87,14 @@ class TestIntegrateControl(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, simulate=True, make_plots=False)
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol = om.CaseReader(sol_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         t_sol = sol.get_val('traj.phase0.timeseries.time')
         t_sim = sim.get_val('traj.phase0.timeseries.time')
@@ -190,8 +197,14 @@ class TestIntegrateControl(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, simulate=True, make_plots=False)
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol = om.CaseReader(sol_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         t_sol = sol.get_val('traj.phase0.timeseries.time')
         t_sim = sim.get_val('traj.phase0.timeseries.time')
@@ -306,8 +319,14 @@ class TestIntegrateControl(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, simulate=True, make_plots=False, simulate_kwargs={'atol': 1.0E-9, 'rtol': 1.0E-9})
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol = om.CaseReader(sol_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         t_sol = sol.get_val('traj.phase0.timeseries.time')
         t_sim = sim.get_val('traj.phase0.timeseries.time')
@@ -429,8 +448,14 @@ class TestIntegratePolynomialControl(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, simulate=True, make_plots=True)
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol = om.CaseReader(sol_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         t_sol = sol.get_val('traj.phase0.timeseries.time')
         t_sim = sim.get_val('traj.phase0.timeseries.time')
@@ -533,8 +558,14 @@ class TestIntegratePolynomialControl(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, simulate=True, make_plots=False)
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol = om.CaseReader(sol_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         t_sol = sol.get_val('traj.phase0.timeseries.time')
         t_sim = sim.get_val('traj.phase0.timeseries.time')
@@ -651,8 +682,14 @@ class TestIntegratePolynomialControl(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, simulate=True, make_plots=True)
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol = om.CaseReader(sol_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         t_sol = sol.get_val('traj.phase0.timeseries.time')
         t_sim = sim.get_val('traj.phase0.timeseries.time')
diff --git a/dymos/examples/brachistochrone/test/test_static_outputs.py b/dymos/examples/brachistochrone/test/test_static_outputs.py
index b67e7802..661da547 100644
--- a/dymos/examples/brachistochrone/test/test_static_outputs.py
+++ b/dymos/examples/brachistochrone/test/test_static_outputs.py
@@ -7,6 +7,7 @@ import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs
 
 import dymos as dm
+from dymos.utils.misc import om_version
 
 
 class BrachODEStaticOutput(om.ExplicitComponent):
@@ -137,8 +138,14 @@ class TestStaticODEOutput(unittest.TestCase):
         # Run the driver to solve the problem
         dm.run_problem(p, simulate=True, make_plots=False)
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sim_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol = om.CaseReader(sim_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         with self.assertRaises(expected_exception=KeyError) as e:
             sol.get_val('traj.phase0.timeseries.foo')
diff --git a/dymos/examples/brachistochrone/test/test_timeseries_units.py b/dymos/examples/brachistochrone/test/test_timeseries_units.py
index e1c15276..52c3744c 100644
--- a/dymos/examples/brachistochrone/test/test_timeseries_units.py
+++ b/dymos/examples/brachistochrone/test/test_timeseries_units.py
@@ -4,6 +4,7 @@ import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 from openmdao.utils.assert_utils import assert_near_equal
 import dymos as dm
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -90,8 +91,14 @@ class TestTimeseriesUnits(unittest.TestCase):
 
         dm.run_problem(p, run_driver=run_driver, simulate=True)
 
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj0.sim_prob.get_outputs_dir() / sim_db
+
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(sol_case.get_val('traj0.phase0.timeseries.xdot', units='degF'),
                           sol_case.get_val('traj0.phase0.timeseries.xdot', units='degF'))
diff --git a/dymos/examples/cannonball/test/test_cannonball_matrix_state.py b/dymos/examples/cannonball/test/test_cannonball_matrix_state.py
index 98361227..8b7af412 100644
--- a/dymos/examples/cannonball/test/test_cannonball_matrix_state.py
+++ b/dymos/examples/cannonball/test/test_cannonball_matrix_state.py
@@ -4,6 +4,7 @@ import dymos as dm
 import unittest
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+from dymos.utils.misc import om_version
 
 
 class ODEComp(om.ExplicitComponent):
@@ -74,7 +75,15 @@ class TestCannonballMatrixState(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -92,7 +101,15 @@ class TestCannonballMatrixState(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -110,7 +127,15 @@ class TestCannonballMatrixState(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -128,7 +153,15 @@ class TestCannonballMatrixState(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -147,7 +180,11 @@ class TestCannonballMatrixState(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -165,7 +202,11 @@ class TestCannonballMatrixState(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -216,7 +257,11 @@ class TestCannonballMatrixStateExplicitShape(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -233,7 +278,11 @@ class TestCannonballMatrixStateExplicitShape(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -251,7 +300,11 @@ class TestCannonballMatrixStateExplicitShape(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -268,7 +321,11 @@ class TestCannonballMatrixStateExplicitShape(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -320,7 +377,11 @@ class TestCannonballMatrixStateNonFlatIndices(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -337,7 +398,11 @@ class TestCannonballMatrixStateNonFlatIndices(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -355,7 +420,11 @@ class TestCannonballMatrixStateNonFlatIndices(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
@@ -372,7 +441,11 @@ class TestCannonballMatrixStateNonFlatIndices(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
-        c = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        c = om.CaseReader(sim_db).get_case('final')
 
         assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
diff --git a/dymos/examples/cannonball/test/test_connect_control_to_parameter.py b/dymos/examples/cannonball/test/test_connect_control_to_parameter.py
index d0c92545..4679e160 100644
--- a/dymos/examples/cannonball/test/test_connect_control_to_parameter.py
+++ b/dymos/examples/cannonball/test/test_connect_control_to_parameter.py
@@ -2,16 +2,12 @@ import unittest
 
 import numpy as np
 
-import openmdao
 import openmdao.api as om
 
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 from dymos.examples.cannonball.cannonball_ode import rho_interp
 
-om_dev_version = openmdao.__version__.endswith('dev')
-om_version = tuple(int(s) for s in openmdao.__version__.split('-')[0].split('.'))
-
 
 class CannonballODEVectorCD(om.ExplicitComponent):
     """
diff --git a/dymos/examples/cannonball/test/test_load_case_missing_phase.py b/dymos/examples/cannonball/test/test_load_case_missing_phase.py
index 61372722..0aa20ec7 100644
--- a/dymos/examples/cannonball/test/test_load_case_missing_phase.py
+++ b/dymos/examples/cannonball/test/test_load_case_missing_phase.py
@@ -5,13 +5,11 @@ from openmdao.utils.assert_utils import assert_near_equal
 
 import numpy as np
 import openmdao.api as om
-import openmdao
 
 import dymos as dm
 from dymos.examples.cannonball.size_comp import CannonballSizeComp
 from dymos.examples.cannonball.cannonball_ode import CannonballODE
-
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -168,8 +166,13 @@ class TestTwoPhaseCannonballLoadCase(unittest.TestCase):
         descent.set_state_val('v', [150, 200])
         descent.set_state_val('gam', [0, -45], units='deg')
 
-        case = om.CaseReader('dymos_solution.db').get_case('final')
-        if om_version < (3, 26, 1):
+        sol_db = 'dymos_solution.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p_ascent_only.get_outputs_dir() / sol_db
+
+        case = om.CaseReader(sol_db).get_case('final')
+
+        if om_version()[0] < (3, 26, 1):
             dm.load_case(p, previous_solution=case)
         else:
             p.load_case(case)
diff --git a/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py b/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py
index 4be60db4..894d674e 100644
--- a/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py
+++ b/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py
@@ -2,16 +2,12 @@ import unittest
 
 import numpy as np
 
-import openmdao
 import openmdao.api as om
 
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 from dymos.examples.cannonball.cannonball_ode import rho_interp
 
-om_dev_version = openmdao.__version__.endswith('dev')
-om_version = tuple(int(s) for s in openmdao.__version__.split('-')[0].split('.'))
-
 
 GRAVITY = 9.80665  # m/s**2
 
diff --git a/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py b/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py
index 0f56813c..c9fffc2d 100644
--- a/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py
+++ b/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py
@@ -8,6 +8,7 @@ from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 import dymos as dm
 from dymos.examples.double_integrator.double_integrator_ode import DoubleIntegratorODE
 from dymos.utils.testing_utils import assert_timeseries_near_equal
+from dymos.utils.misc import om_version
 
 
 @require_pyoptsparse(optimizer='SLSQP')
@@ -76,10 +77,16 @@ class TestDoubleIntegratorExample(unittest.TestCase):
         assert_near_equal(v[-1], 0.0, tolerance=tol)
 
     def test_timeseries_units_gl(self):
-        double_integrator_direct_collocation(dm.GaussLobatto, compressed=True)
+        p = double_integrator_direct_collocation(dm.GaussLobatto, compressed=True)
 
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
 
         t_sol = sol_case.get_val('traj.phase0.timeseries.time')
         t_sim = sim_case.get_val('traj.phase0.timeseries.time')
@@ -91,10 +98,16 @@ class TestDoubleIntegratorExample(unittest.TestCase):
                                          abs_tolerance=0.01)
 
     def test_timeseries_units_radau(self):
-        double_integrator_direct_collocation(dm.Radau, compressed=True)
+        p = double_integrator_direct_collocation(dm.Radau, compressed=True)
+
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
 
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
 
         t_sol = sol_case.get_val('traj.phase0.timeseries.time')
         t_sim = sim_case.get_val('traj.phase0.timeseries.time')
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 d5e102a5..9519e7ca 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
@@ -62,7 +62,7 @@ class TestExampleTwoBurnOrbitRaise(unittest.TestCase):
                               'vt.png', 'pos_y.png', 'theta.png')
 
             for file in expected_files:
-                plotfile = pathlib.Path(_get_reports_dir(p)).joinpath('plots') / file
+                plotfile = pathlib.Path(p.get_reports_dir(p)).joinpath('plots') / file
                 self.assertTrue(plotfile.exists(), msg=f'{plotfile} does not exist!')
 
 
diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py
index 5e45e5e9..2902f24a 100644
--- a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py
+++ b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py
@@ -6,6 +6,7 @@ from openmdao.utils.mpi import MPI
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 from dymos.examples.finite_burn_orbit_raise.finite_burn_orbit_raise_problem import two_burn_orbit_raise_problem
+from dymos.utils.misc import om_version
 
 
 @require_pyoptsparse(optimizer='IPOPT')
@@ -19,12 +20,18 @@ class TestExampleTwoBurnOrbitRaiseMPI(unittest.TestCase):
 
         CONNECTED = False
 
-        two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3,
-                                     compressed=False, optimizer=optimizer, simulate=True,
-                                     connected=CONNECTED, show_output=False)
+        p = two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3,
+                                         compressed=False, optimizer=optimizer, simulate=True,
+                                         connected=CONNECTED, show_output=False)
 
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
 
         # The last phase in this case is run in reverse time if CONNECTED=True,
         # so grab the correct index to test the resulting delta-V.
@@ -38,12 +45,18 @@ class TestExampleTwoBurnOrbitRaiseMPI(unittest.TestCase):
 
         CONNECTED = True
 
-        two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3,
-                                     compressed=False, optimizer=optimizer, simulate=True,
-                                     connected=CONNECTED, show_output=False)
+        p = two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3,
+                                         compressed=False, optimizer=optimizer, simulate=True,
+                                         connected=CONNECTED, show_output=False)
+
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
 
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
 
         # The last phase in this case is run in reverse time if CONNECTED=True,
         # so grab the correct index to test the resulting delta-V.
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 5165cea0..766c8474 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
@@ -7,6 +7,7 @@ from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 from openmdao.utils.mpi import MPI, multi_proc_exception_check
 
 from dymos.examples.finite_burn_orbit_raise.finite_burn_orbit_raise_problem import two_burn_orbit_raise_problem
+from dymos.utils.misc import om_version
 from dymos.utils.testing_utils import assert_cases_equal
 
 
@@ -29,17 +30,27 @@ class TestExampleTwoBurnOrbitRaiseConnectedRestart(unittest.TestCase):
                 assert_near_equal(p.get_val('traj.burn2.states:deltav')[0], 0.3995,
                                   tolerance=4.0E-3)
 
-        sim_case1 = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sim_case1 = om.CaseReader(sim_db).get_case('final')
 
         # Run again without an actual optimizer
         p = two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3,
                                          compressed=False, optimizer=optimizer,
                                          show_output=False, connected=True, run_driver=False,
-                                         restart='dymos_solution.db')
+                                         restart=sol_db)
 
         p.run_model()
 
-        sim_case2 = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sim_case2 = om.CaseReader(sim_db).get_case('final')
 
         # Verify that the second case has the same inputs and outputs
         assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8)
@@ -50,7 +61,13 @@ class TestExampleTwoBurnOrbitRaiseConnectedRestart(unittest.TestCase):
         p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
                                          compressed=False, optimizer=optimizer, show_output=False)
 
-        sim_case1 = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sim_case1 = om.CaseReader(sim_db).get_case('final')
 
         with multi_proc_exception_check(p.comm):
             if p.model.traj.phases.burn2 in p.model.traj.phases._subsystems_myproc:
@@ -60,9 +77,13 @@ class TestExampleTwoBurnOrbitRaiseConnectedRestart(unittest.TestCase):
         # Run again without an actual optimzier
         two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
                                      compressed=False, optimizer=optimizer, run_driver=False,
-                                     show_output=False, restart='dymos_solution.db')
+                                     show_output=False, restart=sol_db)
+
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
 
-        sim_case2 = om.CaseReader('dymos_simulation.db').get_case('final')
+        sim_case2 = om.CaseReader(sim_db).get_case('final')
 
         # Verify that the second case has the same inputs and outputs
         assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8)
@@ -105,18 +126,30 @@ class TestExampleTwoBurnOrbitRaiseConnected(unittest.TestCase):
                 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')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        case1 = om.CaseReader(sol_db).get_case('final')
+        sim_case1 = om.CaseReader(sim_db).get_case('final')
 
         # Run again without an actual optimizer
         p = two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3,
                                          compressed=False, optimizer=optimizer, run_driver=False,
-                                         show_output=False, restart='dymos_solution.db',
+                                         show_output=False, restart=sol_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')
+        sol_db = 'dymos_solution2.db'
+        sim_db = 'dymos_simulation2.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        case2 = om.CaseReader(sol_db).get_case('final')
+        sim_case2 = om.CaseReader(sim_db).get_case('final')
 
         # Verify that the second case has the same inputs and outputs
         assert_cases_equal(case1, case2, tol=1.0E-8)
@@ -135,18 +168,30 @@ class TestExampleTwoBurnOrbitRaiseConnected(unittest.TestCase):
                 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')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        case1 = om.CaseReader(sol_db).get_case('final')
+        sim_case1 = om.CaseReader(sim_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',
+                                         show_output=False, restart=sol_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')
+        sol_db = 'dymos_solution2.db'
+        sim_db = 'dymos_simulation2.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        case2 = om.CaseReader(sol_db).get_case('final')
+        sim_case2 = om.CaseReader(sim_db).get_case('final')
 
         # Verify that the second case has the same inputs and outputs
         assert_cases_equal(case1, case2, tol=1.0E-7)
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 c8addd53..32dfaa1b 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
@@ -1,6 +1,5 @@
 import unittest
 import numpy as np
-import openmdao
 from numpy.polynomial import Polynomial as P
 try:
     import matplotlib.pyplot as plt
@@ -15,11 +14,9 @@ from dymos.utils.introspection import get_promoted_vars
 
 import dymos as dm
 from dymos.examples.min_time_climb.min_time_climb_ode import MinTimeClimbODE
-
+from dymos.utils.misc import om_version
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
-
 
 def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto',
                    transcription_order=3, force_alloc_complex=False, add_rate=False, time_name='time',
@@ -180,14 +177,19 @@ class TestMinTimeClimb(unittest.TestCase):
         ts = {k: v for k, v in output_dict.items() if 'timeseries.' in k}
         self.assertTrue('traj.phase0.timeseries.mach_rate' in ts)
 
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sol_db).get_case('final')
 
         time = case[f'traj.phase0.timeseries.{time_name}'][:, 0]
         mach = case['traj.phase0.timeseries.mach'][:, 0]
         mach_rate = case['traj.phase0.timeseries.mach_rate'][:, 0]
 
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
-
         sim_time = sim_case[f'traj.phase0.timeseries.{time_name}'][:, 0]
         sim_mach = sim_case['traj.phase0.timeseries.mach'][:, 0]
         sim_mach_rate = sim_case['traj.phase0.timeseries.mach_rate'][:, 0]
@@ -277,7 +279,7 @@ class TestMinTimeClimb(unittest.TestCase):
         self._test_mach_rate(p, plot=False)
 
     @require_pyoptsparse(optimizer='IPOPT')
-    @unittest.skipIf(om_version < (3, 32, 2), 'Test requires OpenMDAO 3.31.2 or later')
+    @unittest.skipIf(om_version()[0] < (3, 32, 2), 'Test requires OpenMDAO 3.32.2 or later')
     def test_results_birkhoff(self):
         NUM_SEG = 1
         ORDER = 30
diff --git a/dymos/examples/shuttle_reentry/test/test_reentry.py b/dymos/examples/shuttle_reentry/test/test_reentry.py
index 05be7863..741133f3 100644
--- a/dymos/examples/shuttle_reentry/test/test_reentry.py
+++ b/dymos/examples/shuttle_reentry/test/test_reentry.py
@@ -8,6 +8,7 @@ from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 from dymos.examples.shuttle_reentry.shuttle_ode import ShuttleODE
+from dymos.utils.misc import om_version
 
 # Expected results from Betts
 expected_results = {'constrained': {'time': 2198.67, 'theta': 30.6255},
@@ -272,8 +273,15 @@ class TestReentry(unittest.TestCase):
 
         dm.run_problem(p, simulate=True)
 
-        sol = om.CaseReader('dymos_solution.db').get_case('final')
-        sim = om.CaseReader('dymos_simulation.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / 'dymos_solution.db'
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db'
+        else:
+            sol_db = 'dymos_solution.db'
+            sim_db = 'dymos_simulation.db'
+
+        sol = om.CaseReader(sol_db).get_case('final')
+        sim = om.CaseReader(sim_db).get_case('final')
 
         from scipy.interpolate import interp1d
 
diff --git a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py
index 8897552a..3eced152 100644
--- a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py
+++ b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py
@@ -1,12 +1,9 @@
 import os
 import unittest
 
-import openmdao
 from openmdao.utils.testing_utils import use_tempdirs
 from openmdao.utils.mpi import MPI
-
-
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -25,7 +22,7 @@ class TestVanderpolForDocs(unittest.TestCase):
 
         dm.run_problem(p, run_driver=False, simulate=True, make_plots=True)
 
-    @unittest.skipIf(om_version < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
+    @unittest.skipIf(om_version()[0] < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
     def test_vanderpol_for_docs_optimize(self):
         import dymos as dm
         from dymos.examples.vanderpol.vanderpol_dymos import vanderpol
@@ -36,7 +33,7 @@ class TestVanderpolForDocs(unittest.TestCase):
 
         dm.run_problem(p, simulate=True, make_plots=True)
 
-    @unittest.skipIf(om_version < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
+    @unittest.skipIf(om_version()[0] < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
     def test_vanderpol_for_docs_optimize_refine(self):
         import dymos as dm
         from dymos.examples.vanderpol.vanderpol_dymos import vanderpol
@@ -62,7 +59,7 @@ class TestVanderpolForDocs(unittest.TestCase):
 class TestVanderpolDelayMPI(unittest.TestCase):
     N_PROCS = 2
 
-    @unittest.skipIf(om_version < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
+    @unittest.skipIf(om_version()[0] < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
     def test_vanderpol_delay_mpi(self):
         import openmdao.api as om
         import dymos as dm
diff --git a/dymos/examples/vanderpol/test/test_vanderpol.py b/dymos/examples/vanderpol/test/test_vanderpol.py
index 35170467..1c7d38a3 100644
--- a/dymos/examples/vanderpol/test/test_vanderpol.py
+++ b/dymos/examples/vanderpol/test/test_vanderpol.py
@@ -2,15 +2,12 @@ import unittest
 from numpy.testing import assert_almost_equal
 import numpy as np
 
-import openmdao
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 from openmdao.utils.mpi import MPI
 
 import dymos as dm
 from dymos.examples.vanderpol.vanderpol_dymos import vanderpol
-
-
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -20,7 +17,7 @@ class TestVanderpolExample(unittest.TestCase):
         p = vanderpol(transcription='gauss-lobatto', num_segments=75)
         p.run_model()
 
-    @unittest.skipIf(om_version < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
+    @unittest.skipIf(om_version()[0] < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
     @require_pyoptsparse(optimizer='IPOPT')
     def test_vanderpol_simulate_true(self):
         p = vanderpol(transcription='radau-ps', num_segments=30, transcription_order=3,
@@ -28,7 +25,7 @@ class TestVanderpolExample(unittest.TestCase):
 
         dm.run_problem(p, run_driver=True, simulate=True)
 
-    @unittest.skipIf(om_version < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
+    @unittest.skipIf(om_version()[0] < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
     def test_vanderpol_optimal(self):
         p = vanderpol(transcription='gauss-lobatto', num_segments=75)
         dm.run_problem(p)  # find optimal control solution to stop oscillation
@@ -39,7 +36,7 @@ class TestVanderpolExample(unittest.TestCase):
         assert_almost_equal(p.get_val('traj.phase0.states:x1')[-1, ...], np.zeros(1))
         assert_almost_equal(p.get_val('traj.phase0.controls:u')[-1, ...], np.zeros(1), decimal=3)
 
-    @unittest.skipIf(om_version < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
+    @unittest.skipIf(om_version()[0] < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
     def test_vanderpol_optimal_grid_refinement(self):
         # enabling grid refinement gives a faster and better solution with fewer segments
         p = vanderpol(transcription='gauss-lobatto', num_segments=15)
@@ -59,7 +56,7 @@ class TestVanderpolExampleMPI(unittest.TestCase):
 
     N_PROCS = 4
 
-    @unittest.skipIf(om_version < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
+    @unittest.skipIf(om_version()[0] < (3, 29, 0), 'Test requires OpenMDAO 3.29.0 or later')
     @require_pyoptsparse(optimizer='IPOPT')
     @unittest.skipUnless(MPI, 'this test requires MPI')
     def test_vanderpol_optimal_mpi(self):
diff --git a/dymos/phase/test/test_analytic_phase.py b/dymos/phase/test/test_analytic_phase.py
index 8be126be..06416dde 100644
--- a/dymos/phase/test/test_analytic_phase.py
+++ b/dymos/phase/test/test_analytic_phase.py
@@ -6,6 +6,7 @@ import dymos as dm
 
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.testing_utils import use_tempdirs
+from dymos.utils.misc import om_version
 
 
 class SimpleIVPSolution(om.ExplicitComponent):
@@ -254,8 +255,12 @@ class TestAnalyticPhaseSimpleResults(unittest.TestCase):
         p.set_val('traj.phase.t_duration', 1.0, units='s')
         p.set_val('traj.phase.parameters:y0', 0.6, units='unitless')
 
+        sol_db = 'dymos_solution.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+
         # Load the previous case and rerun
-        dm.run_problem(p, simulate=False, restart='dymos_solution.db')
+        dm.run_problem(p, simulate=False, restart=sol_db)
 
         t = p.get_val('traj.phase.timeseries.time', units='s')
         y = p.get_val('traj.phase.timeseries.y', units='unitless')
diff --git a/dymos/test/test_load_case.py b/dymos/test/test_load_case.py
index fa810109..66eb9420 100644
--- a/dymos/test/test_load_case.py
+++ b/dymos/test/test_load_case.py
@@ -2,11 +2,9 @@ import os
 import unittest
 
 from openmdao.utils.testing_utils import use_tempdirs
-import openmdao
 import openmdao.api as om
 import dymos as dm
-
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
+from dymos.utils.misc import om_version
 
 
 def setup_problem(trans=dm.GaussLobatto(num_segments=10), polynomial_control=False,
@@ -60,11 +58,6 @@ def setup_problem(trans=dm.GaussLobatto(num_segments=10), polynomial_control=Fal
 @use_tempdirs
 class TestLoadCase(unittest.TestCase):
 
-    def tearDown(self):
-        for filename in ['brachistochrone_solution.db']:
-            if os.path.exists(filename):
-                os.remove(filename)
-
     def test_load_case_unchanged_grid(self):
         import openmdao.api as om
         from openmdao.utils.assert_utils import assert_near_equal
@@ -76,7 +69,12 @@ class TestLoadCase(unittest.TestCase):
         dm.run_problem(p)
 
         # Load the solution
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_file = p.get_outputs_dir() / 'dymos_solution.db'
+        else:
+            sol_file = 'dymos_solution.db'
+
+        case = om.CaseReader(sol_file).get_case('final')
 
         # Initialize the system with values from the case.
         # We unnecessarily call setup again just to make sure we obliterate the previous solution
@@ -104,7 +102,12 @@ class TestLoadCase(unittest.TestCase):
         dm.run_problem(p)
 
         # Load the solution
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_file = p.get_outputs_dir() / 'dymos_solution.db'
+        else:
+            sol_file = 'dymos_solution.db'
+
+        case = om.CaseReader(sol_file).get_case('final')
 
         # Initialize the system with values from the case.
         # We unnecessarily call setup again just to make sure we obliterate the previous solution
@@ -130,7 +133,12 @@ class TestLoadCase(unittest.TestCase):
         dm.run_problem(p)
 
         # Load the solution
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_file = p.get_outputs_dir() / 'dymos_solution.db'
+        else:
+            sol_file = 'dymos_solution.db'
+
+        case = om.CaseReader(sol_file).get_case('final')
 
         # create a problem with a different transcription with a different number of variables
         q = setup_problem(dm.Radau(num_segments=20))
@@ -162,7 +170,12 @@ class TestLoadCase(unittest.TestCase):
         dm.run_problem(p)
 
         # Load the solution
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_file = p.get_outputs_dir() / 'dymos_solution.db'
+        else:
+            sol_file = 'dymos_solution.db'
+
+        case = om.CaseReader(sol_file).get_case('final')
 
         # create a problem with a different transcription with a different number of variables
         q = setup_problem(dm.GaussLobatto(num_segments=50))
@@ -194,7 +207,12 @@ class TestLoadCase(unittest.TestCase):
         dm.run_problem(p)
 
         # Load the solution
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_file = p.get_outputs_dir() / 'dymos_solution.db'
+        else:
+            sol_file = 'dymos_solution.db'
+
+        case = om.CaseReader(sol_file).get_case('final')
 
         # create a problem with a different transcription with a different number of variables
         q = setup_problem(dm.Birkhoff(num_nodes=50))
@@ -264,7 +282,12 @@ class TestLoadCase(unittest.TestCase):
         dm.run_problem(p)
 
         # Load the solution
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_file = p.get_outputs_dir() / 'dymos_solution.db'
+        else:
+            sol_file = 'dymos_solution.db'
+
+        case = om.CaseReader(sol_file).get_case('final')
 
         # create a problem with a different transcription with a different number of variables
         q = setup_problem(dm.GaussLobatto(num_segments=50))
@@ -290,7 +313,12 @@ class TestLoadCase(unittest.TestCase):
         dm.run_problem(p)
 
         # Load the solution
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_file = p.get_outputs_dir() / 'dymos_solution.db'
+        else:
+            sol_file = 'dymos_solution.db'
+
+        case = om.CaseReader(sol_file).get_case('final')
 
         # create a problem with a different transcription with a different number of variables
         q = setup_problem(dm.Radau(num_segments=10), fix_final_state=False, fix_final_control=True)
@@ -311,7 +339,12 @@ class TestLoadCase(unittest.TestCase):
         dm.run_problem(p)
 
         # Load the solution
-        case = om.CaseReader('dymos_solution.db').get_case('final')
+        if om_version()[0] > (3, 34, 2):
+            sol_file = p.get_outputs_dir() / 'dymos_solution.db'
+        else:
+            sol_file = 'dymos_solution.db'
+
+        case = om.CaseReader(sol_file).get_case('final')
 
         # create a problem with a different transcription with a different number of variables
         q = setup_problem(dm.Radau(num_segments=10), polynomial_control=True,
diff --git a/dymos/test/test_run_problem.py b/dymos/test/test_run_problem.py
index 4e0930f8..7ad49482 100755
--- a/dymos/test/test_run_problem.py
+++ b/dymos/test/test_run_problem.py
@@ -11,7 +11,6 @@ try:
 except ImportError:
     matplotlib = None
 
-import openmdao
 import openmdao.api as om
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
@@ -21,12 +20,11 @@ import dymos as dm
 from dymos.examples.hyper_sensitive.hyper_sensitive_ode import HyperSensitiveODE
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
 from dymos.examples.brachistochrone.brachistochrone_vector_states_ode import BrachistochroneVectorStatesODE
+from dymos.utils.misc import om_version
 from dymos.utils.testing_utils import _get_reports_dir
 
 _, optimizer = set_pyoptsparse_opt('IPOPT', fallback=True)
 
-om_version = tuple([int(s) for s in openmdao.__version__.split('-')[0].split('.')])
-
 
 @use_tempdirs
 class TestRunProblem(unittest.TestCase):
@@ -268,10 +266,16 @@ class TestRunProblem(unittest.TestCase):
 
         dm.run_problem(p, simulate=True, simulate_kwargs={'times_per_seg': 100})
 
-        self.assertTrue(os.path.exists('dymos_solution.db'))
         # Assert the results are what we expect.
-        sol_case = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol_case = om.CaseReader(sol_db).get_case('final')
+        sim_case = om.CaseReader(sim_db).get_case('final')
+        self.assertTrue(os.path.exists(sol_db))
 
         sol_t = sol_case.get_val('traj.phase0.timeseries.time')
         sim_t = sim_case.get_val('traj.phase0.timeseries.time')
@@ -383,13 +387,19 @@ class TestRunProblem(unittest.TestCase):
 
         dm.run_problem(p, refine_iteration_limit=20)
 
-        self.assertTrue(os.path.exists('dymos_solution.db'))
+        sol_db = 'dymos_solution.db'
+        driver_db = 'brach_driver_rec.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            driver_db = p.get_outputs_dir() / driver_db
+
+        self.assertTrue(os.path.exists(sol_db))
         # Assert the results are what we expect.
-        cr = om.CaseReader('dymos_solution.db')
+        cr = om.CaseReader(sol_db)
         case = cr.get_case('final')
         assert_almost_equal(case.outputs['traj.phase0.timeseries.time'].max(), 1.8016, decimal=4)
 
-        cr_opt = om.CaseReader('brach_driver_rec.db')
+        cr_opt = om.CaseReader(driver_db)
         cases = cr_opt.list_cases(source='driver', out_stream=None)
 
         for case in cases:
@@ -438,13 +448,18 @@ class TestRunProblem(unittest.TestCase):
 
         dm.run_problem(p, refine_iteration_limit=20, case_prefix='brach_test')
 
-        self.assertTrue(os.path.exists('dymos_solution.db'))
+        sol_db = 'dymos_solution.db'
+        driver_db = 'brach_driver_rec.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            driver_db = p.get_outputs_dir() / driver_db
+
+        self.assertTrue(os.path.exists(sol_db))
         # Assert the results are what we expect.
-        cr = om.CaseReader('dymos_solution.db')
-        case = cr.get_case('brach_test_final')
+        case = om.CaseReader(sol_db).get_case('brach_test_final')
         assert_almost_equal(case.outputs['traj.phase0.timeseries.time'].max(), 1.8016, decimal=4)
 
-        cr_opt = om.CaseReader('brach_driver_rec.db')
+        cr_opt = om.CaseReader(driver_db)
         cases = cr_opt.list_cases(source='driver', out_stream=None)
 
         for case in cases:
@@ -534,10 +549,14 @@ class TestRunProblem(unittest.TestCase):
 
         dm.run_problem(p, simulate=True)
 
-        self.assertTrue(os.path.exists('dymos_solution.db'))
+        sol_db = 'dymos_solution.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+
+        self.assertTrue(os.path.exists(sol_db))
+
         # Assert the results are what we expect.
-        cr = om.CaseReader('dymos_solution.db')
-        case = cr.get_case('final')
+        case = om.CaseReader(sol_db).get_case('final')
         assert_almost_equal(case.outputs['traj.phase0.timeseries.time'].max(), 1.8016, decimal=4)
 
     def test_restart_from_file(self):
@@ -550,19 +569,22 @@ class TestRunProblem(unittest.TestCase):
         p = vanderpol(transcription='gauss-lobatto', num_segments=75)
 
         # Run the problem (simulate only)
-        p.run_model()
+        dm.run_problem(p, run_driver=False, simulate=True)
 
-        # simulate and record
-        p.model.traj.simulate(record_file='vanderpol_simulation.sql')
+        # Run the model
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
 
         # create a new problem for restart to simulate a different command line execution
         q = vanderpol(transcription='gauss-lobatto', num_segments=75)
 
-        # # Run the model
-        run_problem(q, run_driver=False, simulate=False, restart='vanderpol_simulation.sql')
+        run_problem(q, run_driver=False, simulate=False, restart=sim_db)
 
-        # s = q.model.traj.simulate(rtol=1.0E-9, atol=1.0E-9)
-        s = om.CaseReader('vanderpol_simulation.sql').get_case('final')
+        sol_db2 = 'dymos_solution.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db2 = q.get_outputs_dir() / sol_db2
+        s = om.CaseReader(sol_db2).get_case('final')
 
         # get_val returns data for duplicate time points; remove them before interpolating
         tq = q.get_val('traj.phase0.timeseries.time')[:, 0]
@@ -598,17 +620,22 @@ class TestRunProblem(unittest.TestCase):
         p = vanderpol(transcription='gauss-lobatto', num_segments=75)
 
         # Run the problem (simulate only)
-        p.run_model()
-
-        # simulate and record
-        p.model.traj.simulate(record_file='vanderpol_simulation.sql')
+        dm.run_problem(p, run_driver=False, simulate=True)
 
         # create a new problem for restart to simulate a different command line execution
         q = vanderpol(transcription='gauss-lobatto', num_segments=75)
 
-        # # Run the model
-        s = om.CaseReader('vanderpol_simulation.sql').get_case('final')
-        run_problem(q, run_driver=False, simulate=False, restart=s)
+        # Run the model
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        run_problem(q, run_driver=False, simulate=False, restart=sim_db)
+
+        sol_db2 = 'dymos_solution.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db2 = q.get_outputs_dir() / sol_db2
+        s = om.CaseReader(sol_db2).get_case('final')
 
         # get_val returns data for duplicate time points; remove them before interpolating
         tq = q.get_val('traj.phase0.timeseries.time')[:, 0]
@@ -718,16 +745,22 @@ class TestRunProblemPlotting(unittest.TestCase):
             self.assertFalse(plotfile.exists(), msg=f'Unexpectedly found plot file {plotfile}')
 
     def test_run_brachistochrone_problem_set_simulation_record_file(self):
-        simulation_record_file = 'simulation_record_file.db'
-        dm.run_problem(self.p, simulate=True, simulation_record_file=simulation_record_file)
+        sim_db = 'simulation_record_file.db'
+        dm.run_problem(self.p, simulate=True, simulation_record_file=sim_db)
+
+        if om_version()[0] > (3, 34, 2):
+            sim_db = self.p.model.traj.sim_prob.get_outputs_dir() / sim_db
 
-        self.assertTrue(os.path.exists(simulation_record_file))
+        self.assertTrue(os.path.exists(sim_db))
 
     def test_run_brachistochrone_problem_set_solution_record_file(self):
-        solution_record_file = 'solution_record_file.db'
-        dm.run_problem(self.p, solution_record_file=solution_record_file)
+        sol_db = 'solution_record_file.db'
+        dm.run_problem(self.p, solution_record_file=sol_db)
+
+        if om_version()[0] > (3, 34, 2):
+            sol_db = self.p.get_outputs_dir() / sol_db
 
-        self.assertTrue(os.path.exists(solution_record_file))
+        self.assertTrue(os.path.exists(sol_db))
 
     @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
     def test_run_brachistochrone_problem_plot_simulation(self):
@@ -832,8 +865,14 @@ class TestSimulateArrayParam(unittest.TestCase):
         dm.run_problem(p, simulate=True)
 
         # Test the results
-        sol_results = om.CaseReader('dymos_solution.db').get_case('final')
-        sim_results = om.CaseReader('dymos_solution.db').get_case('final')
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db
+
+        sol_results = om.CaseReader(sol_db).get_case('final')
+        sim_results = om.CaseReader(sim_db).get_case('final')
 
         sol = sol_results.get_val('traj.phase0.parameter_vals:array')
         sim = sim_results.get_val('traj.phase0.parameter_vals:array')
diff --git a/dymos/test/test_upgrade_guide.py b/dymos/test/test_upgrade_guide.py
index 8e5820c1..b40b02fc 100644
--- a/dymos/test/test_upgrade_guide.py
+++ b/dymos/test/test_upgrade_guide.py
@@ -7,6 +7,7 @@ import dymos as dm
 
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+from dymos.utils.misc import om_version
 
 
 @use_tempdirs
@@ -649,12 +650,6 @@ class TestUpgrade_0_17_0(unittest.TestCase):
 @use_tempdirs
 class TestUpgrade_0_19_0(unittest.TestCase):
 
-    def tearDown(self):
-        if os.path.exists('dymos_solution.db'):
-            os.remove('dymos_solution.db')
-        if os.path.exists('dymos_simulation.db'):
-            os.remove('dymos_simulation.db')
-
     @require_pyoptsparse(optimizer='SLSQP')
     def _make_problem(self, transcription='gauss-lobatto', num_segments=8, transcription_order=3,
                       compressed=True, optimizer='SLSQP', run_driver=True,
@@ -733,9 +728,15 @@ class TestUpgrade_0_19_0(unittest.TestCase):
 
         return p
 
-    def run_asserts(self):
+    def run_asserts(self, p):
+
+        sol_db = 'dymos_solution.db'
+        sim_db = 'dymos_simulation.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.traj0.sim_prob.get_outputs_dir() / sim_db
 
-        for db in ['dymos_solution.db', 'dymos_simulation.db']:
+        for db in (sol_db, sim_db):
             p = om.CaseReader(db).get_case('final')
 
             t_initial = p.get_val('traj0.phase0.timeseries.time')[0]
@@ -768,9 +769,9 @@ class TestUpgrade_0_19_0(unittest.TestCase):
             assert_near_equal(thetaf, 100.12, tolerance=0.01)
 
     def test_ex_brachistochrone_radau_uncompressed(self):
-        self._make_problem(transcription='radau-ps', compressed=False)
-        self.run_asserts()
+        p = self._make_problem(transcription='radau-ps', compressed=False)
+        self.run_asserts(p)
 
     def test_ex_brachistochrone_gl_uncompressed(self):
-        self._make_problem(transcription='gauss-lobatto', compressed=False)
-        self.run_asserts()
+        p = self._make_problem(transcription='gauss-lobatto', compressed=False)
+        self.run_asserts(p)
diff --git a/dymos/utils/test/test_testing_utils.py b/dymos/utils/test/test_testing_utils.py
index cb32f26a..56eb4bd1 100644
--- a/dymos/utils/test/test_testing_utils.py
+++ b/dymos/utils/test/test_testing_utils.py
@@ -5,6 +5,7 @@ import numpy as np
 import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs
 
+from dymos.utils.misc import om_version
 from dymos.utils.testing_utils import assert_cases_equal, assert_timeseries_near_equal
 
 
@@ -356,8 +357,14 @@ class TestAssertCasesEqual(unittest.TestCase):
         p2.record('final')
         p2.cleanup()
 
-        c1 = om.CaseReader('p1.db').get_case('final')
-        c2 = om.CaseReader('p2.db').get_case('final')
+        p1_db = 'p1.db'
+        p2_db = 'p2.db'
+        if om_version()[0] > (3, 34, 2):
+            p1_db = p1.get_outputs_dir() / p1_db
+            p2_db = p2.get_outputs_dir() / p2_db
+
+        c1 = om.CaseReader(p1_db).get_case('final')
+        c2 = om.CaseReader(p2_db).get_case('final')
 
         with self.assertRaises(AssertionError) as e:
             assert_cases_equal(c1, c2)
@@ -396,8 +403,14 @@ class TestAssertCasesEqual(unittest.TestCase):
         p2.record('final')
         p2.cleanup()
 
-        c1 = om.CaseReader('p1.db').get_case('final')
-        c2 = om.CaseReader('p2.db').get_case('final')
+        p1_db = 'p1.db'
+        p2_db = 'p2.db'
+        if om_version()[0] > (3, 34, 2):
+            p1_db = p1.get_outputs_dir() / p1_db
+            p2_db = p2.get_outputs_dir() / p2_db
+
+        c1 = om.CaseReader(p1_db).get_case('final')
+        c2 = om.CaseReader(p2_db).get_case('final')
 
         assert_cases_equal(c1, c2, require_same_vars=False)
 
@@ -428,8 +441,14 @@ class TestAssertCasesEqual(unittest.TestCase):
         p2.record('final')
         p2.cleanup()
 
-        c1 = om.CaseReader('p1.db').get_case('final')
-        c2 = om.CaseReader('p2.db').get_case('final')
+        p1_db = 'p1.db'
+        p2_db = 'p2.db'
+        if om_version()[0] > (3, 34, 2):
+            p1_db = p1.get_outputs_dir() / p1_db
+            p2_db = p2.get_outputs_dir() / p2_db
+
+        c1 = om.CaseReader(p1_db).get_case('final')
+        c2 = om.CaseReader(p2_db).get_case('final')
 
         with self.assertRaises(AssertionError) as e:
             assert_cases_equal(c1, c2)
@@ -467,8 +486,14 @@ class TestAssertCasesEqual(unittest.TestCase):
         p2.record('final')
         p2.cleanup()
 
-        c1 = om.CaseReader('p1.db').get_case('final')
-        c2 = om.CaseReader('p2.db').get_case('final')
+        p1_db = 'p1.db'
+        p2_db = 'p2.db'
+        if om_version()[0] > (3, 34, 2):
+            p1_db = p1.get_outputs_dir() / p1_db
+            p2_db = p2.get_outputs_dir() / p2_db
+
+        c1 = om.CaseReader(p1_db).get_case('final')
+        c2 = om.CaseReader(p2_db).get_case('final')
 
         expected = "\nThe following variables contain different values:\n" \
                    "var        max error       mean error\n" \
diff --git a/dymos/visualization/test/test_timeseries_plots.py b/dymos/visualization/test/test_timeseries_plots.py
index b14180f2..1dbb487b 100644
--- a/dymos/visualization/test/test_timeseries_plots.py
+++ b/dymos/visualization/test/test_timeseries_plots.py
@@ -16,6 +16,7 @@ from dymos.examples.finite_burn_orbit_raise.finite_burn_eom import FiniteBurnODE
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
 from dymos.examples.battery_multibranch.battery_multibranch_ode import BatteryODE
 from dymos.utils.lgl import lgl
+from dymos.utils.misc import om_version
 from dymos.utils.testing_utils import _get_reports_dir
 from dymos.visualization.timeseries_plots import timeseries_plots
 
@@ -89,7 +90,11 @@ class TestTimeSeriesPlotsBasics(unittest.TestCase):
         with dm.options.temporary(plots='matplotlib'):
             dm.run_problem(self.p, make_plots=False)
 
-            timeseries_plots('dymos_solution.db', problem=self.p)
+            sol_db = 'dymos_solution.db'
+            if om_version()[0] > (3, 34, 2):
+                sol_db = self.p.get_outputs_dir() / sol_db
+
+            timeseries_plots(sol_db, problem=self.p)
             plot_dir = pathlib.Path(_get_reports_dir(self.p)).joinpath('plots')
 
             self.assertTrue(plot_dir.joinpath('x.png').exists())
@@ -107,7 +112,11 @@ class TestTimeSeriesPlotsBasics(unittest.TestCase):
         # records to the default file 'dymos_simulation.db'
         dm.run_problem(self.p, make_plots=False, solution_record_file='solution_record_file.db')
 
-        timeseries_plots('solution_record_file.db', problem=self.p)
+        sol_db = 'solution_record_file.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = self.p.get_outputs_dir() / sol_db
+
+        timeseries_plots(sol_db, problem=self.p)
         plot_dir = pathlib.Path(_get_reports_dir(self.p)).joinpath('plots')
 
         self.assertTrue(plot_dir.joinpath('x.png').exists())
@@ -127,7 +136,15 @@ class TestTimeSeriesPlotsBasics(unittest.TestCase):
         dm.run_problem(self.p, simulate=True, make_plots=False,
                        simulation_record_file='simulation_record_file.db')
 
-        timeseries_plots('dymos_solution.db', simulation_record_file='simulation_record_file.db', problem=self.p)
+        sol_db = 'dymos_solution.db'
+        sim_db = 'simulation_record_file.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = self.p.get_outputs_dir() / sol_db
+            sim_db = self.p.model.traj0.sim_prob.get_outputs_dir() / sim_db
+
+        timeseries_plots(sol_db,
+                         simulation_record_file=sim_db,
+                         problem=self.p)
 
         plot_dir = pathlib.Path(_get_reports_dir(self.p)).joinpath("plots").resolve()
         self.assertTrue(plot_dir.joinpath('x.png').exists())
@@ -147,8 +164,12 @@ class TestTimeSeriesPlotsBasics(unittest.TestCase):
 
         dm.run_problem(self.p, make_plots=False)
 
+        sol_db = 'dymos_solution.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = self.p.get_outputs_dir() / sol_db
+
         plot_dir = pathlib.Path(_get_reports_dir(self.p)).joinpath("test_plot_dir").resolve()
-        timeseries_plots('dymos_solution.db', plot_dir=plot_dir, problem=self.p)
+        timeseries_plots(sol_db, plot_dir=plot_dir)
 
         self.assertTrue(plot_dir.joinpath('x.png').exists())
         self.assertTrue(plot_dir.joinpath('y.png').exists())
@@ -307,7 +328,14 @@ class TestTimeSeriesPlotsMultiPhase(unittest.TestCase):
         dm.run_problem(p, simulate=True, make_plots=False,
                        simulation_record_file='simulation_record_file.db')
 
-        timeseries_plots('dymos_solution.db', simulation_record_file='simulation_record_file.db',
+        sol_db = 'dymos_solution.db'
+        sim_db = 'simulation_record_file.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = self.p.get_outputs_dir() / sol_db
+            sim_db = self.p.model.sim_prob.get_outputs_dir() / sim_db
+
+        timeseries_plots(sol_db,
+                         simulation_record_file=sim_db,
                          problem=p)
         plot_dir = pathlib.Path(_get_reports_dir(p)).joinpath("plots")
 
@@ -408,10 +436,18 @@ class TestTimeSeriesPlotsMultiPhase(unittest.TestCase):
         dm.run_problem(prob, simulate=True, make_plots=False,
                        simulation_record_file='simulation_record_file.db')
 
-        timeseries_plots('dymos_solution.db', simulation_record_file='simulation_record_file.db',
-                         problem=prob)
+        sol_db = 'dymos_solution.db'
+        sim_db = 'simulation_record_file.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = prob.get_outputs_dir() / sol_db
+            sim_db = prob.model.traj.sim_prob.get_outputs_dir() / sim_db
+
         plot_dir = pathlib.Path(_get_reports_dir(prob)).joinpath("plots")
 
+        timeseries_plots(sol_db,
+                         simulation_record_file=sim_db,
+                         problem=prob)
+
         self.assertTrue(plot_dir.joinpath('time_phase.png').exists())
         self.assertTrue(plot_dir.joinpath('state_of_charge.png').exists())
         self.assertTrue(plot_dir.joinpath('state_of_charge.png').exists())
@@ -562,8 +598,13 @@ class TestTimeSeriesPlotsMultiPhase(unittest.TestCase):
         dm.run_problem(p, simulate=True, make_plots=False,
                        simulation_record_file='simulation_record_file.db')
 
-        timeseries_plots('dymos_solution.db', simulation_record_file='simulation_record_file.db',
-                         problem=p)
+        sol_db = 'dymos_solution.db'
+        sim_db = 'simulation_record_file.db'
+        if om_version()[0] > (3, 34, 2):
+            sol_db = p.get_outputs_dir() / sol_db
+            sim_db = p.model.sim_prob.get_outputs_dir() / sim_db
+
+        timeseries_plots(sol_db, sim_db, problem=p)
         plot_dir = pathlib.Path(_get_reports_dir(p)).joinpath("plots")
 
         for varname in ['time_phase', 'r', 'r_dot', 'theta',
diff --git a/joss/test/test_brachistochrone_for_joss.py b/joss/test/test_brachistochrone_for_joss.py
index 08e29551..dd4d3868 100644
--- a/joss/test/test_brachistochrone_for_joss.py
+++ b/joss/test/test_brachistochrone_for_joss.py
@@ -120,3 +120,7 @@ class TestBrachistochroneForJOSS(unittest.TestCase):
         assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-4)
         assert_near_equal(x[-1], 10, tolerance=1.0E-4)
         assert_near_equal(y[-1], 5, tolerance=1.0E-4)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/joss/test/test_cannonball_for_joss.py b/joss/test/test_cannonball_for_joss.py
index 4951fa8b..38ef1718 100644
--- a/joss/test/test_cannonball_for_joss.py
+++ b/joss/test/test_cannonball_for_joss.py
@@ -7,6 +7,7 @@ try:
 except ImportError:
     matplotlib = None
 
+
 @use_tempdirs
 class TestCannonballForJOSS(unittest.TestCase):
 
@@ -31,7 +32,6 @@ class TestCannonballForJOSS(unittest.TestCase):
                               np.array(USatm1976Data.rho * english_to_metric_rho, dtype=complex),
                               kind='linear')
 
-
         class CannonballSize(om.ExplicitComponent):
             """
             Static calculations performed before the dynamic model
@@ -57,7 +57,6 @@ class TestCannonballForJOSS(unittest.TestCase):
                 outputs['mass'] = (4 / 3.) * rho * np.pi * radius ** 3
                 outputs['area'] = np.pi * radius ** 2
 
-
         class CannonballODE(om.ExplicitComponent):
             """
             Cannonball ODE assuming flat earth and accounting for air resistance