src/pystencils/backends/__init__.py

+from .cbackend import generate_c
+
+__all__ = ['generate_c']
+try:
+    from .dot import print_dot  # NOQA
+    __all__.append('print_dot')
+except ImportError:
+    pass

src/pystencils/backends/arm_instruction_sets.py

+from pystencils.typing import CFunction
+
+
+def get_argument_string(function_shortcut, first=''):
+    args = function_shortcut[function_shortcut.index('[') + 1: -1]
+    arg_string = "("
+    if first:
+        arg_string += first + ', '
+    for arg in args.split(","):
+        arg = arg.strip()
+        if not arg:
+            continue
+        if arg in ('0', '1', '2', '3', '4', '5'):
+            arg_string += "{" + arg + "},"
+        else:
+            arg_string += arg + ","
+    arg_string = arg_string[:-1] + ")"
+    return arg_string
+
+
+def get_vector_instruction_set_arm(data_type='double', instruction_set='neon'):
+    if instruction_set not in ['neon', 'sme'] and not instruction_set.startswith('sve'):
+        raise NotImplementedError(instruction_set)
+    if instruction_set in ['sve', 'sve2', 'sme']:
+        cmp = 'cmp'
+    elif instruction_set.startswith('sve2') and instruction_set not in ('sve256', 'sve2048'):
+        cmp = 'cmp'
+        bitwidth = int(instruction_set[4:])
+    elif instruction_set.startswith('sve'):
+        cmp = 'cmp'
+        bitwidth = int(instruction_set[3:])
+    elif instruction_set == 'neon':
+        cmp = 'c'
+        bitwidth = 128
+
+    base_names = {
+        '+': 'add[0, 1]',
+        '-': 'sub[0, 1]',
+        '*': 'mul[0, 1]',
+        '/': 'div[0, 1]',
+        'sqrt': 'sqrt[0]',
+
+        'loadU': 'ld1[0]',
+        'storeU': 'st1[0, 1]',
+
+        'abs': 'abs[0]',
+        '==': f'{cmp}eq[0, 1]',
+        '!=': f'{cmp}eq[0, 1]',
+        '<=': f'{cmp}le[0, 1]',
+        '<': f'{cmp}lt[0, 1]',
+        '>=': f'{cmp}ge[0, 1]',
+        '>': f'{cmp}gt[0, 1]',
+    }
+
+    bits = {'double': 64,
+            'float': 32,
+            'int': 32}
+
+    result = dict()
+
+    if instruction_set in ['sve', 'sve2', 'sme']:
+        width = 'svcntd()' if data_type == 'double' else 'svcntw()'
+        intwidth = 'svcntw()'
+        result['bytes'] = 'svcntb()'
+    else:
+        width = bitwidth // bits[data_type]
+        intwidth = bitwidth // bits['int']
+        result['bytes'] = bitwidth // 8
+    if instruction_set.startswith('sve') or instruction_set == 'sme':
+        base_names['stream'] = 'stnt1[0, 1]'
+        prefix = 'sv'
+        suffix = f'_f{bits[data_type]}'
+    elif instruction_set == 'neon':
+        prefix = 'v'
+        suffix = f'q_f{bits[data_type]}'
+
+    if instruction_set in ['sve', 'sve2', 'sme']:
+        predicate = f'{prefix}whilelt_b{bits[data_type]}_u64({{loop_counter}}, {{loop_stop}})'
+        int_predicate = f'{prefix}whilelt_b{bits["int"]}_u64({{loop_counter}}, {{loop_stop}})'
+    else:
+        predicate = f'{prefix}whilelt_b{bits[data_type]}(0, {width})'
+        int_predicate = f'{prefix}whilelt_b{bits["int"]}(0, {intwidth})'
+
+    for intrinsic_id, function_shortcut in base_names.items():
+        function_shortcut = function_shortcut.strip()
+        name = function_shortcut[:function_shortcut.index('[')]
+
+        arg_string = get_argument_string(function_shortcut, first=predicate if prefix == 'sv' else '')
+        if prefix == 'sv' and not name.startswith('ld') and not name.startswith('st') and not name.startswith(cmp):
+            undef = '_x'
+        else:
+            undef = ''
+
+        result[intrinsic_id] = prefix + name + suffix + undef + arg_string
+
+    if instruction_set in ['sve', 'sve2', 'sme']:
+        result['width'] = CFunction(width, "int")
+        result['intwidth'] = CFunction(intwidth, "int")
+    else:
+        result['width'] = width
+        result['intwidth'] = intwidth
+
+    if instruction_set.startswith('sve') or instruction_set == 'sme':
+        result['makeVecConst'] = f'svdup_f{bits[data_type]}' + '({0})'
+        result['makeVecConstInt'] = f'svdup_s{bits["int"]}' + '({0})'
+        result['makeVecIndex'] = f'svindex_s{bits["int"]}' + '({0}, {1})'
+
+        if instruction_set != 'sme':
+            vindex = f'svindex_u{bits[data_type]}(0, {{0}})'
+            result['storeS'] = f'svst1_scatter_u{bits[data_type]}index_f{bits[data_type]}({predicate}, {{0}}, ' + \
+                               vindex.format("{2}") + ', {1})'
+            result['loadS'] = f'svld1_gather_u{bits[data_type]}index_f{bits[data_type]}({predicate}, {{0}}, ' + \
+                              vindex.format("{1}") + ')'
+        if instruction_set.startswith('sve2') and instruction_set not in ('sve256', 'sve2048'):
+            result['streamS'] = f'svstnt1_scatter_u{bits[data_type]}offset_f{bits[data_type]}({predicate}, {{0}}, ' + \
+                                vindex.format(f"{{2}}*{bits[data_type]//8}") + ', {1})'
+
+        result['+int'] = f"svadd_s{bits['int']}_x({int_predicate}, " + "{0}, {1})"
+
+        result['float'] = f'svfloat{bits["float"]}_{"s" if instruction_set not in ["sve", "sve2", "sme"] else ""}t'
+        result['double'] = f'svfloat{bits["double"]}_{"s" if instruction_set not in ["sve", "sve2", "sme"] else ""}t'
+        result['int'] = f'svint{bits["int"]}_{"s" if instruction_set not in ["sve", "sve2", "sme"] else ""}t'
+        result['bool'] = f'svbool_{"s" if instruction_set not in ["sve", "sve2", "sme"] else ""}t'
+
+        result['headers'] = ['<arm_sve.h>', '<arm_acle.h>', '"arm_neon_helpers.h"']
+
+        result['&'] = f'svand_b_z({predicate},' + ' {0}, {1})'
+        result['|'] = f'svorr_b_z({predicate},' + ' {0}, {1})'
+        result['blendv'] = f'svsel_f{bits[data_type]}' + '({2}, {1}, {0})'
+        result['any'] = f'svptest_any({predicate}, {{0}})'
+        result['all'] = f'svcntp_b{bits[data_type]}({predicate}, {{0}}) == {width}'
+
+        result['maskStoreU'] = result['storeU'].replace(predicate, '{2}')
+        result['maskStream'] = result['stream'].replace(predicate, '{2}')
+        if instruction_set != 'sme':
+            result['maskStoreS'] = result['storeS'].replace(predicate, '{3}')
+            if instruction_set.startswith('sve2') and instruction_set not in ('sve256', 'sve2048'):
+                result['maskStreamS'] = result['streamS'].replace(predicate, '{3}')
+
+        result['streamFence'] = '__dmb(15)'
+
+        if instruction_set == 'sme':
+            result['function_prefix'] = '__attribute__((arm_locally_streaming))'
+        elif instruction_set not in ['sve', 'sve2', 'sme']:
+            result['compile_flags'] = [f'-msve-vector-bits={bitwidth}']
+    else:
+        result['makeVecConst'] = f'vdupq_n_f{bits[data_type]}' + '({0})'
+        result['makeVec'] = f'makeVec_f{bits[data_type]}' + '(' + ", ".join(['{' + str(i) + '}' for i in
+                                                                             range(width)]) + ')'
+        result['makeVecConstInt'] = f'vdupq_n_s{bits["int"]}' + '({0})'
+        result['makeVecInt'] = f'makeVec_s{bits["int"]}' + '({0}, {1}, {2}, {3})'
+
+        result['+int'] = f"vaddq_s{bits['int']}" + "({0}, {1})"
+
+        result[data_type] = f'float{bits[data_type]}x{width}_t'
+        result['int'] = f'int{bits["int"]}x{intwidth}_t'
+        result['bool'] = f'uint{bits[data_type]}x{width}_t'
+
+        result['headers'] = ['<arm_neon.h>', '"arm_neon_helpers.h"']
+
+        result['!='] = f'vmvnq_u{bits[data_type]}({result["=="]})'
+
+        result['&'] = f'vandq_u{bits[data_type]}' + '({0}, {1})'
+        result['|'] = f'vorrq_u{bits[data_type]}' + '({0}, {1})'
+        result['blendv'] = f'vbslq_f{bits[data_type]}' + '({2}, {1}, {0})'
+        result['any'] = f'vaddlvq_u8(vreinterpretq_u8_u{bits[data_type]}({{0}})) > 0'
+        result['all'] = f'vaddlvq_u8(vreinterpretq_u8_u{bits[data_type]}({{0}})) == 16*0xff'
+
+        # SVE has real nontemporal stores, so we only need to zero cachlines on Neon
+        result['cachelineZero'] = 'cachelineZero((void*) {0})'
+
+    result['cachelineSize'] = 'cachelineSize()'
+
+    return result

src/pystencils/backends/cbackend.py

