diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 955c70b142..f6a0bf6ff2 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -53,7 +53,7 @@ jobs:
       # h5py is optional; some versions don't have binaries (yet)
       run: |
         python3 -m pip install ruamel.yaml scons==3.1.2 numpy cython pandas pytest \
-        pytest-github-actions-annotate-failures
+        pytest-github-actions-annotate-failures pint
         python3 -m pip install h5py
     - name: Build Cantera
       run: |
@@ -70,6 +70,13 @@ jobs:
     - name: Test Cantera
       run:
         python3 `which scons` test show_long_tests=yes verbose_tests=yes --debug=time
+    - name: Save the wheel file to install Cantera
+      uses: actions/upload-artifact@v3
+      with:
+        path: build/python/dist/Cantera*.whl
+        retention-days: 2
+        name: cantera-wheel-${{ matrix.python-version }}-${{ matrix.os }}
+        if-no-files-found: error
 
   clang-compiler:
     name: LLVM/Clang with Python 3.8
@@ -96,7 +103,7 @@ jobs:
       run: python3 -m pip install -U pip setuptools wheel
     - name: Install Python dependencies
       run: |
-        python3 -m pip install ruamel.yaml scons numpy cython pandas h5py pytest pytest-github-actions-annotate-failures
+        python3 -m pip install ruamel.yaml scons numpy cython pandas h5py pytest pytest-github-actions-annotate-failures pint
     - name: Build Cantera
       run: python3 `which scons` build env_vars=all
         CXX=clang++-12 CC=clang-12 f90_interface=n extra_lib_dirs=/usr/lib/llvm/lib
@@ -147,7 +154,7 @@ jobs:
     - name: Install Python dependencies
       # h5py is optional; may fail if no wheel is present for a given OS/Python version
       run: |
-        $PYTHON_CMD -m pip install ruamel.yaml numpy cython pandas pytest pytest-github-actions-annotate-failures
+        $PYTHON_CMD -m pip install ruamel.yaml numpy cython pandas pytest pytest-github-actions-annotate-failures pint
         $PYTHON_CMD -m pip install h5py || true
     - name: Install Python dependencies for GH Python
       if: matrix.python-version == '3.11'
@@ -197,7 +204,7 @@ jobs:
     - name: Install Python dependencies
       run: |
         python3 -m pip install ruamel.yaml scons numpy cython pandas scipy pytest h5py \
-        pytest-github-actions-annotate-failures pytest-cov gcovr
+        pytest-github-actions-annotate-failures pytest-cov gcovr pint
     - name: Setup .NET Core SDK
       uses: actions/setup-dotnet@v2
       with:
@@ -293,7 +300,7 @@ jobs:
       - name: Install Python dependencies
         run: |
           python3 -m pip install ruamel.yaml scons numpy cython sphinx\<4.0 jinja2\<3.1.0 \
-          sphinxcontrib-katex sphinxcontrib-matlabdomain sphinxcontrib-doxylink
+          sphinxcontrib-katex sphinxcontrib-matlabdomain sphinxcontrib-doxylink pint
       - name: Build Cantera with documentation
         run: python3 `which scons` build -j2 doxygen_docs=y sphinx_docs=y debug=n optimize=n use_pch=n
       - name: Ensure 'scons help' options work
@@ -339,8 +346,9 @@ jobs:
 
   run-examples:
     name: Run the Python examples using bash
-    runs-on: ubuntu-20.04
+    runs-on: ${{ matrix.os }}
     timeout-minutes: 60
+    needs: ["ubuntu-multiple-pythons"]
     strategy:
       matrix:
         python-version: ['3.8', '3.10', '3.11']
@@ -349,10 +357,10 @@ jobs:
       HDF5_LIBDIR: /usr/lib/x86_64-linux-gnu/hdf5/serial
       HDF5_INCLUDEDIR: /usr/include/hdf5/serial
     steps:
+      # We're not building Cantera here, we only need the checkout for the samples
+      # folder, so no need to do a recursive checkout
       - uses: actions/checkout@v3
         name: Checkout the repository
-        with:
-          submodules: recursive
       - name: Setup Python
         uses: actions/setup-python@v4
         with:
@@ -361,27 +369,18 @@ jobs:
       - name: Install Apt dependencies
         run: |
           sudo apt update
-          sudo apt install libboost-dev gfortran graphviz liblapack-dev libblas-dev \
-          gcc-9 g++-9 libhdf5-dev
+          sudo apt install graphviz
+      - name: Download the wheel artifact
+        uses: actions/download-artifact@v3
+        with:
+          name: cantera-wheel-${{ matrix.python-version }}-${{ matrix.os }}
+          path: dist
       - name: Upgrade pip
         run: python3 -m pip install -U pip setuptools wheel
       - name: Install Python dependencies
         run: |
-          python3 -m pip install ruamel.yaml scons numpy cython pandas matplotlib scipy h5py
-      - name: Build Cantera
-        # compile with GCC 9.4.0 on ubuntu-20.04 as an alternative to the default
-        # (GCC 7.5.0 is both default and oldest supported version)
-        # compile without native HDF5 support
-        run: python3 `which scons` build -j2 debug=n CC=gcc-9 CXX=g++-9
-        if: matrix.python-version != '3.10'
-      - name: Build Cantera (Python 3.10 with HDF)
-        # compile with GCC 9.4.0 on ubuntu-20.04 as an alternative to the default
-        # (GCC 7.5.0 is both default and oldest supported version)
-        # compile with native HDF5 support
-        run: |
-          python3 `which scons` build -j2 debug=n CC=gcc-9 CXX=g++-9 \
-          hdf_libdir=$HDF5_LIBDIR hdf_include=$HDF5_INCLUDEDIR
-        if: matrix.python-version == '3.10'
+            python3 -m pip install numpy ruamel.yaml h5py pandas matplotlib scipy pint
+            python3 -m pip install --pre --no-index --find-links dist cantera
       - name: Run the examples
         # See https://unix.stackexchange.com/a/392973 for an explanation of the -exec part
         run: |
@@ -389,7 +388,6 @@ jobs:
           find samples/python -type f -iname "*.py" \
           -exec sh -c 'for n; do echo "$n" | tee -a results.txt && python3 "$n" >> results.txt || exit 1; done' sh {} +
         env:
-          PYTHONPATH: build/python
           # The ignore setting here is due to a new warning introduced in Matplotlib==3.6.0
           PYTHONWARNINGS: "error,ignore:warn_name_set_on_empty_Forward::pyparsing"
           MPLBACKEND: Agg
@@ -434,7 +432,7 @@ jobs:
         run: |
           conda install -q sundials=${{ matrix.sundials-ver }} scons numpy ruamel.yaml \
           cython boost-cpp fmt eigen yaml-cpp h5py pandas libgomp openblas pytest \
-          highfive
+          highfive pint
       - name: Build Cantera
         run: |
           scons build system_fmt=y system_eigen=y system_yamlcpp=y system_sundials=y \
@@ -495,7 +493,7 @@ jobs:
       # use boost-cpp rather than boost from conda-forge
       # Install SCons >=4.4.0 to make sure that MSVC_TOOLSET_VERSION variable is present
       run: |
-        mamba install -q '"scons>=4.4.0"' numpy cython ruamel.yaml boost-cpp eigen yaml-cpp h5py pandas pytest highfive
+        mamba install -q '"scons>=4.4.0"' numpy cython ruamel.yaml boost-cpp eigen yaml-cpp h5py pandas pytest highfive pint
       shell: pwsh
     - name: Build Cantera
       run: scons build system_eigen=y system_yamlcpp=y system_highfive=y logging=debug
@@ -558,7 +556,7 @@ jobs:
       - name: Install Python dependencies
         run: |
           python -m pip install -U pip setuptools wheel
-          python -m pip install '"scons<4.4.0"' pypiwin32 numpy ruamel.yaml cython pandas pytest pytest-github-actions-annotate-failures
+          python -m pip install '"scons<4.4.0"' pypiwin32 numpy ruamel.yaml cython pandas pytest pytest-github-actions-annotate-failures pint
       - name: Restore Boost cache
         uses: actions/cache@v3
         id: cache-boost
@@ -622,7 +620,7 @@ jobs:
     - name: Install Python dependencies
       run: |
         python3 -m pip install ruamel.yaml scons numpy cython pandas h5py pytest \
-        pytest-github-actions-annotate-failures
+        pytest-github-actions-annotate-failures pint
     - name: Setup Intel oneAPI environment
       run: |
         source /opt/intel/oneapi/setvars.sh
@@ -669,7 +667,7 @@ jobs:
       run: python3 -m pip install -U pip setuptools wheel
     - name: Install Python dependencies
       run: python3 -m pip install ruamel.yaml scons numpy cython h5py pandas pytest
-        pytest-github-actions-annotate-failures
+        pytest-github-actions-annotate-failures pint
     - name: Setup Intel oneAPI environment
       run: |
         source /opt/intel/oneapi/setvars.sh
@@ -701,7 +699,7 @@ jobs:
     - name: Install Python dependencies
       run: |
         python -m pip install -U pip setuptools wheel
-        python -m pip install scons pypiwin32 numpy ruamel.yaml cython h5py pandas pytest pytest-github-actions-annotate-failures
+        python -m pip install scons pypiwin32 numpy ruamel.yaml cython h5py pandas pytest pytest-github-actions-annotate-failures pint
     - name: Restore Boost cache
       uses: actions/cache@v3
       id: cache-boost
diff --git a/AUTHORS b/AUTHORS
index 52ee1adfb5..83362a284c 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -6,6 +6,7 @@ tracker.
 
 Rounak Agarwal (@agarwalrounak)
 David Akinpelu (@DavidAkinpelu), Louisiana State University
+Halla Ali (@hallaali)
 Emil Atz (@EmilAtz)
 Jongyoon Bae (@jongyoonbae), Brown University
 Philip Berndt
diff --git a/SConstruct b/SConstruct
index 69c770f1d1..11ed47d564 100644
--- a/SConstruct
+++ b/SConstruct
@@ -813,7 +813,7 @@ else:
     config.add(config_options)
     toolchain = ["default"]
 
