Merge branch 'threadSafety' into 'master'

Thread safety

See merge request !391

src/pystencils/boundaries/boundaryhandling.py

@@ -450,5 +450,6 @@ def create_boundary_kernel(field, index_field, stencil, boundary_functor, target
     dir_symbol = TypedSymbol("dir", np.int32)
     elements += [SympyAssignment(dir_symbol, index_field[0]('dir'))]
     elements += boundary_functor(field, direction_symbol=dir_symbol, index_field=index_field)
-    config = CreateKernelConfig(index_fields=[index_field], target=target, **kernel_creation_args)
+    config = CreateKernelConfig(index_fields=[index_field], target=target, skip_independence_check=True,
+                                **kernel_creation_args)
     return create_kernel(elements, config=config)

src/pystencils/config.py

@@ -135,8 +135,9 @@ class CreateKernelConfig:
     """
     skip_independence_check: bool = False
     """
-    Don't check that loop iterations are independent. This is needed e.g. for
-    periodicity kernel, that access the field outside the iteration bounds. Use with care!
+    By default the assignment list is checked for read/write independence. This means fields are only written at
+    locations where they are read. Doing so guarantees thread safety. In some cases e.g. for
+    periodicity kernel, this can not be assured and does the check needs to be deactivated. Use with care!
     """
 
     class DataTypeFactory:

src/pystencils/kernel_contrains_check.py

@@ -38,6 +38,7 @@ class KernelConstraintsCheck:
 
     def __init__(self, check_independence_condition=True, check_double_write_condition=True):
         self.scopes = NestedScopes()
+        self.field_reads = defaultdict(set)
         self.field_writes = defaultdict(set)
         self.fields_read = set()
         self.check_independence_condition = check_independence_condition
@@ -111,6 +112,13 @@ class KernelConstraintsCheck:
             if self.check_double_write_condition and len(self.field_writes[fai]) > 1:
                 raise ValueError(
                     f"Field {lhs.field.name} is written at two different locations")
+
+            if fai in self.field_reads:
+                reads = tuple(self.field_reads[fai])
+                if len(reads) > 1 or lhs.offsets != reads[0]:
+                    if self.check_independence_condition:
+                        raise ValueError(f"Field {lhs.field.name} is written at different location than it was read. "
+                                         f"This means the resulting kernel would not be thread safe")
         elif isinstance(lhs, sp.Symbol):
             if self.scopes.is_defined_locally(lhs):
                 raise ValueError(f"Assignments not in SSA form, multiple assignments to {lhs.name}")
@@ -120,8 +128,9 @@ class KernelConstraintsCheck:
 
     def update_accesses_rhs(self, rhs):
         if isinstance(rhs, Field.Access) and self.check_independence_condition:
-            writes = self.field_writes[self.FieldAndIndex(
-                rhs.field, rhs.index)]
+            fai = self.FieldAndIndex(rhs.field, rhs.index)
+            writes = self.field_writes[fai]
+            self.field_reads[fai].add(rhs.offsets)
             for write_offset in writes:
                 assert len(writes) == 1
                 if write_offset != rhs.offsets:

src/pystencils/kernelcreation.py

@@ -128,6 +128,7 @@ def create_domain_kernel(assignments: NodeCollection, *, config: CreateKernelCon
     # --- check constrains
     check = KernelConstraintsCheck(check_independence_condition=not config.skip_independence_check,
                                    check_double_write_condition=not config.allow_double_writes)
+
     check.visit(assignments)
 
     assignments.bound_fields = check.fields_written

tests/test_vectorization_specific.py

@@ -226,7 +226,7 @@ def test_issue62(dtype, instruction_set, padding):
     dy = sp.Symbol("dy")
     src, dst, rhs = ps.fields(f"src, src_tmp, rhs: {dtype}[2D]", layout='fzyx')
 
-    up = ps.Assignment(src[0, 0], ((dy ** 2 * (src[1, 0] + src[-1, 0]))
+    up = ps.Assignment(dst[0, 0], ((dy ** 2 * (src[1, 0] + src[-1, 0]))
                                    + (dx ** 2 * (src[0, 1] + src[0, -1]))
                                    - (rhs[0, 0] * dx ** 2 * dy ** 2)) / (2 * (dx ** 2 + dy ** 2)))