+import re
+from collections import namedtuple
+import hashlib
+from typing import Set
+
+import numpy as np
+import sympy as sp
+from sympy.core import S
+from sympy.logic.boolalg import BooleanFalse, BooleanTrue
+from sympy.functions.elementary.trigonometric import TrigonometricFunction, InverseTrigonometricFunction
+from sympy.functions.elementary.hyperbolic import HyperbolicFunction
+
+from pystencils.astnodes import KernelFunction, LoopOverCoordinate, Node
+from pystencils.cpu.vectorization import vec_all, vec_any, CachelineSize
+from pystencils.typing import (
+    PointerType, VectorType, CastFunc, create_type, get_type_of_expression,
+    ReinterpretCastFunc, VectorMemoryAccess, BasicType, TypedSymbol, CFunction)
+from pystencils.enums import Backend
+from pystencils.fast_approximation import fast_division, fast_inv_sqrt, fast_sqrt
+from pystencils.functions import DivFunc, AddressOf
+from pystencils.integer_functions import (
+    bit_shift_left, bit_shift_right, bitwise_and, bitwise_or, bitwise_xor,
+    int_div, int_power_of_2, modulo_ceil)
+
+try:
+    from sympy.printing.c import C99CodePrinter as CCodePrinter  # for sympy versions > 1.6
+except ImportError:
+    from sympy.printing.ccode import C99CodePrinter as CCodePrinter
+
+__all__ = ['generate_c', 'CustomCodeNode', 'PrintNode', 'get_headers', 'CustomSympyPrinter']
+
+
+HEADER_REGEX = re.compile(r'^[<"].*[">]$')
+
+
+def generate_c(ast_node: Node,
+               signature_only: bool = False,
+               dialect: Backend = Backend.C,
+               custom_backend=None,
+               with_globals=True) -> str:
+    """Prints an abstract syntax tree node as C or CUDA code.
+
+    This function does not need to distinguish for most AST nodes between C, C++ or CUDA code, it just prints 'C-like'
+    code as encoded in the abstract syntax tree (AST). The AST is built differently for C or CUDA by calling different
+    create_kernel functions.
+
+    Args:
+        ast_node: ast representation of kernel
+        signature_only: generate signature without function body
+        dialect: `Backend`: 'C' or 'CUDA'
+        custom_backend: use own custom printer for code generation
+        with_globals: enable usage of global variables
+    Returns:
+        C-like code for the ast node and its descendants
+    """
+    global_declarations = get_global_declarations(ast_node)
+    for d in global_declarations:
+        if hasattr(ast_node, "global_variables"):
+            ast_node.global_variables.update(d.symbols_defined)
+        else:
+            ast_node.global_variables = d.symbols_defined
+    if custom_backend:
+        printer = custom_backend
+    elif dialect == Backend.C:
+        try:
+            # TODO Vectorization Revamp: instruction_set should not be just slapped on ast
+            instruction_set = ast_node.instruction_set
+        except Exception:
+            instruction_set = None
+        printer = CBackend(signature_only=signature_only,
+                           vector_instruction_set=instruction_set)
+    elif dialect == Backend.CUDA:
+        from pystencils.backends.cuda_backend import CudaBackend
+        printer = CudaBackend(signature_only=signature_only)
+    else:
+        raise ValueError(f'Unknown {dialect=}')
+    code = printer(ast_node)
+    if not signature_only and isinstance(ast_node, KernelFunction):
+        if with_globals and global_declarations:
+            code = "\n" + code
+            for declaration in global_declarations:
+                code = printer(declaration) + "\n" + code
+
+    return code
+
+
+def get_global_declarations(ast):
+    global_declarations = []
+
+    def visit_node(sub_ast):
+        nonlocal global_declarations
+        if hasattr(sub_ast, "required_global_declarations"):
+            global_declarations += sub_ast.required_global_declarations
+
+        if hasattr(sub_ast, "args"):
+            for node in sub_ast.args:
+                visit_node(node)
+
+    visit_node(ast)
+
+    return sorted(set(global_declarations), key=str)
+
+
+def get_headers(ast_node: Node) -> Set[str]:
+    """Return a set of header files, necessary to compile the printed C-like code."""
+    headers = set()
+
+    if isinstance(ast_node, KernelFunction) and ast_node.instruction_set:
+        headers.update(ast_node.instruction_set['headers'])
+
+    if hasattr(ast_node, 'headers'):
+        headers.update(ast_node.headers)
+    for a in ast_node.args:
+        if isinstance(a, (sp.Expr, Node)):
+            headers.update(get_headers(a))
+
+    for g in get_global_declarations(ast_node):
+        if isinstance(g, Node):
+            headers.update(get_headers(g))
+
+    for h in headers:
+        assert HEADER_REGEX.match(h), f'header /{h}/ does not follow the pattern /"..."/ or /<...>/'
+
+    return headers
+
+
+# --------------------------------------- Backend Specific Nodes -------------------------------------------------------
+
+# TODO future CustomCodeNode should not be backend specific move it elsewhere
+class CustomCodeNode(Node):
+    def __init__(self, code, symbols_read, symbols_defined, parent=None):
+        super(CustomCodeNode, self).__init__(parent=parent)
+        self._code = "\n" + code
+        self._symbols_read = set(symbols_read)
+        self._symbols_defined = set(symbols_defined)
+        self.headers = []
+
+    def get_code(self, dialect, vector_instruction_set, print_arg):
+        return self._code
+
+    @property
+    def args(self):
+        return []
+
+    @property
+    def symbols_defined(self):
+        return self._symbols_defined
+
+    @property
+    def undefined_symbols(self):
+        return self._symbols_read - self._symbols_defined
+
+    def __eq__(self, other):
+        return type(self) is type(other) and self._code == other._code
+
+    def __hash__(self):
+        return hash(self._code)
+
+
+class PrintNode(CustomCodeNode):
+    # noinspection SpellCheckingInspection
+    def __init__(self, symbol_to_print):
+        code = f'\nstd::cout << "{symbol_to_print.name}  =  " << {symbol_to_print.name} << std::endl; \n'
+        super(PrintNode, self).__init__(code, symbols_read=[symbol_to_print], symbols_defined=set())
+        self.headers.append("<iostream>")
+
+
+# ------------------------------------------- Printer ------------------------------------------------------------------
+
+
+# noinspection PyPep8Naming
+class CBackend:
+
+    def __init__(self, sympy_printer=None, signature_only=False, vector_instruction_set=None, dialect=Backend.C):
+        if sympy_printer is None:
+            if vector_instruction_set is not None:
+                self.sympy_printer = VectorizedCustomSympyPrinter(vector_instruction_set)
+            else:
+                self.sympy_printer = CustomSympyPrinter()
+        else:
+            self.sympy_printer = sympy_printer
+
+        self._vector_instruction_set = vector_instruction_set
+        self._indent = "   "
+        self._dialect = dialect
+        self._signatureOnly = signature_only
+        self._kwargs = {}
+        self.sympy_printer._kwargs = self._kwargs
+
+    def __call__(self, node):
+        prev_is = VectorType.instruction_set
+        VectorType.instruction_set = self._vector_instruction_set
+        result = str(self._print(node))
+        VectorType.instruction_set = prev_is
+        return result
+
+    def _print(self, node):
+        if isinstance(node, str):
+            return node
+        for cls in type(node).__mro__:
+            method_name = f"_print_{cls.__name__}"
+            if hasattr(self, method_name):
+                return getattr(self, method_name)(node)
+        raise NotImplementedError(f"{self.__class__.__name__} does not support node of type {node.__class__.__name__}")
+
+    def _print_AbstractType(self, node):
+        return str(node)
+
+    def _print_KernelFunction(self, node):
+        function_arguments = [f"{self._print(s.symbol.dtype)} {s.symbol.name}" for s in node.get_parameters()
+                              if not type(s.symbol) is CFunction]
+        launch_bounds = ""
+        if self._dialect == Backend.CUDA:
+            max_threads = node.indexing.max_threads_per_block()
+            if max_threads:
+                launch_bounds = f"__launch_bounds__({max_threads}) "
+        func_declaration = "FUNC_PREFIX %svoid %s(%s)" % (launch_bounds, node.function_name,
+                                                          ", ".join(function_arguments))
+        if self._signatureOnly:
+            return func_declaration
+
+        body = self._print(node.body)
+        return func_declaration + "\n" + body
+
+    def _print_Block(self, node):
+        block_contents = "\n".join([self._print(child) for child in node.args])
+        return "{\n%s\n}" % (self._indent + self._indent.join(block_contents.splitlines(True)))
+
+    def _print_PragmaBlock(self, node):
+        return f"{node.pragma_line}\n{self._print_Block(node)}"
+
+    def _print_LoopOverCoordinate(self, node):
+        counter_name = node.loop_counter_name
+        counter_dtype = node.loop_counter_symbol.dtype.c_name
+        start = f"{counter_dtype} {counter_name} = {self.sympy_printer.doprint(node.start)}"
+        condition = f"{counter_name} < {self.sympy_printer.doprint(node.stop)}"
+        update = f"{counter_name} += {self.sympy_printer.doprint(node.step)}"
+        loop_str = f"for ({start}; {condition}; {update})"
+        self._kwargs['loop_counter'] = counter_name
+        self._kwargs['loop_stop'] = node.stop
+
+        prefix = "\n".join(node.prefix_lines)
+        if prefix:
+            prefix += "\n"
+        return f"{prefix}{loop_str}\n{self._print(node.body)}"
+
+    def _print_SympyAssignment(self, node):
+        printed_lhs = self.sympy_printer.doprint(node.lhs)
+        printed_rhs = self.sympy_printer.doprint(node.rhs)
+
+        if node.is_declaration:
+            if node.use_auto:
+                data_type = 'auto'
+            else:
+                data_type = self._print(node.lhs.dtype).replace(' const', '')
+                if node.is_const:
+                    data_type = f'const {data_type}'
+            return f"{data_type} {printed_lhs} = {printed_rhs};"
+        else:
+            lhs_type = get_type_of_expression(node.lhs)  # TOOD: this should have been typed
+            printed_mask = ""
+            if type(lhs_type) is VectorType and isinstance(node.lhs, CastFunc):
+                arg, data_type, aligned, nontemporal, mask, stride = node.lhs.args
+                instr = 'storeU'
+                if nontemporal and 'storeA' not in self._vector_instruction_set and \
+                        'stream' in self._vector_instruction_set:
+                    instr = 'stream'
+                elif aligned:
+                    instr = 'stream' if nontemporal and 'stream' in self._vector_instruction_set else 'storeA'
+                if mask != True:  # NOQA
+                    instr = 'maskStream' if nontemporal and 'maskStream' in self._vector_instruction_set else \
+                            'maskStoreA' if aligned else 'maskStoreU'
+                    if instr not in self._vector_instruction_set:
+                        if instr == 'maskStream' and 'stream' in self._vector_instruction_set:
+                            store, load = 'stream', 'loadA'
+                        elif (instr in ('maskStream', 'maskStoreA')) and 'storeA' in self._vector_instruction_set:
+                            store, load = 'storeA', 'loadA'
+                        else:
+                            store, load = 'storeU', 'loadU'
+                        load = load if load in self._vector_instruction_set else 'loadU'
+                        self._vector_instruction_set[instr] = self._vector_instruction_set[store].format(
+                            '{0}', self._vector_instruction_set['blendv'].format(
+                                self._vector_instruction_set[load].format('{0}', **self._kwargs),
+                                '{1}', '{2}', **self._kwargs), **self._kwargs)
+                    printed_mask = self.sympy_printer.doprint(mask)
+                    if data_type.base_type.c_name == 'double':
+                        if self._vector_instruction_set['double'] == '__m256d':
+                            printed_mask = f"_mm256_castpd_si256({printed_mask})"
+                        elif self._vector_instruction_set['double'] == '__m128d':
+                            printed_mask = f"_mm_castpd_si128({printed_mask})"
+                    elif data_type.base_type.c_name == 'float':
+                        if self._vector_instruction_set['float'] == '__m256':
+                            printed_mask = f"_mm256_castps_si256({printed_mask})"
+                        elif self._vector_instruction_set['float'] == '__m128':
+                            printed_mask = f"_mm_castps_si128({printed_mask})"
+
+                rhs_type = get_type_of_expression(node.rhs)
+                if type(rhs_type) is not VectorType:
+                    raise ValueError(f'Cannot vectorize {node.rhs} of type {rhs_type} inside of the pretty printer! '
+                                     f'This should have happen earlier!')
+                    # rhs = CastFunc(node.rhs, VectorType(rhs_type)) # Unknown width
+                else:
+                    rhs = node.rhs
+
+                ptr = "&" + self.sympy_printer.doprint(node.lhs.args[0])
+
+                if stride != 1:
+                    instr = ('maskStreamS' if nontemporal and 'maskStreamS' in self._vector_instruction_set else
+                             'maskStoreS') if mask != True else \
+                            ('streamS' if nontemporal and 'streamS' in self._vector_instruction_set else 'storeS')  # NOQA
+                    return self._vector_instruction_set[instr].format(ptr, self.sympy_printer.doprint(rhs),
+                                                                      stride, printed_mask, **self._kwargs) + ';'
+
+                pre_code = ''
+                if nontemporal and 'cachelineZero' in self._vector_instruction_set and mask == True:  # NOQA
+                    first_cond = f"((uintptr_t) {ptr} & {CachelineSize.mask_symbol}) == 0"
+                    offset = sp.Add(*[sp.Symbol(LoopOverCoordinate.get_loop_counter_name(i))
+                                      * node.lhs.args[0].field.spatial_strides[i] for i in
+                                      range(len(node.lhs.args[0].field.spatial_strides))])
+                    if stride == 1:
+                        offset = offset.subs({node.lhs.args[0].field.spatial_strides[0]: 1})
+                    size = sp.Mul(*node.lhs.args[0].field.spatial_shape)
+                    element_size = 8 if data_type.base_type.c_name == 'double' else 4
+                    size_cond = f"({offset} + {CachelineSize.symbol/element_size}) < {size}"
+                    pre_code = f"if ({first_cond} && {size_cond}) " + "{\n\t" + \
+                        self._vector_instruction_set['cachelineZero'].format(ptr, **self._kwargs) + ';\n}\n'
+
+                code = self._vector_instruction_set[instr].format(ptr, self.sympy_printer.doprint(rhs),
+                                                                  printed_mask, **self._kwargs) + ';'
+                flushcond = f"((uintptr_t) {ptr} & {CachelineSize.mask_symbol}) == {CachelineSize.last_symbol}"
+                if nontemporal and 'flushCacheline' in self._vector_instruction_set:
+                    code2 = self._vector_instruction_set['flushCacheline'].format(
+                        ptr, self.sympy_printer.doprint(rhs), **self._kwargs) + ';'
+                    code = f"{code}\nif ({flushcond}) {{\n\t{code2}\n}}"
+                elif aligned and nontemporal and 'storeAAndFlushCacheline' in self._vector_instruction_set:
+                    lhs_hash = hashlib.sha1(self.sympy_printer.doprint(node.lhs).encode('ascii')).hexdigest()[:8]
+                    rhs_hash = hashlib.sha1(self.sympy_printer.doprint(rhs).encode('ascii')).hexdigest()[:8]
+                    tmpvar = f'_tmp_{lhs_hash}_{rhs_hash}'
+                    code = 'const ' + self._print(node.lhs.dtype).replace(' const', '') + ' ' + tmpvar + ' = ' \
+                        + self.sympy_printer.doprint(rhs) + ';'
+                    code1 = self._vector_instruction_set[instr].format(ptr, tmpvar, printed_mask, **self._kwargs) + ';'
+                    maskStore, store, load = 'maskStoreAAndFlushCacheline', 'storeAAndFlushCacheline', 'loadA'
+                    instr2 = maskStore if mask != True else store  # NOQA
+                    if instr2 not in self._vector_instruction_set:
+                        self._vector_instruction_set[maskStore] = self._vector_instruction_set[store].format(
+                            '{0}', self._vector_instruction_set['blendv'].format(
+                                self._vector_instruction_set[load].format('{0}', **self._kwargs),
+                                '{1}', '{2}', **self._kwargs),
+                            **self._kwargs)
+                    code2 = self._vector_instruction_set[instr2].format(ptr, tmpvar, printed_mask, **self._kwargs) + ';'
+                    code += f"\nif ({flushcond}) {{\n\t{code2}\n}} else {{\n\t{code1}\n}}"
+                return pre_code + code
+            else:
+                return f"{printed_lhs} = {printed_rhs};"
+
+    def _print_NontemporalFence(self, _):
+        if 'streamFence' in self._vector_instruction_set:
+            return self._vector_instruction_set['streamFence'] + ';'
+        else:
+            return ''
+
+    def _print_CachelineSize(self, node):
+        if 'cachelineSize' in self._vector_instruction_set:
+            code = f'const size_t {node.symbol} = {self._vector_instruction_set["cachelineSize"]};\n'
+            code += f'const size_t {node.mask_symbol} = {node.symbol} - 1;\n'
+            vectorsize = self._vector_instruction_set['bytes']
+            code += f'const size_t {node.last_symbol} = {node.symbol} - {vectorsize};\n'
+            return code
+        else:
+            return ''
+
+    def _print_TemporaryMemoryAllocation(self, node):
+        if self._vector_instruction_set:
+            align = self._vector_instruction_set['bytes']
+        else:
+            align = node.symbol.dtype.base_type.numpy_dtype.itemsize
+
+        np_dtype = node.symbol.dtype.base_type.numpy_dtype
+        required_size = np_dtype.itemsize * node.size + align
+        size = modulo_ceil(required_size, align)
+        code = "#if defined(_MSC_VER)\n"
+        code += "{dtype} {name}=({dtype})_aligned_malloc({size}, {align}) + {offset};\n"
+        code += "#elif __cplusplus >= 201703L || __STDC_VERSION__ >= 201112L\n"
+        code += "{dtype} {name}=({dtype})aligned_alloc({align}, {size}) + {offset};\n"
+        code += "#else\n"
+        code += "{dtype} {name};\n"
+        code += "posix_memalign((void**) &{name}, {align}, {size});\n"
+        code += "{name} += {offset};\n"
+        code += "#endif"
+        return code.format(dtype=node.symbol.dtype,
+                           name=self.sympy_printer.doprint(node.symbol.name),
+                           size=self.sympy_printer.doprint(size),
+                           offset=int(node.offset(align)),
+                           align=align)
+
+    def _print_TemporaryMemoryFree(self, node):
+        if self._vector_instruction_set:
+            align = self._vector_instruction_set['bytes']
+        else:
+            align = node.symbol.dtype.base_type.numpy_dtype.itemsize
+
+        code = "#if defined(_MSC_VER)\n"
+        code += "_aligned_free(%s - %d);\n" % (self.sympy_printer.doprint(node.symbol.name), node.offset(align))
+        code += "#else\n"
+        code += "free(%s - %d);\n" % (self.sympy_printer.doprint(node.symbol.name), node.offset(align))
+        code += "#endif"
+        return code
+
+    def _print_SkipIteration(self, _):
+        return "continue;"
+
+    def _print_CustomCodeNode(self, node):
+        return node.get_code(self._dialect, self._vector_instruction_set, print_arg=self.sympy_printer._print)
+
+    def _print_SourceCodeComment(self, node):
+        return f"/* {node.text } */"
+
+    def _print_EmptyLine(self, node):
+        return ""
+
+    def _print_Conditional(self, node):
+        if type(node.condition_expr) is BooleanTrue:
+            return self._print_Block(node.true_block)
+        elif type(node.condition_expr) is BooleanFalse:
+            return self._print_Block(node.false_block)
+        cond_type = get_type_of_expression(node.condition_expr)
+        if isinstance(cond_type, VectorType):
+            raise ValueError("Problem with Conditional inside vectorized loop - use vec_any or vec_all")
+        condition_expr = self.sympy_printer.doprint(node.condition_expr)
+        true_block = self._print_Block(node.true_block)
+        result = f"if ({condition_expr})\n{true_block} "
+        if node.false_block:
+            false_block = self._print_Block(node.false_block)
+            result += f"else {false_block}"
+        return result
+
+
+# ------------------------------------------ Helper function & classes -------------------------------------------------
+
+
+# noinspection PyPep8Naming
+class CustomSympyPrinter(CCodePrinter):
+
+    def __init__(self):
+        super(CustomSympyPrinter, self).__init__()
+
+    def _print_Pow(self, expr):
+        """Don't use std::pow function, for small integer exponents, write as multiplication"""
+        # Ideally the printer has as little logic as possible. Therefore,
+        # powers should be rewritten as `DivFunc`s / unevaluated `Mul`s before
+        # printing. `NodeCollection` offers a convenience function to do just
+        # that. However, `cut_loops` rewrites unevaluated multiplications as
+        # `Pow`s again. Neither `deepcopy` nor `func(*args)` are suited to
+        # rebuild unevaluated expressions. Therefore, as long as we stick with
+        # SymPy, this is the only way to avoid printing `pow`s.
+        exp = expr.exp.expr if isinstance(expr.exp, CastFunc) else expr.exp
+        one_type = expr.base.dtype if hasattr(expr.base, "dtype") else get_type_of_expression(expr.base)
+
+        if exp.is_integer and exp.is_number and (0 < exp <= 8):
+            return f"({self._print(sp.Mul(*[expr.base] * exp, evaluate=False))})"
+        elif exp.is_integer and exp.is_number and (-8 <= exp < 0):
+            return f"{self._typed_number(1, one_type)} / ({self._print(sp.Mul(*([expr.base] * -exp), evaluate=False))})"
+        else:
+            return super(CustomSympyPrinter, self)._print_Pow(expr)
+
+    # TODO don't print ones in sp.Mul
+
+    def _print_Rational(self, expr):
+        """Evaluate all rationals i.e. print 0.25 instead of 1.0/4.0"""
+        res = str(expr.evalf(17))
+        return res
+
+    def _print_Equality(self, expr):
+        """Equality operator is not printable in default printer"""
+        return '((' + self._print(expr.lhs) + ") == (" + self._print(expr.rhs) + '))'
+
+    def _print_Piecewise(self, expr):
+        """Print piecewise in one line (remove newlines)"""
+        result = super(CustomSympyPrinter, self)._print_Piecewise(expr)
+        return result.replace("\n", "")
+
+    def _print_Abs(self, expr):
+        if expr.args[0].is_integer:
+            return f'abs({self._print(expr.args[0])})'
+        else:
+            return f'fabs({self._print(expr.args[0])})'
+
+    def _print_AbstractType(self, node):
+        return str(node)
+
+    def _print_Function(self, expr):
+        infix_functions = {
+            bitwise_xor: '^',
+            bit_shift_right: '>>',
+            bit_shift_left: '<<',
+            bitwise_or: '|',
+            bitwise_and: '&',
+        }
+        if hasattr(expr, 'to_c'):
+            return expr.to_c(self._print)
+        if isinstance(expr, ReinterpretCastFunc):
+            arg, data_type = expr.args
+            if isinstance(data_type, PointerType):
+                const_str = "const" if data_type.const else ""
+                return f"(({const_str} {self._print(data_type.base_type)} *)(& {self._print(arg)}))"
+            else:
+                return f"*(({self._print(PointerType(data_type, restrict=False))})(& {self._print(arg)}))"
+        elif isinstance(expr, AddressOf):
+            assert len(expr.args) == 1, "address_of must only have one argument"
+            return f"&({self._print(expr.args[0])})"
+        elif isinstance(expr, CastFunc):
+            cast = "(({data_type})({code}))"
+            arg, data_type = expr.args
+            if arg.is_Number and not isinstance(arg, (sp.core.numbers.Infinity, sp.core.numbers.NegativeInfinity)):
+                return self._typed_number(arg, data_type)
+            elif isinstance(arg, (InverseTrigonometricFunction, TrigonometricFunction, HyperbolicFunction)) \
+                    and data_type == BasicType('float32'):
+                known = self.known_functions[arg.__class__.__name__.lower()]
+                code = self._print(arg)
+                return code.replace(known, f"{known}f")
+            elif isinstance(arg, (sp.Pow, sp.exp)) and data_type == BasicType('float32'):
+                known = ['sqrt', 'cbrt', 'pow', 'exp']
+                code = self._print(arg)
+                for k in known:
+                    if k in code:
+                        return code.replace(k, f'{k}f')
+                # Powers of small integers are printed as divisions/multiplications.
+                if '/' in code or '*' in code:
+                    return cast.format(data_type=data_type, code=code)
+                raise ValueError(f"{code} doesn't give {known=} function back.")
+            else:
+                return cast.format(data_type=data_type, code=self._print(arg))
+        elif isinstance(expr, fast_division):
+            raise ValueError("fast_division is only supported for Taget.GPU")
+        elif isinstance(expr, fast_sqrt):
+            raise ValueError("fast_sqrt is only supported for Taget.GPU")
+        elif isinstance(expr, fast_inv_sqrt):
+            raise ValueError("fast_inv_sqrt is only supported for Taget.GPU")
+        elif isinstance(expr, vec_any) or isinstance(expr, vec_all):
+            return self._print(expr.args[0])
+        elif isinstance(expr, sp.Abs):
+            return f"abs({self._print(expr.args[0])})"
+        elif isinstance(expr, sp.Mod):
+            if expr.args[0].is_integer and expr.args[1].is_integer:
+                return f"({self._print(expr.args[0])} % {self._print(expr.args[1])})"
+            else:
+                return f"fmod({self._print(expr.args[0])}, {self._print(expr.args[1])})"
+        elif expr.func in infix_functions:
+            return f"({self._print(expr.args[0])} {infix_functions[expr.func]} {self._print(expr.args[1])})"
+        elif expr.func == int_power_of_2:
+            return f"(1 << ({self._print(expr.args[0])}))"
+        elif expr.func == int_div:
+            return f"(({self._print(expr.args[0])}) / ({self._print(expr.args[1])}))"
+        elif expr.func == DivFunc:
+            return f'(({self._print(expr.divisor)}) / ({self._print(expr.dividend)}))'
+        else:
+            name = expr.name if hasattr(expr, 'name') else expr.__class__.__name__
+            arg_str = ', '.join(self._print(a) for a in expr.args)
+            return f'{name}({arg_str})'
+
+    def _typed_number(self, number, dtype):
+        res = self._print(number)
+        if dtype.numpy_dtype == np.float32:
+            return res + '.0f' if '.' not in res else res + 'f'
+        elif dtype.numpy_dtype == np.float64:
+            return res + '.0' if '.' not in res else res
+        elif dtype.is_int():
+            tokens = res.split('.')
+            if len(tokens) == 1: 
+                return res
+            elif int(tokens[1]) != 0:
+                raise ValueError(f"Cannot print non-integer number {res} as an integer.")
+            else:
+                return tokens[0]
+        else:
+            return res
+
+    def _print_ConditionalFieldAccess(self, node):
+        return self._print(sp.Piecewise((node.outofbounds_value, node.outofbounds_condition), (node.access, True)))
+
+    def _print_Max(self, expr):
+        def inner_print_max(args):
+            if len(args) == 1:
+                return self._print(args[0])
+            half = len(args) // 2
+            a = inner_print_max(args[:half])
+            b = inner_print_max(args[half:])
+            return f"(({a} > {b}) ? {a} : {b})"
+        return inner_print_max(expr.args)
+
+    def _print_Min(self, expr):
+        def inner_print_min(args):
+            if len(args) == 1:
+                return self._print(args[0])
+            half = len(args) // 2
+            a = inner_print_min(args[:half])
+            b = inner_print_min(args[half:])
+            return f"(({a} < {b}) ? {a} : {b})"
+        return inner_print_min(expr.args)
+
+
+# noinspection PyPep8Naming
+class VectorizedCustomSympyPrinter(CustomSympyPrinter):
+    SummandInfo = namedtuple("SummandInfo", ['sign', 'term'])
+
+    def __init__(self, instruction_set):
+        super(VectorizedCustomSympyPrinter, self).__init__()
+        self.instruction_set = instruction_set
+
+    def _scalarFallback(self, func_name, expr, *args, **kwargs):
+        expr_type = get_type_of_expression(expr)
+        if type(expr_type) is not VectorType:
+            return getattr(super(VectorizedCustomSympyPrinter, self), func_name)(expr, *args, **kwargs)
+        else:
+            assert self.instruction_set['width'] == expr_type.width
+            return None
+
+    def _print_Abs(self, expr):
+        if isinstance(get_type_of_expression(expr), (VectorType, VectorMemoryAccess)):
+            return self.instruction_set['abs'].format(self._print(expr.args[0]), **self._kwargs)
+        return super()._print_Abs(expr)
+
+    def _typed_vectorized_number(self, expr, data_type):
+        basic_data_type = data_type.base_type
+        number = self._typed_number(expr, basic_data_type)
+        instruction = 'makeVecConst'
+        if basic_data_type.is_bool():
+            instruction = 'makeVecConstBool'
+        # TODO Vectorization Revamp: is int, or sint, or uint (my guess is sint)
+        elif basic_data_type.is_int():
+            instruction = 'makeVecConstInt'
+        return self.instruction_set[instruction].format(number, **self._kwargs)
+
+    def _typed_vectorized_symbol(self, expr, data_type):
+        if not isinstance(expr, TypedSymbol):
+            raise ValueError(f'{expr} is not a TypeSymbol. It is {expr.type=}')
+        basic_data_type = data_type.base_type
+        symbol = self._print(expr)
+        if basic_data_type != expr.dtype:
+            symbol = f'(({basic_data_type})({symbol}))'
+
+        instruction = 'makeVecConst'
+        if basic_data_type.is_bool():
+            instruction = 'makeVecConstBool'
+        # TODO Vectorization Revamp: is int, or sint, or uint (my guess is sint)
+        elif basic_data_type.is_int():
+            instruction = 'makeVecConstInt'
+        return self.instruction_set[instruction].format(symbol, **self._kwargs)
+
+    def _print_CastFunc(self, expr):
+        arg, data_type = expr.args
+        if type(data_type) is VectorType:
+            base_type = data_type.base_type
+            # vector_memory_access is a cast_func itself so it should't be directly inside a cast_func
+            assert not isinstance(arg, VectorMemoryAccess)
+            if isinstance(arg, sp.Tuple):
+                is_boolean = get_type_of_expression(arg[0]) == create_type("bool")
+                is_integer = get_type_of_expression(arg[0]) == create_type("int")
+                printed_args = [self._print(a) for a in arg]
+                instruction = 'makeVecBool' if is_boolean else 'makeVecInt' if is_integer else 'makeVec'
+                if instruction == 'makeVecInt' and 'makeVecIndex' in self.instruction_set:
+                    increments = np.array(arg)[1:] - np.array(arg)[:-1]
+                    if len(set(increments)) == 1:
+                        return self.instruction_set['makeVecIndex'].format(printed_args[0], increments[0],
+                                                                           **self._kwargs)
+                return self.instruction_set[instruction].format(*printed_args, **self._kwargs)
+            else:
+                if arg.is_Number and not isinstance(arg, (sp.core.numbers.Infinity, sp.core.numbers.NegativeInfinity)):
+                    return self._typed_vectorized_number(arg, data_type)
+                elif isinstance(arg, TypedSymbol):
+                    return self._typed_vectorized_symbol(arg, data_type)
+                elif isinstance(arg, (InverseTrigonometricFunction, TrigonometricFunction, HyperbolicFunction)) \
+                        and base_type == BasicType('float32'):
+                    raise NotImplementedError('Vectorizer is not tested for trigonometric functions yet')
+                    # known = self.known_functions[arg.__class__.__name__.lower()]
+                    # code = self._print(arg)
+                    # return code.replace(known, f"{known}f")
+                elif isinstance(arg, sp.Pow):
+                    if base_type == BasicType('float32') or base_type == BasicType('float64'):
+                        return self._print_Pow(arg)
+                    else:
+                        raise NotImplementedError('Integer Pow is not implemented')
+                elif isinstance(arg, sp.UnevaluatedExpr):
+                    return self._print(arg.args[0])
+                else:
+                    raise NotImplementedError('Vectorizer cannot cast between different datatypes')
+                    # to_type = self.instruction_set['suffix'][data_type.base_type.c_name]
+                    # from_type = self.instruction_set['suffix'][get_type_of_expression(arg).base_type.c_name]
+                    # return self.instruction_set['cast'].format(from_type, to_type, self._print(arg))
+        else:
+            return self._scalarFallback('_print_Function', expr)
+            # raise ValueError(f'Non VectorType cast "{data_type}" in vectorized code.')
+
+    def _print_Function(self, expr):
+        if isinstance(expr, VectorMemoryAccess):
+            arg, data_type, aligned, _, mask, stride = expr.args
+            if stride != 1:
+                return self.instruction_set['loadS'].format(f"& {self._print(arg)}", stride, **self._kwargs)
+            instruction = self.instruction_set['loadA'] if aligned else self.instruction_set['loadU']
+            return instruction.format(f"& {self._print(arg)}", **self._kwargs)
+        elif expr.func == DivFunc:
+            result = self._scalarFallback('_print_Function', expr)
+            if not result:
+                result = self.instruction_set['/'].format(self._print(expr.divisor), self._print(expr.dividend),
+                                                          **self._kwargs)
+            return result
+        elif isinstance(expr, fast_division):
+            raise ValueError("fast_division is only supported for Taget.GPU")
+        elif isinstance(expr, fast_sqrt):
+            raise ValueError("fast_sqrt is only supported for Taget.GPU")
+        elif isinstance(expr, fast_inv_sqrt):
+            raise ValueError("fast_inv_sqrt is only supported for Taget.GPU")
+        elif isinstance(expr, vec_any) or isinstance(expr, vec_all):
+            instr = 'any' if isinstance(expr, vec_any) else 'all'
+            expr_type = get_type_of_expression(expr.args[0])
+            if type(expr_type) is not VectorType:
+                return self._print(expr.args[0])
+            else:
+                if isinstance(expr.args[0], sp.Rel):
+                    op = expr.args[0].rel_op
+                    if (instr, op) in self.instruction_set:
+                        return self.instruction_set[(instr, op)].format(*[self._print(a) for a in expr.args[0].args],
+                                                                        **self._kwargs)
+                return self.instruction_set[instr].format(self._print(expr.args[0]), **self._kwargs)
+
+        return super(VectorizedCustomSympyPrinter, self)._print_Function(expr)
+
+    def _print_And(self, expr):
+        result = self._scalarFallback('_print_And', expr)
+        if result:
+            return result
+
+        arg_strings = [self._print(a) for a in expr.args]
+        assert len(arg_strings) > 0
+        result = arg_strings[0]
+        for item in arg_strings[1:]:
+            result = self.instruction_set['&'].format(result, item, **self._kwargs)
+        return result
+
+    def _print_Or(self, expr):
+        result = self._scalarFallback('_print_Or', expr)
+        if result:
+            return result
+
+        arg_strings = [self._print(a) for a in expr.args]
+        assert len(arg_strings) > 0
+        result = arg_strings[0]
+        for item in arg_strings[1:]:
+            result = self.instruction_set['|'].format(result, item, **self._kwargs)
+        return result
+
+    def _print_Add(self, expr, order=None):
+        try:
+            result = self._scalarFallback('_print_Add', expr)
+        except Exception:
+            result = None
+        if result:
+            return result
+        args = expr.args
+
+        # special treatment for all-integer args, for loop index arithmetic until we have proper int vectorization
+        suffix = ""
+        if all([(type(e) is CastFunc and str(e.dtype) == self.instruction_set['int']) or isinstance(e, sp.Integer)
+                or (type(e) is TypedSymbol and isinstance(e.dtype, BasicType) and e.dtype.is_int()) for e in args]):
+            dtype = set([e.dtype for e in args if type(e) is CastFunc])
+            assert len(dtype) == 1
+            dtype = dtype.pop()
+            args = [CastFunc(e, dtype) if (isinstance(e, sp.Integer) or isinstance(e, TypedSymbol)) else e
+                    for e in args]
+            suffix = "int"
+
+        summands = []
+        for term in args:
+            if term.func == sp.Mul:
+                sign, t = self._print_Mul(term, inside_add=True)
+            else:
+                t = self._print(term)
+                sign = 1
+            summands.append(self.SummandInfo(sign, t))
+        # Use positive terms first
+        summands.sort(key=lambda e: e.sign, reverse=True)
+        # if no positive term exists, prepend a zero
+        if summands[0].sign == -1:
+            summands.insert(0, self.SummandInfo(1, "0"))
+
+        assert len(summands) >= 2
+        processed = summands[0].term
+        for summand in summands[1:]:
+            func = self.instruction_set['-' + suffix] if summand.sign == -1 else self.instruction_set['+' + suffix]
+            processed = func.format(processed, summand.term, **self._kwargs)
+        return processed
+
+    def _print_Pow(self, expr):
+        # Due to loop cutting sp.Mul is evaluated again.
+
+        try:
+            result = self._scalarFallback('_print_Pow', expr)
+        except ValueError:
+            result = None
+        if result:
+            return result
+
+        one = self.instruction_set['makeVecConst'].format(1.0, **self._kwargs)
+        root = self.instruction_set['sqrt'].format(self._print(expr.base), **self._kwargs)
+
+        if isinstance(expr.exp, CastFunc) and expr.exp.args[0].is_number:
+            exp = expr.exp.args[0]
+        else:
+            exp = expr.exp
+
+        # TODO the printer should not have any intelligence like this.
+        # TODO To remove all of these cases the vectoriser needs to be reworked. See loop cutting
+        if exp.is_integer and exp.is_number and 0 < exp < 8:
+            return self._print(sp.Mul(*[expr.base] * exp, evaluate=False))
+        elif exp == 0.5:
+            return root
+        elif exp == -0.5:
+            return self.instruction_set['/'].format(one, root, **self._kwargs)
+        else:
+            raise ValueError("Generic exponential not supported: " + str(expr))
+
+    def _print_Mul(self, expr, inside_add=False):
+        # noinspection PyProtectedMember
+        from sympy.core.mul import _keep_coeff
+
+        if not inside_add:
+            result = self._scalarFallback('_print_Mul', expr)
+        else:
+            result = None
+        if result:
+            return result
+
+        c, e = expr.as_coeff_Mul()
+        if c < 0:
+            expr = _keep_coeff(-c, e)
+            sign = -1
+        else:
+            sign = 1
+
+        a = []  # items in the numerator
+        b = []  # items that are in the denominator (if any)
+
+        # Gather args for numerator/denominator
+        for item in expr.as_ordered_factors():
+            if item.is_commutative and item.is_Pow and item.exp.is_Rational and item.exp.is_negative:
+                if item.exp != -1:
+                    b.append(sp.Pow(item.base, -item.exp, evaluate=False))
+                else:
+                    b.append(sp.Pow(item.base, -item.exp))
+            else:
+                a.append(item)
+
+        a = a or [S.One]
+
+        a_str = [self._print(x) for x in a]
+        b_str = [self._print(x) for x in b]
+
+        result = a_str[0]
+        for item in a_str[1:]:
+            result = self.instruction_set['*'].format(result, item, **self._kwargs)
+
+        if len(b) > 0:
+            denominator_str = b_str[0]
+            for item in b_str[1:]:
+                denominator_str = self.instruction_set['*'].format(denominator_str, item, **self._kwargs)
+            result = self.instruction_set['/'].format(result, denominator_str, **self._kwargs)
+
+        if inside_add:
+            return sign, result
+        else:
+            if sign < 0:
+                return self.instruction_set['*'].format(self._print(S.NegativeOne), result, **self._kwargs)
+            else:
+                return result
+
+    def _print_Relational(self, expr):
+        result = self._scalarFallback('_print_Relational', expr)
+        if result:
+            return result
+        return self.instruction_set[expr.rel_op].format(self._print(expr.lhs), self._print(expr.rhs), **self._kwargs)
+
+    def _print_Equality(self, expr):
+        result = self._scalarFallback('_print_Equality', expr)
+        if result:
+            return result
+        return self.instruction_set['=='].format(self._print(expr.lhs), self._print(expr.rhs), **self._kwargs)
+
+    def _print_Piecewise(self, expr):
+        result = self._scalarFallback('_print_Piecewise', expr)
+        if result:
+            return result
+
+        if expr.args[-1].cond.args[0] is not sp.sympify(True):
+            # We need the last conditional to be a True, otherwise the resulting
+            # function may not return a result.
+            raise ValueError("All Piecewise expressions must contain an "
+                             "(expr, True) statement to be used as a default "
+                             "condition. Without one, the generated "
+                             "expression may not evaluate to anything under "
+                             "some condition.")
+
+        result = self._print(expr.args[-1][0])
+        for true_expr, condition in reversed(expr.args[:-1]):
+            if isinstance(condition, CastFunc) and get_type_of_expression(condition.args[0]) == create_type("bool"):
+                result = "(({}) ? ({}) : ({}))".format(self._print(condition.args[0]), self._print(true_expr),
+                                                       result, **self._kwargs)
+            else:
+                # noinspection SpellCheckingInspection
+                result = self.instruction_set['blendv'].format(result, self._print(true_expr), self._print(condition),
+                                                               **self._kwargs)
+        return result