-env = Environment(tools=toolchain+["textfile", "subst", "recursiveInstall", "wix", "gch"],
+env = Environment(tools=toolchain+["textfile", "subst", "recursiveInstall", "UnitsInterfaceBuilder", "wix", "gch"],
                   ENV={"PATH": os.environ["PATH"]},
                   toolchain=toolchain,
                   **extraEnvArgs)
@@ -1737,6 +1737,7 @@ elif env['python_package'] == 'n':
 
 env['install_python_action'] = ''
 env['python_module_loc'] = ''
+env["python_module"] = None
 env["ct_pyscriptdir"] = "<not installed>"
 
 def check_module(name):
diff --git a/doc/example-keywords.txt b/doc/example-keywords.txt
index 22a82c264d..8cc74559f0 100644
--- a/doc/example-keywords.txt
+++ b/doc/example-keywords.txt
@@ -41,5 +41,6 @@ thermodynamic cycle
 thermodynamics
 transport
 tutorial
+units
 user-defined model
 well-stirred reactor
diff --git a/doc/sphinx/conf.py b/doc/sphinx/conf.py
index 05eeb6923d..cf1625b3a1 100644
--- a/doc/sphinx/conf.py
+++ b/doc/sphinx/conf.py
@@ -84,18 +84,18 @@ def escape_splats(app, what, name, obj, options, lines):
                     lines[i] = l.replace("*", r"\*")
     app.connect('autodoc-process-docstring', escape_splats)
 
-
 autoclass_content = 'both'
 
 doxylink = {
-        'ct': (os.path.abspath('../../build/docs/Cantera.tag'),
-               '../../doxygen/html/')
+    'ct': (os.path.abspath('../../build/docs/Cantera.tag'),
+            '../../doxygen/html/')
 }
 
 intersphinx_mapping = {
     'python': ('https://docs.python.org/3', None),
     'pandas': ('https://pandas.pydata.org/pandas-docs/stable/', None),
     'numpy': ('https://numpy.org/doc/stable/', None),
+    'pint': ('https://pint.readthedocs.io/en/stable/', None),
 }
 
 # Ensure that the primary domain is the Python domain, since we've added the
diff --git a/doc/sphinx/cython/index.rst b/doc/sphinx/cython/index.rst
index e3ddbf046a..25490f21ea 100644
--- a/doc/sphinx/cython/index.rst
+++ b/doc/sphinx/cython/index.rst
@@ -15,4 +15,5 @@ Contents:
    zerodim
    onedim
    constants
+   units
    utilities
diff --git a/doc/sphinx/cython/units.rst b/doc/sphinx/cython/units.rst
new file mode 100644
index 0000000000..f5036671a1
--- /dev/null
+++ b/doc/sphinx/cython/units.rst
@@ -0,0 +1,49 @@
+.. py:currentmodule:: cantera.with_units
+
+Python Interface With Units
+===========================
+
+This interface allows users to specify physical units associated with quantities.
+To do so, this interface leverages the `pint <https://pint.readthedocs.io/en/stable/>`__
+library to provide consistent unit conversion. Several examples of this interface can
+be found in the ``samples/python`` folder in the 
+`root of the repository <https://github.com/Cantera/cantera/tree/main/samples/python>`__.
+Examples that use this interface are suffixed with `_units`.
+
+The overall goal is to provide a compatible implementation of the `cantera.Solution` and
+`cantera.PureFluid` interfaces. Please see those pages for further documentation of the
+available properties.
+
+Solution with Units
+-------------------
+
+.. autoclass:: Solution
+   :no-members:
+
+PureFluid Phases With Units
+---------------------------
+
+The following convenience classes are available to create `PureFluid <PureFluid>`
+objects with the indicated equation of state:
+
+.. autoclass:: CarbonDioxide
+   :no-members:
+.. autoclass:: Heptane
+   :no-members:
+.. autoclass:: Hfc134a
+   :no-members:
+.. autoclass:: Hydrogen
+   :no-members:
+.. autoclass:: Methane
+   :no-members:
+.. autoclass:: Nitrogen
+   :no-members:
+.. autoclass:: Oxygen
+   :no-members:
+.. autoclass:: Water
+   :no-members:
+
+The full documentation for the `PureFluid <PureFluid>` class and its properties is here:
+
+.. autoclass:: PureFluid
+   :no-members:
diff --git a/interfaces/cython/SConscript b/interfaces/cython/SConscript
index debe8dd9b1..528b7d86bb 100644
--- a/interfaces/cython/SConscript
+++ b/interfaces/cython/SConscript
@@ -3,6 +3,8 @@ import re
 from pathlib import Path
 from buildutils import *
 
+from string import Template
+
 Import('env', 'build', 'install')
 
 localenv = env.Clone()
@@ -56,7 +58,7 @@ for pyxfile in multi_glob(localenv, "cantera", "pyx"):
     cython_obj.append(obj)
 
 copy_header = localenv.Command('#src/extensions/delegator.h', '#build/python/cantera/delegator.h',
-                          Copy('$TARGET', '$SOURCE'))
+                               Copy('$TARGET', '$SOURCE'))
 ext_manager = localenv.SharedObject("#src/extensions/PythonExtensionManager.cpp")
 
 cython_obj.extend(ext_manager)
@@ -107,6 +109,12 @@ for f in (multi_glob(localenv, 'cantera', 'py') +
           multi_glob(localenv, 'cantera/*', 'py')):
     localenv.Depends(mod, f)
 
+units = localenv.UnitsInterfaceBuilder(
+    "cantera/with_units/solution.py",
+    "cantera/with_units/solution.py.in",
+)
+localenv.Depends(mod, units)
+
 # Determine installation path and install the Python module
 install_cmd = ["$python_cmd_esc", "-m", "pip", "install"]
 user_install = False
diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx
index 0354b101aa..a56eaeaa59 100644
--- a/interfaces/cython/cantera/thermo.pyx
+++ b/interfaces/cython/cantera/thermo.pyx
@@ -1960,7 +1960,7 @@ cdef class PureFluid(ThermoPhase):
         def __set__(self, Q):
             if (self.P >= self.critical_pressure or
                 abs(self.P-self.P_sat)/self.P > 1e-4):
