@@ -869,8 +869,15 @@ def is_storage_related(self, dims=None):
         """
         for d in self.findices:
             if d._defines & set(as_tuple(dims)):  # noqa: SIM102
-                if any(i.is_NonlinearDerived for i in d._defines) or \
-                   self.is_const(d):
+                # An injective Dimension maps distinct iterations to distinct
+                # locations by construction -- e.g. a factor-based
+                # ConditionalDimension, whose guard admits only
+                # `parent % factor == 0` and whose `i -> i//factor` is
+                # one-to-one on those. Unlike, say, a SteppingDimension, where
+                # `time % 2` deliberately revisits the same slot
+                nonlinear = [i for i in d._defines
+                             if i.is_NonlinearDerived and not i.is_Injective]
+                if nonlinear or self.is_const(d):
                     return True
         return False
 
 
@@ -116,6 +116,7 @@ class Dimension(ArgProvider):
     is_Custom = False
     is_Derived = False
     is_NonlinearDerived = False
+    is_Injective = False
     is_AbstractSub = False
     is_Sub = False
     is_MultiSub = False
@@ -1090,6 +1091,19 @@ def __init_finalize__(self, name, parent=None, factor=None, condition=None,
         self._indirect = indirect
         self._relation = relation
 
+    @property
+    def is_Injective(self):
+        """
+        True if distinct iterations of the parent that pass the guard map to
+        distinct indices of `self`, so that no two iterations can access the
+        same memory location through it.
+
+        A pure factor subsampling is injective: the guard admits only
+        `parent % factor == 0`, and `i -> i//factor` is one-to-one on those.
+        An arbitrary `condition` carries no such promise.
+        """
+        return self._factor is not None and self._condition is None
+
     @property
     def uses_symbolic_factor(self):
         return self._factor is not None
 
@@ -5,7 +5,7 @@
 from sympy import Or
 
 from devito import (
-    CondEq, ConditionalDimension, Constant, Eq, Function, Grid, Operator, SubGrid,
+    CondEq, ConditionalDimension, Constant, Eq, Function, Grid, Inc, Operator, SubGrid,
     TimeFunction, switchconfig
 )
 from devito.ir.iet import Conditional, Expression, FindNodes, Iteration
@@ -915,3 +915,132 @@ def test_time_subsampling(self, mode):
 
         for i in range((nt - 2)//factor + 1):
             assert np.allclose(vsave.data[i], i*factor)
+
+
+class TestIncReduction:
+
+    """
+    `Inc` into a factor-subsampled Function must not be treated as a
+    reduction. The write index `i -> i//factor` is one-to-one on the guarded
+    points, so no two iterations touch the same location and the loop is
+    plainly parallel.
+
+    Treating it as one emitted an OpenMP array reduction over the whole
+    Function; with several such Functions the private copies overflowed the
+    stack and the kernel segfaulted.
+    """
+
+    @pytest.mark.parametrize('factor', [1, 2, 3])
+    def test_no_reduction_clause(self, factor):
+        n = 12
+        grid = Grid(shape=(n, n))
+        x, y = grid.dimensions
+        cd = ConditionalDimension(name='xs', parent=x, factor=factor)
+
+        f = Function(name='f', shape=((n + factor - 1)//factor, n),
+                     dimensions=(cd, y), grid=grid, space_order=0)
+
+        with switchconfig(language='openmp'):
+            op = Operator(Inc(f, 1))
+
+        assert 'reduction' not in str(op.ccode)
+
+    @pytest.mark.parametrize('factor', [1, 2, 3])
+    def test_dimensions_are_parallel(self, factor):
+        """
+        The properties must match what the same `Inc` gets on a plain
+        Function, i.e. PARALLEL rather than PARALLEL_IF_ATOMIC.
+        """
+        n = 12
+        grid = Grid(shape=(n, n))
+        x, y = grid.dimensions
+        cd = ConditionalDimension(name='xs', parent=x, factor=factor)
+
+        f = Function(name='f', shape=((n + factor - 1)//factor, n),
+                     dimensions=(cd, y), grid=grid, space_order=0)
+
+        op = Operator(Inc(f, 1))
+
+        for i in FindNodes(Iteration).visit(op):
+            assert i.is_Parallel
+
+    @pytest.mark.parametrize('factor', [1, 2, 3])
+    def test_accumulates(self, factor):
+        """Dropping the reduction must not drop the accumulation."""
+        n = 12
+        grid = Grid(shape=(n, n))
+        x, y = grid.dimensions
+        cd = ConditionalDimension(name='xs', parent=x, factor=factor)
+
+        f = Function(name='f', shape=((n + factor - 1)//factor, n),
+                     dimensions=(cd, y), grid=grid, space_order=0)
+        f.data[:] = 0.
+
+        op = Operator(Inc(f, 1))
+        op.apply()
+        op.apply()
+
+        assert np.all(f.data == 2.)
+
+    def test_many_functions(self):
+        """
+        The original failure: nine subsampled Functions incremented together.
+        The spurious array reductions needed ~13 MiB of thread-private copies
+        and overflowed the stack.
+        """
+        n = 40
+        grid = Grid(shape=(n, n, n))
+        dims = tuple(ConditionalDimension(name=f'{d.name}s', parent=d, factor=1)
+                     for d in grid.dimensions)
+
+        fs = [Function(name=f'og_{i}', grid=grid, shape=(n, n, n),
+                       dimensions=dims) for i in range(9)]
+        for f in fs:
+            f.data[:] = 0.
+
+        with switchconfig(language='openmp'):
+            op = Operator([Inc(f, 1) for f in fs])
+            assert 'reduction' not in str(op.ccode)
+            op.apply()
+
+        for f in fs:
+            assert np.all(f.data == 1.)
+
+
+class TestInjectivity:
+
+    """`is_Injective` must be claimed only where it holds."""
+
+    def test_factor_is_injective(self):
+        grid = Grid(shape=(8, 8))
+        x, _ = grid.dimensions
+
+        assert ConditionalDimension(name='a', parent=x, factor=2).is_Injective
+
+    def test_stepping_is_not(self):
+        """`time % 2` deliberately revisits the same slot."""
+        grid = Grid(shape=(8, 8))
+
+        assert not grid.stepping_dim.is_Injective
+
+    def test_condition_is_not(self):
+        """An arbitrary condition carries no one-to-one promise."""
+        grid = Grid(shape=(8, 8))
+        x, _ = grid.dimensions
+        cd = ConditionalDimension(name='b', parent=x, condition=CondEq(x, 0))
+
+        assert not cd.is_Injective
+
+    def test_factor_with_condition_is_not(self):
+        grid = Grid(shape=(8, 8))
+        x, _ = grid.dimensions
+        cd = ConditionalDimension(name='c', parent=x, factor=2,
+                                  condition=CondEq(x % 4, 0))
+
+        assert not cd.is_Injective
+
+    def test_plain_dimension_is_not(self):
+        grid = Grid(shape=(8, 8))
+        x, _ = grid.dimensions
+
+        assert not x.is_Injective