src/pystencils/backends/cuda_backend.py

+from pystencils.astnodes import Node
+from pystencils.backends.cbackend import CBackend, CustomSympyPrinter, generate_c
+from pystencils.enums import Backend
+from pystencils.fast_approximation import fast_division, fast_inv_sqrt, fast_sqrt
+
+
+def generate_cuda(ast_node: Node, signature_only: bool = False, custom_backend=None, with_globals=True) -> str:
+    """Prints an abstract syntax tree node as CUDA code.
+
+    Args:
+        ast_node: ast representation of kernel
+        signature_only: generate signature without function body
+        custom_backend: use own custom printer for code generation
+        with_globals: enable usage of global variables
+
+    Returns:
+        CUDA code for the ast node and its descendants
+    """
+    return generate_c(ast_node, signature_only, dialect=Backend.CUDA,
+                      custom_backend=custom_backend, with_globals=with_globals)
+
+
+class CudaBackend(CBackend):
+
+    def __init__(self, sympy_printer=None,
+                 signature_only=False):
+        if not sympy_printer:
+            sympy_printer = CudaSympyPrinter()
+
+        super().__init__(sympy_printer, signature_only, dialect=Backend.CUDA)
+
+    def _print_SharedMemoryAllocation(self, node):
+        dtype = node.symbol.dtype
+        name = self.sympy_printer.doprint(node.symbol.name)
+        num_elements = '*'.join([str(s) for s in node.shared_mem.shape])
+        code = f"__shared__ {dtype} {name}[{num_elements}];"
+        return code
+
+    @staticmethod
+    def _print_ThreadBlockSynchronization(_):
+        return "__synchtreads();"
+
+    def _print_TextureDeclaration(self, node):
+        cond = node.texture.field.dtype.numpy_dtype.itemsize > 4
+        return f'texture<{"fp_tex_" if cond else ""}{str(node.texture.field.dtype)}, ' \
+               f'cudaTextureType{node.texture.field.spacial_dimensions}D, cudaReadModeElementType> {node.texture};'
+
+    def _print_SkipIteration(self, _):
+        return "return;"
+
+
+class CudaSympyPrinter(CustomSympyPrinter):
+    language = "CUDA"
+
+    def __init__(self):
+        super(CudaSympyPrinter, self).__init__()
+
+    def _print_Function(self, expr):
+        if isinstance(expr, fast_division):
+            assert len(expr.args) == 2, f"__fdividef has two arguments, but {len(expr.args)} where given"
+            return f"__fdividef({self._print(expr.args[0])}, {self._print(expr.args[1])})"
+        elif isinstance(expr, fast_sqrt):
+            assert len(expr.args) == 1, f"__fsqrt_rn has one argument, but {len(expr.args)} where given"
+            return f"__fsqrt_rn({self._print(expr.args[0])})"
+        elif isinstance(expr, fast_inv_sqrt):
+            assert len(expr.args) == 1, f"__frsqrt_rn has one argument, but {len(expr.args)} where given"
+            return f"__frsqrt_rn({self._print(expr.args[0])})"
+        return super()._print_Function(expr)