-                raise ValueError('Cannot set vapor quality outside the'
+                raise ValueError('Cannot set vapor quality outside the '
                                  'two-phase region')
             self.thermo.setState_Psat(self.P, Q)
 
diff --git a/interfaces/cython/cantera/with_units/__init__.py b/interfaces/cython/cantera/with_units/__init__.py
new file mode 100644
index 0000000000..d3d6aa7a47
--- /dev/null
+++ b/interfaces/cython/cantera/with_units/__init__.py
@@ -0,0 +1,12 @@
+# This file is part of Cantera. See License.txt in the top-level directory or
+# at https://cantera.org/license.txt for license and copyright information.
+
+# This code to set the application registry has to come before any code that wants to
+# use the registry is imported. In particular, it has to come before any of our code
+# that uses units!
+from pint import UnitRegistry, set_application_registry
+cantera_units_registry = UnitRegistry()
+set_application_registry(cantera_units_registry)
+
+# Now we can import our code
+from .solution import *
diff --git a/interfaces/cython/cantera/with_units/solution.py.in b/interfaces/cython/cantera/with_units/solution.py.in
new file mode 100644
index 0000000000..d6aa2993ce
--- /dev/null
+++ b/interfaces/cython/cantera/with_units/solution.py.in
@@ -0,0 +1,195 @@
+# This file is part of Cantera. See License.txt in the top-level directory or
+# at https://cantera.org/license.txt for license and copyright information.
+from textwrap import dedent
+
+from .. import Solution as _Solution, PureFluid as _PureFluid, CanteraError
+from .. import (Heptane as _Heptane, Water as _Water, Hfc134a as _Hfc134a,
+                CarbonDioxide as _CarbonDioxide, Hydrogen as _Hydrogen,
+                Methane as _Methane, Nitrogen as _Nitrogen, Oxygen as _Oxygen)
+from pint import get_application_registry
+
+__all__ = ("units", "Q_", "Solution", "PureFluid", "Heptane", "CarbonDioxide",
+           "Hfc134a", "Hydrogen", "Methane", "Nitrogen", "Oxygen", "Water",
+           "CanteraError")
+
+units = get_application_registry()
+Q_ = units.Quantity
+
+
+def copy_doc(method):
+    """Decorator to copy docstrings from related methods in upstream classes.
+
+    This decorator will copy the docstring from the same named method in the upstream
+    class, either `Solution` or `PureFluid`. The docstring in the method being
+    decorated is appended to the upstream documentation.
+    """
+    doc = getattr(method, "__doc__", None) or ""
+    if isinstance(method, property):
+        method = method.fget
+        if not doc:
+            doc = getattr(method, "__doc__", None) or ""
+    original_klass = method.__qualname__.split(".")[0]
+    klass = {"Solution": _Solution, "PureFluid": _PureFluid}[original_klass]
+    original_method = getattr(klass, method.__name__)
+    original_doc = dedent(getattr(original_method, "__doc__", ""))
+    method.__doc__ = f"{original_doc}\n{doc}"
+    return method
+
+
+class Solution:
+    """
+    This implementation of `Solution <cantera.with_units.Solution>` operates with
+    units by using the `pint` library to convert between unit systems. All properties
+    are assigned units in the standard MKS system that Cantera uses, substituting kmol
+    instead of mol. Each property is an instance of the `pint.Quantity` class.
+
+    Similarly, properties must be instances of `pint.Quantity` classes when they are
+    used for assignment to set the state. The properties may have any units, so long
+    as the dimensions for the quantity are consistent. For example, temperatures can
+    be provided in K, degC, degF, or degR; conversion will be done internally to
+    Cantera's consistent unit system.
+
+    See the `pint documentation <https://pint.readthedocs.io>`__ for more information
+    about using pint's ``Quantity`` classes.
+
+    **Note:** This class is experimental. It only implements methods from `ThermoPhase`.
+    Methods from other classes are not yet supported. If you are interested in contributing
+    to this feature, please chime in on our enhancements issue:
+    `<https://github.com/Cantera/enhancements/issues/174>`__.
+    """
+    def __init__(self, infile="", name="", *, yaml=None):
+        self.__dict__["_phase"] = _Solution(infile, name, yaml=yaml)
+
+@common_properties@
+
+@thermophase_properties@
+
+Solution.__doc__ = f"{Solution.__doc__}\n{_Solution.__doc__}"
+
+class PureFluid:
+    """
+    This implementation of `PureFluid <cantera.with_units.PureFluid>` operates with
+    units by using the `pint` library to convert between unit systems. All properties
+    are assigned units in the standard MKS system that Cantera uses, substituting kmol
+    instead of mol. Each property is an instance of the `pint.Quantity` class.
+
+    Similarly, properties must be instances of `pint.Quantity` classes when they are
+    used for assignment to set the state. The properties may have any units, so long
+    as the dimensions for the quantity are consistent. For example, temperatures can
+    be provided in K, degC, degF, or degR; conversion will be done internally to
+    Cantera's consistent unit system.
+
+    See the `pint documentation <https://pint.readthedocs.io>`__ for more information
+    about using pint's ``Quantity`` classes.
+    """
+    def __init__(self, infile, name="", *, yaml=None, **kwargs):
+        self.__dict__["_phase"] = _PureFluid(infile, name, yaml=yaml, **kwargs)
+
+@common_properties@
+
+    @property
+    @copy_doc
+    def Q(self):
+        """Must be set using a quantity with dimensionless units."""
+        Q = self._phase.Q
+        return Q_(Q, "dimensionless")
+
+    @Q.setter
+    def Q(self, value):
+        if value is not None:
+            try:
+                Q = value.to("dimensionless").magnitude
+            except AttributeError as e:
+                if "'to'" in str(e):
+                    raise CanteraError(
+                        f"Value {value!r} must be an instance of a pint.Quantity class"
+                    ) from None
+                else:
+                    raise
+        else:
+            Q = self.Q.magnitude
+        self._phase.Q = Q
+
+@thermophase_properties@
+
+    @property
+    @copy_doc
+    def TPQ(self):
+        T, P, Q = self._phase.TPQ
+        return Q_(T, "K"), Q_(P, "Pa"), Q_(Q, "dimensionless")
+
+    @TPQ.setter
+    def TPQ(self, value):
+        T = value[0] if value[0] is not None else self.T
+        P = value[1] if value[1] is not None else self.P
+        Q = value[2] if value[2] is not None else self.Q
+        for val, unit in ((T, "K"), (P, "Pa"), (Q, "dimensionless")):
+            try:
+                val.ito(unit)
+            except AttributeError as e:
+                if "'ito'" in str(e):
+                    raise CanteraError(
+                        f"Value {val!r} must be an instance of a pint.Quantity class"
+                    ) from None
+                else:
+                    raise
+        self._phase.TPQ = T.magnitude, P.magnitude, Q.magnitude
+
+@purefluid_properties@
+
+
+PureFluid.__doc__ = f"{PureFluid.__doc__}\n{_PureFluid.__doc__}"
+
+
+def Heptane():
+    return PureFluid("liquidvapor.yaml", "heptane")
+
+Heptane.__doc__ = _Heptane.__doc__
+
+
+def CarbonDioxide():
+    return PureFluid("liquidvapor.yaml", "carbon-dioxide")
+
+
+CarbonDioxide.__doc__ = _CarbonDioxide.__doc__
+
+
+def Hfc134a():
+    return PureFluid("liquidvapor.yaml", "HFC-134a")
+
+
+Hfc134a.__doc__ = _Hfc134a.__doc__
+
+
+def Hydrogen():
+    return PureFluid("liquidvapor.yaml", "hydrogen")
+
+
+Hydrogen.__doc__ = _Hydrogen.__doc__
+
+
+def Methane():
+    return PureFluid("liquidvapor.yaml", "methane")
+
+
+Methane.__doc__ = _Methane.__doc__
+
+
+def Nitrogen():
+    return PureFluid("liquidvapor.yaml", "nitrogen")
+
+
+Nitrogen.__doc__ = _Nitrogen.__doc__
+
+
+def Oxygen():
+    return PureFluid("liquidvapor.yaml", "oxygen")
+
+
+Oxygen.__doc__ = _Oxygen.__doc__
+
+
+def Water(backend="Reynolds"):
+    return PureFluid("liquidvapor.yaml", "water", backend=backend)
+
+Water.__doc__ = _Water.__doc__
diff --git a/interfaces/cython/setup.cfg.in b/interfaces/cython/setup.cfg.in
index 0373a73817..1a0358a77d 100644
--- a/interfaces/cython/setup.cfg.in
+++ b/interfaces/cython/setup.cfg.in
@@ -38,7 +38,6 @@ project_urls =
 
 [options]
 zip_safe = False
-include_package_data = True
 install_requires =
     numpy >= 1.12.0
     ruamel.yaml >= 0.15.34
@@ -49,6 +48,7 @@ packages =
     cantera.data
     cantera.test
     cantera.examples
+    cantera.with_units
 
 [options.package_data]
 # The module extension needs to be here since we don't want setuptools to compile
@@ -61,10 +61,11 @@ cantera.examples = *.txt
 
 [options.extras_require]
 pandas = pandas
+units = pint
 
 [options.entry_points]
 console_scripts =
     ck2yaml = cantera.ck2yaml:script_entry_point
     cti2yaml = cantera.cti2yaml:main
     ctml2yaml = cantera.ctml2yaml:main
-	yaml2ck = cantera.yaml2ck:main
+    yaml2ck = cantera.yaml2ck:main
diff --git a/interfaces/python_minimal/setup.cfg.in b/interfaces/python_minimal/setup.cfg.in
index 96e0c6995c..2d96c75eda 100644
--- a/interfaces/python_minimal/setup.cfg.in
+++ b/interfaces/python_minimal/setup.cfg.in
@@ -45,4 +45,4 @@ console_scripts =
     ck2yaml = cantera.ck2yaml:script_entry_point
     cti2yaml = cantera.cti2yaml:main
     ctml2yaml = cantera.ctml2yaml:main
-	yaml2ck = cantera.yaml2ck:main
+    yaml2ck = cantera.yaml2ck:main
diff --git a/interfaces/python_sdist/MANIFEST.in b/interfaces/python_sdist/MANIFEST.in
index dc998757e1..9b1ae04f48 100644
--- a/interfaces/python_sdist/MANIFEST.in
+++ b/interfaces/python_sdist/MANIFEST.in
@@ -3,7 +3,6 @@
 # listed in setup.cfg:options.exclude_package_data.
 
 include cantera/_cantera.cpp
-include cantera/_cantera.h
 recursive-include cantera *.pyx
 include sundials_config.h.in
 include config.h.in
diff --git a/interfaces/python_sdist/SConscript b/interfaces/python_sdist/SConscript
index 1bf55332df..0aaf71b7e3 100644
--- a/interfaces/python_sdist/SConscript
+++ b/interfaces/python_sdist/SConscript
@@ -38,6 +38,8 @@ def copy_ext_src(target, source, env):
     if not FMT_SOURCE.is_dir():
         raise ValueError("Missing fmt submodule. " + error_message)
     for cc_file in (FMT_SOURCE / "src").glob("*.cc"):
+        if cc_file.name == "fmt.cc":
+            continue
         shutil.copy2(cc_file, FMT_ROOT)
     shutil.copytree(FMT_SOURCE / "include" / "fmt", FMT_ROOT / "fmt")
 
@@ -159,6 +161,10 @@ config_h_in_target = sdist(localenv.Command(
     config_builder,
 ))
 
+sdist(localenv.UnitsInterfaceBuilder(
+    "cantera/with_units/solution.py",
+    "#interfaces/cython/cantera/with_units/solution.py.in",
+))
 # Use RecursiveInstall to make sure that files are not overwritten during the copy.
 # A normal Copy Action would fail because of the existing directories.
 sdist(localenv.RecursiveInstall("cantera",
@@ -166,8 +172,6 @@ sdist(localenv.RecursiveInstall("cantera",
                                 exclude=["__pycache__"]))
 sdist(localenv.RecursiveInstall("cantera/data",
                                 "#build/data"))
-sdist(localenv.RecursiveInstall("cantera/test/data",
-                                "#test/data"))
 
 # Copy the minimal Sundials configuration template into the sdist so that
 # it can be filled in at compile time on the user's machine
diff --git a/interfaces/python_sdist/pyproject.toml b/interfaces/python_sdist/pyproject.toml
index a272bb9561..1bb9493c62 100644
--- a/interfaces/python_sdist/pyproject.toml
+++ b/interfaces/python_sdist/pyproject.toml
@@ -1,3 +1,10 @@
 [build-system]
-requires = ["setuptools>=43.0.0", "wheel", "oldest-supported-numpy", "Cython>=0.29.12"]
+# These versions are pinned to the latest versions at the time of this commit.
+# Feel free to update as required.
+requires = [
+  "setuptools==67.7.1",
+  "wheel",
+  "oldest-supported-numpy",
+  "Cython==0.29.34",
+]
 build-backend = "setuptools.build_meta"
diff --git a/interfaces/python_sdist/setup.cfg.in b/interfaces/python_sdist/setup.cfg.in
index 41d122a6db..77aa5557b7 100644
--- a/interfaces/python_sdist/setup.cfg.in
+++ b/interfaces/python_sdist/setup.cfg.in
@@ -45,15 +45,19 @@ install_requires =
     packaging
 python_requires @py_requires_ver_str@
 packages =
-    cantera
     cantera.data
+    cantera.examples
+    cantera.test
+    cantera.with_units
 
 # These options include data in the sdist and wheel if the files are also listed in
 # MANIFEST.in. Note that only files that are inside the "cantera" packages will be
 # included in the wheel.
 [options.package_data]
-cantera.data = *.*, */*.*
-cantera = *.pxd, test/*.txt, examples/*.txt
+cantera = *.pxd
+cantera.data = *.*
+cantera.examples = *.txt
+cantera.test = *.txt
 
 # These options exclude data from the wheel file but not from the sdist
 [options.exclude_package_data]
@@ -61,10 +65,11 @@ cantera = *.cpp, *.h, *.pyx
 
 [options.extras_require]
 pandas = pandas
+units = pint
 
 [options.entry_points]
 console_scripts =
     ck2yaml = cantera.ck2yaml:script_entry_point
     cti2yaml = cantera.cti2yaml:main
     ctml2yaml = cantera.ctml2yaml:main
-	yaml2ck = cantera.yaml2ck:main
+    yaml2ck = cantera.yaml2ck:main
diff --git a/interfaces/python_sdist/setup.py b/interfaces/python_sdist/setup.py
index 705d2608f0..127297a178 100644
--- a/interfaces/python_sdist/setup.py
+++ b/interfaces/python_sdist/setup.py
@@ -8,14 +8,13 @@
 import numpy
 import shutil
 
-HERE = Path(__file__).parent
-CT_SRC = HERE / "src"
-EXT_SRC = HERE / "ext"
-CT_INCLUDE = HERE / "include"
+CT_SRC = Path("src")
+EXT_SRC = Path("ext")
+CT_INCLUDE = Path("include")
 BOOST_INCLUDE = None
 FORCE_CYTHON_COMPILE = False
 
-CYTHON_BUILT_FILES = [HERE / "cantera" / f"_cantera.{ext}" for ext in ("cpp", "h")]
+CYTHON_BUILT_FILES = [Path("cantera") / f"_cantera.{ext}" for ext in ("cpp", "h")]
 
 
 class CanteraOptionsMixin:
@@ -114,9 +113,6 @@ def configure_build():
 
     config_h = {}
 
-    def create_config(key, value):
-        config_h[key] = f"#define {key} {value}"
-
     if not boost_version:
         raise ValueError(
             "Could not find Boost headers. Please set an environment variable called "
@@ -134,7 +130,7 @@ def create_config(key, value):
         )
 
     if sys.platform != "win32":
-        extra_compile_flags = ["-std=c++14", "-g0"]
+        extra_compile_flags = ["-std=c++17", "-g0"]
         sundials_configh = {
             "SUNDIALS_USE_GENERIC_MATH": "#define SUNDIALS_USE_GENERIC_MATH 1",
             "SUNDIALS_BLAS_LAPACK": "/* #undef SUNDIALS_BLAS_LAPACK */"
@@ -150,14 +146,14 @@ def create_config(key, value):
         }
         sundials_cflags = []
 
-    sun_config_h_in = (HERE / "sundials_config.h.in").read_text()
-    sun_config_h = HERE / "sundials_config.h"
+    sun_config_h_in = Path("sundials_config.h.in").read_text()
+    sun_config_h = Path("sundials_config.h")
     sun_config_h.write_text(sun_config_h_in.format_map(sundials_configh))
     shutil.copy2(sun_config_h, EXT_SRC / "sundials" / "sundials")
     shutil.copy2(sun_config_h, CT_INCLUDE / "cantera" / "ext" / "sundials")
 
-    config_h_in = (HERE / "config.h.in").read_text()
-    ct_config_h = HERE / "include" / "cantera" / "base" / "config.h"
+    config_h_in = Path("config.h.in").read_text()
+    ct_config_h = Path("include") / "cantera" / "base" / "config.h"
     ct_config_h.write_text(config_h_in.format_map(config_h))
 
     return extra_compile_flags, sundials_cflags, sundials_macros
diff --git a/samples/python/thermo/isentropic_units.py b/samples/python/thermo/isentropic_units.py
new file mode 100644
index 0000000000..a5a1a4f655
--- /dev/null
+++ b/samples/python/thermo/isentropic_units.py
@@ -0,0 +1,71 @@
+"""
+Isentropic, adiabatic flow example - calculate area ratio vs. Mach number curve.
+Uses the pint library to include customized units in the calculation.
+
+
+Requires: Cantera >= 3.0.0, pint
+Keywords: thermodynamics, compressible flow, units
+"""
+
+import cantera.with_units as ctu
+import numpy as np
+
+# This sets the default output format of the units to have 2 significant digits
+# and the units are printed with a Unicode font. See:
+# https://pint.readthedocs.io/en/stable/user/formatting.html
+ctu.units.default_format = ".2F~P"
+
+
+def soundspeed(gas):
+    """The speed of sound. Assumes an ideal gas."""
+
+    gamma = gas.cp / gas.cv
+    specific_gas_constant = ctu.units.molar_gas_constant / gas.mean_molecular_weight
+    return np.sqrt(gamma * specific_gas_constant * gas.T).to("m/s")
+
+
+def isentropic(gas=None):
+    """
+    In this example, the area ratio vs. Mach number curve is computed. If a gas
+    object is supplied, it will be used for the calculations, with the
+    stagnation state given by the input gas state. Otherwise, the calculations
+    will be done for a 10:1 hydrogen/nitrogen mixture with stagnation T0 = 1700.33
+    degrees Fahrenheit, P0 = 10 atm.
+    """
+    if gas is None:
+        gas = ctu.Solution('gri30.yaml')
+        gas.TPX = 2160 * ctu.units.degR, 10.0 * ctu.units.atm, 'H2:1,N2:0.1'
+
+    # get the stagnation state parameters
+    s0 = gas.s
+    h0 = gas.h
+    p0 = gas.P
+
+    mdot = 1 * ctu.units.kg / ctu.units.s  # arbitrary
+    amin = 1.e14 * ctu.units.m**2
+
+    data = []
+
+    # compute values for a range of pressure ratios
+    p_range = np.logspace(-3, 0, 10) * p0
+    for p in p_range:
+
+        # set the state using (p,s0)
+        gas.SP = s0, p
+
+        v = np.sqrt(2.0*(h0 - gas.h)).to("m/s")     # h + V^2/2 = h0
+        area = mdot/(gas.density*v)    # rho*v*A = constant
+        amin = min(amin, area)
+        data.append([area, v/soundspeed(gas), gas.T, p/p0])
+
+    return data, amin
+
+
+if __name__ == "__main__":
+    print(__doc__)
+    data, amin = isentropic()
+    label_string = "area ratio\tMach number\ttemperature\tpressure ratio"
+    output_string = "{0:.2E~P}\t{1}            {2}\t{3:.2E~P}"
+    print(label_string)
+    for row in data:
+        print(output_string.format(row[0] / amin, row[1], row[2], row[3]))
diff --git a/samples/python/thermo/rankine_units.py b/samples/python/thermo/rankine_units.py
new file mode 100644
index 0000000000..cae7d93bc9
--- /dev/null
+++ b/samples/python/thermo/rankine_units.py
@@ -0,0 +1,79 @@
+"""
+Calculate the efficiency of a Rankine vapor power cycle using a pure fluid model
+for water. Includes the units of quantities in the calculations.
+
+Requires: Cantera >= 3.0.0, pint
+Keywords: thermodynamics, thermodynamic cycle, non-ideal fluid, units
+"""
+
+import cantera.with_units as ctu
+
+# parameters
+eta_pump = 0.6 * ctu.units.dimensionless  # pump isentropic efficiency
+eta_turbine = 0.8 * ctu.units.dimensionless # turbine isentropic efficiency
+p_max = 116.03 * ctu.units.psi  # maximum pressure
+
+
+def pump(fluid, p_final, eta):
+    """Adiabatically pump a fluid to pressure p_final, using
+    a pump with isentropic efficiency eta."""
+    h0 = fluid.h
+    s0 = fluid.s
+    fluid.SP = s0, p_final
+    h1s = fluid.h
+    isentropic_work = h1s - h0
+    actual_work = isentropic_work / eta
+    h1 = h0 + actual_work
+    fluid.HP = h1, p_final
+    return actual_work
+
+
+def expand(fluid, p_final, eta):
+    """Adiabatically expand a fluid to pressure p_final, using
+    a turbine with isentropic efficiency eta."""
+    h0 = fluid.h
+    s0 = fluid.s
+    fluid.SP =s0, p_final
+    h1s = fluid.h
+    isentropic_work = h0 - h1s
+    actual_work = isentropic_work * eta
+    h1 = h0 - actual_work
+    fluid.HP = h1, p_final
+    return actual_work
+
+
+def print_state(n, fluid):
+    print('\n***************** State {0} ******************'.format(n))
+    print(fluid.report())
+
+
+if __name__ == '__main__':
+    # create an object representing water
+    w = ctu.Water()
+
+    # start with saturated liquid water at 80.33 degrees Fahrenheit
+    w.TQ = ctu.Q_(80.33, "degF"), 0.0 * ctu.units.dimensionless
+    h1 = w.h
+    p1 = w.P
+    print_state(1, w)
+
+    # pump it adiabatically to p_max
+    pump_work = pump(w, p_max, eta_pump)
+    h2 = w.h
+    print_state(2, w)
+
+    # heat it at constant pressure until it reaches the saturated vapor state
+    # at this pressure
+    w.PQ = p_max, 1.0 * ctu.units.dimensionless
+    h3 = w.h
+    heat_added = h3 - h2
+    print_state(3, w)
+
+    # expand back to p1
+    turbine_work = expand(w, p1, eta_turbine)
+    print_state(4, w)
+
+    # efficiency
+    eff = (turbine_work - pump_work)/heat_added
+
+    print('efficiency = ', eff)
diff --git a/samples/python/thermo/sound_speed.py b/samples/python/thermo/sound_speed.py
index 3a152fc2e4..5a93ba726f 100644
--- a/samples/python/thermo/sound_speed.py
+++ b/samples/python/thermo/sound_speed.py
@@ -7,6 +7,7 @@
 
 import cantera as ct
 import math
+import numpy as np
 
 
 def equilSoundSpeeds(gas, rtol=1.0e-6, max_iter=5000):
@@ -53,7 +54,7 @@ def equilSoundSpeeds(gas, rtol=1.0e-6, max_iter=5000):
 if __name__ == "__main__":
     gas = ct.Solution('gri30.yaml')
     gas.X = 'CH4:1.00, O2:2.0, N2:7.52'
-    for n in range(27):
-        T = 300.0 + 100.0 * n
+    T_range = np.arange(300, 2901, 100)
+    for T in T_range:
         gas.TP = T, ct.one_atm
         print(T, equilSoundSpeeds(gas))
diff --git a/samples/python/thermo/sound_speed_units.py b/samples/python/thermo/sound_speed_units.py
new file mode 100644
index 0000000000..7500ae6a3b
--- /dev/null
+++ b/samples/python/thermo/sound_speed_units.py
@@ -0,0 +1,65 @@
+"""
+Compute the "equilibrium" and "frozen" sound speeds for a gas. Uses the pint library to
+include customized units in the calculation.
+
+Requires: Cantera >= 3.0.0, pint
+Keywords: thermodynamics, equilibrium, units
+"""
+
+import cantera.with_units as ctu
+import numpy as np
+
+# This sets the default output format of the units to have 2 significant digits
+# and the units are printed with a Unicode font. See:
+# https://pint.readthedocs.io/en/stable/formatting.html#unit-format-types
+ctu.units.default_format = ".2F~P"
+
+def equilibrium_sound_speeds(gas, rtol=1.0e-6, max_iter=5000):
+    """
+    Returns a tuple containing the equilibrium and frozen sound speeds for a
+    gas with an equilibrium composition.  The gas is first set to an
+    equilibrium state at the temperature and pressure of the gas, since
+    otherwise the equilibrium sound speed is not defined.
+    """
+
+    # set the gas to equilibrium at its current T and P
+    gas.equilibrate('TP', rtol=rtol, max_iter=max_iter)
+
+    # save properties
+    s0 = gas.s
+    p0 = gas.P
+    r0 = gas.density
+
+    # perturb the pressure
+    p1 = p0*1.0001
+
+    # set the gas to a state with the same entropy and composition but
+    # the perturbed pressure
+    gas.SP = s0, p1
+
+    # frozen sound speed
+    afrozen = np.sqrt((p1 - p0)/(gas.density - r0)).to("ft/s")
+
+    # now equilibrate the gas holding S and P constant
+    gas.equilibrate('SP', rtol=rtol, max_iter=max_iter)
+
+    # equilibrium sound speed
+    aequil = np.sqrt((p1 - p0)/(gas.density - r0)).to("ft/s")
+
+    # compute the frozen sound speed using the ideal gas expression as a check
+    gamma = gas.cp/gas.cv
+    gamma * ctu.units.molar_gas_constant
+    afrozen2 = np.sqrt(gamma * ctu.units.molar_gas_constant * gas.T /
+                         gas.mean_molecular_weight).to("ft/s")
+
+    return aequil, afrozen, afrozen2
+
+# test program
+if __name__ == "__main__":
+    gas = ctu.Solution('gri30.yaml')
+    gas.X = 'CH4:1.00, O2:2.0, N2:7.52'
+    T_range = np.linspace(80, 4880, 25) * ctu.units.degF
+    print("Temperature      Equilibrium Sound Speed     Frozen Sound Speed      Frozen Sound Speed Check")
+    for T in T_range:
+        gas.TP = T, 1.0 * ctu.units.atm
+        print(T.to("degF"), *equilibrium_sound_speeds(gas), sep = "               ")
diff --git a/site_scons/buildutils.py b/site_scons/buildutils.py
index 1b560e8334..601653338a 100644
--- a/site_scons/buildutils.py
+++ b/site_scons/buildutils.py
@@ -2,7 +2,6 @@
 import json
 import os
 import sys
-import platform
 import textwrap
 import re
 import subprocess
diff --git a/site_scons/site_tools/UnitsInterfaceBuilder.py b/site_scons/site_tools/UnitsInterfaceBuilder.py
new file mode 100644
index 0000000000..0c9af9a51b
--- /dev/null
+++ b/site_scons/site_tools/UnitsInterfaceBuilder.py
@@ -0,0 +1,248 @@
+from textwrap import dedent, indent
+from string import Template
+
+def UnitsInterfaceBuilder(env, target, source):
+    """This builder creates the cantera.with_units interface for the Python package.
+
+    This builder is meant to be called like
+    ``env.UnitsInterfaceBuilder(target_file, template_source_file)``
+
+    The builder will create string templates for all the ``ThermoPhase`` methods and
+    fill them in with the appropriate SI + kmol units to be converted with pint.
+
+    The return value is an item of type ``env.SubstFile`` which can be put into a
+    ``Depends`` call, so that this builder will be called for a particular module.
+    """
+    UNITS = {
+        "cp_mass": '"J/kg/K"', "cp_mole": '"J/kmol/K"', "cv_mass": '"J/kg/K"',
+        "cv_mole": '"J/kmol/K"', "density_mass": '"kg/m**3"', "density_mole": '"kmol/m**3"',
+        "enthalpy_mass": '"J/kg"', "enthalpy_mole": '"J/kmol"', "entropy_mass": '"J/kg/K"',
+        "entropy_mole": '"J/kmol/K"', "gibbs_mass": '"J/kg"', "gibbs_mole": '"J/kmol"',
+        "int_energy_mass": '"J/kg"', "int_energy_mole": '"J/kmol"',
+        "volume_mass": '"m**3/kg"', "volume_mole": '"m**3/kmol"', "T": '"K"', "P": '"Pa"',
+        "X": '"dimensionless"', "Y": '"dimensionless"', "Q": '"dimensionless"',
+        "cp": '"J/K/" + self.basis_units', "cv": '"J/K/" + self.basis_units',
+        "density": 'self.basis_units + "/m**3"', "h": '"J/" + self.basis_units',
+        "s": '"J/K/" + self.basis_units', "g": '"J/" + self.basis_units',
+        "u": '"J/" + self.basis_units', "v": '"m**3/" + self.basis_units',
+        "H": '"J/" + self.basis_units', "V": '"m**3/" + self.basis_units',
+        "S": '"J/K/" + self.basis_units', "D": 'self.basis_units + "/m**3"',
+        "U": '"J/" + self.basis_units', "P_sat": '"Pa"', "T_sat": '"K"',
+        "atomic_weight": '"kg/kmol"', "chemical_potentials": '"J/kmol"',
+        "concentrations": '"kmol/m**3"', "critical_pressure": '"Pa"',
+        "critical_temperature": '"K"', "critical_density": 'self.basis_units + "/m**3"',
+        "electric_potential": '"V"', "electrochemical_potentials": '"J/kmol"',
+        "isothermal_compressibility": '"1/Pa"', "max_temp": '"K"',
+        "mean_molecular_weight": '"kg/kmol"', "min_temp": '"K"',
+        "molecular_weights": '"kg/kmol"', "partial_molar_cp": '"J/kmol/K"',
+        "partial_molar_enthalpies": '"J/kmol"', "partial_molar_entropies": '"J/kmol/K"',
+        "partial_molar_int_energies": '"J/kmol"', "partial_molar_volumes": '"m**3/kmol"',
+        "reference_pressure": '"Pa"', "thermal_expansion_coeff": '"1/K"'
+    }
+
+    SYMBOL = {
+        "T": "T", "P": "P", "D": "density", "H": "h", "S": "s",
+        "V": "v", "U": "u", "Q": "Q", "X": "X", "Y": "Y"
+    }
+
+    getter_properties = [
+        "density_mass", "density_mole", "enthalpy_mass", "enthalpy_mole", "entropy_mass",
+        "entropy_mole", "int_energy_mass", "int_energy_mole", "volume_mass", "volume_mole",
+        "gibbs_mass", "gibbs_mole", "cp_mass", "cp_mole", "cv_mass", "cv_mole", "P",
+        "P_sat", "T", "T_sat", "atomic_weight", "chemical_potentials", "concentrations",
+        "critical_pressure", "critical_temperature", "critical_density",
+        "electric_potential", "electrochemical_potentials", "isothermal_compressibility",
+        "max_temp", "mean_molecular_weight", "min_temp", "molecular_weights",
+        "partial_molar_cp", "partial_molar_enthalpies", "partial_molar_entropies",
+        "partial_molar_int_energies", "partial_molar_volumes", "reference_pressure",
+        "thermal_expansion_coeff", "cp", "cv", "density", "h", "s", "g", "u", "v",
+    ]
+
+    getter_template = Template(dedent("""
+        @property
+        @copy_doc
+        def ${name}(self):
+            return Q_(self._phase.${name}, ${units})
+    """))
+
+    thermophase_getters = []
+    for name in getter_properties:
+        thermophase_getters.append(getter_template.substitute(name=name, units=UNITS[name]))
+
+    setter_2_template = Template(dedent("""
+        @property
+        @copy_doc
+        def ${name}(self):
+            ${n0}, ${n1} = self._phase.${name}
+            return Q_(${n0}, ${u0}), Q_(${n1}, ${u1})
+
+        @${name}.setter
+        def ${name}(self, value):
+            ${n0} = value[0] if value[0] is not None else self.${n0}
+            ${n1} = value[1] if value[1] is not None else self.${n1}
+            for val, unit in ((${n0}, ${u0}), (${n1}, ${u1})):
+                try:
+                    val.ito(unit)
+                except AttributeError as e:
+                    if "'ito'" in str(e):
+                        raise CanteraError(
+                            f"Value {val!r} must be an instance of a pint.Quantity class"
+                        ) from None
+                    else:
+                        raise
+            self._phase.${name} = ${n0}.magnitude, ${n1}.magnitude
+    """))
+
+    tp_setter_2_properties = ["TP", "DP", "HP", "SP", "SV", "TD", "UV"]
+    pf_setter_2_properties = ["PQ", "TQ", "PV", "SH", "ST", "TH", "TV", "UP", "VH"]
+
+    thermophase_2_setters = []
+    for name in tp_setter_2_properties:
+        d = dict(name=name, n0=SYMBOL[name[0]], u0=UNITS[name[0]], n1=SYMBOL[name[1]],
+                u1=UNITS[name[1]])
+        thermophase_2_setters.append(setter_2_template.substitute(d))
+
+    purefluid_2_setters = []
+    for name in pf_setter_2_properties:
+        d = dict(name=name, n0=SYMBOL[name[0]], u0=UNITS[name[0]], n1=SYMBOL[name[1]],
+                u1=UNITS[name[1]])
+        purefluid_2_setters.append(setter_2_template.substitute(d))
+
+    setter_3_template = Template(dedent("""
+        @property
+        @copy_doc
+        def ${name}(self):
+            ${n0}, ${n1}, ${n2} = self._phase.${name}
+            return Q_(${n0}, ${u0}), Q_(${n1}, ${u1}), Q_(${n2}, ${u2})
+
+        @${name}.setter
+        def ${name}(self, value):
+            ${n0} = value[0] if value[0] is not None else self.${n0}
+            ${n1} = value[1] if value[1] is not None else self.${n1}
+            for val, unit in ((${n0}, ${u0}), (${n1}, ${u1})):
+                try:
+                    val.ito(unit)
+                except AttributeError as e:
+                    if "'ito'" in str(e):
+                        raise CanteraError(
+                            f"Value {val!r} must be an instance of a pint.Quantity class"
+                        ) from None
+                    else:
+                        raise
+            if value[2] is not None:
+                try:
+                    ${n2} = value[2].to(${u2}).magnitude
+                except AttributeError:
+                    ${n2} = value[2]
+            else:
+                ${n2} = self.${n2}.magnitude
+            self._phase.${name} = ${n0}.magnitude, ${n1}.magnitude, ${n2}
+    """))
+
+    tp_setter_3_properties = [
+        "TPX", "TPY", "DPX", "DPY", "HPX", "HPY", "SPX", "SPY", "SVX", "SVY", "TDX", "TDY",
+        "UVX", "UVY"
+    ]
+
+    thermophase_3_setters = []
+    for name in tp_setter_3_properties:
+        d = dict(name=name, n0=SYMBOL[name[0]], u0=UNITS[name[0]], n1=SYMBOL[name[1]],
+                u1=UNITS[name[1]], n2=SYMBOL[name[2]], u2=UNITS[name[2]])
+        thermophase_3_setters.append(setter_3_template.substitute(d))
+
+    getter_3_template = Template(dedent("""
+        @property
+        @copy_doc
+        def ${name}(self):
+            ${n0}, ${n1}, ${n2} = self._phase.${name}
+            return Q_(${n0}, ${u0}), Q_(${n1}, ${u1}), Q_(${n2}, ${u2})
+    """))
+
+    pf_getter_3_properties = ["DPQ", "HPQ", "SPQ", "SVQ", "TDQ", "UVQ"]
+
+    purefluid_3_getters = []
+    for name in pf_getter_3_properties:
+        d = dict(name=name, n0=name[0], u0=UNITS[name[0]], n1=name[1], u1=UNITS[name[1]],
+                n2=name[2], u2=UNITS[name[2]])
+        purefluid_3_getters.append(getter_3_template.substitute(d))
+
+
+    def recursive_join(*args, joiner=""):
+        result = ""
+        for arg in args:
+            result = result + joiner.join(arg)
+        return result
+
+
+    thermophase_properties = recursive_join(thermophase_getters, thermophase_2_setters,
+                                            thermophase_3_setters)
+    purefluid_properties = recursive_join(purefluid_2_setters, purefluid_3_getters)
+
+    common_properties = dedent("""
+        def __getattr__(self, name):
+            return getattr(self._phase, name)
+
+        def __setattr__(self, name, value):
+            if name in dir(self):
+                object.__setattr__(self, name, value)
+            else:
+                setattr(self._phase, name, value)
+
+        @property
+        def basis_units(self):
+            \"\"\"The units associated with the mass/molar basis of this phase.\"\"\"
+            if self._phase.basis == "mass":
+                return "kg"
+            else:
+                return "kmol"
+
+        @property
+        @copy_doc
+        def X(self):
+            \"\"\"If an array is used for setting, the units must be dimensionless.\"\"\"
+            X = self._phase.X
+            return Q_(X, "dimensionless")
+
+        @X.setter
+        def X(self, value):
+            if value is not None:
+                try:
+                    X = value.to("dimensionless").magnitude
+                except AttributeError:
+                    X = value
+            else:
+                X = self.X.magnitude
+            self._phase.X = X
+
+        @property
+        @copy_doc
+        def Y(self):
+            \"\"\"If an array is used for setting, the units must be dimensionless.\"\"\"
+            Y = self._phase.Y
+            return Q_(Y, "dimensionless")
+
+        @Y.setter
+        def Y(self, value):
+            if value is not None:
+                try:
+                    Y = value.to("dimensionless").magnitude
+                except AttributeError:
+                    Y = value
+            else:
+                Y = self.Y.magnitude
+            self._phase.Y = Y
+    """)
+
+    env["common_properties"] = indent(common_properties.strip(), " "*4)
+    env["thermophase_properties"] = indent(thermophase_properties.strip(), " "*4)
+    env["purefluid_properties"] = indent(purefluid_properties.strip(), " "*4)
+    units = env.SubstFile(target, source)
+    return units
+
+
+def generate(env, **kw):
+    env.AddMethod(UnitsInterfaceBuilder)
+
+
+def exists(env):
+    return True
diff --git a/test/python/conftest.py b/test/python/conftest.py
new file mode 100644
index 0000000000..debb72318b
--- /dev/null
+++ b/test/python/conftest.py
@@ -0,0 +1,3 @@
+import pytest
+
+pytest.register_assert_rewrite("pint.testing")
diff --git a/test/python/coverage.ini b/test/python/coverage.ini
index 14a8e48bdd..80e95de18d 100644
--- a/test/python/coverage.ini
+++ b/test/python/coverage.ini
@@ -4,3 +4,4 @@ omit =
     *.cti
     stringsource
 source = ../../build/python/cantera
+branch = true
diff --git a/test/python/test_units.py b/test/python/test_units.py
new file mode 100644
index 0000000000..67a8db610f
--- /dev/null
+++ b/test/python/test_units.py
@@ -0,0 +1,573 @@
+from contextlib import nullcontext
+from dataclasses import dataclass
+from typing import Optional, Tuple, Dict
+import sys
+
+import pytest
+import cantera.with_units as ctu
+import cantera as ct
+try:
+    from pint.testing import assert_allclose
+except ModuleNotFoundError:
+    # pint.testing was introduced in pint 0.20
+    from pint.testsuite.helpers import assert_quantity_almost_equal as assert_allclose
+
+
+@pytest.fixture(scope="function")
+def ideal_gas():
+    return ctu.Solution("h2o2.yaml")
+
+
+@pytest.fixture(scope="function")
+def pure_fluid():
+    return ctu.Water()
+
+
+@pytest.fixture(params=["ideal_gas", "pure_fluid"])
+def generic_phase(request):
+    return request.getfixturevalue(request.param)
+
+
+def test_setting_basis_units_fails(generic_phase):
+    with pytest.raises(AttributeError):
+        generic_phase.basis_units = "some random string"
+
+
+def test_mass_basis(generic_phase):
+    """Check that mass basis units have kg and the generic getter returns the same
+    value as the mass-specific getter."""
+    generic_phase.basis = "mass"
+    assert generic_phase.basis_units == "kg"
+    assert_allclose(generic_phase.density_mass, generic_phase.density)
+    assert_allclose(generic_phase.enthalpy_mass, generic_phase.h)
+    assert_allclose(generic_phase.entropy_mass, generic_phase.s)
+    assert_allclose(generic_phase.int_energy_mass, generic_phase.u)
+    assert_allclose(generic_phase.volume_mass, generic_phase.v)
+    assert_allclose(generic_phase.gibbs_mass, generic_phase.g)
+    assert_allclose(generic_phase.cp_mass, generic_phase.cp)
+    assert_allclose(generic_phase.cv_mass, generic_phase.cv)
+
+
+def test_molar_basis(generic_phase):
+    """Check that molar basis units have kmol and the generic getter returns the
+    same value as the molar-specific getter."""
+    generic_phase.basis = "molar"
+    assert generic_phase.basis_units == "kmol"
+    assert_allclose(generic_phase.density_mole, generic_phase.density)
+    assert_allclose(generic_phase.enthalpy_mole, generic_phase.h)
+    assert_allclose(generic_phase.entropy_mole, generic_phase.s)
+    assert_allclose(generic_phase.int_energy_mole, generic_phase.u)
+    assert_allclose(generic_phase.volume_mole, generic_phase.v)
+    assert_allclose(generic_phase.gibbs_mole, generic_phase.g)
+    assert_allclose(generic_phase.cp_mole, generic_phase.cp)
+    assert_allclose(generic_phase.cv_mole, generic_phase.cv)
+
+
+@dataclass(frozen=True)
+class Dimensions:
+    name: str
+    mass: Optional[float] = None
+    length: Optional[float] = None
+    time: Optional[float] = None
+    substance: Optional[float] = None
+    temperature: Optional[float] = None
+    current: Optional[float] = None
+    dimensions: Tuple[str, ...] = (
+        "mass",
+        "length",
+        "time",
+        "substance",
+        "temperature",
+        "current",
+    )
+
+    def __str__(self):
+        return self.name
+
+    def inverse(self, name: Optional[str] = None) -> "Dimensions":
+        dimensionality = {}
+        for dimension in self.dimensions:
+            value = getattr(self, dimension)
+            if value is not None:
+                dimensionality[dimension] = -1 * value
+        new_name = name if name is not None else f"inverse_{self.name}"
+        return Dimensions(name=new_name, **dimensionality)
+
+    def as_dict(self) -> Dict[str, float]:
+        """Add the square brackets around the dimension for comparison with pint"""
+        dimensionality = {}
+        for dimension in self.dimensions:
+            value = getattr(self, dimension)
+            if value is not None:
+                dimensionality[f"[{dimension}]"] = value
+        return dimensionality
+
+
+# Create instances of Dimensions that correspond to all the combinations of dimensions
+# that are implemented in the with_units interface.
+temperature = Dimensions(temperature=1, name="temperature")
+pressure = Dimensions(mass=1, length=-1, time=-2, name="pressure")
+isothermal_compressiblity = pressure.inverse()
+inverse_temperature = temperature.inverse()
+atomic_molecular_weights = Dimensions("atomic_molecular_weights", mass=1, substance=-1)
+mole_mass_fractions = Dimensions(name="mole_mass_fractions")
+chemical_potential = Dimensions(
+    name="chemical_potential", mass=1, length=2, time=-2, substance=-1
+)
+electric_potential = Dimensions(
+    name="electric_potential", mass=1, length=2, time=-3, current=-1
+)
+concentrations_like = Dimensions(name="concentrations_like", substance=1, length=-3)
+molar_volume = Dimensions(name="volume_mole", substance=-1, length=3)
+volume_mass = Dimensions(name="volume_mass", mass=-1, length=3)
+density_mass = volume_mass.inverse(name="density_mass")
+mass_basis_energy_like = Dimensions(name="mass_basis_energy_like", length=2, time=-2)
+mass_basis_entropy_like = Dimensions(
+    name="mass_basis_entropy_like", length=2, time=-2, temperature=-1
+)
+molar_basis_energy_like = Dimensions(
+    name="molar_basis_energy_like", mass=1, length=2, time=-2, substance=-1
+)
+molar_basis_entropy_like = Dimensions(
+    name="molar_basis_entropy_like",
+    mass=1,
+    length=2,
+    time=-2,
+    substance=-1,
+    temperature=-1,
+)
+
+
+def yield_dimensions():
+    """Yield pytest.param instances with the dimensions"""
+    # This dictionary maps the dimensions to the relevant property names
+    dims: Dict[Dimensions, Tuple[str, ...]] = {
+        temperature: ("T", "max_temp", "min_temp"),
+        inverse_temperature: ("thermal_expansion_coeff",),
+        pressure: ("P", "reference_pressure"),
+        isothermal_compressiblity: ("isothermal_compressibility",),
+        atomic_molecular_weights: (
+            "atomic_weight",
+            "molecular_weights",
+            "mean_molecular_weight",
+        ),
+        mole_mass_fractions: ("X", "Y"),
+        chemical_potential: ("chemical_potentials", "electrochemical_potentials"),
+        electric_potential: ("electric_potential",),
+        concentrations_like: ("concentrations", "density_mole"),
+        molar_volume: ("volume_mole", "partial_molar_volumes"),
+        volume_mass: ("volume_mass",),
+        density_mass: ("density_mass",),
+        mass_basis_energy_like: ("enthalpy_mass", "int_energy_mass", "gibbs_mass"),
+        mass_basis_entropy_like: ("entropy_mass", "cp_mass", "cv_mass"),
+        molar_basis_energy_like: (
+            "enthalpy_mole",
+            "int_energy_mole",
+            "gibbs_mole",
+            "partial_molar_enthalpies",
+            "partial_molar_int_energies",
+        ),
+        molar_basis_entropy_like: (
+            "entropy_mole",
+            "cp_mole",
+            "cv_mole",
+            "partial_molar_cp",
+            "partial_molar_entropies",
+        ),
+    }
+    for dimension, props in dims.items():
+        for prop in props:
+            yield pytest.param(prop, dimension.as_dict(), id=f"{dimension}-{prop}")
+
+
+@pytest.mark.parametrize("prop,dimensions", yield_dimensions())
+def test_dimensions(generic_phase, prop, dimensions):
+    """Test that the dimensions returned for a property are correct.
+
+    Arguments
+    =========
+
+    generic_phase: A phase definition created from cantera.with_units objects.
+        Currently, one of `Solution` or `PureFluid`. Created by the generic_phase
+        fixture.
+    prop: A string of the property name
+    dimensions: The known dimensions for this property
+
+    The latter two arguments are supplied by the parametrize on this test. That
+    parametrize is effectively a loop over all the implemented properties on the
+    classes. The loop is implemented in the `yield_dimensions()` function. The
+    parametrize call is kinda like calling ``list(generator_function())`` to
+    discover all the values in the generator, except pytest does that automatically
+    for us and fills in the arguments.
+    """
+    pint_dim = dict(getattr(generic_phase, prop).dimensionality)
+    assert pint_dim == dimensions
+
+
+@pytest.mark.parametrize(
+    "phase",
+    (
+        pytest.param(ctu.Solution("liquidvapor.yaml", "heptane"), id="Solution"),
+        pytest.param(ctu.Water(), id="PureFluid"),
+    ),
+)
+def test_purefluid_dimensions(phase):
+    # Test some dimensions that weren't tested as part of the Solution tests
+    # Create and test a liquidvapor phase in a Solution object, since an ideal gas phase
+    # doesn't implement saturation or critical properties.
+    assert dict(phase.T_sat.dimensionality) == temperature.as_dict()
+    assert dict(phase.critical_temperature.dimensionality) == temperature.as_dict()
+    assert dict(phase.P_sat.dimensionality) == pressure.as_dict()
+    assert dict(phase.critical_pressure.dimensionality) == pressure.as_dict()
+    assert dict(phase.critical_density.dimensionality) == density_mass.as_dict()
+
+
+def yield_prop_pairs():
+    pairs = [
+        pytest.param(("TP", "T", "P"), id="TP"),
+        pytest.param(("SP", "s", "P"), id="SP"),
+        pytest.param(("UV", "u", "v"), id="UV"),
+        pytest.param(("DP", "density", "P"), id="DP"),
+        pytest.param(("HP", "h", "P"), id="HP"),
+        pytest.param(("SV", "s", "v"), id="SV"),
+        pytest.param(("TD", "T", "density"), id="TD"),
+    ]
+    yield from pairs
+
+
+def yield_prop_triples():
+    for pair in yield_prop_pairs():
+        values = pair.values[0]
+        yield pytest.param(
+            (values[0] + "X", *values[1:], "X"),
+            id=pair.id + "X",
+        )
+        yield pytest.param(
+            (values[0] + "Y", *values[1:], "Y"),
+            id=pair.id + "Y",
+        )
+
+
+def yield_prop_pairs_and_triples():
+    yield from yield_prop_pairs()
+    yield from yield_prop_triples()
+
+
+@pytest.fixture
+def TD_in_the_right_basis(request):
+    if request.param == "mass":
+        T = ctu.Q_(500, "K")
+        rho = ctu.Q_(1.5, "kg/m**3")
+    elif request.param == "molar":
+        T = ctu.Q_(750, "K")
+        rho = ctu.Q_(0.02, "kmol/m**3")
+    return (T, rho)
+
+
+@pytest.fixture
+def initial_TDY(request, generic_phase):
+    generic_phase.basis = request.param
+    return generic_phase.TDY
+
+
+@pytest.fixture
+def some_setters_arent_implemented_for_purefluid(request):
+    pair_or_triple = request.getfixturevalue("props")[0]
+    is_pure_fluid = isinstance(request.getfixturevalue("generic_phase"), ctu.PureFluid)
+    if is_pure_fluid and pair_or_triple.startswith("DP"):
+        request.applymarker(
+            pytest.mark.xfail(
+                raises=NotImplementedError,
+                reason=f"The {pair_or_triple} method isn't implemented",
+            )
+        )
+
+
+# The parameterization is done here with the indirect kwarg to make sure that the same
+# basis is passed to both fixtures. The alternative is to use the params kwarg to the
+# fixture decorator, which would give us (mass, molar) basis pairs, and that doesn't
+# make sense.
+@pytest.mark.parametrize(
+    "TD_in_the_right_basis,initial_TDY",
+    [
+        pytest.param("mass", "mass", id="mass"),
+        pytest.param("molar", "molar", id="molar"),
+    ],
+    indirect=True,
+)
+@pytest.mark.parametrize("props", yield_prop_pairs_and_triples())
+@pytest.mark.usefixtures("some_setters_arent_implemented_for_purefluid")
+def test_setters(generic_phase, TD_in_the_right_basis, initial_TDY, props):
+    pair_or_triple = props[0]
+    if isinstance(generic_phase, ctu.PureFluid):
+        Y_1 = ctu.Q_([1.0], "dimensionless")
+    else:
+        Y_1 = ctu.Q_(
+            [0.1, 0.0, 0.0, 0.1, 0.4, 0.2, 0.0, 0.0, 0.2, 0.0], "dimensionless"
+        )
+    generic_phase.TDY = *TD_in_the_right_basis, Y_1
+
+    # Use TDY setting to get the properties at the modified state
+    new_props = getattr(generic_phase, pair_or_triple)
+
+    # Reset to the initial state so that the next state setting actually has to do
+    # something.
+    generic_phase.TDY = initial_TDY
+
+    # If we're only setting a pair of properties, reset the mass fractions to the
+    # expected state before using the pair to set.
+    if len(pair_or_triple) == 2:
+        generic_phase.Y = Y_1
+
+    # Use the test pair or triple to set the state and assert that the
+    # natural properties are equal to the modified state
+    setattr(generic_phase, pair_or_triple, new_props)
+    T_1, rho_1 = TD_in_the_right_basis
+    assert_allclose(generic_phase.T, T_1)
+    assert_allclose(generic_phase.density, rho_1)
+    assert_allclose(generic_phase.Y, Y_1)
+
+
+@pytest.mark.parametrize("props", yield_prop_triples())
+@pytest.mark.usefixtures("some_setters_arent_implemented_for_purefluid")
+def test_setters_hold_constant(generic_phase, props):
+    triple, first, second, third = props
+
+    # Set an arbitrary initial state
+    if generic_phase.n_species == 1:
+        generic_phase.X = ctu.Q_([1.0], "dimensionless")
+        composition = ctu.Q_([1.0], "dimensionless")
+    else:
+        generic_phase.X = "H2O:0.1, O2:0.95, AR:3.0"
+        composition = "H2:0.1, O2:1.0, AR:3.0"
+    generic_phase.TD = ctu.Q_(1000, "K"), ctu.Q_(1.5, "kg/m**3")
+    property_3 = getattr(generic_phase, third)
+
+    # Change to another arbitrary state and store values to compare when a property
+    # isn't changed
+    reset_state = (ctu.Q_(500, "K"), ctu.Q_(2.5, "kg/m**3"), composition)
+    generic_phase.TDX = reset_state
+    first_val, second_val, third_val = getattr(generic_phase, triple)
+
+    setattr(generic_phase, triple, (None, None, property_3))
+    assert_allclose(getattr(generic_phase, first), first_val)
+    assert_allclose(getattr(generic_phase, second), second_val)
+    assert_allclose(getattr(generic_phase, third), property_3)
+
+    generic_phase.TDX = reset_state
+    setattr(generic_phase, triple, (None, None, None))
+    assert_allclose(getattr(generic_phase, first), first_val)
+    assert_allclose(getattr(generic_phase, second), second_val)
+    assert_allclose(getattr(generic_phase, third), third_val)
+
+
+@pytest.mark.parametrize("props", yield_prop_pairs_and_triples())
+@pytest.mark.parametrize(
+    "basis,rho_0",
+    [
+        pytest.param("mass", ctu.Q_(0.7, "kg/m**3"), id="mass"),
+        pytest.param("molar", ctu.Q_(0.01, "kmol/m**3"), id="molar"),
+    ],
+)
+def test_multi_prop_getters_are_equal_to_single(generic_phase, props, basis, rho_0):
+    pair_or_triple, first, second, *third = props
+    generic_phase.basis = basis
+    if generic_phase.n_species != 1:
+        generic_phase.Y = "H2:0.1, H2O2:0.1, AR:0.8"
+    generic_phase.TD = ctu.Q_(350.0, "K"), rho_0
+    # This test is equivalent to
+    # T, P, X = solution.TPX
+    # assert isclose(T, solution.T)
+    # assert isclose(P, solution.P)
+    # assert all(isclose(X, solution.X))
+    # Where T, P, and X loop through all the valid property pairs and triples, for
+    # both mass and molar basis units.
+    first_value, second_value, *third_value = getattr(generic_phase, pair_or_triple)
+    assert_allclose(getattr(generic_phase, first), first_value)
+    assert_allclose(getattr(generic_phase, second), second_value)
+    if third:
+        assert_allclose(getattr(generic_phase, third[0]), third_value[0])
+
+
+@pytest.mark.parametrize("pair", yield_prop_pairs())
+def test_set_pair_without_units_is_an_error(generic_phase, pair):
+    with pytest.raises(ctu.CanteraError, match="an instance of a pint"):
+        setattr(generic_phase, pair[0], [300, None])
+
+
+
+@pytest.mark.parametrize("triple", yield_prop_triples())
+def test_set_triple_without_units_is_an_error(generic_phase, triple):
+    value_1 = [300, None, [1] + [0] * (generic_phase.n_species - 1)]
+    with pytest.raises(ctu.CanteraError, match="an instance of a pint"):
+        setattr(generic_phase, triple[0], value_1)
+
+@pytest.mark.parametrize("props", yield_prop_triples())
+@pytest.mark.usefixtures("some_setters_arent_implemented_for_purefluid")
+def test_set_triple_with_no_units_on_composition_succeeds(
+    generic_phase,
+    props,
+):
+    value_3 = [None, None, [1] + [0] * (generic_phase.n_species - 1)]
+    setattr(generic_phase, props[0], value_3)
+    assert_allclose(getattr(generic_phase, props[0][2]), value_3[2])
+
+
+@pytest.fixture
+def pure_fluid_in_vapordome(pure_fluid):
+    pure_fluid.TQ = None, ctu.Q_(0.5, "dimensionless")
+    return pure_fluid
+
+
+def test_set_Q(pure_fluid_in_vapordome):
+    P = pure_fluid_in_vapordome.P
+    T = pure_fluid_in_vapordome.T
+    pure_fluid_in_vapordome.Q = ctu.Q_(0.6, "dimensionless")
+    assert_allclose(pure_fluid_in_vapordome.Q, 0.6 * ctu.units.dimensionless)
+    assert_allclose(pure_fluid_in_vapordome.T, T)
+    assert_allclose(pure_fluid_in_vapordome.P, P)
+
+    pure_fluid_in_vapordome.Q = None
+    assert_allclose(pure_fluid_in_vapordome.Q, 0.6 * ctu.units.dimensionless)
+    assert_allclose(pure_fluid_in_vapordome.T, T)
+    assert_allclose(pure_fluid_in_vapordome.P, P)
+
+    with pytest.raises(ctu.CanteraError, match="an instance of a pint"):
+        pure_fluid_in_vapordome.Q = 0.5
+
+
+def yield_purefluid_only_setters():
+    props = [
+        pytest.param(("TPQ", "T", "P", "Q"), id="TPQ"),
+        pytest.param(("TQ", "T", "Q"), id="TQ"),
+        pytest.param(("PQ", "P", "Q"), id="PQ"),
+        pytest.param(("PV", "P", "v"), id="PV"),
+        pytest.param(("SH", "s", "h"), id="SH"),
+        pytest.param(("ST", "s", "T"), id="ST"),
+        pytest.param(("TH", "T", "h"), id="TH"),
+        pytest.param(("TV", "T", "v"), id="TV"),
+        pytest.param(("VH", "v", "h"), id="VH"),
+        pytest.param(("UP", "u", "P"), id="UP"),
+    ]
+    yield from props
+
+
+def yield_purefluid_only_getters():
+    props = [
+        pytest.param(("DPQ", "density", "P", "Q"), id="DPQ"),
+        pytest.param(("HPQ", "h", "P", "Q"), id="HPQ"),
+        pytest.param(("SPQ", "s", "P", "Q"), id="SPQ"),
+        pytest.param(("SVQ", "s", "v", "Q"), id="SVQ"),
+        pytest.param(("TDQ", "T", "density", "Q"), id="TDQ"),
+        pytest.param(("UVQ", "u", "v", "Q"), id="UVQ"),
+    ]
+    yield from props
+
+
+def yield_all_purefluid_only_props():
+    yield from yield_purefluid_only_setters()
+    yield from yield_purefluid_only_getters()
+
+
+@pytest.mark.parametrize("prop", yield_purefluid_only_setters())
+def test_set_without_units_is_error_purefluid(prop, pure_fluid):
+    value = [None] * (len(prop[0]) - 1) + [0.5]
+    with pytest.raises(ctu.CanteraError, match="an instance of a pint"):
+        setattr(pure_fluid, prop[0], value)
+
+
+@pytest.mark.parametrize("props", yield_all_purefluid_only_props())
+def test_multi_prop_getters_purefluid(pure_fluid, props):
+    # This test is equivalent to
+    # T, P, X = pure_fluid.TPX
+    # assert isclose(T, pure_fluid.T)
+    # assert isclose(P, pure_fluid.P)
+    # assert all(isclose(X, pure_fluid.X))
+    # Where T, P, and X loop through all the valid property pairs and triples.
+    pair_or_triple, first, second, *third = props
+    first_value, second_value, *third_value = getattr(pure_fluid, pair_or_triple)
+    assert_allclose(getattr(pure_fluid, first), first_value)
+    assert_allclose(getattr(pure_fluid, second), second_value)
+    if third:
+        assert_allclose(getattr(pure_fluid, third[0]), third_value[0])
+
+
+@pytest.mark.parametrize("props", yield_purefluid_only_setters())
+def test_setters_purefluid(props, pure_fluid):
+    # Only need to run this for a single pure fluid
+    initial_TD = pure_fluid.TD
+    pair_or_triple = props[0]
+
+    T_1 = ctu.Q_(500, "K")
+    if pair_or_triple in ("SH", "TH"):
+        # This state is able to converge for these setters, whereas the state below
+        # does not converge
+        rho_1 = ctu.Q_(1000, "kg/m**3")
+    else:
+        # This state is located inside the vapor dome to be able to test the
+        # TQ and PQ setters
+        rho_1 = ctu.Q_(25.93245092697775, "kg/m**3")
+
+    # Use TD setting to get the properties at the modified state
+    pure_fluid.TD = T_1, rho_1
+    new_props = getattr(pure_fluid, pair_or_triple)
+
+    # Reset to the initial state so that the next state setting actually has to do
+    # something.
+    pure_fluid.TD = initial_TD
+
+    # Use the test pair or triple to set the state and assert that the
+    # natural properties are equal to the modified state
+    setattr(pure_fluid, pair_or_triple, new_props)
+    assert_allclose(pure_fluid.T, T_1)
+    assert_allclose(pure_fluid.density, rho_1)
+
+
+@pytest.mark.parametrize("prop", ("X", "Y"))
+def test_X_Y_setters_with_none(generic_phase, prop):
+    comparison = getattr(generic_phase, prop)
+    setattr(generic_phase, prop, None)
+    # Assert that the value hasn't changed
+    assert_allclose(comparison, getattr(generic_phase, prop))
+
+
+@pytest.mark.parametrize("prop", ("X", "Y"))
+def test_X_Y_setters_without_units_works(generic_phase, prop):
+    composition = f"{generic_phase.species_names[0]}:1"
+    setattr(generic_phase, prop, composition)
+    assert_allclose(getattr(generic_phase, prop)[0], ctu.Q_([1], "dimensionless"))
+
+
+def test_thermophase_properties_exist(ideal_gas):
+    # Since the Solution class in the with_units subpackage only implements
+    # the ThermoPhase interface for now, instantiate a regular ThermoPhase
+    # to compare the attributes and make sure all of them exist on the with_units
+    # object
+    tp = ct.ThermoPhase("h2o2.yaml")
+
+    for attr in dir(tp):
+        if attr.startswith("_"):
+            continue
+
+        try:
+            getattr(tp, attr)
+        except (NotImplementedError, ct.ThermoModelMethodError):
+            continue
+
+        assert hasattr(ideal_gas, attr)
+
+
+def test_purefluid_properties_exist(pure_fluid):
+    # Test that all the properties implemented on the "upstream" PureFluid class
+    # are also implemented for the "with_units" variety.
+    pf = ct.PureFluid("liquidvapor.yaml", "water")
+    for attr in dir(pf):
+        if attr.startswith("_"):
+            continue
+
+        try:
+            getattr(pf, attr)
+        except (NotImplementedError, ct.ThermoModelMethodError):
+            continue
+
+        assert hasattr(pure_fluid, attr)