diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py
index 39bf5692..d5082199 100644
--- a/src/underworld3/systems/ddt.py
+++ b/src/underworld3/systems/ddt.py
@@ -1765,11 +1765,28 @@ def update_pre_solve(
 
             except Exception:
                 # Fallback to projection solver for expressions that can't be directly evaluated
-                # (e.g., containing derivatives)
-                self._psi_star_projection_solver.uw_function = self.psi_fn
+                # (e.g., containing derivatives — true for the NS viscous flux every step).
+                # Route via _build_projection_source so the (1, Nc) row-matrix flattening
+                # required by SNES_MultiComponent_Projection is applied for tensor vtypes.
+                # Without this, a (dim, dim) tensor function meets a (1, Nc) solver field
+                # and SymPy raises "Matrix size mismatch: (1, Nc) + (dim, dim)" (issue #180).
+                self._psi_star_projection_solver.uw_function = self._build_projection_source(
+                    self.psi_fn
+                )
                 self._psi_star_projection_solver.smoothing = 0.0
                 self._psi_star_projection_solver.solve(verbose=verbose)
 
+                # For tensor vtypes the projection writes into the flat (1, Nc) variable,
+                # so we must fan it back out to psi_star[0] — otherwise subsequent
+                # history operations read a stale tensor. Mirrors the same fan-out in
+                # initialise_history() (~line 1540).
+                if getattr(self, '_psi_star_use_multicomponent', False):
+                    for k, (i, j) in enumerate(self._psi_star_indep_indices):
+                        vals = self._psi_star_flat_var.array[:, 0, k]
+                        self.psi_star[0].array[:, i, j] = vals
+                        if i != j:
+                            self.psi_star[0].array[:, j, i] = vals
+
         # 3. Compute the upstream values from the psi_fn
 
         # We use the u_star variable as a working value here so we have to work backwards
diff --git a/tests/test_0610_navier_stokes_slcn_projection.py b/tests/test_0610_navier_stokes_slcn_projection.py
new file mode 100644
index 00000000..a849ebe6
--- /dev/null
+++ b/tests/test_0610_navier_stokes_slcn_projection.py
@@ -0,0 +1,54 @@
+"""
+Regression test for the NavierStokesSLCN DFDt projection fallback.
+
+When ``SNES_MultiComponent_Projection`` was wired into the ``SemiLagrangian``
+DDt path, the ``psi_fn`` setter and ``_setup_projections`` were updated to
+flatten the source tensor to a ``(1, Nc)`` row matrix via
+``_build_projection_source``. The fallback path inside ``update_pre_solve``
+(taken when ``uw.function.evaluate`` raises on expressions containing
+derivatives — true for the NavierStokes viscous flux every step) missed the
+migration and assigned ``self.psi_fn`` directly, producing:
+
+    sympy.matrices.exceptions.ShapeError:
+        Matrix size mismatch: (1, 3) + (2, 2).
+
+See: https://github.com/underworldcode/underworld3/issues/180
+"""
+import pytest
+import sympy
+import underworld3 as uw
+
+
+@pytest.mark.level_2
+@pytest.mark.tier_a
+def test_navier_stokes_slcn_solve_does_not_raise_shape_error():
+    """First solve completes without ShapeError from the DDt projection fallback."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0),
+        maxCoords=(1.0, 1.0),
+        cellSize=1 / 8,
+        regular=False,
+    )
+
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, continuous=True)
+
+    ns = uw.systems.NavierStokes(
+        mesh,
+        velocityField=v,
+        pressureField=p,
+        rho=1.0,
+        order=2,
+    )
+    ns.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    ns.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+
+    ns.bodyforce = sympy.Matrix([0, 0]).T
+    ns.add_dirichlet_bc((1.0, 0.0), "Top")
+    ns.add_dirichlet_bc((0.0, 0.0), "Bottom")
+    ns.add_dirichlet_bc((0.0, 0.0), "Left")
+    ns.add_dirichlet_bc((0.0, 0.0), "Right")
+
+    # Take a single step. Pre-fix this raised ShapeError("(1, 3) + (2, 2)")
+    # from the DDt projection fallback. Post-fix it just converges.
+    ns.solve(timestep=0.01, zero_init_guess=True)