src/pystencils/backends/dot.py

+import graphviz
+try:
+    from graphviz import Digraph
+    import graphviz.quoting as quote
+except ImportError:
+    from graphviz import Digraph
+    import graphviz.lang as quote
+from sympy.printing.printer import Printer
+
+
+# noinspection PyPep8Naming
+class DotPrinter(Printer):
+    """
+    A printer which converts ast to DOT (graph description language).
+    """
+    def __init__(self, node_to_str_function, **kwargs):
+        super(DotPrinter, self).__init__()
+        self._node_to_str_function = node_to_str_function
+        self.dot = Digraph(**kwargs)
+        self.dot.quote_edge = quote.quote
+
+    def _print_KernelFunction(self, func):
+        self.dot.node(str(id(func)), style='filled', fillcolor='#a056db', label=self._node_to_str_function(func))
+        self._print(func.body)
+        self.dot.edge(str(id(func)), str(id(func.body)))
+
+    def _print_LoopOverCoordinate(self, loop):
+        self.dot.node(str(id(loop)), style='filled', fillcolor='#3498db', label=self._node_to_str_function(loop))
+        self._print(loop.body)
+        self.dot.edge(str(id(loop)), str(id(loop.body)))
+
+    def _print_Block(self, block):
+        for node in block.args:
+            self._print(node)
+
+        self.dot.node(str(id(block)), style='filled', fillcolor='#dbc256', label=repr(block))
+        for node in block.args:
+            self.dot.edge(str(id(block)), str(id(node)))
+
+    def _print_SympyAssignment(self, assignment):
+        self.dot.node(str(id(assignment)), style='filled', fillcolor='#56db7f',
+                      label=self._node_to_str_function(assignment))
+
+    def _print_Conditional(self, expr):
+        self.dot.node(str(id(expr)), style='filled', fillcolor='#56bd7f', label=self._node_to_str_function(expr))
+        self._print(expr.true_block)
+        self.dot.edge(str(id(expr)), str(id(expr.true_block)))
+        if expr.false_block:
+            self._print(expr.false_block)
+            self.dot.edge(str(id(expr)), str(id(expr.false_block)))
+
+    def doprint(self, expr):
+        self._print(expr)
+        return self.dot.source
+
+
+def __shortened(node):
+    from pystencils.astnodes import LoopOverCoordinate, KernelFunction, SympyAssignment, Conditional
+    if isinstance(node, LoopOverCoordinate):
+        return "Loop over dim %d" % (node.coordinate_to_loop_over,)
+    elif isinstance(node, KernelFunction):
+        params = node.get_parameters()
+        param_names = [p.field_name for p in params if p.is_field_pointer]
+        param_names += [p.symbol.name for p in params if not p.is_field_parameter]
+        return f"Func: {node.function_name} ({','.join(param_names)})"
+    elif isinstance(node, SympyAssignment):
+        return repr(node.lhs)
+    elif isinstance(node, Conditional):
+        return repr(node)
+    else:
+        raise NotImplementedError(f"Cannot handle node type {type(node)}")
+
+
+def print_dot(node, view=False, short=False, **kwargs):
+    """
+    Returns a string which can be used to generate a DOT-graph
+    :param node: The ast which should be generated
+    :param view: Boolean, if rendering of the image directly should occur.
+    :param short: Uses the __shortened output
+    :param kwargs: is directly passed to the DotPrinter class: http://graphviz.readthedocs.io/en/latest/api.html#digraph
+    :return: string in DOT format
+    """
+    node_to_str_function = repr
+    if short:
+        node_to_str_function = __shortened
+    printer = DotPrinter(node_to_str_function, **kwargs)
+    dot = printer.doprint(node)
+    if view:
+        return graphviz.Source(dot)
+    return dot

src/pystencils/backends/json.py

+# -*- coding: utf-8 -*-
+#
+# Copyright © 2019 Stephan Seitz <stephan.seitz@fau.de>
+#
+# Distributed under terms of the GPLv3 license.
+
+"""
+
+"""
+import json
+
+from pystencils.astnodes import NodeOrExpr
+from pystencils.backends.cbackend import CustomSympyPrinter, generate_c
+
+try:
+    import yaml
+except ImportError:
+    raise ImportError('yaml not installed')
+
+
+def expr_to_dict(expr_or_node: NodeOrExpr, with_c_code=True, full_class_names=False):
+    """Converts a SymPy expression to a serializable dict (mainly for debugging purposes)
+
+    The dict recursively contains all args of the expression as ``dict``s
+
+    See :func:`.write_json`
+
+    Args:
+        expr_or_node (NodeOrExpr): a SymPy expression or a :class:`pystencils.astnodes.Node`
+        with_c_code (bool, optional): include C representation of the nodes
+        full_class_names (bool, optional): use full class names (type(object) instead of ``type(object).__name__``
+    """
+
+    self = {'str': str(expr_or_node)}
+    if with_c_code:
+        try:
+            self.update({'c': generate_c(expr_or_node)})
+        except Exception:
+            try:
+                self.update({'c': CustomSympyPrinter().doprint(expr_or_node)})
+            except Exception:
+                pass
+    for a in expr_or_node.args:
+        self.update({str(a.__class__ if full_class_names else a.__class__.__name__): expr_to_dict(a)})
+
+    return self
+
+
+def print_json(expr_or_node: NodeOrExpr):
+    """Print debug JSON of an AST to string
+
+    Args:
+        expr_or_node (NodeOrExpr): a SymPy expression or a :class:`pystencils.astnodes.Node`
+
+    Returns:
+        str: JSON representation of AST
+    """
+    expr_or_node_dict = expr_to_dict(expr_or_node)
+    return json.dumps(expr_or_node_dict, indent=4)
+
+
+def write_json(filename: str, expr_or_node: NodeOrExpr):
+    """Writes debug JSON represenation of AST to file
+
+    Args:
+        filename (str): Path for the file to write
+        expr_or_node (NodeOrExpr): a SymPy expression or a :class:`pystencils.astnodes.Node`
+    """
+    expr_or_node_dict = expr_to_dict(expr_or_node)
+    with open(filename, 'w') as f:
+        json.dump(expr_or_node_dict, f, indent=4)
+
+
+def print_yaml(expr_or_node):
+    expr_or_node_dict = expr_to_dict(expr_or_node, full_class_names=False)
+    return yaml.dump(expr_or_node_dict)
+
+
+def write_yaml(filename, expr_or_node):
+    expr_or_node_dict = expr_to_dict(expr_or_node)
+    with open(filename, 'w') as f:
+        yaml.dump(expr_or_node_dict, f)

src/pystencils/backends/ppc_instruction_sets.py

+def get_argument_string(function_shortcut):
+    args = function_shortcut[function_shortcut.index('[') + 1: -1]
+    arg_string = "("
+    for arg in args.split(","):
+        arg = arg.strip()
+        if not arg:
+            continue
+        if arg in ('0', '1', '2', '3', '4', '5'):
+            arg_string += "{" + arg + "},"
+        else:
+            arg_string += arg + ","
+    arg_string = arg_string[:-1] + ")"
+    return arg_string
+
+
+def get_vector_instruction_set_ppc(data_type='double', instruction_set='vsx'):
+    if instruction_set != 'vsx':
+        raise NotImplementedError(instruction_set)
+
+    base_names = {
+        '+': 'add[0, 1]',
+        '-': 'sub[0, 1]',
+        '*': 'mul[0, 1]',
+        '/': 'div[0, 1]',
+        'sqrt': 'sqrt[0]',
+        'rsqrt': 'rsqrte[0]',  # rsqrt is available too, but not on Clang
+
+        'loadU': 'xl[0x0, 0]',
+        'loadA': 'ld[0x0, 0]',
+        'storeU': 'xst[1, 0x0, 0]',
+        'storeA': 'st[1, 0x0, 0]',
+        'storeAAndFlushCacheline': 'stl[1, 0x0, 0]',
+
+        'abs': 'abs[0]',
+        '==': 'cmpeq[0, 1]',
+        '!=': 'cmpne[0, 1]',
+        '<=': 'cmple[0, 1]',
+        '<': 'cmplt[0, 1]',
+        '>=': 'cmpge[0, 1]',
+        '>': 'cmpgt[0, 1]',
+        '&': 'and[0, 1]',
+        '|': 'or[0, 1]',
+        'blendv': 'sel[0, 1, 2]',
+
+        ('any', '=='): 'any_eq[0, 1]',
+        ('any', '!='): 'any_ne[0, 1]',
+        ('any', '<='): 'any_le[0, 1]',
+        ('any', '<'): 'any_lt[0, 1]',
+        ('any', '>='): 'any_ge[0, 1]',
+        ('any', '>'): 'any_gt[0, 1]',
+        ('all', '=='): 'all_eq[0, 1]',
+        ('all', '!='): 'all_ne[0, 1]',
+        ('all', '<='): 'all_le[0, 1]',
+        ('all', '<'): 'all_lt[0, 1]',
+        ('all', '>='): 'all_ge[0, 1]',
+        ('all', '>'): 'all_gt[0, 1]',
+    }
+
+    bits = {'double': 64,
+            'float': 32,
+            'int': 32}
+
+    width = 128 // bits[data_type]
+    intwidth = 128 // bits['int']
+
+    result = dict()
+    result['bytes'] = 16
+
+    for intrinsic_id, function_shortcut in base_names.items():
+        function_shortcut = function_shortcut.strip()
+        name = function_shortcut[:function_shortcut.index('[')]
+
+        arg_string = get_argument_string(function_shortcut)
+
+        result[intrinsic_id] = 'vec_' + name + arg_string
+
+    if data_type == 'double':
+        # Clang and XL C++ are missing these for doubles
+        result['loadA'] = '(__vector double)' + result['loadA'].format('(float*) {0}')
+        result['storeA'] = result['storeA'].format('(float*) {0}', '(__vector float) {1}')
+        result['storeAAndFlushCacheline'] = result['storeAAndFlushCacheline'].format('(float*) {0}',
+                                                                                     '(__vector float) {1}')
+
+    result['+int'] = "vec_add({0}, {1})"
+
+    result['width'] = width
+    result['intwidth'] = intwidth
+    result[data_type] = f'__vector {data_type}'
+    result['int'] = '__vector int'
+    result['bool'] = f'__vector __bool {"long long" if data_type == "double" else "int"}'
+    result['headers'] = ['<altivec.h>', '"ppc_altivec_helpers.h"']
+
+    result['makeVecConst'] = '((' + result[data_type] + '){{' + \
+        ", ".join(['(' + data_type + ') {0}' for _ in range(width)]) + '}})'
+    result['makeVec'] = '((' + result[data_type] + '){{' + \
+        ", ".join(['{' + data_type + '} {' + str(i) + '}' for i in range(width)]) + '}})'
+    result['makeVecConstInt'] = '((' + result['int'] + '){{' + ", ".join(['(int) {0}' for _ in range(intwidth)]) + '}})'
+    result['makeVecInt'] = '((' + result['int'] + '){{(int) {0}, (int) {1}, (int) {2}, (int) {3}}})'
+
+    result['any'] = 'vec_any_ne({0}, ((' + result['bool'] + ') {{' + ", ".join(['0'] * width) + '}}))'
+    result['all'] = 'vec_all_ne({0}, ((' + result['bool'] + ') {{' + ", ".join(['0'] * width) + '}}))'
+
+    result['cachelineSize'] = 'cachelineSize()'
+    result['cachelineZero'] = 'cachelineZero((void*) {0})'
+
+    return result

src/pystencils/backends/riscv_instruction_sets.py

+from pystencils.typing import CFunction
+
+
+def get_argument_string(function_shortcut, last=''):
+    args = function_shortcut[function_shortcut.index('[') + 1: -1]
+    arg_string = "("
+    for arg in args.split(","):
+        arg = arg.strip()
+        if not arg:
+            continue
+        if arg in ('0', '1', '2', '3', '4', '5'):
+            arg_string += "{" + arg + "},"
+        else:
+            arg_string += arg + ","
+    if last:
+        arg_string += last + ','
+    arg_string = arg_string[:-1] + ")"
+    return arg_string
+
+
+def get_vector_instruction_set_riscv(data_type='double', instruction_set='rvv'):
+    assert instruction_set == 'rvv'
+
+    bits = {'double': 64,
+            'float': 32,
+            'int': 32}
+
+    base_names = {
+        '+': 'fadd_vv[0, 1]',
+        '-': 'fsub_vv[0, 1]',
+        '*': 'fmul_vv[0, 1]',
+        '/': 'fdiv_vv[0, 1]',
+        'sqrt': 'fsqrt_v[0]',
+
+        'loadU': f'le{bits[data_type]}_v[0]',
+        'storeU': f'se{bits[data_type]}_v[0, 1]',
+        'maskStoreU': f'se{bits[data_type]}_v[2, 0, 1]',
+        'loadS': f'lse{bits[data_type]}_v[0, 1]',
+        'storeS': f'sse{bits[data_type]}_v[0, 2, 1]',
+        'maskStoreS': f'sse{bits[data_type]}_v[3, 0, 2, 1]',
+
+        'abs': 'fabs_v[0]',
+        '==': 'mfeq_vv[0, 1]',
+        '!=': 'mfne_vv[0, 1]',
+        '<=': 'mfle_vv[0, 1]',
+        '<': 'mflt_vv[0, 1]',
+        '>=': 'mfge_vv[0, 1]',
+        '>': 'mfgt_vv[0, 1]',
+        '&': 'mand_mm[0, 1]',
+        '|': 'mor_mm[0, 1]',
+
+        'blendv': 'merge_vvm[2, 0, 1]',
+        'any': 'cpop_m[0]',
+        'all': 'cpop_m[0]',
+    }
+
+    result = dict()
+
+    width = f'vsetvlmax_e{bits[data_type]}m1()'
+    intwidth = 'vsetvlmax_e{bits["int"]}m1()'
+    result['bytes'] = 'vsetvlmax_e8m1()'
+    prefix = 'v'
+    suffix = f'_f{bits[data_type]}m1'
+
+    vl = '{loop_stop} - {loop_counter}'
+    int_vl = f'({vl})*{bits[data_type]//bits["int"]}'
+
+    for intrinsic_id, function_shortcut in base_names.items():
+        function_shortcut = function_shortcut.strip()
+        name = function_shortcut[:function_shortcut.index('[')]
+        if name.startswith('mf'):
+            suffix2 = suffix + f'_b{bits[data_type]}'
+        elif name.endswith('_mm') or name.endswith('_m'):
+            suffix2 = f'_b{bits[data_type]}'
+        elif intrinsic_id.startswith('mask'):
+            suffix2 = suffix + '_m'
+        else:
+            suffix2 = suffix
+
+        arg_string = get_argument_string(function_shortcut, last=vl)
+
+        result[intrinsic_id] = prefix + name + suffix2 + arg_string
+
+    result['width'] = CFunction(width, "int")
+    result['intwidth'] = CFunction(intwidth, "int")
+
+    result['makeVecConst'] = f'vfmv_v_f_f{bits[data_type]}m1({{0}}, {vl})'
+    result['makeVecConstInt'] = f'vmv_v_x_i{bits["int"]}m1({{0}}, {int_vl})'
+    result['makeVecIndex'] = f'vmacc_vx_i{bits["int"]}m1({result["makeVecConstInt"]}, {{1}}, ' + \
+                             f'vid_v_i{bits["int"]}m1({int_vl}), {int_vl})'
+
+    result['storeS'] = result['storeS'].replace('{2}', f'{{2}}*{bits[data_type]//8}')
+    result['loadS'] = result['loadS'].replace('{1}', f'{{1}}*{bits[data_type]//8}')
+    result['maskStoreS'] = result['maskStoreS'].replace('{2}', f'{{2}}*{bits[data_type]//8}')
+
+    result['+int'] = f"vadd_vv_i{bits['int']}m1({{0}}, {{1}}, {int_vl})"
+
+    result['float'] = f'vfloat{bits["float"]}m1_t'
+    result['double'] = f'vfloat{bits["double"]}m1_t'
+    result['int'] = f'vint{bits["int"]}m1_t'
+    result['bool'] = f'vbool{bits[data_type]}_t'
+
+    result['headers'] = ['<riscv_vector.h>', '"riscv_v_helpers.h"']
+
+    result['any'] += ' > 0x0'
+    result['all'] += f' == vsetvl_e{bits[data_type]}m1({vl})'
+
+    result['cachelineSize'] = 'cachelineSize()'
+    result['cachelineZero'] = 'cachelineZero((void*) {0})'
+
+    return result

src/pystencils/backends/simd_instruction_sets.py

