List LBM: Boundary conditions & new mapping class

boundaries/boundaryhandling.py

@@ -395,9 +395,9 @@ class BoundaryOffsetInfo(CustomCppCode):
 
     @staticmethod
     def _offset_symbols(dim):
-        return [TypedSymbol("c_%d" % (d,), create_type(np.int64)) for d in range(dim)]
+        return [TypedSymbol("c%s" % (d,), create_type(np.int64)) for d in ['x', 'y', 'z'][:dim]]
 
-    INV_DIR_SYMBOL = TypedSymbol("inv_dir", "int")
+    INV_DIR_SYMBOL = TypedSymbol("invdir", "int")
 
 
 def create_boundary_kernel(field, index_field, stencil, boundary_functor, target='cpu', openmp=True):

integer_functions.py

@@ -55,7 +55,7 @@ class modulo_ceil(sp.Function):
         >>> from pystencils import TypedSymbol
         >>> a, b = TypedSymbol("a", "int64"), TypedSymbol("b", "int32")
         >>> modulo_ceil(a, b).to_c(str)
-        '(a) % (b) == 0 ? a : ((int64_t)((a) / (b))+1) * (b)'
+        '((a) % (b) == 0 ? a : ((int64_t)((a) / (b))+1) * (b))'
     """
     nargs = 2
 

transformations.py

@@ -72,6 +72,20 @@ def get_common_shape(field_set):
     return shape
 
 
+def get_field_accesses(expr, result=set()):
+    if isinstance(expr, Field.Access):
+        result.add(expr)
+        for o in expr.offsets:
+            get_field_accesses(o, result)
+        for i in expr.index:
+            get_field_accesses(i, result)
+    elif hasattr(expr, 'atoms'):
+        new_accesses = expr.atoms(Field.Access)
+        result.update(new_accesses)
+        for a in new_accesses:
+            get_field_accesses(a, result)
+
+
 def make_loop_over_domain(body, function_name, iteration_slice=None, ghost_layers=None, loop_order=None):
     """Uses :class:`pystencils.field.Field.Access` to create (multiple) loops around given AST.
 
@@ -88,7 +102,10 @@ def make_loop_over_domain(body, function_name, iteration_slice=None, ghost_layer
         :class:`LoopOverCoordinate` instance with nested loops, ordered according to field layouts
     """
     # find correct ordering by inspecting participating FieldAccesses
-    field_accesses = body.atoms(Field.Access)
+    field_accesses = set()
+    get_field_accesses(body, field_accesses)
+    field_accesses = {e for e in field_accesses if not e.is_absolute_access}
+
     # exclude accesses to buffers from field_list, because buffers are treated separately
     field_list = [e.field for e in field_accesses if not FieldType.is_buffer(e.field)]
     fields = set(field_list)