compiler: Guard device loops over SubDimensions against empty iteration spaces

devito/passes/iet/parpragma.py

@@ -580,10 +580,13 @@ def _make_guard(self, parregion, *args):
             return super()._make_guard(parregion, *args)
 
         cond = []
-        # There must be at least one iteration or potential crash
-        if not parregion.is_Affine:
-            trees = retrieve_iteration_tree(parregion.root)
-            tree = trees[0][:parregion.ncollapsed]
+        # There must be at least one iteration or potential crash. The bounds of
+        # a SubDimension depend on runtime thicknesses, so its iteration space may
+        # be empty, possibly with a negative size, even if the region is affine;
+        # some compilers (e.g., nvc) then launch a huge number of empty gangs
+        trees = retrieve_iteration_tree(parregion.root)
+        tree = trees[0][:parregion.ncollapsed]
+        if not parregion.is_Affine or any(i.dim.is_AbstractSub for i in tree):
             cond.extend([i.symbolic_size > 0 for i in tree])
 
         # SparseFunctions may occasionally degenerate to zero-size arrays. In such

tests/test_gpu_common.py

@@ -1102,7 +1102,9 @@ def test_composite_buffering_tasking(self):
         op1 = Operator(eqns, opt=('buffering', 'tasking', 'orchestrate'))
 
         # Check generated code -- thanks to buffering only expect 1 lock!
-        assert len(retrieve_iteration_tree(op0)) == 2
+        # In `op0`, the SubDomain loop nest is guarded against an empty iteration
+        # space, so the time loop also hosts a Conditional (an extra section)
+        assert len(retrieve_iteration_tree(op0)) == 3
         assert len(retrieve_iteration_tree(op1)) == 3
         symbols = FindSymbols().visit(op1)
         assert len([i for i in symbols if isinstance(i, Lock)]) == 1
@@ -1141,7 +1143,9 @@ def test_composite_buffering_tasking_multi_output(self):
                                    'npthreads': 1}))
 
         # Check generated code -- thanks to buffering only expect 1 lock!
-        assert len(retrieve_iteration_tree(op0)) == 2
+        # In `op0`, the SubDomain loop nest is guarded against an empty iteration
+        # space, so the time loop also hosts a Conditional (an extra section)
+        assert len(retrieve_iteration_tree(op0)) == 3
         assert len(retrieve_iteration_tree(op1)) == 4
         assert len(retrieve_iteration_tree(op2)) == 3
         symbols = FindSymbols().visit(op1)

tests/test_subdomains.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 import pytest
-from sympy import sin, tan
+from sympy import And, sin, tan
 
 from conftest import assert_structure, opts_tiling
 from devito import (
@@ -12,7 +12,8 @@
 )
 from devito.exceptions import InvalidArgument
 from devito.ir import (
-    Expression, FindNodes, FindSymbols, Iteration, SymbolRegistry, retrieve_iteration_tree
+    Conditional, Expression, FindNodes, FindSymbols, Iteration, SymbolRegistry,
+    retrieve_iteration_tree
 )
 from devito.tools import timed_region
 
@@ -932,6 +933,30 @@ def test_subdomainset_w_sparse(self):
         check[0, 0] = 1
         assert np.all(np.isclose(f.data, check))
 
+    @pytest.mark.parametrize('language', ['openacc', 'openmp'])
+    def test_device_guard_empty_ispace(self, language):
+        """
+        The thicknesses of a SubDomainSet are runtime values, so its iteration
+        space may be empty, even with a negative size (e.g., in the MPI regions
+        outside of a subdomain). A device parallel loop must not be entered
+        in that case, as some compilers (e.g., nvc) may otherwise launch a huge
+        number of empty gangs. Generated code only, so it runs without a GPU.
+        """
+        grid = Grid(shape=(11, 11))
+        f = TimeFunction(name='f', grid=grid)
+
+        sds = SubDomainSet(N=1, bounds=(2, 2, 2, 2), grid=grid)
+
+        op = Operator(Eq(f.forward, f + 1, subdomain=sds),
+                      platform='nvidiaX', language=language)
+
+        iterations = [i for i in FindNodes(Iteration).visit(op) if i.dim.is_MultiSub]
+        assert len(iterations) == 2
+        conds = [c for c in FindNodes(Conditional).visit(op)
+                 if iterations[0] in FindNodes(Iteration).visit(c.then_body)]
+        assert len(conds) == 1
+        assert conds[0].condition == And(*[i.symbolic_size > 0 for i in iterations])
+
 
 class TestBorder:
     # Note: This class is partially covered by doctests