+import os
+import platform
+from ctypes import CDLL, c_int, c_size_t, sizeof, byref
+from warnings import warn
+
+import numpy as np
+
+from pystencils.backends.x86_instruction_sets import get_vector_instruction_set_x86
+from pystencils.backends.arm_instruction_sets import get_vector_instruction_set_arm
+from pystencils.backends.ppc_instruction_sets import get_vector_instruction_set_ppc
+from pystencils.backends.riscv_instruction_sets import get_vector_instruction_set_riscv
+from pystencils.cache import memorycache
+from pystencils.typing import numpy_name_to_c
+
+
+def get_vector_instruction_set(data_type='double', instruction_set='avx'):
+    if data_type == 'float':
+        warn(f"Ambiguous input for data_type: {data_type}. For single precision please use float32. "
+             f"For more information please take numpy.dtype as a reference. This input will not be supported in future "
+             f"releases")
+        data_type = 'float64'
+
+    type_name = numpy_name_to_c(np.dtype(data_type).name)
+
+    if instruction_set in ['neon', 'sme'] or instruction_set.startswith('sve'):
+        return get_vector_instruction_set_arm(type_name, instruction_set)
+    elif instruction_set in ['vsx']:
+        return get_vector_instruction_set_ppc(type_name, instruction_set)
+    elif instruction_set in ['rvv']:
+        return get_vector_instruction_set_riscv(type_name, instruction_set)
+    else:
+        return get_vector_instruction_set_x86(type_name, instruction_set)
+
+
+@memorycache
+def get_supported_instruction_sets():
+    """List of supported instruction sets on current hardware, or None if query failed."""
+    if 'PYSTENCILS_SIMD' in os.environ:
+        return os.environ['PYSTENCILS_SIMD'].split(',')
+    if platform.system() == 'Darwin' and platform.machine() == 'arm64':
+        result = ['neon']
+        libc = CDLL('/usr/lib/libc.dylib')
+        value = c_int(0)
+        size = c_size_t(sizeof(value))
+        status = libc.sysctlbyname(b"hw.optional.arm.FEAT_SME", byref(value), byref(size), None, 0)
+        if status == 0 and value.value == 1:
+            result.insert(0, "sme")
+        return result
+    elif platform.system() == 'Windows' and platform.machine() == 'ARM64':
+        return ['neon']
+    elif platform.system() == 'Linux' and platform.machine() == 'aarch64':
+        result = ['neon']  # Neon is mandatory on 64-bit ARM
+        libc = CDLL('libc.so.6')
+        hwcap = libc.getauxval(16)  # AT_HWCAP
+        hwcap2 = libc.getauxval(26)  # AT_HWCAP2
+        if hwcap & (1 << 22):  # HWCAP_SVE
+            if hwcap2 & (1 << 1):  # HWCAP2_SVE2
+                name = 'sve2'
+            else:
+                name = 'sve'
+            length = 8 * libc.prctl(51, 0, 0, 0, 0)  # PR_SVE_GET_VL
+            if length < 0:
+                raise OSError("SVE length query failed")
+            while length >= 128:
+                result.append(f"{name}{length}")
+                length //= 2
+            result.append(name)
+        if hwcap2 & (1 << 23):  # HWCAP2_SME
+            result.insert(0, "sme")  # prepend to list so it is not automatically chosen as best instruction set
+        return result
+    elif platform.system() == 'Linux' and platform.machine().startswith('riscv'):
+        libc = CDLL('libc.so.6')
+        hwcap = libc.getauxval(16)  # AT_HWCAP
+        hwcap_isa_v = 1 << (ord('V') - ord('A'))  # COMPAT_HWCAP_ISA_V
+        return ['rvv'] if hwcap & hwcap_isa_v else []
+    elif platform.system() == 'Linux' and platform.machine().startswith('ppc64'):
+        libc = CDLL('libc.so.6')
+        hwcap = libc.getauxval(16)  # AT_HWCAP
+        return ['vsx'] if hwcap & 0x00000080 else []  # PPC_FEATURE_HAS_VSX
+    elif platform.machine() in ['x86_64', 'x86', 'AMD64', 'i386']:
+        try:
+            from cpuinfo import get_cpu_info
+        except ImportError:
+            return None
+
+        result = []
+        required_sse_flags = {'sse', 'sse2', 'ssse3', 'sse4_1', 'sse4_2'}
+        required_avx_flags = {'avx', 'avx2'}
+        required_avx512_flags = {'avx512f'}
+        possible_avx512vl_flags = {'avx512vl', 'avx10_1'}
+        flags = set(get_cpu_info()['flags'])
+        if flags.issuperset(required_sse_flags):
+            result.append("sse")
+        if flags.issuperset(required_avx_flags):
+            result.append("avx")
+        if flags.issuperset(required_avx512_flags):
+            result.append("avx512")
+        if not flags.isdisjoint(possible_avx512vl_flags):
+            result.append("avx512vl")
+        return result
+    else:
+        raise NotImplementedError('Instruction set detection for %s on %s is not implemented' %
+                                  (platform.system(), platform.machine()))
+
+
+@memorycache
+def get_cacheline_size(instruction_set):
+    """Get the size (in bytes) of a cache block that can be zeroed without memory access.
+       Usually, this is identical to the cache line size."""
+    
+    instruction_sets = get_vector_instruction_set('double', instruction_set)
+    if 'cachelineSize' not in instruction_sets:
+        return None
+    
+    import pystencils as ps
+    from pystencils.astnodes import SympyAssignment
+    import numpy as np
+    from pystencils.cpu.vectorization import CachelineSize
+    
+    arr = np.zeros((1, 1), dtype=np.float32)
+    f = ps.Field.create_from_numpy_array('f', arr, index_dimensions=0)
+    ass = [CachelineSize(), SympyAssignment(f.center, CachelineSize.symbol)]
+    ast = ps.create_kernel(ass, cpu_vectorize_info={'instruction_set': instruction_set})
+    kernel = ast.compile()
+    kernel(**{f.name: arr, CachelineSize.symbol.name: 0})
+    return int(arr[0, 0])

src/pystencils/backends/x86_instruction_sets.py

