diff --git a/firedrake/adjoint_utils/constant.py b/firedrake/adjoint_utils/constant.py
index fed5978a6d..ae39df14f6 100644
--- a/firedrake/adjoint_utils/constant.py
+++ b/firedrake/adjoint_utils/constant.py
@@ -108,7 +108,7 @@ def _ad_copy(self):
         return self._constant_from_values()
 
     def _ad_dim(self):
-        return self.dat.cdim
+        return self.dat.data_ro.size
 
     def _ad_imul(self, other):
         self.assign(self._constant_from_values(self.dat.data_ro.reshape(-1) * other))
diff --git a/firedrake/assemble.py b/firedrake/assemble.py
index 1060abb52e..35889fe6af 100644
--- a/firedrake/assemble.py
+++ b/firedrake/assemble.py
@@ -1323,8 +1323,8 @@ def _get_mat_type(mat_type, sub_mat_type, arguments):
     if sub_mat_type is None:
         sub_mat_type = parameters.parameters["default_sub_matrix_type"]
 
-    if has_real_subspace and mat_type != "nest":
-        raise ValueError
+    if has_real_subspace and mat_type not in ["nest", "matfree"]:
+        raise ValueError("Matrices containing real space arguments must have type 'nest' or 'matfree'")
     if sub_mat_type not in {"aij", "baij"}:
         raise ValueError(
             f"Invalid submatrix type, '{sub_mat_type}' (not 'aij' or 'baij')"
diff --git a/firedrake/assign.py b/firedrake/assign.py
index d221b143cb..463ef771c5 100644
--- a/firedrake/assign.py
+++ b/firedrake/assign.py
@@ -5,6 +5,7 @@
 import finat.ufl
 import numpy as np
 import pyop3 as op3
+from pyop3.exceptions import DataValueError
 import pytools
 from pyadjoint.tape import annotate_tape
 from pyop2.utils import cached_property
@@ -256,8 +257,9 @@ def _assign_single_dat(self, lhs, subset, rvalue, assign_to_halos):
         if isinstance(rvalue, numbers.Number) or rvalue.shape in {(1,), assignee.shape}:
             assignee[...] = rvalue
         else:
-            cdim = self._assignee.function_space()._cdim
-            assert rvalue.shape == (cdim,)
+            cdim = self._assignee.function_space().value_size
+            if rvalue.shape != (cdim,):
+                raise DataValueError("Assignee and assignment values are different shapes")
             assignee.reshape((-1, cdim))[...] = rvalue
 
     def _compute_rvalue(self, func_data):
diff --git a/firedrake/constant.py b/firedrake/constant.py
index 244fa8d0db..15e2ada542 100644
--- a/firedrake/constant.py
+++ b/firedrake/constant.py
@@ -3,8 +3,8 @@
 import finat.ufl
 
 from tsfc.ufl_utils import TSFCConstantMixin
-from pyop2.exceptions import DataTypeError, DataValueError
 import pyop3 as op3
+from pyop3.exceptions import DataValueError
 from firedrake.petsc import PETSc
 from firedrake.utils import ScalarType
 from ufl.classes import all_ufl_classes, ufl_classes, terminal_classes
@@ -29,15 +29,12 @@ def _create_const(value, comm):
     shape = data.shape
     rank = len(shape)
 
-    if comm is not None:
-        raise NotImplementedError("Won't be a back door for real space here, do elsewhere")
-
     if rank == 0:
-        axes = op3.AxisTree(op3.Axis(1))
+        axes = op3.AxisTree()
     else:
-        axes = op3.AxisTree(op3.Axis(shape[0]))
-        for size in shape[1:]:
-            axes = axes.add_axis(op3.Axis(size), *axes.leaf)
+        axes = op3.AxisTree(op3.Axis({"XXX": shape[0]}, label="dim0"))
+        for i, s in enumerate(shape[1:]):
+            axes = axes.add_axis(op3.Axis({"XXX": s}, label=f"dim{i+1}"), *axes.leaf)
     dat = op3.HierarchicalArray(axes, data=data.flatten())
     return dat, rank, shape
 
@@ -198,11 +195,10 @@ def assign(self, value):
         self
 
         """
+        if self.ufl_shape() and np.array(value).shape != self.ufl_shape():
+            raise DataValueError("Cannot assign to constant, value has incorrect shape")
         self.dat.data_wo[...] = value
         return self
-        # TODO pyop3
-        # except (DataTypeError, DataValueError) as e:
-        #     raise ValueError(e)
 
     def __iadd__(self, o):
         raise NotImplementedError("Augmented assignment to Constant not implemented")
diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py
index d2faa0ea08..30c2107cdd 100644
--- a/firedrake/functionspaceimpl.py
+++ b/firedrake/functionspaceimpl.py
@@ -675,6 +675,10 @@ def _local_ises(self):
     @utils.cached_property
     def local_section(self):
         section = PETSc.Section().create(comm=self.comm)
+        if self._ufl_function_space.ufl_element().family() == "Real":
+            # If real we don't need to populate the section
+            return section
+
         points = self._mesh.points
         section.setChart(0, points.size)
         perm = PETSc.IS().createGeneral(points.numbering.data_ro, comm=self.comm)
diff --git a/firedrake/logging.py b/firedrake/logging.py
index 7b524fbe9c..8c279ec496 100644
--- a/firedrake/logging.py
+++ b/firedrake/logging.py
@@ -15,7 +15,7 @@
            "RED", "GREEN", "BLUE")
 
 
-packages = ("pyop2", "tsfc", "firedrake", "UFL")
+packages = ("pyop2", "pyop3", "tsfc", "firedrake", "UFL")
 
 
 logger = logging.getLogger("firedrake")
diff --git a/firedrake/parloops.py b/firedrake/parloops.py
index b5ab7540d9..f33950b99b 100644
--- a/firedrake/parloops.py
+++ b/firedrake/parloops.py
@@ -377,6 +377,9 @@ def _(
 ):
     plex = V.mesh().topology
 
+    if V.ufl_element().family() == "Real":
+        return array
+
     if integral_type == "cell":
         # TODO ideally the FIAT permutation would not need to be known
         # about by the mesh topology and instead be handled here. This
@@ -518,30 +521,24 @@ def _(
 def _cell_integral_pack_indices(V: WithGeometry, cell: op3.LoopIndex) -> op3.IndexTree:
     plex = V.mesh().topology
 
-    if V.ufl_element().family() == "Real":
-        indices = op3.IndexTree(op3.Slice("dof", [op3.AffineSliceComponent("XXX")]))
-    else:
-        indices = op3.IndexTree.from_nest({
-            plex._fiat_closure(cell): [
-                op3.Slice("dof", [op3.AffineSliceComponent("XXX")])
-                for _ in range(plex.dimension+1)
-            ]
-        })
+    indices = op3.IndexTree.from_nest({
+        plex._fiat_closure(cell): [
+            op3.Slice("dof", [op3.AffineSliceComponent("XXX")])
+            for _ in range(plex.dimension+1)
+        ]
+    })
     return _with_shape_indices(V, indices)
 
 
 def _facet_integral_pack_indices(V: WithGeometry, facet: op3.LoopIndex) -> op3.IndexTree:
     plex = V.ufl_domain().topology
 
-    if V.ufl_element().family() == "Real":
-        indices = op3.IndexTree(op3.ScalarIndex(plex.name, "XXX", 0))
-    else:
-        indices = op3.IndexTree.from_nest({
-            plex._fiat_closure(plex.support(facet)): [
-                op3.Slice("dof", [op3.AffineSliceComponent("XXX")])
-                for _ in range(plex.dimension+1)
-            ]
-        })
+    indices = op3.IndexTree.from_nest({
+        plex._fiat_closure(plex.support(facet)): [
+            op3.Slice("dof", [op3.AffineSliceComponent("XXX")])
+            for _ in range(plex.dimension+1)
+        ]
+    })
     # don't add support as an extra axis here, done already
     return _with_shape_indices(V, indices, and_support=False)
 
@@ -624,13 +621,19 @@ def _with_shape_axes(V, axes, target_paths, index_exprs, integral_type):
     trees_ = []
     for space, tree in zip(spaces, trees):
         if space.shape:
-            for leaf in tree.leaves:
-                for i, dim in enumerate(space.shape):
-                    label = f"dim{i}"
-                    subaxis = op3.Axis({"XXX": dim}, label)
-                    tree = tree.add_axis(subaxis, *leaf)
-                    new_target_paths[subaxis.id, "XXX"] = pmap({label: "XXX"})
-                    new_index_exprs[subaxis.id, "XXX"] = pmap({label: op3.AxisVariable(label)})
+            for parent, component in tree.leaves:
+                axis_list = [
+                    op3.Axis({"XXX": dim}, f"dim{ii}")
+                    for ii, dim in enumerate(space.shape)
+                ]
+                tree = tree.add_subtree(
+                    op3.AxisTree.from_iterable(axis_list),
+                    parent=parent,
+                    component=component
+                )
+                for axis in axis_list:
+                    new_target_paths[axis.id, "XXX"] = pmap({axis.label: "XXX"})
+                    new_index_exprs[axis.id, "XXX"] = pmap({axis.label: op3.AxisVariable(axis.label)})
 
         trees_.append(tree)
     trees = tuple(trees_)
diff --git a/tests/regression/test_appctx_cleanup.py b/tests/regression/test_appctx_cleanup.py
index 39c33aaa56..7021dcd1fc 100644
--- a/tests/regression/test_appctx_cleanup.py
+++ b/tests/regression/test_appctx_cleanup.py
@@ -1,5 +1,9 @@
 import numpy
 from firedrake import *
+import pytest
+
+
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
 
 
 class NonePC(PCBase):
diff --git a/tests/regression/test_bcs.py b/tests/regression/test_bcs.py
index dd9934491f..08c40391ad 100644
--- a/tests/regression/test_bcs.py
+++ b/tests/regression/test_bcs.py
@@ -50,7 +50,6 @@ def test_init_bcs_illegal(mesh, v):
         DirichletBC(FunctionSpace(mesh, "CG", 1), v, 0)
 
 
-@pytest.mark.skip(reason="pyop3 TODO")
 @pytest.mark.parametrize('measure', [dx, ds])
 def test_assemble_bcs_wrong_fs(V, measure):
     "Assemble a Matrix with a DirichletBC on an incompatible FunctionSpace."
@@ -241,7 +240,6 @@ def test_preassembly_doesnt_modify_assembled_rhs(V, f):
     assert np.allclose(b_vals, b.dat.data_ro)
 
 
-@pytest.mark.skip(reason="pyop3 TODO")
 def test_preassembly_bcs_caching(V):
     bc1 = DirichletBC(V, 0, 1)
     bc2 = DirichletBC(V, 1, 2)
@@ -268,6 +266,7 @@ def test_preassembly_bcs_caching(V):
     assert not any(Aneither.M.values.diagonal() == 0)
 
 
+@pytest.mark.skip(reason="pyop3 TODO")
 def test_assemble_mass_bcs_2d(V):
     if V.value_size > 1:
         pytest.skip(reason="pyop3 TODO")
@@ -295,6 +294,7 @@ def test_assemble_mass_bcs_2d(V):
     assert assemble(inner((w - f), (w - f))*dx) < 1e-12
 
 
+@pytest.mark.skip(reason="pyop3 TODO")
 @pytest.mark.parametrize("quad",
                          [False, True],
                          ids=["triangle", "quad"])
diff --git a/tests/regression/test_constant.py b/tests/regression/test_constant.py
index 99c349b066..04e59a2d97 100644
--- a/tests/regression/test_constant.py
+++ b/tests/regression/test_constant.py
@@ -5,9 +5,6 @@
 import pytest
 
 
-pytest.skip(allow_module_level=True, reason="pyop3 TODO")
-
-
 def test_scalar_constant():
     for m in [UnitIntervalMesh(5), UnitSquareMesh(2, 2), UnitCubeMesh(2, 2, 2)]:
         c = Constant(1, domain=m)
@@ -122,8 +119,8 @@ def test_constant_vector_assign_works():
 
     f.assign(c)
 
-    assert np.allclose(f.dat.data_ro[:, 0], 10)
-    assert np.allclose(f.dat.data_ro[:, 1], 11)
+    assert np.allclose(f.sub(0).dat.data_ro, 10)
+    assert np.allclose(f.sub(1).dat.data_ro, 11)
 
 
 def test_constant_vector_assign_to_scalar_error():
@@ -162,9 +159,11 @@ def test_constant_assign_to_mixed():
     f.sub(0).assign(c)
     f.sub(1).assign(c)
 
-    for d in f.dat.data_ro:
-        assert np.allclose(d[:, 0], 10)
-        assert np.allclose(d[:, 1], 11)
+
+    assert np.allclose(f.sub(0).sub(0).dat.data_ro, 10)
+    assert np.allclose(f.sub(0).sub(1).dat.data_ro, 11)
+    assert np.allclose(f.sub(1).sub(0).dat.data_ro, 10)
+    assert np.allclose(f.sub(0).sub(1).dat.data_ro, 11)
 
 
 def test_constant_multiplies_function():
@@ -207,6 +206,7 @@ def test_constant_names_are_not_used_in_generated_code():
 
 
 @pytest.mark.skipcomplex
+@pytest.mark.xfail(reason="requires matnest")
 def test_correct_constants_are_used_in_split_form():
     # see https://github.com/firedrakeproject/firedrake/issues/3091
     mesh = UnitSquareMesh(3, 3)
diff --git a/tests/regression/test_custom_callbacks.py b/tests/regression/test_custom_callbacks.py
index 26c164c96b..3b6c635aa9 100644
--- a/tests/regression/test_custom_callbacks.py
+++ b/tests/regression/test_custom_callbacks.py
@@ -1,6 +1,10 @@
 from firedrake import *
 from firedrake.utils import ScalarType
 import numpy as np
+import pytest
+
+
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
 
 
 def test_callbacks():
diff --git a/tests/regression/test_facets.py b/tests/regression/test_facets.py
index a01457c255..a1d18cf99b 100644
--- a/tests/regression/test_facets.py
+++ b/tests/regression/test_facets.py
@@ -183,6 +183,7 @@ def test_internal_integral_unit_tet():
     assert abs(assemble(u('+') * dS)) < 1.0e-14
 
 
+@pytest.mark.xfail(reason="pyop3 TODO")
 def test_facet_map_no_reshape():
     m = UnitSquareMesh(1, 1)
     V = FunctionSpace(m, "DG", 0)
@@ -190,6 +191,7 @@ def test_facet_map_no_reshape():
     assert efnm.values_with_halo.shape == (4, 1)
 
 
+@pytest.mark.skip(reason="pyop3 TODO")
 def test_mesh_with_no_facet_markers():
     mesh = UnitTriangleMesh()
     mesh.init()
diff --git a/tests/regression/test_interior_facets.py b/tests/regression/test_interior_facets.py
index 257109d778..b39cf0a87c 100644
--- a/tests/regression/test_interior_facets.py
+++ b/tests/regression/test_interior_facets.py
@@ -7,6 +7,9 @@
 PETSc.Sys.popErrorHandler()
 
 
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
+
+
 def run_test():
     # mesh = UnitSquareMesh(10, 10)
     mesh = UnitSquareMesh(2, 2)
diff --git a/tests/regression/test_interp_dual.py b/tests/regression/test_interp_dual.py
index 5928d81e4e..46e6bd5644 100644
--- a/tests/regression/test_interp_dual.py
+++ b/tests/regression/test_interp_dual.py
@@ -5,6 +5,9 @@
 import ufl
 
 
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
+
+
 @pytest.fixture(scope='module')
 def mesh():
     return UnitSquareMesh(5, 5)
diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py
index c20f5d34b9..5e2cfab84c 100644
--- a/tests/regression/test_interpolate.py
+++ b/tests/regression/test_interpolate.py
@@ -7,6 +7,9 @@
 cwd = abspath(dirname(__file__))
 
 
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
+
+
 def test_constant():
     cg1 = FunctionSpace(UnitSquareMesh(5, 5), "CG", 1)
     f = assemble(interpolate(Constant(1.0), cg1))
diff --git a/tests/regression/test_multiple_domains.py b/tests/regression/test_multiple_domains.py
index bb9391f3e9..1a0aa60dd8 100644
--- a/tests/regression/test_multiple_domains.py
+++ b/tests/regression/test_multiple_domains.py
@@ -48,7 +48,6 @@ def test_mismatching_meshes_indexed_function(mesh1, mesh3):
         assemble(inner(d1, TestFunction(V2))*dx(domain=mesh1))
 
 
-@pytest.mark.skip(reason="pyop3 TODO")
 def test_mismatching_meshes_constant(mesh1, mesh3):
     V2 = FunctionSpace(mesh3, "CG", 1)
 
diff --git a/tests/regression/test_point_eval_api.py b/tests/regression/test_point_eval_api.py
index b129472dec..ab0ffd603d 100644
--- a/tests/regression/test_point_eval_api.py
+++ b/tests/regression/test_point_eval_api.py
@@ -7,7 +7,7 @@
 cwd = abspath(dirname(__file__))
 
 
-pytest.mark.skip(allow_module_level=True, reason="pyop3 point location")
+pytest.skip(allow_module_level=True, reason="pyop3 point location")
 
 
 def test_1d_args():
diff --git a/tests/regression/test_point_eval_cells.py b/tests/regression/test_point_eval_cells.py
index db7c601ada..5fa5f4a5d3 100644
--- a/tests/regression/test_point_eval_cells.py
+++ b/tests/regression/test_point_eval_cells.py
@@ -7,7 +7,7 @@
 cwd = abspath(dirname(__file__))
 
 
-pytest.mark.skip(allow_module_level=True, reason="pyop3 point location")
+pytest.skip(allow_module_level=True, reason="pyop3 point location")
 
 
 @pytest.fixture(params=[False, True])
diff --git a/tests/regression/test_point_eval_fs.py b/tests/regression/test_point_eval_fs.py
index 37cc66c255..c594bc1b34 100644
--- a/tests/regression/test_point_eval_fs.py
+++ b/tests/regression/test_point_eval_fs.py
@@ -7,7 +7,7 @@
 cwd = abspath(dirname(__file__))
 
 
-pytest.mark.skip(allow_module_level=True, reason="pyop3 point location")
+pytest.skip(allow_module_level=True, reason="pyop3 point location")
 
 
 @pytest.fixture
diff --git a/tests/regression/test_real_space.py b/tests/regression/test_real_space.py
index f92cf0dbc9..93f1d3b6c3 100644
--- a/tests/regression/test_real_space.py
+++ b/tests/regression/test_real_space.py
@@ -5,6 +5,9 @@
 from firedrake.__future__ import *
 
 
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
+
+
 @pytest.mark.skipcomplex
 def test_real_assembly():
     mesh = UnitIntervalMesh(3)
diff --git a/tests/regression/test_serendipity_biharmonic.py b/tests/regression/test_serendipity_biharmonic.py
index d66d9b0ef6..701c561e68 100644
--- a/tests/regression/test_serendipity_biharmonic.py
+++ b/tests/regression/test_serendipity_biharmonic.py
@@ -1,5 +1,9 @@
 from firedrake import *
 import numpy
+import pytest
+
+
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
 
 
 def test_serendipity_biharmonic():
diff --git a/tests/regression/test_vfs_component_bcs.py b/tests/regression/test_vfs_component_bcs.py
index c670534cf3..4e16394cbd 100644
--- a/tests/regression/test_vfs_component_bcs.py
+++ b/tests/regression/test_vfs_component_bcs.py
@@ -3,6 +3,9 @@
 import numpy as np
 
 
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
+
+
 @pytest.fixture
 def m():
     return UnitSquareMesh(4, 4)
diff --git a/tests/regression/test_zero_forms.py b/tests/regression/test_zero_forms.py
index 0fe3878948..d65a2be067 100644
--- a/tests/regression/test_zero_forms.py
+++ b/tests/regression/test_zero_forms.py
@@ -4,9 +4,6 @@
 from firedrake import *
 
 
-pytest.skip(allow_module_level=True, reason="pyop3 TODO")
-
-
 @pytest.fixture(scope='module', params=[False, True])
 def mesh(request):
     quadrilateral = request.param
@@ -20,12 +17,10 @@ def mesh(request):
            (1, 2, 3, 4)]
 
 
-@pytest.mark.skip(reason="pyop3 TODO")
 def test_ds_dx(mesh):
     assert np.allclose(assemble(1*dx(domain=mesh) + 1*ds(domain=mesh)), 5.0)
 
 
-@pytest.mark.skip(reason="pyop3 TODO")
 @pytest.mark.parametrize('domains', domains)
 def test_dsn(mesh, domains):
 
@@ -38,7 +33,7 @@ def test_dsn(mesh, domains):
     assert np.allclose(assemble(form), len(domains))
 
 
-@pytest.mark.skip(reason="pyop3 TODO")
+@pytest.mark.skip(reason="pyop3 parallel")
 @pytest.mark.parallel
 def test_dsn_parallel(mesh):
     for d in domains:
@@ -67,9 +62,6 @@ def test_dsn_parallel(mesh):
                                            ['function', 'constant'],
                                            ['scalar', 'vector', 'tensor']))
 def test_math_functions(mesh, expr, value, typ, fs_type):
-    if mesh.ufl_cell().cellname == "quadrilateral":
-        pytest.skip("pyop3 TODO")
-
     if typ == 'function':
         family, degree = 'CG', 1
     elif typ == 'constant':
diff --git a/tests/slate/test_local_logging.py b/tests/slate/test_local_logging.py
index ca48fe7424..d279507ff4 100644
--- a/tests/slate/test_local_logging.py
+++ b/tests/slate/test_local_logging.py
@@ -1,4 +1,8 @@
 import os
+import pytest
+
+
+pytest.skip(allow_module_level=True, reason="pyop3 TODO")
 
 
 def test_slate_logging():
diff --git a/tests/slate/test_slate_hybridization_extr.py b/tests/slate/test_slate_hybridization_extr.py
index f66bc5a7b6..2b8d1a631e 100644
--- a/tests/slate/test_slate_hybridization_extr.py
+++ b/tests/slate/test_slate_hybridization_extr.py
@@ -65,4 +65,4 @@ def test_hybrid_extr_helmholtz(quad):
     u_err = errornorm(u_h, nh_u)
 
     assert sigma_err < 5e-8
-    assert u_err < 1e-
+    assert u_err < 1e-8