diff --git a/scripts/test_levels.sh b/scripts/test_levels.sh
index bb376884..0207c6f5 100755
--- a/scripts/test_levels.sh
+++ b/scripts/test_levels.sh
@@ -146,26 +146,34 @@ run_tests() {
 }
 
 # Functions for each test level using pytest markers
+#
+# A level selection EXCLUDES the levels above it, and has to: pytest MERGES
+# marks rather than overriding them, so a file whose module declares
+# `pytestmark = pytest.mark.level_1` and whose heavy test then carries
+# `@pytest.mark.level_2` leaves that test marked BOTH. A plain `-m level_1`
+# selects it, and the demotion the author wrote does nothing. Nine files rely
+# on that demotion; the heaviest of their tests is a 96-second homotopy solve
+# that was running in the "quick" tier because of it.
 run_level_1() {
     echo "⚡ Running LEVEL 1: Quick Tests (Core Functionality)"
-    echo "Using pytest marker: -m level_1"
+    echo "Using pytest marker: -m 'level_1 and not level_2 and not level_3'"
     echo "Expected runtime: ~2 minutes"
     echo ""
 
-    # Run all tests marked with level_1
+    # Tests marked level_1 and NOT demoted to a higher level (see above)
     run_tests "Level 1 tests (quick core functionality)" \
-        tests/ -m level_1
+        tests/ -m "level_1 and not level_2 and not level_3"
 }
 
 run_level_2() {
     echo "🔧 Running LEVEL 2: Intermediate Tests"
-    echo "Using pytest marker: -m level_2"
+    echo "Using pytest marker: -m 'level_2 and not level_3'"
     echo "Expected runtime: ~5 minutes"
     echo ""
 
-    # Run all tests marked with level_2
+    # Tests marked level_2 and NOT demoted to level_3 (see above)
     run_tests "Level 2 tests (units, integration, projections)" \
-        tests/ -m level_2
+        tests/ -m "level_2 and not level_3"
 
     # Parallel tests for global statistics (requires MPI)
     if [ $RUN_PARALLEL -eq 1 ]; then
diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index 03fe527f..fa3f4b4d 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -143,9 +143,20 @@ def _from_gmsh(filename, comm=None, markVertices=False, useRegions=True, useMult
             plex_0.setName("uw_mesh")
             plex_0.markBoundaryFaces("All_Boundaries", 1001)
 
-            viewer = PETSc.ViewerHDF5().create(filename + ".h5", "w", comm=PETSc.COMM_SELF)
+            # Write aside and rename, so ``filename + ".h5"`` never names a
+            # half-written file. The name comes from the mesh PARAMETERS, so a
+            # second process building the same geometry in this directory picks
+            # the same one and would otherwise read what this one is still
+            # writing (issue #563).
+            # Imported here, not at module scope: the meshing package imports
+            # this module, so a top-level import would close a cycle.
+            from underworld3.meshing._mesh_files import _scratch_name
+
+            scratch = _scratch_name(filename + ".h5")
+            viewer = PETSc.ViewerHDF5().create(str(scratch), "w", comm=PETSc.COMM_SELF)
             viewer(plex_0)
             viewer.destroy()
+            os.replace(scratch, filename + ".h5")
     finally:
         # The gmsh import options are import-time scratch — meaningful only for
         # the createFromFile above. Clear the whole namespace so a value set by
@@ -155,7 +166,10 @@ def _from_gmsh(filename, comm=None, markVertices=False, useRegions=True, useMult
         # read as 2-D). Runs on success or failure.
         _clear_gmsh_import_options()
 
-    # Now we have an h5 file and we can hand this to _from_plexh5
+    # Now we have an h5 file and we can hand this to _from_plexh5. The barrier
+    # is what the atomic write above makes necessary AND sufficient: the other
+    # ranks must not look for the file before rank 0 has renamed it into place.
+    uw.mpi.barrier()
 
     return _from_plexh5(filename + ".h5", comm, return_sf=True)
 
diff --git a/src/underworld3/meshing/_mesh_files.py b/src/underworld3/meshing/_mesh_files.py
new file mode 100644
index 00000000..30676741
--- /dev/null
+++ b/src/underworld3/meshing/_mesh_files.py
@@ -0,0 +1,77 @@
+r"""Where generated mesh files are written, and how they are written.
+
+Building a mesh hands gmsh a file to write and then reads it back through
+PETSc, so the pair — ``<name>.msh`` and the ``<name>.msh.h5`` PETSc converts
+it to — is scratch shared between those two steps. It is not a cache: every
+construction regenerates both, and nothing checks whether they already exist.
+
+The name is derived from the mesh PARAMETERS, so two processes building the
+same geometry in one working directory choose the same name, and one can read
+a file the other is still writing (issue #563). Identical geometry is exactly
+the colliding case, which is why a parameter sweep or a parallel test run hits
+it and ordinary use does not.
+
+Two mechanisms make that safe, and both are needed:
+
+* the directory is settable per process through ``UW_MESH_CACHE_DIR``, so
+  independent jobs can be given somewhere of their own;
+* every write lands atomically, so a reader that shares a directory anyway
+  sees a complete file or no file, never a half-written one.
+"""
+import os
+from pathlib import Path
+
+import underworld3 as uw
+
+DEFAULT_MESH_FILE_DIR = ".meshes"
+
+
+def mesh_file_dir():
+    """The directory generated mesh files are written to.
+
+    ``UW_MESH_CACHE_DIR`` overrides the default ``.meshes``. Every rank of one
+    job must agree on it, so this reads the environment — inherited identically
+    by every rank — and never anything process-local such as the pid.
+    """
+    return Path(os.environ.get("UW_MESH_CACHE_DIR", DEFAULT_MESH_FILE_DIR))
+
+
+def mesh_file_path(basename):
+    """Full path for a generated mesh file, with its directory created.
+
+    Parameters
+    ----------
+    basename : str
+        The file's name, conventionally ``uw_<generator>_<parameters>.msh``.
+    """
+    directory = mesh_file_dir()
+    if uw.mpi.rank == 0:
+        directory.mkdir(parents=True, exist_ok=True)
+    return str(directory / basename)
+
+
+def _scratch_name(final):
+    """A process-unique sibling of ``final`` KEEPING ITS EXTENSION.
+
+    The extension has to survive: gmsh chooses its output format from it, so
+    writing to ``mesh.msh.1234.tmp`` would silently produce something that is
+    not a gmsh mesh.
+    """
+    final = Path(final)
+    return final.with_name(f"{final.stem}.{os.getpid()}.tmp{final.suffix}")
+
+
+def write_gmsh(filename):
+    """``gmsh.write``, landing atomically at ``filename``.
+
+    gmsh writes in place, so a concurrent reader can open a file that is still
+    being filled. Writing under a process-unique name and renaming makes the
+    appearance of the final name atomic — :func:`os.replace` is atomic within a
+    filesystem — so a reader sees either the previous complete file or the new
+    one.
+    """
+    import gmsh
+
+    scratch = _scratch_name(filename)
+    gmsh.write(str(scratch))
+    os.replace(scratch, filename)
diff --git a/src/underworld3/meshing/annulus.py b/src/underworld3/meshing/annulus.py
index 01fa3cbc..bbc6b41c 100644
--- a/src/underworld3/meshing/annulus.py
+++ b/src/underworld3/meshing/annulus.py
@@ -15,6 +15,7 @@
 import math
 
 import underworld3 as uw
+from underworld3.meshing._mesh_files import mesh_file_dir, write_gmsh
 from underworld3.discretisation import Mesh
 from underworld3 import VarType
 from underworld3.coordinates import CoordinateSystemType
@@ -129,9 +130,9 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f"uw_QuarterAnnulus_ro{radiusOuter}_ri{radiusInner}_csize{cellSize}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_QuarterAnnulus_ro{radiusOuter}_ri{radiusInner}_csize{cellSize}.msh"
     else:
         uw_filename = filename
 
@@ -237,7 +238,7 @@ class boundaries(Enum):
 
         print("generate")
 
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     new_mesh = Mesh(
@@ -397,9 +398,9 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f".meshes/uw_annulus_ro{radiusOuter}_ri{radiusInner}_csize{cellSize}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_annulus_ro{radiusOuter}_ri{radiusInner}_csize{cellSize}.msh"
     else:
         uw_filename = filename
 
@@ -466,7 +467,7 @@ class boundaries(Enum):
         gmsh.model.geo.synchronize()
 
         gmsh.model.mesh.generate(2)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     # Ensure boundaries conform (if refined)
@@ -625,9 +626,9 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f"uw_SegmentOfAnnulus_ro{radiusOuter}_ri{radiusInner}_extent{angleExtent}_csize{cellSize}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_SegmentOfAnnulus_ro{radiusOuter}_ri{radiusInner}_extent{angleExtent}_csize{cellSize}.msh"
     else:
         uw_filename = filename
 
@@ -720,7 +721,7 @@ class boundaries(Enum):
         gmsh.model.geo.synchronize()
 
         gmsh.model.mesh.generate(2)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     # Ensure boundaries conform (if refined)
@@ -896,9 +897,9 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f".meshes/uw_annulus_ro{radiusOuter}_ri{radiusInner}_csize{cellSizeOuter}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_annulus_ro{radiusOuter}_ri{radiusInner}_csize{cellSizeOuter}.msh"
     else:
         uw_filename = filename
 
@@ -1004,7 +1005,7 @@ class boundaries(Enum):
 
         gmsh.model.mesh.generate(2)
 
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
 
         # We need to build the plex here in order to make some changes
         # before the mesh gets built
@@ -1249,9 +1250,9 @@ class regions(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f".meshes/uw_annulus_internalBoundary_rO{radiusOuter}rInt{radiusInternal}_rI{radiusInner}_csize{cellSize}_csizefs{cellSize_Outer}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_annulus_internalBoundary_rO{radiusOuter}rInt{radiusInternal}_rI{radiusInner}_csize{cellSize}_csizefs{cellSize_Outer}.msh"
     else:
         uw_filename = filename
 
@@ -1345,7 +1346,7 @@ class regions(Enum):
         gmsh.model.geo.synchronize()
 
         gmsh.model.mesh.generate(2)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     ## This is the same as the simple annulus
@@ -1558,9 +1559,9 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f".meshes/uw_disc_internalBoundaries_rO{radiusUpper}rInt{radiusInternal}_rI{radiusLower}_csize{cellSize}_csizefs{cellSize_Upper}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_disc_internalBoundaries_rO{radiusUpper}rInt{radiusInternal}_rI{radiusLower}_csize{cellSize}_csizefs{cellSize_Upper}.msh"
     else:
         uw_filename = filename
 
@@ -1643,7 +1644,7 @@ class boundaries(Enum):
         gmsh.model.geo.synchronize()
 
         gmsh.model.mesh.generate(2)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     ## This is the same as the simple annulus
diff --git a/src/underworld3/meshing/cartesian.py b/src/underworld3/meshing/cartesian.py
index 20579204..3e4d3ca2 100644
--- a/src/underworld3/meshing/cartesian.py
+++ b/src/underworld3/meshing/cartesian.py
@@ -16,6 +16,7 @@
 import math
 
 import underworld3 as uw
+from underworld3.meshing._mesh_files import mesh_file_dir, write_gmsh
 from underworld3.discretisation import Mesh
 from underworld3 import VarType
 from underworld3.coordinates import CoordinateSystemType
@@ -73,7 +74,8 @@ def UnstructuredSimplexBox(
         Currently only works for 2D meshes.
     filename : str, optional
         Path to save the mesh file. If None, generates a unique name
-        in the ``.meshes/`` directory based on mesh parameters.
+        in the mesh-file directory (``.meshes/`` by default, or
+        ``UW_MESH_CACHE_DIR``) based on mesh parameters.
     refinement : int, optional
         Number of uniform refinement levels to apply after mesh
         generation. Each level approximately quadruples element count.
@@ -197,9 +199,9 @@ class boundary_normals_3D(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f".meshes/uw_simplexbox_minC{minCoords}_maxC{maxCoords}_csize{cellSize}_reg{regular}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_simplexbox_minC{minCoords}_maxC{maxCoords}_csize{cellSize}_reg{regular}.msh"
     else:
         uw_filename = filename
 
@@ -309,7 +311,7 @@ class boundary_normals_3D(Enum):
 
         # Generate Mesh
         gmsh.model.mesh.generate(dim)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     def box_return_coords_to_bounds(coords):
@@ -413,7 +415,8 @@ def BoxInternalBoundary(
     qdegree : int, default=2
         Quadrature degree for numerical integration.
     filename : str, optional
-        Path to save the mesh file. If None, auto-generates in ``.meshes/``.
+        Path to save the mesh file. If None, auto-generates in the mesh-file
+        directory (``.meshes/`` by default, or ``UW_MESH_CACHE_DIR``).
     refinement : int, optional
         Number of uniform refinement levels to apply.
     gmsh_verbosity : int, default=0
@@ -540,12 +543,12 @@ class boundary_normals_3D(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
         if not simplex:
             # structuredQuadBoxIB
-            uw_filename = f".meshes/uw_sqbIB_minC{minCoords}_maxC{maxCoords}.msh"
+            uw_filename = f"{mesh_file_dir()}/uw_sqbIB_minC{minCoords}_maxC{maxCoords}.msh"
         else:
-            uw_filename = f".meshes/uw_usbIB_minC{minCoords}_maxC{maxCoords}.msh"
+            uw_filename = f"{mesh_file_dir()}/uw_usbIB_minC{minCoords}_maxC{maxCoords}.msh"
     else:
         uw_filename = filename
 
@@ -646,7 +649,7 @@ class boundary_normals_3D(Enum):
                 gmsh.model.mesh.set_recombine(2, surface2)
 
             gmsh.model.mesh.generate(dim)
-            gmsh.write(uw_filename)
+            write_gmsh(uw_filename)
             gmsh.finalize()
 
         if dim == 3:
@@ -883,7 +886,7 @@ class boundary_normals_3D(Enum):
                 gmsh.model.mesh.set_recombine(3, volume_b)
 
             gmsh.model.mesh.generate(dim)
-            gmsh.write(uw_filename)
+            write_gmsh(uw_filename)
             gmsh.finalize()
 
     def box_return_coords_to_bounds(coords):
@@ -1132,11 +1135,11 @@ def _one_patch(entry, default_name):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
         grade = ("" if patch_cellSize is None
                  else f"_pcs{patch_cellSize}_gd{grading_distance}")
         tagline = "_".join(f"{nm}{len(p)}" for nm, p, _f in patches)
-        uw_filename = (f".meshes/uw_boxpatch_minC{minCoords}_"
+        uw_filename = (f"{mesh_file_dir()}/uw_boxpatch_minC{minCoords}_"
                        f"maxC{maxCoords}_csize{cellSize}{grade}_"
                        f"{tagline}.msh")
     else:
@@ -1224,7 +1227,7 @@ def _one_patch(entry, default_name):
         # message. Fault-session follow-up.
         try:
             gmsh.model.mesh.generate(3)
-            gmsh.write(uw_filename)
+            write_gmsh(uw_filename)
         finally:
             # A mesher failure must not leave the gmsh session
             # initialized: a poisoned session makes the NEXT mesh
@@ -1298,7 +1301,8 @@ def StructuredQuadBox(
     qdegree : int, default=2
         Quadrature degree for numerical integration.
     filename : str, optional
-        Path to save the mesh file. If None, auto-generates in ``.meshes/``.
+        Path to save the mesh file. If None, auto-generates in the mesh-file
+        directory (``.meshes/`` by default, or ``UW_MESH_CACHE_DIR``).
     refinement : int, optional
         Number of uniform refinement levels to apply.
     gmsh_verbosity : int, default=0
@@ -1421,9 +1425,9 @@ class boundary_normals_3D(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f".meshes/uw_structuredQuadBox_minC{minCoords}_maxC{maxCoords}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_structuredQuadBox_minC{minCoords}_maxC{maxCoords}.msh"
     else:
         uw_filename = filename
 
@@ -1616,7 +1620,7 @@ class boundary_normals_3D(Enum):
 
         # Generate Mesh
         gmsh.model.mesh.generate(dim)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     def box_return_coords_to_bounds(coords):
diff --git a/src/underworld3/meshing/geographic.py b/src/underworld3/meshing/geographic.py
index 29a8ddbd..6a2aafba 100644
--- a/src/underworld3/meshing/geographic.py
+++ b/src/underworld3/meshing/geographic.py
@@ -15,6 +15,7 @@
 import math
 
 import underworld3 as uw
+from underworld3.meshing._mesh_files import mesh_file_dir, write_gmsh
 from underworld3.discretisation import Mesh
 from underworld3 import VarType
 from underworld3.coordinates import CoordinateSystemType
@@ -206,8 +207,8 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
-        uw_filename = f".meshes/uw_cubed_spherical_shell_ro{radiusOuter}_ri{radiusInner}_elts{numElementsDepth}_plex{simplex}.msh"
+            os.makedirs(mesh_file_dir(), exist_ok=True)
+        uw_filename = f"{mesh_file_dir()}/uw_cubed_spherical_shell_ro{radiusOuter}_ri{radiusInner}_elts{numElementsDepth}_plex{simplex}.msh"
     else:
         uw_filename = filename
 
@@ -326,7 +327,7 @@ class boundaries(Enum):
 
         # Generate Mesh
         gmsh.model.mesh.generate(3)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     def spherical_mesh_refinement_callback(dm):
@@ -618,9 +619,9 @@ class boundaries(Enum):
     # Generate mesh filename if not provided
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
         uw_filename = (
-            f".meshes/uw_geographic_{ellipsoid_dict['planet']}_"
+            f"{mesh_file_dir()}/uw_geographic_{ellipsoid_dict['planet']}_"
             f"lon{lon_min:.1f}_{lon_max:.1f}_"
             f"lat{lat_min:.1f}_{lat_max:.1f}_"
             f"d{depth_min:.0f}_{depth_max:.0f}_"
@@ -764,7 +765,7 @@ class boundaries(Enum):
 
         # Generate mesh
         gmsh.model.mesh.generate(3)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     def geographic_return_coords_to_bounds(coords):
diff --git a/src/underworld3/meshing/segmented.py b/src/underworld3/meshing/segmented.py
index fdf9433b..e71c86a8 100644
--- a/src/underworld3/meshing/segmented.py
+++ b/src/underworld3/meshing/segmented.py
@@ -16,6 +16,7 @@
 import math
 
 import underworld3 as uw
+from underworld3.meshing._mesh_files import mesh_file_dir, write_gmsh
 from underworld3.discretisation import Mesh
 from underworld3 import VarType
 from underworld3.coordinates import CoordinateSystemType
@@ -97,8 +98,8 @@ def SegmentedSphericalSurface2D(
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
-        uw_filename = f".meshes/uw_segmented_spherical_surface_r{radius}_csize{cellSize}_segs{num_segments}.msh"
+            os.makedirs(mesh_file_dir(), exist_ok=True)
+        uw_filename = f"{mesh_file_dir()}/uw_segmented_spherical_surface_r{radius}_csize{cellSize}_segs{num_segments}.msh"
     else:
         uw_filename = filename
 
@@ -182,11 +183,11 @@ def SegmentedSphericalSurface2D(
 
         # Generate Mesh
         gmsh.model.mesh.generate(2)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
 
         # xyz coordinates of the mesh
         xyz = gmsh.model.mesh.get_nodes()[1].reshape(-1, 3)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
         plex_0 = gmsh2dmplex(
@@ -338,8 +339,8 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
-        uw_filename = f".meshes/uw_segmented_sphere_ro{radiusOuter}_ri{radiusInner}_csize{cellSize}_segs{num_segments}.msh"
+            os.makedirs(mesh_file_dir(), exist_ok=True)
+        uw_filename = f"{mesh_file_dir()}/uw_segmented_sphere_ro{radiusOuter}_ri{radiusInner}_csize{cellSize}_segs{num_segments}.msh"
     else:
         uw_filename = filename
 
@@ -541,7 +542,7 @@ class boundaries(Enum):
 
         gmsh.model.mesh.generate(3)
 
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
         # We need to build the plex here in order to make some changes
@@ -786,8 +787,8 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
-        uw_filename = f".meshes/uw_segmented_ball_ro{radius}_csize{cellSize}_segs{num_segments}.msh"
+            os.makedirs(mesh_file_dir(), exist_ok=True)
+        uw_filename = f"{mesh_file_dir()}/uw_segmented_ball_ro{radius}_csize{cellSize}_segs{num_segments}.msh"
     else:
         uw_filename = filename
 
@@ -961,7 +962,7 @@ class boundaries(Enum):
 
         gmsh.model.mesh.generate(3)
 
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
         # We need to build the plex here in order to make some changes
diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py
index 60bfc952..0f79a4d4 100644
--- a/src/underworld3/meshing/spherical.py
+++ b/src/underworld3/meshing/spherical.py
@@ -16,6 +16,7 @@
 import math
 
 import underworld3 as uw
+from underworld3.meshing._mesh_files import mesh_file_dir, write_gmsh
 from underworld3.discretisation import Mesh
 from underworld3 import VarType
 from underworld3.coordinates import CoordinateSystemType
@@ -137,10 +138,10 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
         uw_filename = (
-            f".meshes/uw_spherical_shell_ro{radiusOuter}_ri{radiusInner}_csize{cellSize}.msh"
+            f"{mesh_file_dir()}/uw_spherical_shell_ro{radiusOuter}_ri{radiusInner}_csize{cellSize}.msh"
         )
     else:
         uw_filename = filename
@@ -210,7 +211,7 @@ class boundaries(Enum):
         gmsh.model.occ.synchronize()
 
         gmsh.model.mesh.generate(3)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     # Ensure boundaries conform (if refined)
@@ -382,9 +383,9 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
         uw_filename = (
-            f".meshes/uw_spherical_manifold_r{radius}_csize{cellSize}.msh"
+            f"{mesh_file_dir()}/uw_spherical_manifold_r{radius}_csize{cellSize}.msh"
         )
     else:
         uw_filename = filename
@@ -414,7 +415,7 @@ class boundaries(Enum):
 
         gmsh.option.setNumber("Mesh.CharacteristicLengthMax", cellSize)
         gmsh.model.mesh.generate(2)  # 2-D mesh in 3-D space
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     # Force the PETSc gmsh reader to preserve the 3-D embedding when
@@ -572,9 +573,9 @@ class regions(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f".meshes/uw_spherical_shell_ro{radiusOuter}_rint{radiusInternal}_ri{radiusInner}_csize{cellSize}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_spherical_shell_ro{radiusOuter}_rint{radiusInternal}_ri{radiusInner}_csize{cellSize}.msh"
     else:
         uw_filename = filename
 
@@ -692,7 +693,7 @@ def bbox_radius(dimtag):
         gmsh.model.addPhysicalGroup(shell_vol[0], [shell_vol[1]], 99999, "Elements")
 
         gmsh.model.mesh.generate(3)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     # Ensure boundaries conform (if refined)
@@ -876,9 +877,9 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
+            os.makedirs(mesh_file_dir(), exist_ok=True)
 
-        uw_filename = f".meshes/uw_segmentofsphere_ro{radiusOuter}_ri{radiusInner}_longext{longitudeExtent}_latext{latitudeExtent}_csize{cellSize}.msh"
+        uw_filename = f"{mesh_file_dir()}/uw_segmentofsphere_ro{radiusOuter}_ri{radiusInner}_longext{longitudeExtent}_latext{latitudeExtent}_csize{cellSize}.msh"
     else:
         uw_filename = filename
 
@@ -997,7 +998,7 @@ def getSphericalXYZ(point):
         gmsh.model.occ.synchronize()
 
         gmsh.model.mesh.generate(3)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     # Ensure boundaries conform (if refined)
@@ -1156,8 +1157,8 @@ class boundaries(Enum):
 
     if filename is None:
         if uw.mpi.rank == 0:
-            os.makedirs(".meshes", exist_ok=True)
-        uw_filename = f".meshes/uw_cubed_spherical_shell_ro{radiusOuter}_ri{radiusInner}_elts{numElements}_plex{simplex}.msh"
+            os.makedirs(mesh_file_dir(), exist_ok=True)
+        uw_filename = f"{mesh_file_dir()}/uw_cubed_spherical_shell_ro{radiusOuter}_ri{radiusInner}_elts{numElements}_plex{simplex}.msh"
     else:
         uw_filename = filename
 
@@ -1266,7 +1267,7 @@ class boundaries(Enum):
 
         # Generate Mesh
         gmsh.model.mesh.generate(3)
-        gmsh.write(uw_filename)
+        write_gmsh(uw_filename)
         gmsh.finalize()
 
     def sphere_return_coords_to_bounds(coords):
diff --git a/tests/conftest.py b/tests/conftest.py
index 18cadd3e..2762656c 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,3 +1,5 @@
+import os
+
 # ==============================================================================
 # VISUALIZATION BACKENDS - must be set before any visualization imports
 # ==============================================================================
@@ -29,6 +31,22 @@
 import pytest
 
 
+# ==============================================================================
+# MESH FILES - one directory per xdist worker
+# ==============================================================================
+# Generated mesh files are named from the mesh PARAMETERS, so two workers
+# building the same geometry choose the same name and one can read what the
+# other is still writing (issue #563). The writes are atomic, which makes that
+# safe; giving each worker its own directory also stops them doing the identical
+# work twice. Set at import, before any test builds a mesh.
+#
+# `PYTEST_XDIST_WORKER` is absent in a serial run, which correctly leaves the
+# default `.meshes/` in place.
+_xdist_worker = os.environ.get("PYTEST_XDIST_WORKER")
+if _xdist_worker:
+    os.environ.setdefault("UW_MESH_CACHE_DIR", f".meshes/{_xdist_worker}")
+
+
 @pytest.fixture(scope="function", autouse=True)
 def isolate_test_state(request):
     """
diff --git a/tests/pytest.ini b/tests/pytest.ini
index 2dccd2db..8b6cf6b2 100644
--- a/tests/pytest.ini
+++ b/tests/pytest.ini
@@ -25,6 +25,10 @@ markers =
     tier_c: Experimental tests (development only, not for automation)
 
     # Complexity levels (what kind of test, independent of number prefix)
+    # Select a level by EXCLUDING the ones above it — pytest merges marks, so a
+    # module-level level_1 plus a per-test @pytest.mark.level_2 leaves the test
+    # marked BOTH, and a plain `-m level_1` still selects it:
+    #     pytest -m "level_1 and not level_2 and not level_3"
     level_1: Quick core tests - imports, basic setup, no solving (~seconds)
     level_2: Intermediate tests - integration, units, regression (~minutes)
     level_3: Physics tests - solvers, time-stepping, coupled systems (~minutes to hours)
diff --git a/tests/test_0650_recursion_prevention_regression.py b/tests/test_0650_recursion_prevention_regression.py
index b2b2fbb7..aa698b28 100644
--- a/tests/test_0650_recursion_prevention_regression.py
+++ b/tests/test_0650_recursion_prevention_regression.py
@@ -21,6 +21,25 @@
 import sys
 import os
 
+
+def _headroom(frames):
+    """A recursion limit ``frames`` above the CURRENT stack depth.
+
+    These tests mean "this operation does not recurse without bound", and an
+    absolute ``setrecursionlimit(50)`` does not say that: it also assumes the
+    stack is nearly empty when the test starts. Run under pytest-xdist, whose
+    worker adds its own frames, the budget is spent before the test body
+    begins and the test fails for a reason that has nothing to do with
+    recursion. Measuring from where we actually are keeps the assertion about
+    the operation.
+    """
+    depth = 0
+    frame = sys._getframe()
+    while frame is not None:
+        depth += 1
+        frame = frame.f_back
+    return depth + frames
+
 # Add src to path for testing
 # REMOVED: sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "src"))
 
@@ -38,7 +57,7 @@ def test_uwquantity_atoms_no_recursion(self):
 
         # Set recursion limit to catch infinite recursion quickly
         old_limit = sys.getrecursionlimit()
-        sys.setrecursionlimit(100)  # Low limit to catch recursion fast
+        sys.setrecursionlimit(_headroom(100))  # Low limit to catch recursion fast
 
         try:
             # This was causing infinite recursion before the fix
@@ -86,7 +105,7 @@ def test_mathematical_object_chain_safety(self):
 
         # Create compound expressions (these should not cause recursion)
         old_limit = sys.getrecursionlimit()
-        sys.setrecursionlimit(50)
+        sys.setrecursionlimit(_headroom(50))
 
         try:
             # Mathematical operations should not trigger recursion
@@ -129,7 +148,7 @@ def test_advection_diffusion_parameter_evaluation(self):
 
         # Set recursion limit to catch the issue
         old_limit = sys.getrecursionlimit()
-        sys.setrecursionlimit(50)
+        sys.setrecursionlimit(_headroom(50))
 
         try:
             # This function evaluation was causing recursion in estimate_dt()
@@ -158,7 +177,7 @@ def test_sympy_function_calls_with_uwexpressions(self):
         expr = uw.function.expression(r"func_test", sym=0.5)
 
         old_limit = sys.getrecursionlimit()
-        sys.setrecursionlimit(100)
+        sys.setrecursionlimit(_headroom(100))
 
         try:
             # SymPy functions should not cause recursion when applied to UWexpressions
@@ -183,7 +202,7 @@ def test_sympy_substitution_no_recursion(self):
         expr = uw.function.expression(r"sub_test", sym=sympy.Symbol("x"))
 
         old_limit = sys.getrecursionlimit()
-        sys.setrecursionlimit(100)
+        sys.setrecursionlimit(_headroom(100))
 
         try:
             # Substitution operations should not cause recursion
@@ -204,7 +223,7 @@ def test_sympy_differentiation_no_recursion(self):
         expr = uw.function.expression(r"diff_test", sym=x**2 + 2 * x + 1)
 
         old_limit = sys.getrecursionlimit()
-        sys.setrecursionlimit(100)
+        sys.setrecursionlimit(_headroom(100))
 
         try:
             # Differentiation should not cause recursion
@@ -244,7 +263,7 @@ def test_estimate_dt_no_recursion(self):
         old_limit = sys.getrecursionlimit()
         # Set limit high enough for SymPy tree traversal but low enough to catch infinite loops
         # Original bug (UWQuantity._sympify_() returning self) would hit even high limits
-        sys.setrecursionlimit(300)
+        sys.setrecursionlimit(_headroom(300))
 
         try:
             # This was the specific call that failed with the original recursion bug
@@ -284,7 +303,7 @@ def test_constitutive_model_parameter_access_no_recursion(self):
 
         old_limit = sys.getrecursionlimit()
         # Set reasonable limit to catch infinite recursion but allow normal operations
-        sys.setrecursionlimit(300)
+        sys.setrecursionlimit(_headroom(300))
 
         try:
             # Accessing parameters should not cause recursion
@@ -313,7 +332,7 @@ def recursive_function(n):
             return recursive_function(n - 1)
 
         old_limit = sys.getrecursionlimit()
-        sys.setrecursionlimit(50)
+        sys.setrecursionlimit(_headroom(50))
 
         try:
             # This should hit recursion limit
@@ -346,7 +365,7 @@ def atoms(self, *types):
         safe_obj = SafeObject(sympy.Symbol("x"))
 
         old_limit = sys.getrecursionlimit()
-        sys.setrecursionlimit(50)
+        sys.setrecursionlimit(_headroom(50))
 
         try:
             atoms = safe_obj.atoms(sympy.Symbol)
@@ -369,7 +388,7 @@ def check_for_recursion_risk(obj):
             # The real test: can we call atoms() without infinite recursion?
             import sys
             old_limit = sys.getrecursionlimit()
-            sys.setrecursionlimit(100)
+            sys.setrecursionlimit(_headroom(100))
             try:
                 result = obj.atoms(sympy.Symbol)
                 return False  # No risk - it worked