+def get_argument_string(intrinsic_id, width, function_shortcut):
+    if intrinsic_id == 'makeVecConst' or intrinsic_id == 'makeVecConstInt':
+        arg_string = f"({','.join(['{0}'] * width)})"
+    elif intrinsic_id == 'makeVec' or intrinsic_id == 'makeVecInt':
+        params = ["{" + str(i) + "}" for i in reversed(range(width))]
+        arg_string = f"({','.join(params)})"
+    elif intrinsic_id == 'makeVecBool':
+        params = [f"(({{{i}}} ? -1.0 : 0.0)" for i in reversed(range(width))]
+        arg_string = f"({','.join(params)})"
+    elif intrinsic_id == 'makeVecConstBool':
+        params = ["(({0}) ? -1.0 : 0.0)" for _ in range(width)]
+        arg_string = f"({','.join(params)})"
+    else:
+        args = function_shortcut[function_shortcut.index('[') + 1: -1]
+        arg_string = "("
+        for arg in args.split(","):
+            arg = arg.strip()
+            if not arg:
+                continue
+            if arg in ('0', '1', '2', '3', '4', '5'):
+                arg_string += "{" + arg + "},"
+            else:
+                arg_string += arg + ","
+        arg_string = arg_string[:-1] + ")"
+    return arg_string
+
+
+def get_vector_instruction_set_x86(data_type='double', instruction_set='avx'):
+    comparisons = {
+        '==': '_CMP_EQ_UQ',
+        '!=': '_CMP_NEQ_UQ',
+        '>=': '_CMP_GE_OQ',
+        '<=': '_CMP_LE_OQ',
+        '<': '_CMP_NGE_UQ',
+        '>': '_CMP_NLE_UQ',
+    }
+    base_names = {
+        '+': 'add[0, 1]',
+        '-': 'sub[0, 1]',
+        '*': 'mul[0, 1]',
+        '/': 'div[0, 1]',
+        '&': 'and[0, 1]',
+        '|': 'or[0, 1]',
+        'blendv': 'blendv[0, 1, 2]',
+
+        'sqrt': 'sqrt[0]',
+
+        'makeVecConst': 'set[]',
+        'makeVec': 'set[]',
+        'makeVecBool': 'set[]',
+        'makeVecConstBool': 'set[]',
+        'makeVecInt': 'set[]',
+        'makeVecConstInt': 'set[]',
+
+        'loadU': 'loadu[0]',
+        'loadA': 'load[0]',
+        'storeU': 'storeu[0,1]',
+        'storeA': 'store[0,1]',
+        'stream': 'stream[0,1]',
+        'maskStoreA': 'mask_store[0, 2, 1]' if instruction_set.startswith('avx512') else 'maskstore[0, 2, 1]',
+        'maskStoreU': 'mask_storeu[0, 2, 1]' if instruction_set.startswith('avx512') else 'maskstore[0, 2, 1]',
+    }
+
+    for comparison_op, constant in comparisons.items():
+        base_names[comparison_op] = f'cmp[0, 1, {constant}]'
+
+    headers = {
+        'avx512': ['<immintrin.h>'],
+        'avx512vl': ['<immintrin.h>'],
+        'avx': ['<immintrin.h>'],
+        'sse': ['<immintrin.h>', '<xmmintrin.h>', '<emmintrin.h>', '<pmmintrin.h>',
+                '<tmmintrin.h>', '<smmintrin.h>', '<nmmintrin.h>']
+    }
+
+    suffix = {
+        'double': 'pd',
+        'float': 'ps',
+        'int': 'epi32'
+    }
+    prefix = {
+        'sse': '_mm',
+        'avx': '_mm256',
+        'avx512vl': '_mm256',
+        'avx512': '_mm512',
+    }
+
+    width = {
+        ("double", "sse"): 2,
+        ("float", "sse"): 4,
+        ("int", "sse"): 4,
+        ("double", "avx"): 4,
+        ("float", "avx"): 8,
+        ("int", "avx"): 8,
+        ("double", "avx512vl"): 4,
+        ("float", "avx512vl"): 8,
+        ("int", "avx512vl"): 8,
+        ("double", "avx512"): 8,
+        ("float", "avx512"): 16,
+        ("int", "avx512"): 16,
+    }
+    result = {
+        'width': width[(data_type, instruction_set)],
+        'intwidth': width[('int', instruction_set)],
+        'bytes': 4 * width[("float", instruction_set)]
+    }
+    pre = prefix[instruction_set]
+    for intrinsic_id, function_shortcut in base_names.items():
+        function_shortcut = function_shortcut.strip()
+        name = function_shortcut[:function_shortcut.index('[')]
+
+        if 'Int' in intrinsic_id:
+            suf = suffix['int']
+            arg_string = get_argument_string(intrinsic_id, result['intwidth'], function_shortcut)
+        else:
+            suf = suffix[data_type]
+            arg_string = get_argument_string(intrinsic_id, result['width'], function_shortcut)
+
+        mask_suffix = '_mask' if instruction_set.startswith('avx512') and intrinsic_id in comparisons.keys() else ''
+        result[intrinsic_id] = pre + "_" + name + "_" + suf + mask_suffix + arg_string
+
+    bit_width = result['width'] * (64 if data_type == 'double' else 32)
+    result['double'] = f"__m{bit_width}d"
+    result['float'] = f"__m{bit_width}"
+    result['int'] = f"__m{bit_width}i"
+    result['bool'] = result[data_type]
+
+    result['headers'] = headers[instruction_set]
+    result['any'] = f"{pre}_movemask_{suf}({{0}}) > 0"
+    result['all'] = f"{pre}_movemask_{suf}({{0}}) == {hex(2**result['width']-1)}"
+
+    setsuf = "x" if bit_width < 512 and bit_width // result['width'] == 64 else ""
+
+    if instruction_set.startswith('avx512'):
+        size = result['width']
+        masksize = max(size, 8)
+        result['&'] = f'_kand_mask{masksize}({{0}}, {{1}})'
+        result['|'] = f'_kor_mask{masksize}({{0}}, {{1}})'
+        result['any'] = f'!_ktestz_mask{masksize}_u8({{0}}, {{0}})'
+        result['all'] = f'_kortestc_mask{masksize}_u8({{0}}, {{0}})'
+        result['blendv'] = f'{pre}_mask_blend_{suf}({{2}}, {{0}}, {{1}})'
+        result['rsqrt'] = f"{pre}_rsqrt14_{suf}({{0}})"
+        result['bool'] = f"__mmask{masksize}"
+
+        params = " | ".join(["({{{i}}} ? {power} : 0)".format(i=i, power=2 ** i) for i in range(8)])
+        result['makeVecBool'] = f"__mmask8(({params}) )"
+        params = " | ".join(["({{0}} ? {power} : 0)".format(power=2 ** i) for i in range(8)])
+        result['makeVecConstBool'] = f"__mmask8(({params}) )"
+
+        vindex = f'{pre}_set_epi{bit_width//size}{setsuf}(' + \
+                 ', '.join([str(i) for i in range(result['width'])][::-1]) + ')'
+        vindex = f'{pre}_mullo_epi{bit_width//size}({vindex}, {pre}_set1_epi{bit_width//size}{setsuf}({{0}}))'
+        scale = bit_width // size // 8
+        result['storeS'] = f'{pre}_i{bit_width//size}scatter_{suf}({{0}}, ' + vindex.format("{2}") + \
+                           f', {{1}}, {scale})'
+        result['maskStoreS'] = f'{pre}_mask_i{bit_width//size}scatter_{suf}({{0}}, {{3}}, ' + vindex.format("{2}") + \
+                               f', {{1}}, {scale})'
+        if bit_width == 512:
+            result['loadS'] = f'{pre}_i{bit_width//size}gather_{suf}(' + vindex.format("{1}") + f', {{0}}, {scale})'
+        else:
+            result['loadS'] = f'{pre}_i{bit_width//size}gather_{suf}({{0}}, ' + vindex.format("{1}") + f', {scale})'
+
+    # abs intrinsic exists in 512 bits, but expands to a sequence. We generate that same sequence for 128 and 256 bits
+    if instruction_set == 'avx512':
+        result['abs'] = f"{pre}_abs_{suf}({{0}})"
+    else:
+        result['abs'] = f"{pre}_castsi{bit_width}_{suf}({pre}_and_si{bit_width}(" + \
+                        f"{pre}_set1_epi{bit_width // result['width']}{setsuf}(0x7" + \
+                        'f' * (bit_width // result['width'] // 4 - 1) + "), " + \
+                        f"{pre}_cast{suf}_si{bit_width}({{0}})))"
+
+    if instruction_set == 'avx' and data_type == 'float':
+        result['rsqrt'] = f"{pre}_rsqrt_{suf}({{0}})"
+
+    result['+int'] = f"{pre}_add_{suffix['int']}({{0}}, {{1}})"
+
+    result['streamFence'] = '_mm_mfence()'
+
+    return result

src/pystencils/bit_masks.py

+import sympy as sp
+# from pystencils.typing import get_type_of_expression
+
+
+# noinspection PyPep8Naming
+class flag_cond(sp.Function):
+    """Evaluates a flag condition on a bit mask, and returns the value of one of two expressions,
+    depending on whether the flag is set. 
+
+    Three argument version:
+    ```
+        flag_cond(flag_bit, mask, expr) = expr if (flag_bit is set in mask) else 0
+    ```
+
+    Four argument version:
+    ```
+        flag_cond(flag_bit, mask, expr_then, expr_else) = expr_then if (flag_bit is set in mask) else expr_else
+    ```
+    """
+
+    nargs = (3, 4)
+
+    def __new__(cls, flag_bit, mask_expression, *expressions):
+
+        # TODO Jan reintroduce checking
+        # flag_dtype = get_type_of_expression(flag_bit)
+        # if not flag_dtype.is_int():
+        #     raise ValueError('Argument flag_bit must be of integer type.')
+        #
+        # mask_dtype = get_type_of_expression(mask_expression)
+        # if not mask_dtype.is_int():
+        #     raise ValueError('Argument mask_expression must be of integer type.')
+
+        return super().__new__(cls, flag_bit, mask_expression, *expressions)
+
+    def to_c(self, print_func):
+        flag_bit = self.args[0]
+        mask = self.args[1]
+
+        then_expression = self.args[2]
+
+        flag_bit_code = print_func(flag_bit)
+        mask_code = print_func(mask)
+        then_code = print_func(then_expression)
+
+        code = f"(({mask_code}) >> ({flag_bit_code}) & 1) * ({then_code})"
+
+        if len(self.args) > 3:
+            else_expression = self.args[3]
+            else_code = print_func(else_expression)
+            code += f" + (({mask_code}) >> ({flag_bit_code}) ^ 1) * ({else_code})"
+
+        return code

src/pystencils/boundaries/__init__.py

+from pystencils.boundaries.boundaryconditions import Dirichlet, Neumann
+from pystencils.boundaries.boundaryhandling import BoundaryHandling
+from pystencils.boundaries.inkernel import add_neumann_boundary
+
+__all__ = ['BoundaryHandling', 'Neumann', 'Dirichlet', 'add_neumann_boundary']

src/pystencils/boundaries/boundaryconditions.py

+from typing import Any, List, Tuple
+
+from pystencils.astnodes import SympyAssignment
+from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
+from pystencils.typing import create_type
+
+
+class Boundary:
+    """Base class all boundaries should derive from"""
+
+    inner_or_boundary = True
+    single_link = False
+
+    def __init__(self, name=None):
+        self._name = name
+
+    def __call__(self, field, direction_symbol, index_field) -> List[SympyAssignment]:
+        """Defines the boundary behavior and must therefore be implemented by all boundaries.
+
+        Here the boundary is defined as a list of sympy assignments, from which a boundary kernel is generated.
+
+        Args:
+            field: pystencils field where boundary condition should be applied.
+                   The current cell is cell next to the boundary, which is influenced by the boundary
+                   cell i.e. has a link from the boundary cell to itself.
+            direction_symbol: a sympy symbol that can be used as index to the pdf_field. It describes
+                              the direction pointing from the fluid to the boundary cell
+            index_field: the boundary index field that can be used to retrieve and update boundary data
+        """
+        raise NotImplementedError("Boundary class has to overwrite __call__")
+
+    @property
+    def additional_data(self) -> Tuple[str, Any]:
+        """Return a list of (name, type) tuples for additional data items required in this boundary
+        These data items can either be initialized in separate kernel see additional_data_kernel_init or by
+        Python callbacks - see additional_data_callback """
+        return ()
+
+    @property
+    def additional_data_init_callback(self):
+        """Return a callback function called with a boundary data setter object and returning a dict of
+        data-name to data for each element that should be initialized"""
+        return None
+
+    @property
+    def name(self):
+        if self._name:
+            return self._name
+        else:
+            return type(self).__name__
+
+    @name.setter
+    def name(self, new_value):
+        self._name = new_value
+
+
+class Neumann(Boundary):
+
+    inner_or_boundary = False
+    single_link = True
+
+    def __call__(self, field, direction_symbol, **kwargs):
+
+        neighbor = BoundaryOffsetInfo.offset_from_dir(direction_symbol, field.spatial_dimensions)
+        if field.index_dimensions == 0:
+            return [SympyAssignment(field.center, field[neighbor])]
+        else:
+            from itertools import product
+            if not field.has_fixed_index_shape:
+                raise NotImplementedError("Neumann boundary works only for fields with fixed index shape")
+            index_iter = product(*(range(i) for i in field.index_shape))
+            return [SympyAssignment(field(*idx), field[neighbor](*idx)) for idx in index_iter]
+
+    def __hash__(self):
+        # All boundaries of these class behave equal -> should also be equal
+        return hash("Neumann")
+
+    def __eq__(self, other):
+        return type(other) is Neumann
+
+
+class Dirichlet(Boundary):
+
+    inner_or_boundary = False
+    single_link = True
+
+    def __init__(self, value, name=None):
+        super().__init__(name)
+        self._value = value
+
+    @property
+    def additional_data(self):
+        if callable(self._value):
+            return [('value', create_type("double"))]
+        else:
+            return []
+
+    @property
+    def additional_data_init_callback(self):
+        if callable(self._value):
+            return self._value
+
+    def __call__(self, field, direction_symbol, index_field, **kwargs):
+
+        if field.index_dimensions == 0:
+            return [SympyAssignment(field.center, index_field("value") if self.additional_data else self._value)]
+        elif field.index_dimensions == 1:
+            assert not self.additional_data
+            if not field.has_fixed_index_shape:
+                raise NotImplementedError("Field needs fixed index shape")
+            assert len(self._value) == field.index_shape[0], "Dirichlet value does not match index shape of field"
+            return [SympyAssignment(field(i), self._value[i]) for i in range(field.index_shape[0])]
+        raise NotImplementedError("Dirichlet boundary not implemented for fields with more than one index dimension")

src/pystencils/boundaries/boundaryhandling.py

+from functools import lru_cache
+
+import numpy as np
+import sympy as sp
+
+from pystencils import create_kernel, CreateKernelConfig, Target
+from pystencils.astnodes import SympyAssignment
+from pystencils.backends.cbackend import CustomCodeNode
+from pystencils.boundaries.createindexlist import (
+    create_boundary_index_array, numpy_data_type_for_boundary_object)
+from pystencils.typing import TypedSymbol, create_type
+from pystencils.gpu.gpu_array_handler import GPUArrayHandler
+from pystencils.field import Field
+from pystencils.typing.typed_sympy import FieldPointerSymbol
+
+try:
+    # noinspection PyPep8Naming
+    import waLBerla as wlb
+    if wlb.cpp_available:
+        from pystencils.datahandling.parallel_datahandling import ParallelDataHandling
+        import cupy.cuda.runtime
+    else:
+        ParallelDataHandling = None
+except ImportError:
+    ParallelDataHandling = None
+
+DEFAULT_FLAG_TYPE = np.uint32
+
+
+class FlagInterface:
+    """Manages the reservation of bits (i.e. flags) in an array of unsigned integers.
+
+    Examples:
+        >>> from pystencils import create_data_handling
+        >>> dh = create_data_handling((4, 5))
+        >>> fi = FlagInterface(dh, 'flag_field', np.uint8)
+        >>> assert dh.has_data('flag_field')
+        >>> int(fi.reserve_next_flag())
+        2
+        >>> int(fi.reserve_flag(4))
+        4
+        >>> int(fi.reserve_next_flag())
+        8
+    """
+
+    def __init__(self, data_handling, flag_field_name, dtype=DEFAULT_FLAG_TYPE):
+        self.flag_field_name = flag_field_name
+        self.domain_flag = dtype(1 << 0)
+        self._used_flags = {self.domain_flag}
+        self.data_handling = data_handling
+        self.dtype = dtype
+        self.max_bits = self.dtype().itemsize * 8
+
+        # Add flag field to data handling if it does not yet exist
+        if data_handling.has_data(self.flag_field_name):
+            raise ValueError("There is already a boundary handling registered at the data handling."
+                             "If you want to add multiple handling objects, choose a different name.")
+
+        self.flag_field = data_handling.add_array(self.flag_field_name, dtype=self.dtype, cpu=True, gpu=False)
+        ff_ghost_layers = data_handling.ghost_layers_of_field(self.flag_field_name)
+        for b in data_handling.iterate(ghost_layers=ff_ghost_layers):
+            b[self.flag_field_name].fill(self.domain_flag)
+
+    def reserve_next_flag(self):
+        for i in range(1, self.max_bits):
+            flag = self.dtype(1 << i)
+            if flag not in self._used_flags:
+                self._used_flags.add(flag)
+                assert self._is_power_of_2(flag)
+                return flag
+        raise ValueError(f"All available {self.max_bits} flags are reserved")
+
+    def reserve_flag(self, flag):
+        assert self._is_power_of_2(flag)
+        flag = self.dtype(flag)
+        if flag in self._used_flags:
+            raise ValueError(f"The flag {flag} is already reserved")
+        self._used_flags.add(flag)
+        return flag
+
+    @staticmethod
+    def _is_power_of_2(num):
+        return num != 0 and ((num & (num - 1)) == 0)
+
+
+class BoundaryHandling:
+
+    def __init__(self, data_handling, field_name, stencil, name="boundary_handling", flag_interface=None,
+                 target: Target = Target.CPU, openmp=True):
+        assert data_handling.has_data(field_name)
+        assert data_handling.dim == len(stencil[0]), "Dimension of stencil and data handling do not match"
+        self._data_handling = data_handling
+        self._field_name = field_name
+        self._index_array_name = name + "IndexArrays"
+        self._target = target
+        self._openmp = openmp
+        self._boundary_object_to_boundary_info = {}
+        self.stencil = stencil
+        self._dirty = True
+        fi = flag_interface
+        self.flag_interface = fi if fi is not None else FlagInterface(data_handling, name + "Flags")
+
+        if ParallelDataHandling and isinstance(self.data_handling, ParallelDataHandling):
+            array_handler = GPUArrayHandler(cupy.cuda.runtime.getDevice())
+        else:
+            array_handler = self.data_handling.array_handler
+
+        def to_cpu(gpu_version, cpu_version):
+            gpu_version = gpu_version.boundary_object_to_index_list
+            cpu_version = cpu_version.boundary_object_to_index_list
+            for obj, cpu_arr in cpu_version.items():
+                array_handler.download(gpu_version[obj], cpu_arr)
+
+        def to_gpu(gpu_version, cpu_version):
+            gpu_version = gpu_version.boundary_object_to_index_list
+            cpu_version = cpu_version.boundary_object_to_index_list
+
+            for obj, cpu_arr in cpu_version.items():
+                if obj not in gpu_version or gpu_version[obj].shape != cpu_arr.shape:
+                    gpu_version[obj] = array_handler.empty(cpu_arr.shape, cpu_arr.dtype)
+                    array_handler.upload(gpu_version[obj], cpu_arr)
+                else:
+                    array_handler.upload(gpu_version[obj], cpu_arr)
+
+        class_ = self.IndexFieldBlockData
+        class_.to_cpu = to_cpu
+        class_.to_gpu = to_gpu
+        gpu = self._target in data_handling._GPU_LIKE_TARGETS
+        data_handling.add_custom_class(self._index_array_name, class_, cpu=True, gpu=gpu)
+
+    @property
+    def data_handling(self):
+        return self._data_handling
+
+    def get_flag(self, boundary_obj):
+        return self._boundary_object_to_boundary_info[boundary_obj].flag
+
+    @property
+    def shape(self):
+        return self._data_handling.shape
+
+    @property
+    def dim(self):
+        return self._data_handling.dim
+
+    @property
+    def boundary_objects(self):
+        return tuple(self._boundary_object_to_boundary_info.keys())
+
+    @property
+    def flag_array_name(self):
+        return self.flag_interface.flag_field_name
+
+    def get_mask(self, slice_obj, boundary_obj, inverse=False):
+        if isinstance(boundary_obj, str) and boundary_obj.lower() == 'domain':
+            flag = self.flag_interface.domain_flag
+        else:
+            flag = self._boundary_object_to_boundary_info[boundary_obj].flag
+
+        arr = self.data_handling.gather_array(self.flag_array_name, slice_obj)
+        if arr is None:
+            return None
+        else:
+            result = np.bitwise_and(arr, flag)
+            if inverse:
+                result = np.logical_not(result)
+            return result
+
+    def set_boundary(self, boundary_obj, slice_obj=None, mask_callback=None,
+                     ghost_layers=True, inner_ghost_layers=True, replace=True, force_flag_value=None):
+        """Sets boundary using either a rectangular slice, a boolean mask or a combination of both.
+
+        Args:
+            boundary_obj: instance of a boundary object that should be set
+            slice_obj: a slice object (can be created with make_slice[]) that selects a part of the domain where
+                       the boundary should be set. If none, the complete domain is selected which makes only sense
+                       if a mask_callback is passed. The slice can have ':' placeholders, which are interpreted
+                       depending on the 'inner_ghost_layers' parameter i.e. if it is True, the slice extends
+                       into the ghost layers
+            mask_callback: callback function getting x,y (z) parameters of the cell midpoints and returning a
+                          boolean mask with True entries where boundary cells should be set.
+                          The x, y, z arrays have 2D/3D shape such that they can be used directly
+                          to create the boolean return array. i.e return x < 10 sets boundaries in cells with
+                          midpoint x coordinate smaller than 10.
+            ghost_layers: see DataHandling.iterate()
+            inner_ghost_layers: see DataHandling.iterate()
+            replace: by default all other flags are erased in the cells where the boundary is set. To add a
+                     boundary condition, set this replace flag to False
+            force_flag_value: flag that should be reserved for this boundary. Has to be an integer that is a power of 2
+                              and was not reserved before for another boundary.
+        """
+        if isinstance(boundary_obj, str) and boundary_obj.lower() == 'domain':
+            flag = self.flag_interface.domain_flag
+        else:
+            if force_flag_value:
+                self.flag_interface.reserve_flag(force_flag_value)
+            flag = self._add_boundary(boundary_obj, force_flag_value)
+
+        for b in self._data_handling.iterate(slice_obj, ghost_layers=ghost_layers,
+                                             inner_ghost_layers=inner_ghost_layers):
+            flag_arr = b[self.flag_interface.flag_field_name]
+            if mask_callback is not None:
+                mask = mask_callback(*b.midpoint_arrays)
+                if replace:
+                    flag_arr[mask] = flag
+                else:
+                    np.bitwise_or(flag_arr, flag, where=mask, out=flag_arr)
+                    np.bitwise_and(flag_arr, ~self.flag_interface.domain_flag, where=mask, out=flag_arr)
+            else:
+                if replace:
+                    flag_arr.fill(flag)
+                else:
+                    np.bitwise_or(flag_arr, flag, out=flag_arr)
+                    np.bitwise_and(flag_arr, ~self.flag_interface.domain_flag, out=flag_arr)
+
+        self._dirty = True
+
+        return flag
+
+    def set_boundary_where_flag_is_set(self, boundary_obj, flag):
+        """Adds an (additional) boundary to all cells that have been previously marked with the passed flag."""
+        self._add_boundary(boundary_obj, flag)
+        self._dirty = True
+        return flag
+
+    def prepare(self):
+        if not self._dirty:
+            return
+        self._create_index_fields()
+        self._dirty = False
+
+    def trigger_reinitialization_of_boundary_data(self, **kwargs):
+        if self._dirty:
+            self.prepare()
+        else:
+            ff_ghost_layers = self._data_handling.ghost_layers_of_field(self.flag_interface.flag_field_name)
+            for b in self._data_handling.iterate(ghost_layers=ff_ghost_layers):
+                for b_obj, setter in b[self._index_array_name].boundary_object_to_data_setter.items():
+                    self._boundary_data_initialization(b_obj, setter, **kwargs)
+
+    def __call__(self, **kwargs):
+        if self._dirty:
+            self.prepare()
+
+        for b in self._data_handling.iterate(gpu=self._target in self._data_handling._GPU_LIKE_TARGETS):
+            for b_obj, idx_arr in b[self._index_array_name].boundary_object_to_index_list.items():
+                kwargs[self._field_name] = b[self._field_name]
+                kwargs['indexField'] = idx_arr
+                data_used_in_kernel = (p.fields[0].name
+                                       for p in self._boundary_object_to_boundary_info[b_obj].kernel.parameters
+                                       if isinstance(p.symbol, FieldPointerSymbol) and p.fields[0].name not in kwargs)
+                kwargs.update({name: b[name] for name in data_used_in_kernel})
+
+                self._boundary_object_to_boundary_info[b_obj].kernel(**kwargs)
+
+    def add_fixed_steps(self, fixed_loop, **kwargs):
+        if self._dirty:
+            self.prepare()
+
+        for b in self._data_handling.iterate(gpu=self._target in self._data_handling._GPU_LIKE_TARGETS):
+            for b_obj, idx_arr in b[self._index_array_name].boundary_object_to_index_list.items():
+                arguments = kwargs.copy()
+                arguments[self._field_name] = b[self._field_name]
+                arguments['indexField'] = idx_arr
+                data_used_in_kernel = (p.fields[0].name
+                                       for p in self._boundary_object_to_boundary_info[b_obj].kernel.parameters
+                                       if isinstance(p.symbol, FieldPointerSymbol) and p.field_name not in arguments)
+                arguments.update({name: b[name] for name in data_used_in_kernel if name not in arguments})
+
+                kernel = self._boundary_object_to_boundary_info[b_obj].kernel
+                fixed_loop.add_call(kernel, arguments)
+
+    def geometry_to_vtk(self, file_name='geometry', boundaries='all', ghost_layers=False):
+        """
+        Writes a VTK field where each cell with the given boundary is marked with 1, other cells are 0
+        This can be used to display the simulation geometry in Paraview
+
+        Params:
+            file_name: vtk filename
+            boundaries: boundary object, or special string 'domain' for domain cells or special string 'all' for all
+                      boundary conditions.
+                      can also  be a sequence, to write multiple boundaries to VTK file
+            ghost_layers: number of ghost layers to write, or True for all, False for none
+        """
+        if boundaries == 'all':
+            boundaries = list(self._boundary_object_to_boundary_info.keys()) + ['domain']
+        elif not hasattr(boundaries, "__len__"):
+            boundaries = [boundaries]
+
+        masks_to_name = {}
+        for b in boundaries:
+            if b == 'domain':
+                masks_to_name[self.flag_interface.domain_flag] = 'domain'
+            else:
+                flag = self._boundary_object_to_boundary_info[b].flag
+                masks_to_name[flag] = b.name
+
+        writer = self.data_handling.create_vtk_writer_for_flag_array(file_name, self.flag_interface.flag_field_name,
+                                                                     masks_to_name, ghost_layers=ghost_layers)
+        writer(1)
+
+    # ------------------------------ Implementation Details ------------------------------------------------------------
+
+    def _add_boundary(self, boundary_obj, flag=None):
+        if boundary_obj not in self._boundary_object_to_boundary_info:
+            sym_index_field = Field.create_generic('indexField', spatial_dimensions=1,
+                                                   dtype=numpy_data_type_for_boundary_object(boundary_obj, self.dim))
+            ast = self._create_boundary_kernel(self._data_handling.fields[self._field_name],
+                                               sym_index_field, boundary_obj)
+            if flag is None:
+                flag = self.flag_interface.reserve_next_flag()
+            boundary_info = self.BoundaryInfo(boundary_obj, flag=flag, kernel=ast.compile())
+            self._boundary_object_to_boundary_info[boundary_obj] = boundary_info
+        return self._boundary_object_to_boundary_info[boundary_obj].flag
+
+    def _create_boundary_kernel(self, symbolic_field, symbolic_index_field, boundary_obj):
+        return create_boundary_kernel(symbolic_field, symbolic_index_field, self.stencil, boundary_obj,
+                                      target=self._target, cpu_openmp=self._openmp)
+
+    def _create_index_fields(self):
+        dh = self._data_handling
+        ff_ghost_layers = dh.ghost_layers_of_field(self.flag_interface.flag_field_name)
+        for b in dh.iterate(ghost_layers=ff_ghost_layers):
+            flag_arr = b[self.flag_interface.flag_field_name]
+            pdf_arr = b[self._field_name]
+            index_array_bd = b[self._index_array_name]
+            index_array_bd.clear()
+            for b_info in self._boundary_object_to_boundary_info.values():
+                boundary_obj = b_info.boundary_object
+                idx_arr = create_boundary_index_array(flag_arr, self.stencil, b_info.flag,
+                                                      self.flag_interface.domain_flag, boundary_obj,
+                                                      ff_ghost_layers, boundary_obj.inner_or_boundary,
+                                                      boundary_obj.single_link)
+                if idx_arr.size == 0:
+                    continue
+
+                boundary_data_setter = BoundaryDataSetter(idx_arr, b.offset, self.stencil, ff_ghost_layers, pdf_arr)
+                index_array_bd.boundary_object_to_index_list[b_info.boundary_object] = idx_arr
+                index_array_bd.boundary_object_to_data_setter[b_info.boundary_object] = boundary_data_setter
+                self._boundary_data_initialization(b_info.boundary_object, boundary_data_setter)
+
+    def _boundary_data_initialization(self, boundary_obj, boundary_data_setter, **kwargs):
+        if boundary_obj.additional_data_init_callback:
+            boundary_obj.additional_data_init_callback(boundary_data_setter, **kwargs)
+        if self._target in self._data_handling._GPU_LIKE_TARGETS:
+            self._data_handling.to_gpu(self._index_array_name)
+
+    class BoundaryInfo(object):
+        def __init__(self, boundary_obj, flag, kernel):
+            self.boundary_object = boundary_obj
+            self.flag = flag
+            self.kernel = kernel
+
+    class IndexFieldBlockData:
+        def __init__(self, *args, **kwargs):
+            self.boundary_object_to_index_list = {}
+            self.boundary_object_to_data_setter = {}
+
+        def clear(self):
+            self.boundary_object_to_index_list.clear()
+            self.boundary_object_to_data_setter.clear()
+
+
+class BoundaryDataSetter:
+
+    def __init__(self, index_array, offset, stencil, ghost_layers, pdf_array):
+        self.index_array = index_array
+        self.offset = offset
+        self.stencil = np.array(stencil)
+        self.pdf_array = pdf_array.view()
+        self.pdf_array.flags.writeable = False
+
+        arr_field_names = index_array.dtype.names
+        self.dim = 3 if 'z' in arr_field_names else 2
+        assert 'x' in arr_field_names and 'y' in arr_field_names and 'dir' in arr_field_names, str(arr_field_names)
+        self.boundary_data_names = set(self.index_array.dtype.names) - {'x', 'y', 'z', 'dir'}
+        self.coord_map = {0: 'x', 1: 'y', 2: 'z'}
+        self.ghost_layers = ghost_layers
+
+    def non_boundary_cell_positions(self, coord):
+        assert coord < self.dim
+        return self.index_array[self.coord_map[coord]] + self.offset[coord] - self.ghost_layers + 0.5
+
+    @lru_cache()
+    def link_offsets(self):
+        return self.stencil[self.index_array['dir']]
+
+    @lru_cache()
+    def link_positions(self, coord):
+        return self.non_boundary_cell_positions(coord) + 0.5 * self.link_offsets()[:, coord]
+
+    @lru_cache()
+    def boundary_cell_positions(self, coord):
+        return self.non_boundary_cell_positions(coord) + self.link_offsets()[:, coord]
+
+    def __setitem__(self, key, value):
+        if key not in self.boundary_data_names:
+            raise KeyError(f"Invalid boundary data name {key}. Allowed are {self.boundary_data_names}")
+        self.index_array[key] = value
+
+    def __getitem__(self, item):
+        if item not in self.boundary_data_names:
+            raise KeyError(f"Invalid boundary data name {item}. Allowed are {self.boundary_data_names}")
+        return self.index_array[item]
+
+
+class BoundaryOffsetInfo(CustomCodeNode):
+
+    # --------------------------- Functions to be used by boundaries --------------------------
+
+    @staticmethod
+    def offset_from_dir(dir_idx, dim):
+        return tuple([sp.IndexedBase(symbol, shape=(1,))[dir_idx]
+                      for symbol in BoundaryOffsetInfo._offset_symbols(dim)])
+
+    @staticmethod
+    def inv_dir(dir_idx):
+        return sp.IndexedBase(BoundaryOffsetInfo.INV_DIR_SYMBOL, shape=(1,))[dir_idx]
+
+    # ---------------------------------- Internal ---------------------------------------------
+
+    def __init__(self, stencil):
+        dim = len(stencil[0])
+
+        offset_sym = BoundaryOffsetInfo._offset_symbols(dim)
+        code = "\n"
+        for i in range(dim):
+            offset_str = ", ".join([str(d[i]) for d in stencil])
+            code += "const int32_t %s [] = { %s };\n" % (offset_sym[i].name, offset_str)
+
+        inv_dirs = []
+        for direction in stencil:
+            inverse_dir = tuple([-i for i in direction])
+            inv_dirs.append(str(stencil.index(inverse_dir)))
+
+        code += "const int32_t %s [] = { %s };\n" % (self.INV_DIR_SYMBOL.name, ", ".join(inv_dirs))
+        offset_symbols = BoundaryOffsetInfo._offset_symbols(dim)
+        super(BoundaryOffsetInfo, self).__init__(code, symbols_read=set(),
+                                                 symbols_defined=set(offset_symbols + [self.INV_DIR_SYMBOL]))
+
+    @staticmethod
+    def _offset_symbols(dim):
+        return [TypedSymbol(f"c{d}", create_type(np.int32)) for d in ['x', 'y', 'z'][:dim]]
+
+    INV_DIR_SYMBOL = TypedSymbol("invdir", np.int32)
+
+
+def create_boundary_kernel(field, index_field, stencil, boundary_functor, target=Target.CPU, **kernel_creation_args):
+    elements = [BoundaryOffsetInfo(stencil)]
+    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, skip_independence_check=True,
+                                **kernel_creation_args)
+    return create_kernel(elements, config=config)

src/pystencils/boundaries/createindexlist.py

+import warnings
+
+import numpy as np
+
+
+try:
+    import pyximport
+
+    pyximport.install(language_level=3)
+    cython_funcs_available = True
+except ImportError:
+    cython_funcs_available = False
+
+if cython_funcs_available:
+    from pystencils.boundaries.createindexlistcython import (
+        create_boundary_neighbor_index_list_2d,
+        create_boundary_neighbor_index_list_3d,
+        create_boundary_cell_index_list_2d,
+        create_boundary_cell_index_list_3d,
+    )
+
+boundary_index_array_coordinate_names = ["x", "y", "z"]
+direction_member_name = "dir"
+default_index_array_dtype = np.int32
+
+
+def numpy_data_type_for_boundary_object(boundary_object, dim):
+    coordinate_names = boundary_index_array_coordinate_names[:dim]
+    return np.dtype(
+        [(name, default_index_array_dtype) for name in coordinate_names]
+        + [(direction_member_name, default_index_array_dtype)]
+        + [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data],
+        align=True,
+    )
+
+
+def _create_index_list_python(
+    flag_field_arr,
+    boundary_mask,
+    fluid_mask,
+    stencil,
+    single_link,
+    inner_or_boundary=False,
+    nr_of_ghost_layers=None,
+):
+    if inner_or_boundary and nr_of_ghost_layers is None:
+        raise ValueError(
+            "If inner_or_boundary is set True the number of ghost layers "
+            "around the inner domain has to be specified"
+        )
+
+    if nr_of_ghost_layers is None:
+        nr_of_ghost_layers = 0
+
+    coordinate_names = boundary_index_array_coordinate_names[
+        : len(flag_field_arr.shape)
+    ]
+    index_arr_dtype = np.dtype(
+        [(name, default_index_array_dtype) for name in coordinate_names]
+        + [(direction_member_name, default_index_array_dtype)]
+    )
+
+    # boundary cells are extracted via np.where. To ensure continous memory access in the compute kernel these cells
+    # have to be sorted.
+    boundary_cells = np.transpose(np.nonzero(flag_field_arr == boundary_mask))
+    for i in range(len(flag_field_arr.shape)):
+        boundary_cells = boundary_cells[boundary_cells[:, i].argsort(kind="mergesort")]
+
+    # First a set is created to save all fluid cells which are near boundary
+    fluid_cells = set()
+    for cell in boundary_cells:
+        cell = tuple(cell)
+        for dir_idx, direction in enumerate(stencil):
+            neighbor_cell = tuple(
+                [cell_i + dir_i for cell_i, dir_i in zip(cell, direction)]
+            )
+            # prevent out ouf bounds access. If boundary cell is at the border, some stencil directions would be out.
+            if any(
+                not 0 + nr_of_ghost_layers <= e < upper - nr_of_ghost_layers
+                for e, upper in zip(neighbor_cell, flag_field_arr.shape)
+            ):
+                continue
+            if flag_field_arr[neighbor_cell] & fluid_mask:
+                fluid_cells.add(neighbor_cell)
+
+    # then this is set is transformed to a list to make it sortable. This ensures continoous memory access later.
+    fluid_cells = list(fluid_cells)
+    if len(flag_field_arr.shape) == 3:
+        fluid_cells.sort(key=lambda tup: (tup[-1], tup[-2], tup[0]))
+    else:
+        fluid_cells.sort(key=lambda tup: (tup[-1], tup[0]))
+
+    cells_to_iterate = fluid_cells if inner_or_boundary else boundary_cells
+    checkmask = boundary_mask if inner_or_boundary else fluid_mask
+
+    result = []
+    for cell in cells_to_iterate:
+        cell = tuple(cell)
+        sum_cells = np.zeros(len(cell))
+        for dir_idx, direction in enumerate(stencil):
+            neighbor_cell = tuple(
+                [cell_i + dir_i for cell_i, dir_i in zip(cell, direction)]
+            )
+            # prevent out ouf bounds access. If boundary cell is at the border, some stencil directions would be out.
+            if any(
+                not 0 <= e < upper
+                for e, upper in zip(neighbor_cell, flag_field_arr.shape)
+            ):
+                continue
+            if flag_field_arr[neighbor_cell] & checkmask:
+                if single_link:
+                    sum_cells += np.array(direction)
+                else:
+                    result.append(tuple(cell) + (dir_idx,))
+
+        # the discrete normal direction is the one which gives the maximum inner product to the stencil direction
+        if single_link and any(sum_cells != 0):
+            idx = np.argmax(np.inner(sum_cells, stencil))
+            result.append(tuple(cell) + (idx,))
+
+    return np.array(result, dtype=index_arr_dtype)
+
+
+def create_boundary_index_list(
+    flag_field,
+    stencil,
+    boundary_mask,
+    fluid_mask,
+    nr_of_ghost_layers=1,
+    inner_or_boundary=True,
+    single_link=False,
+):
+    """Creates a numpy array storing links (connections) between domain cells and boundary cells.
+
+    Args:
+        flag_field: flag integer array where boundary and domain cells are marked (interpreted as bit vector)
+        stencil: list of directions, for possible links. When single_link is set to true the order matters, because
+                 then only the first link is added to the list
+        boundary_mask: cells where (cell & mask) is true are considered boundary cells
+        fluid_mask: cells where (cell & mask) is true are considered fluid/inner cells cells
+        nr_of_ghost_layers: only relevant if neighbors is True
+        inner_or_boundary: if true, the result contains the cell coordinates of the domain cells -
+                    if false the boundary cells are listed
+        single_link: if true only the link in normal direction to this cell is reported
+
+    """
+    dim = len(flag_field.shape)
+    coordinate_names = boundary_index_array_coordinate_names[:dim]
+    index_arr_dtype = np.dtype(
+        [(name, default_index_array_dtype) for name in coordinate_names]
+        + [(direction_member_name, default_index_array_dtype)]
+    )
+
+    stencil = np.array(stencil, dtype=default_index_array_dtype)
+    args = (
+        flag_field,
+        nr_of_ghost_layers,
+        boundary_mask,
+        fluid_mask,
+        stencil,
+        single_link,
+    )
+    args_no_gl = (flag_field, boundary_mask, fluid_mask, stencil, single_link)
+
+    if cython_funcs_available:
+        if dim == 2:
+            if inner_or_boundary:
+                idx_list = create_boundary_neighbor_index_list_2d(*args)
+            else:
+                idx_list = create_boundary_cell_index_list_2d(*args_no_gl)
+        elif dim == 3:
+            if inner_or_boundary:
+                idx_list = create_boundary_neighbor_index_list_3d(*args)
+            else:
+                idx_list = create_boundary_cell_index_list_3d(*args_no_gl)
+        else:
+            raise ValueError("Flag field has to be a 2 or 3 dimensional numpy array")
+        return np.array(idx_list, dtype=index_arr_dtype)
+    else:
+        if flag_field.size > 1e6:
+            warnings.warn(
+                "Boundary setup may take very long! Consider installing cython to speed it up"
+            )
+        return _create_index_list_python(
+            *args_no_gl,
+            inner_or_boundary=inner_or_boundary,
+            nr_of_ghost_layers=nr_of_ghost_layers,
+        )
+
+
+def create_boundary_index_array(
+    flag_field,
+    stencil,
+    boundary_mask,
+    fluid_mask,
+    boundary_object,
+    nr_of_ghost_layers=1,
+    inner_or_boundary=True,
+    single_link=False,
+):
+    idx_array = create_boundary_index_list(
+        flag_field,
+        stencil,
+        boundary_mask,
+        fluid_mask,
+        nr_of_ghost_layers,
+        inner_or_boundary,
+        single_link,
+    )
+    dim = len(flag_field.shape)
+
+    if boundary_object.additional_data:
+        coordinate_names = boundary_index_array_coordinate_names[:dim]
+        index_arr_dtype = numpy_data_type_for_boundary_object(boundary_object, dim)
+        extended_idx_field = np.empty(len(idx_array), dtype=index_arr_dtype)
+        for prop in coordinate_names + ["dir"]:
+            extended_idx_field[prop] = idx_array[prop]
+
+        idx_array = extended_idx_field
+
+    return idx_array

src/pystencils/boundaries/createindexlistcython.pyx

+# cython: language_level=3str
+
+import cython
+
+ctypedef fused IntegerType:
+    short
+    int
+    long
+    long long
+    unsigned short
+    unsigned int
+    unsigned long
+
+@cython.boundscheck(False) # turn off bounds-checking for entire function
+@cython.wraparound(False)  # turn off negative index wrapping for entire function
+def create_boundary_neighbor_index_list_2d(object[IntegerType, ndim=2] flag_field,
+                                           int nr_of_ghost_layers, IntegerType boundary_mask, IntegerType fluid_mask,
+                                           object[int, ndim=2] stencil, int single_link):
+    cdef int xs, ys, x, y
+    cdef int dirIdx, num_directions, dx, dy
+
+    cdef int sum_x, sum_y
+    cdef float dot, maxn
+    cdef int calculated_idx
+
+    xs, ys = flag_field.shape
+    boundary_index_list = []
+    num_directions = stencil.shape[0]
+
+    for y in range(nr_of_ghost_layers, ys - nr_of_ghost_layers):
+        for x in range(nr_of_ghost_layers, xs - nr_of_ghost_layers):
+            sum_x = 0; sum_y = 0;
+            if flag_field[x, y] & fluid_mask:
+                for dirIdx in range(num_directions):
+                    dx = stencil[dirIdx,0]; dy = stencil[dirIdx,1]
+                    if flag_field[x + dx, y + dy] & boundary_mask:
+                        if single_link:
+                            sum_x += dx; sum_y += dy;
+                        else:
+                            boundary_index_list.append((x, y, dirIdx))
+
+            dot = 0; maxn = 0; calculated_idx = 0
+            if single_link and (sum_x != 0 or sum_y != 0):
+                for dirIdx in range(num_directions):
+                    dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1];
+
+                    dot = dx * sum_x + dy * sum_y
+                    if dot > maxn:
+                        maxn = dot
+                        calculated_idx = dirIdx
+
+                boundary_index_list.append((x, y, calculated_idx))
+    return boundary_index_list
+
+
+@cython.boundscheck(False) # turn off bounds-checking for entire function
+@cython.wraparound(False)  # turn off negative index wrapping for entire function
+def create_boundary_neighbor_index_list_3d(object[IntegerType, ndim=3] flag_field,
+                                           int nr_of_ghost_layers, IntegerType boundary_mask, IntegerType fluid_mask,
+                                           object[int, ndim=2] stencil, int single_link):
+    cdef int xs, ys, zs, x, y, z
+    cdef int dirIdx, num_directions, dx, dy, dz
+
+    cdef int sum_x, sum_y, sum_z
+    cdef float dot, maxn
+    cdef int calculated_idx
+
+    xs, ys, zs = flag_field.shape
+    boundary_index_list = []
+    num_directions = stencil.shape[0]
+
+    for z in range(nr_of_ghost_layers, zs - nr_of_ghost_layers):
+        for y in range(nr_of_ghost_layers, ys - nr_of_ghost_layers):
+            for x in range(nr_of_ghost_layers, xs - nr_of_ghost_layers):
+                sum_x = 0; sum_y = 0; sum_z = 0
+                if flag_field[x, y, z] & fluid_mask:
+                    for dirIdx in range(num_directions):
+                        dx = stencil[dirIdx,0]; dy = stencil[dirIdx,1]; dz = stencil[dirIdx,2]
+                        if flag_field[x + dx, y + dy, z + dz] & boundary_mask:
+                            if single_link:
+                                sum_x += dx; sum_y += dy; sum_z += dz
+                            else:
+                                boundary_index_list.append((x, y, z, dirIdx))
+
+                dot = 0; maxn = 0; calculated_idx = 0
+                if single_link and (sum_x != 0 or sum_y != 0 or sum_z != 0):
+                    for dirIdx in range(num_directions):
+                        dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1]; dz = stencil[dirIdx, 2]
+
+                        dot = dx * sum_x + dy * sum_y + dz * sum_z
+                        if dot > maxn:
+                            maxn = dot
+                            calculated_idx = dirIdx
+
+                    boundary_index_list.append((x, y, z, calculated_idx))
+    return boundary_index_list
+
+
+
+@cython.boundscheck(False) # turn off bounds-checking for entire function
+@cython.wraparound(False)  # turn off negative index wrapping for entire function
+def create_boundary_cell_index_list_2d(object[IntegerType, ndim=2] flag_field,
+                                       IntegerType boundary_mask, IntegerType fluid_mask,
+                                       object[int, ndim=2] stencil, int single_link):
+    cdef int xs, ys, x, y
+    cdef int dirIdx, num_directions, dx, dy
+
+    cdef int sum_x, sum_y
+    cdef float dot, maxn
+    cdef int calculated_idx
+
+    xs, ys = flag_field.shape
+    boundary_index_list = []
+    num_directions = stencil.shape[0]
+
+    for y in range(0, ys):
+        for x in range(0, xs):
+            sum_x = 0; sum_y = 0;
+            if flag_field[x, y] & boundary_mask:
+                for dirIdx in range(num_directions):
+                    dx = stencil[dirIdx,0]; dy = stencil[dirIdx,1]
+                    if 0 <= x + dx < xs and 0 <= y + dy < ys:
+                        if flag_field[x + dx, y + dy] & fluid_mask:
+                            if single_link:
+                                sum_x += dx; sum_y += dy
+                            else:
+                                boundary_index_list.append((x, y, dirIdx))
+
+            dot = 0; maxn = 0; calculated_idx = 0
+            if single_link and (sum_x != 0 or sum_y != 0):
+                for dirIdx in range(num_directions):
+                    dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1]
+
+                    dot = dx * sum_x + dy * sum_y
+                    if dot > maxn:
+                        maxn = dot
+                        calculated_idx = dirIdx
+
+                boundary_index_list.append((x, y, calculated_idx))
+
+    return boundary_index_list
+
+
+@cython.boundscheck(False) # turn off bounds-checking for entire function
+@cython.wraparound(False)  # turn off negative index wrapping for entire function
+def create_boundary_cell_index_list_3d(object[IntegerType, ndim=3] flag_field,
+                                       IntegerType boundary_mask, IntegerType fluid_mask,
+                                       object[int, ndim=2] stencil, int single_link):
+    cdef int xs, ys, zs, x, y, z
+    cdef int dirIdx, num_directions, dx, dy, dz
+
+    cdef int sum_x, sum_y, sum_z
+    cdef float dot, maxn
+    cdef int calculated_idx
+
+    xs, ys, zs = flag_field.shape
+    boundary_index_list = []
+    num_directions = stencil.shape[0]
+
+    for z in range(0, zs):
+        for y in range(0, ys):
+            for x in range(0, xs):
+                sum_x = 0; sum_y = 0; sum_z = 0
+                if flag_field[x, y, z] & boundary_mask:
+                    for dirIdx in range(num_directions):
+                        dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1]; dz = stencil[dirIdx, 2]
+                        if 0 <= x + dx < xs and 0 <= y + dy < ys and 0 <= z + dz < zs:
+                            if flag_field[x + dx, y + dy, z + dz] & fluid_mask:
+                                if single_link:
+                                    sum_x += dx; sum_y += dy; sum_z += dz
+                                else:
+                                    boundary_index_list.append((x, y, z, dirIdx))
+
+                dot = 0; maxn = 0; calculated_idx=0
+                if single_link and (sum_x != 0 or sum_y !=0 or sum_z !=0):
+                    for dirIdx in range(num_directions):
+                        dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1]; dz = stencil[dirIdx, 2]
+
+                        dot = dx*sum_x + dy*sum_y + dz*sum_z
+                        if dot > maxn:
+                            maxn = dot
+                            calculated_idx = dirIdx
+
+                    boundary_index_list.append((x, y, z, calculated_idx))
+
+    return boundary_index_list
\ No newline at end of file

src/pystencils/boundaries/inkernel.py

+import sympy as sp
+
+from pystencils.boundaries.boundaryhandling import DEFAULT_FLAG_TYPE
+from pystencils.typing import TypedSymbol, create_type
+from pystencils.field import Field
+from pystencils.integer_functions import bitwise_and
+
+
+def add_neumann_boundary(eqs, fields, flag_field, boundary_flag="neumann_flag", inverse_flag=False):
+    """
+    Replaces all neighbor accesses by flag field guarded accesses.
+    If flag in neighboring cell is set, the center value is used instead
+
+    Args:
+        eqs: list of equations containing field accesses to direct neighbors
+        fields: fields for which the Neumann boundary should be applied
+        flag_field: integer field marking boundary cells
+        boundary_flag: if flag field has value 'boundary_flag' (no bit operations yet)
+                       the cell is assumed to be boundary
+        inverse_flag: if true, boundary cells are where flag field has not the value of boundary_flag
+
+    Returns:
+        list of equations with guarded field accesses
+    """
+    if not hasattr(fields, "__len__"):
+        fields = [fields]
+    fields = set(fields)
+
+    if type(boundary_flag) is str:
+        boundary_flag = TypedSymbol(boundary_flag, dtype=create_type(DEFAULT_FLAG_TYPE))
+
+    substitutions = {}
+    for eq in eqs:
+        for fa in eq.atoms(Field.Access):
+            if fa.field not in fields:
+                continue
+            if not all(offset in (-1, 0, 1) for offset in fa.offsets):
+                raise ValueError("Works only for single neighborhood stencils")
+            if all(offset == 0 for offset in fa.offsets):
+                continue
+
+            if inverse_flag:
+                condition = sp.Eq(bitwise_and(flag_field[tuple(fa.offsets)], boundary_flag), 0)
+            else:
+                condition = sp.Ne(bitwise_and(flag_field[tuple(fa.offsets)], boundary_flag), 0)
+
+            center = fa.field(*fa.index)
+            substitutions[fa] = sp.Piecewise((center, condition), (fa, True))
+    return [eq.subs(substitutions) for eq in eqs]

src/pystencils/cache.py

+import os
+from collections.abc import Hashable
+from functools import partial, wraps
+from itertools import chain
+
+from functools import lru_cache as memorycache
+
+from joblib import Memory
+from appdirs import user_cache_dir
+
+if 'PYSTENCILS_CACHE_DIR' in os.environ:
+    cache_dir = os.environ['PYSTENCILS_CACHE_DIR']
+else:
+    cache_dir = user_cache_dir('pystencils')
+disk_cache = Memory(cache_dir, verbose=False).cache
+disk_cache_no_fallback = disk_cache
+
+
+def _wrapper(wrapped_func, cached_func, *args, **kwargs):
+    if all(isinstance(a, Hashable) for a in chain(args, kwargs.values())):
+        return cached_func(*args, **kwargs)
+    else:
+        return wrapped_func(*args, **kwargs)
+
+
+def memorycache_if_hashable(maxsize=128, typed=False):
+    def wrapper(func):
+        return partial(_wrapper, func, memorycache(maxsize, typed)(func))
+
+    return wrapper
+
+
+def sharedmethodcache(cache_id: str):
+    """Decorator for memoization of instance methods, allowing multiple methods to use the same cache.
+
+    This decorator caches results of instance methods per instantiated object of the surrounding class.
+    It allows multiple methods to use the same cache, by passing them the same `cache_id` string.
+    Cached values are stored in a dictionary, which is added as a member `self.<cache_id>` to the 
+    `self` object instance. Make sure that this doesn't cause any naming conflicts with other members!
+    Of course, for this to be useful, said methods must have the same signature (up to additional kwargs)
+    and must return the same result when called with the same arguments."""
+    def _decorator(user_method):
+        @wraps(user_method)
+        def _decorated_func(self, *args, **kwargs):
+            objdict = self.__dict__
+            cache = objdict.setdefault(cache_id, dict())
+
+            key = args
+            for item in kwargs.items():
+                key += item
+
+            if key not in cache:
+                result = user_method(self, *args, **kwargs)
+                cache[key] = result
+                return result
+            else:
+                return cache[key]
+        return _decorated_func
+    return _decorator
+
+
+def clear_cache():
+    """
+    Clears the pystencils cache created by joblib.
+    """
+    memory = Memory(cache_dir, verbose=0)
+    memory.clear(warn=False)
+
+
+# Disable memory cache:
+# disk_cache = lambda o: o
+# disk_cache_no_fallback = lambda o: o

src/pystencils/config.py

+from copy import copy
+from collections import defaultdict
+from dataclasses import dataclass, field
+from types import MappingProxyType
+from typing import Union, Tuple, List, Dict, Callable, Any, DefaultDict, Iterable
+
+from pystencils import Target, Backend, Field
+from pystencils.typing.typed_sympy import BasicType
+from pystencils.typing.utilities import collate_types
+
+import numpy as np
+
+# TODO: There exists DTypeLike in NumPy which would be better than type for type hinting, to new at the moment
+# from numpy.typing import DTypeLike
+
+
+# TODO: CreateKernelConfig is bloated think of more classes better usage, factory whatever ...
+# Proposition: CreateKernelConfigs Classes for different targets?
+@dataclass
+class CreateKernelConfig:
+    """
+    **Below all parameters for the CreateKernelConfig are explained**
+    """
+    target: Target = Target.CPU
+    """
+    All targets are defined in :class:`pystencils.enums.Target`
+    """
+    backend: Backend = None
+    """
+    All backends are defined in :class:`pystencils.enums.Backend`
+    """
+    function_name: str = 'kernel'
+    """
+    Name of the generated function - only important if generated code is written out
+    """
+    data_type: Union[type, str, DefaultDict[str, BasicType], Dict[str, BasicType]] = np.float64
+    """
+    Data type used for all untyped symbols (i.e. non-fields), can also be a dict from symbol name to type.
+    If specified as a dict ideally a defaultdict is used to define a default value for symbols not listed in the
+    dict. If a plain dict is provided it will be transformed into a defaultdict internally. The default value 
+    will then be specified via type collation then.
+    """
+    default_number_float: Union[type, str, BasicType] = None
+    """
+    Data type used for all untyped floating point numbers (i.e. 0.5). By default the value of data_type is used.
+    If data_type is given as a defaultdict its default_factory is used.
+    """
+    default_number_int: Union[type, str, BasicType] = np.int64
+    """
+    Data type used for all untyped integer numbers (i.e. 1)
+    """
+    iteration_slice: Tuple = None
+    """
+    Rectangular subset to iterate over, if not specified the complete non-ghost layer part of the field is iterated over
+    """
+    ghost_layers: Union[bool, int, List[Tuple[int]]] = None
+    """
+    A single integer specifies the ghost layer count at all borders, can also be a sequence of
+    pairs ``[(x_lower_gl, x_upper_gl), .... ]``. These layers are excluded from the iteration.
+    If left to default, the number of ghost layers is determined automatically from the assignments.
+    """
+    cpu_openmp: Union[bool, int] = False
+    """
+    `True` or number of threads for OpenMP parallelization, `False` for no OpenMP. If set to `True`, the maximum number
+    of available threads will be chosen.
+    """
+    cpu_vectorize_info: Dict = None
+    """
+    A dictionary with keys, 'vector_instruction_set', 'assume_aligned' and 'nontemporal'
+    for documentation of these parameters see vectorize function. Example:
+    '{'instruction_set': 'avx512', 'assume_aligned': True, 'nontemporal':True}'
+    """
+    cpu_blocking: Tuple[int] = None
+    """
+    A tuple of block sizes or `None` if no blocking should be applied
+    """
+    omp_single_loop: bool = True
+    """
+    If OpenMP is active: whether multiple outer loops are permitted
+    """
+    base_pointer_specification: Union[List[Iterable[str]], List[Iterable[int]]] = None
+    """
+    Specification of how many and which intermediate pointers are created for a field access.
+    For example [ (0), (2,3,)]  creates on base pointer for coordinates 2 and 3 and writes the offset for coordinate
+    zero directly in the field access. These specifications are defined dependent on the loop ordering.
+    This function translates more readable version into the specification above.
+    
+    For more information see: `pystencils.transformations.create_intermediate_base_pointer`
+    """
+    gpu_indexing: str = 'block'
+    """
+    Either 'block' or 'line' , or custom indexing class, see `pystencils.gpu.AbstractIndexing`
+    """
+    gpu_indexing_params: MappingProxyType = field(default_factory=lambda: MappingProxyType({}))
+    """
+    Dict with indexing parameters (constructor parameters of indexing class)
+    e.g. for 'block' one can specify '{'block_size': (20, 20, 10) }'.
+    """
+    # TODO Markus rework this docstring
+    default_assignment_simplifications: bool = False
+    """
+    If `True` default simplifications are first performed on the Assignments. If problems occur during the
+    simplification a warning will be thrown.
+    Furthermore, it is essential to know that this is a two-stage process. The first stage of the process acts
+    on the level of the `pystencils.AssignmentCollection`.  In this part,
+    `pystencil.simp.create_simplification_strategy` from pystencils.simplificationfactory will be used to
+    apply optimisations like insertion of constants to
+    remove pressure from the registers. Thus the first part of the optimisations can only be executed if
+    an `AssignmentCollection` is passed. The second part of the optimisation acts on the level of each Assignment
+    individually. In this stage, all optimisations from `sympy.codegen.rewriting.optims_c99` are applied
+    to each Assignment. Thus this stage can also be applied if a list of Assignments is passed.
+    """
+    cpu_prepend_optimizations: List[Callable] = field(default_factory=list)
+    """
+    List of extra optimizations to perform first on the AST.
+    """
+    use_auto_for_assignments: bool = False
+    """
+    If set to `True`, auto can be used in the generated code for data types. This makes the type system more robust.
+    """
+    index_fields: List[Field] = None
+    """
+    List of index fields, i.e. 1D fields with struct data type. If not `None`, `create_index_kernel`
+    instead of `create_domain_kernel` is used.
+    """
+    coordinate_names: Tuple[str, Any] = ('x', 'y', 'z')
+    """
+    Name of the coordinate fields in the struct data type.
+    """
+    allow_double_writes: bool = False
+    """
+    If True, don't check if every field is only written at a single location. This is required
+    for example for kernels that are compiled with loop step sizes > 1, that handle multiple
+    cells at once. Use with care!
+    """
+    skip_independence_check: bool = False
+    """
+    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:
+        """Because of pickle, we need to have a nested class, instead of a lambda in __post_init__"""
+        def __init__(self, dt):
+            self.dt = dt
+
+        def __call__(self):
+            return BasicType(self.dt)
+
+    def _check_type(self, dtype_to_check):
+        if isinstance(dtype_to_check, str) and (dtype_to_check == 'float' or dtype_to_check == 'int'):
+            self._typing_error()
+
+        if isinstance(dtype_to_check, type) and not hasattr(dtype_to_check, "dtype"):
+            # NumPy-types are also of type 'type'. However, they have more properties
+            self._typing_error()
+
+    @staticmethod
+    def _typing_error():
+        raise ValueError("It is not possible to use python types (float, int) for datatypes because these "
+                         "types are ambiguous. For example float will map to double. "
+                         "Also the string version like 'float' is not allowed, e.g. use 'float64' instead")
+
+    def __post_init__(self):
+        # ----  Legacy parameters
+        if not isinstance(self.target, Target):
+            raise ValueError("target must be provided by the 'Target' enum")
+
+        # ---- Auto Backend
+        if not self.backend:
+            if self.target == Target.CPU:
+                self.backend = Backend.C
+            elif self.target == Target.GPU:
+                self.backend = Backend.CUDA
+            else:
+                raise NotImplementedError(f'Target {self.target} has no default backend')
+
+        if not isinstance(self.backend, Backend):
+            raise ValueError("backend must be provided by the 'Backend' enum")
+
+        # Normalise data types
+        for dtype in [self.data_type, self.default_number_float, self.default_number_int]:
+            self._check_type(dtype)
+
+        if not isinstance(self.data_type, dict):
+            dt = copy(self.data_type)  # The copy is necessary because BasicType has sympy shinanigans
+            self.data_type = defaultdict(self.DataTypeFactory(dt))
+
+        if isinstance(self.data_type, dict) and not isinstance(self.data_type, defaultdict):
+            for dtype in self.data_type.values():
+                self._check_type(dtype)
+
+            dt = collate_types([BasicType(dtype) for dtype in self.data_type.values()])
+            dtype_dict = self.data_type
+            self.data_type = defaultdict(self.DataTypeFactory(dt), dtype_dict)
+
+        assert isinstance(self.data_type, defaultdict), "At this point data_type must be a defaultdict!"
+        for dtype in self.data_type.values():
+            self._check_type(dtype)
+        self._check_type(self.data_type.default_factory())
+
+        if self.default_number_float is None:
+            self.default_number_float = self.data_type.default_factory()
+
+        if not isinstance(self.default_number_float, BasicType):
+            self.default_number_float = BasicType(self.default_number_float)
+        if not isinstance(self.default_number_int, BasicType):
+            self.default_number_int = BasicType(self.default_number_int)

src/pystencils/cpu/__init__.py

+from pystencils.cpu.cpujit import make_python_function
+from pystencils.cpu.kernelcreation import add_openmp, create_indexed_kernel, create_kernel, add_pragmas
+
+__all__ = ['create_kernel', 'create_indexed_kernel', 'add_openmp', 'add_pragmas', 'make_python_function']