Merge branch 'kohl/moving-loop-strategy-to-kernel-type' into 'main'

Multiple kernels per method

See merge request !14

generate_all_operators.py

@@ -53,10 +53,10 @@ from hog.function_space import (
 )
 from hog.logger import get_logger, TimedLogger
 from hog.operator_generation.kernel_types import (
-    Apply,
-    Assemble,
-    AssembleDiagonal,
-    KernelType,
+    ApplyWrapper,
+    AssembleWrapper,
+    AssembleDiagonalWrapper,
+    KernelWrapperType,
 )
 from hog.operator_generation.loop_strategies import (
     LoopStrategy,
@@ -354,7 +354,7 @@ def main():
             args.quad_degree if args.quad_rule is None else args.quad_rule,
         )
     }
-    
+
     enabled_geometries: Set[TriangleElement | TetrahedronElement] = set()
     if 2 in args.dimensions:
         enabled_geometries.add(TriangleElement())
@@ -488,7 +488,7 @@ class OperatorInfo:
         sp.Matrix,
     ]
     type_descriptor: HOGType
-    kernel_types: List[KernelType] = None  # type: ignore[assignment] # will definitely be initialized in __post_init__
+    kernel_types: List[KernelWrapperType] = None  # type: ignore[assignment] # will definitely be initialized in __post_init__
     geometries: Sequence[ElementGeometry] = field(
         default_factory=lambda: [TriangleElement(), TetrahedronElement()]
     )
@@ -507,7 +507,7 @@ class OperatorInfo:
         if self.kernel_types is None:
             dims = [g.dimensions for g in self.geometries]
             self.kernel_types = [
-                Apply(
+                ApplyWrapper(
                     self.test_space,
                     self.trial_space,
                     type_descriptor=self.type_descriptor,
@@ -518,7 +518,7 @@ class OperatorInfo:
             all_opts = set().union(*[o for (o, _, _) in self.opts])
             if not ({Opts.VECTORIZE, Opts.VECTORIZE512}.intersection(all_opts)):
                 self.kernel_types.append(
-                    Assemble(
+                    AssembleWrapper(
                         self.test_space,
                         self.trial_space,
                         type_descriptor=self.type_descriptor,
@@ -528,8 +528,10 @@ class OperatorInfo:
 
             if self.test_space == self.trial_space:
                 self.kernel_types.append(
-                    AssembleDiagonal(
-                        self.test_space, type_descriptor=self.type_descriptor, dims=dims
+                    AssembleDiagonalWrapper(
+                        self.test_space,
+                        type_descriptor=self.type_descriptor,
+                        dims=dims,
                     )
                 )
 
@@ -591,7 +593,7 @@ def all_operators(
     ops.append(OperatorInfo(mapping="P2", name="SUPGDiffusion", trial_space=P2, test_space=P2, 
                             form=partial(supg_diffusion, velocity_function_space=P2, diffusivityXdelta_function_space=P2), 
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
-    
+
     # fmt: on
 
     p2vec_epsilon = partial(
@@ -703,7 +705,7 @@ def generate_elementwise_op(
         name,
         symbolizer,
         opts=optimizations,
-        kernel_types=op_info.kernel_types,
+        kernel_wrapper_types=op_info.kernel_types,
         type_descriptor=type_descriptor,
     )
 
@@ -725,19 +727,20 @@ def generate_elementwise_op(
             blending=blending,  # type: ignore[call-arg] # kw-args are not supported by Callable
         )
 
-        operator.set_element_matrix(
+        operator.add_integral(
+            name="".join(name.split()),
             dim=geometry.dimensions,
             geometry=geometry,
             integration_domain=MacroIntegrationDomain.VOLUME,
             quad=quad,
             blending=blending,
             form=form,
+            loop_strategy=loop_strategy,
         )
 
     dir_path = os.path.join(args.output, op_info.name.split("_")[0])
     operator.generate_class_code(
         dir_path,
-        loop_strategy=loop_strategy,
         clang_format_binary=args.clang_format_binary,
     )
 

hog/operator_generation/kernel_types.py

@@ -42,16 +42,54 @@ from hog.operator_generation.function_space_impls import FunctionSpaceImpl
 from hog.operator_generation.indexing import FaceType, CellType
 from hog.operator_generation.pystencils_extensions import create_generic_fields
 from hog.operator_generation.types import HOGPrecision, HOGType, hyteg_type
+from hog.operator_generation.loop_strategies import LoopStrategy, SAWTOOTH
 
 
 class KernelType(ABC):
-    name: str
-    src_fields: List[FunctionSpaceImpl]
-    dst_fields: List[FunctionSpaceImpl]
-    _template: Template
-
-    """Type to specify the type of kernel.
-    E.g. a GEMv kernel defines multiple source fields and scalar parameters."""
+    """
+    A HyTeG operator may implement multiple types of methods that perform some kind of operation on a function that
+    is inferred from the bilinear form. For instance, a (matrix-free) matrix-vector multiplication or the assembly
+    of its diagonal.
+
+    Certain operations such as setting boundary data to zero or communication have to be executed before the actual
+    kernel can be executed.
+
+    E.g.:
+
+    ```
+    void apply() {
+        communication()       // "pre-processing"
+        kernel()              // actual kernel
+        post_communication()  // "post-processing"
+    }
+    ```
+
+    For some applications, it might be necessary or at least comfortable to execute multiple kernels inside such a
+    method. For instance, if boundary conditions are applied:
+
+    ```
+    // form is sum of volume and boundary integral:
+    //     ∫ ... dΩ + ∫ ... dS
+    void assemble() {
+        communication()       // "pre-processing"
+        kernel()              // volume kernel (∫ ... dΩ)
+        kernel_boundary()     // boundary kernel (∫ ... dS | additive execution)
+        post_communication()  // "post-processing"
+    }
+    ```
+
+    This class (KernelType) describes the "action" of the kernel (matvec, assembly, ...).
+    Another class (KernelWrapperType) then describes what kind of pre- and post-processing is required.
+
+    ```
+    void assemble() {         // from KernelTypeWrapper
+        communication()       // from KernelTypeWrapper
+        kernel()              // from KernelType + IntegrationInfo 1
+        kernel_boundary()     // from KernelType + IntegrationInfo 2
+        post_communication()  // from KernelTypeWrapper
+    }
+    ```
+    """
 
     @abstractmethod
     def kernel_operation(
@@ -90,27 +128,243 @@ class KernelType(ABC):
         """
         ...
 
+
+class Apply(KernelType):
+    def __init__(self):
+        self.result_prefix = "elMatVec_"
+
+    def kernel_operation(
+        self,
+        src_vecs: List[sp.MatrixBase],
+        dst_vecs: List[sp.MatrixBase],
+        mat: sp.MatrixBase,
+        rows: int,
+    ) -> List[SympyAssignment]:
+        kernel_ops = mat * src_vecs[0]
+
+        tmp_symbols = sp.numbered_symbols(self.result_prefix)
+
+        kernel_op_assignments = [
+            SympyAssignment(tmp, kernel_op)
+            for tmp, kernel_op in zip(tmp_symbols, kernel_ops)
+        ]
+
+        return kernel_op_assignments
+
+    def kernel_post_operation(
+        self,
+        geometry: ElementGeometry,
+        element_index: Tuple[int, int, int],
+        element_type: Union[FaceType, CellType],
+        src_vecs_accesses: List[List[Field.Access]],
+        dst_vecs_accesses: List[List[Field.Access]],
+    ) -> List[ps.astnodes.Node]:
+        tmp_symbols = sp.numbered_symbols(self.result_prefix)
+
+        # Add and store result to destination.
+        store_dst_vecs = [
+            SympyAssignment(a, a + s) for a, s in zip(dst_vecs_accesses[0], tmp_symbols)
+        ]
+        return store_dst_vecs
+
+
+class AssembleDiagonal(KernelType):
+    def __init__(self):
+        self.result_prefix = "elMatDiag_"
+
+    def kernel_operation(
+        self,
+        src_vecs: List[sp.MatrixBase],
+        dst_vecs: List[sp.MatrixBase],
+        mat: sp.MatrixBase,
+        rows: int,
+    ) -> List[SympyAssignment]:
+        kernel_ops = mat.diagonal()
+
+        tmp_symbols = sp.numbered_symbols(self.result_prefix)
+
+        kernel_op_assignments = [
+            SympyAssignment(tmp, kernel_op)
+            for tmp, kernel_op in zip(tmp_symbols, kernel_ops)
+        ]
+
+        return kernel_op_assignments
+
+    def kernel_post_operation(
+        self,
+        geometry: ElementGeometry,
+        element_index: Tuple[int, int, int],
+        element_type: Union[FaceType, CellType],
+        src_vecs_accesses: List[List[Field.Access]],
+        dst_vecs_accesses: List[List[Field.Access]],
+    ) -> List[ps.astnodes.Node]:
+        tmp_symbols = sp.numbered_symbols(self.result_prefix)
+
+        # Add and store result to destination.
+        store_dst_vecs = [
+            SympyAssignment(a, a + s) for a, s in zip(dst_vecs_accesses[0], tmp_symbols)
+        ]
+        return store_dst_vecs
+
+
+class Assemble(KernelType):
+    def __init__(
+        self,
+        src_space: FunctionSpace,
+        dst_space: FunctionSpace,
+    ):
+        self.result_prefix = "elMat_"
+        idx_t = HOGType("idx_t", np.int64)
+
+        self.src: FunctionSpaceImpl = FunctionSpaceImpl.create_impl(
+            src_space, "src", idx_t
+        )
+        self.dst: FunctionSpaceImpl = FunctionSpaceImpl.create_impl(
+            dst_space, "dst", idx_t
+        )
+
+    def kernel_operation(
+        self,
+        src_vecs: List[sp.MatrixBase],
+        dst_vecs: List[sp.MatrixBase],
+        mat: sp.MatrixBase,
+        rows: int,
+    ) -> List[SympyAssignment]:
+        return [
+            SympyAssignment(sp.Symbol(f"{self.result_prefix}{r}_{c}"), mat[r, c])
+            for r in range(mat.shape[0])
+            for c in range(mat.shape[1])
+        ]
+
+    def kernel_post_operation(
+        self,
+        geometry: ElementGeometry,
+        element_index: Tuple[int, int, int],
+        element_type: Union[FaceType, CellType],
+        src_vecs_accesses: List[List[Field.Access]],
+        dst_vecs_accesses: List[List[Field.Access]],
+    ) -> List[ps.astnodes.Node]:
+        src, dst = dst_vecs_accesses
+        el_mat = sp.Matrix(
+            [
+                [sp.Symbol(f"{self.result_prefix}{r}_{c}") for c in range(len(src))]
+                for r in range(len(dst))
+            ]
+        )
+
+        # apply basis/dof transformations
+        transform_src_code, transform_src_mat = self.src.dof_transformation(
+            geometry, element_index, element_type
+        )
+        transform_dst_code, transform_dst_mat = self.dst.dof_transformation(
+            geometry, element_index, element_type
+        )
+        transformed_el_mat = transform_dst_mat.T * el_mat * transform_src_mat
+
+        nr = len(dst)
+        nc = len(src)
+        mat_size = nr * nc
+
+        rowIdx, colIdx = create_generic_fields(["rowIdx", "colIdx"], dtype=np.uint64)
+        # NOTE: The type is 'hyteg_type().pystencils_type', i.e. `real_t`
+        #       (not `self._type_descriptor.pystencils_type`)
+        #       because this function is implemented manually in HyTeG with
+        #       this signature. Passing `np.float64` is not ideal (if `real_t !=
+        #       double`) but it makes sure that casts are inserted if necessary
+        #       (though some might be superfluous).
+        mat = create_generic_fields(
+            ["mat"], dtype=hyteg_type(HOGPrecision.REAL_T).pystencils_type
+        )[0]
+        rowIdxSymb = FieldPointerSymbol(rowIdx.name, rowIdx.dtype, False)
+        colIdxSymb = FieldPointerSymbol(colIdx.name, colIdx.dtype, False)
+        matSymb = FieldPointerSymbol(mat.name, mat.dtype, False)
+
+        body: List[ps.astnodes.Node] = [
+            CustomCodeNode(
+                f"std::vector< uint_t > {rowIdxSymb.name}( {nr} );\n"
+                f"std::vector< uint_t > {colIdxSymb.name}( {nc} );\n"
+                f"std::vector< {mat.dtype} > {matSymb.name}( {mat_size} );",
+                [],
+                [rowIdxSymb, colIdxSymb, matSymb],
+            ),
+            ps.astnodes.EmptyLine(),
+        ]
+        body += [
+            SympyAssignment(
+                rowIdx.absolute_access((k,), (0,)), CastFunc(dst_access, np.uint64)
+            )
+            for k, dst_access in enumerate(dst)
+        ]
+        body += [
+            SympyAssignment(
+                colIdx.absolute_access((k,), (0,)), CastFunc(src_access, np.uint64)
+            )
+            for k, src_access in enumerate(src)
+        ]
+
+        body += [
+            ps.astnodes.EmptyLine(),
+            ps.astnodes.SourceCodeComment("Apply basis transformation"),
+            *sorted({transform_src_code, transform_dst_code}, key=lambda n: n._code),
+            ps.astnodes.EmptyLine(),
+        ]
+
+        body += [
+            SympyAssignment(
+                mat.absolute_access((r * nc + c,), (0,)),
+                CastFunc(transformed_el_mat[r, c], mat.dtype),
+            )
+            for r in range(nr)
+            for c in range(nc)
+        ]
+
+        body += [
+            ps.astnodes.EmptyLine(),
+            CustomCodeNode(
+                f"mat->addValues( {rowIdxSymb.name}, {colIdxSymb.name}, {matSymb.name} );",
+                [
+                    TypedSymbol("mat", "std::shared_ptr< SparseMatrixProxy >"),
+                    rowIdxSymb,
+                    colIdxSymb,
+                    matSymb,
+                ],
+                [],
+            ),
+        ]
+
+        return body
+
+
+class KernelWrapperType(ABC):
+    name: str
+    src_fields: List[FunctionSpaceImpl]
+    dst_fields: List[FunctionSpaceImpl]
+    _template: Template
+    """
+    See documentation of class KernelType. 
+    """
+
+    @property
     @abstractmethod
-    def includes(self) -> Set[str]:
-        ...
+    def kernel_type(self) -> KernelType: ...
 
     @abstractmethod
-    def base_classes(self) -> List[str]:
-        ...
+    def includes(self) -> Set[str]: ...
 
     @abstractmethod
-    def wrapper_methods(self) -> List[CppMethod]:
-        ...
+    def base_classes(self) -> List[str]: ...
 
     @abstractmethod
-    def member_variables(self) -> List[CppMemberVariable]:
-        ...
+    def wrapper_methods(self) -> List[CppMethod]: ...
+
+    @abstractmethod
+    def member_variables(self) -> List[CppMemberVariable]: ...
 
     def substitute(self, subs: Mapping[str, object]) -> None:
         self._template = Template(self._template.safe_substitute(subs))
 
 
-class Apply(KernelType):
+class ApplyWrapper(KernelWrapperType):
     def __init__(
         self,
         src_space: FunctionSpace,
@@ -224,39 +478,9 @@ class Apply(KernelType):
             f'this->stopTiming( "{self.name}" );'
         )
 
-    def kernel_operation(
-        self,
-        src_vecs: List[sp.MatrixBase],
-        dst_vecs: List[sp.MatrixBase],
-        mat: sp.MatrixBase,
-        rows: int,
-    ) -> List[SympyAssignment]:
-        kernel_ops = mat * src_vecs[0]
-
-        tmp_symbols = sp.numbered_symbols(self.result_prefix)
-
-        kernel_op_assignments = [
-            SympyAssignment(tmp, kernel_op)
-            for tmp, kernel_op in zip(tmp_symbols, kernel_ops)
-        ]
-
-        return kernel_op_assignments
-
-    def kernel_post_operation(
-        self,
-        geometry: ElementGeometry,
-        element_index: Tuple[int, int, int],
-        element_type: Union[FaceType, CellType],
-        src_vecs_accesses: List[List[Field.Access]],
-        dst_vecs_accesses: List[List[Field.Access]],
-    ) -> List[ps.astnodes.Node]:
-        tmp_symbols = sp.numbered_symbols(self.result_prefix)
-
-        # Add and store result to destination.
-        store_dst_vecs = [
-            SympyAssignment(a, a + s) for a, s in zip(dst_vecs_accesses[0], tmp_symbols)
-        ]
-        return store_dst_vecs
+    @property
+    def kernel_type(self) -> KernelType:
+        return Apply()
 
     def includes(self) -> Set[str]:
         return (
@@ -302,7 +526,7 @@ class Apply(KernelType):
         return []
 
 
-class GEMV(KernelType):
+class GEMVWrapper(KernelWrapperType):
     def __init__(
         self,
         type_descriptor: HOGType,
@@ -325,7 +549,7 @@ class GEMV(KernelType):
         raise HOGException("Not implemented")
 
 
-class AssembleDiagonal(KernelType):
+class AssembleDiagonalWrapper(KernelWrapperType):
     def __init__(
         self,
         fe_space: FunctionSpace,
@@ -340,7 +564,6 @@ class AssembleDiagonal(KernelType):
         self.src_fields = []
         self.dst_fields = [self.dst]
         self.dims = dims
-        self.result_prefix = "elMatDiag_"
 
         def macro_loop(dim: int) -> str:
             Macro = {2: "Face", 3: "Cell"}[dim]
@@ -413,39 +636,9 @@ class AssembleDiagonal(KernelType):
             f'this->stopTiming( "{self.name}" );'
         )
 
-    def kernel_operation(
-        self,
-        src_vecs: List[sp.MatrixBase],
-        dst_vecs: List[sp.MatrixBase],
-        mat: sp.MatrixBase,
-        rows: int,
-    ) -> List[SympyAssignment]:
-        kernel_ops = mat.diagonal()
-
-        tmp_symbols = sp.numbered_symbols(self.result_prefix)
-
-        kernel_op_assignments = [
-            SympyAssignment(tmp, kernel_op)
-            for tmp, kernel_op in zip(tmp_symbols, kernel_ops)
-        ]
-
-        return kernel_op_assignments
-
-    def kernel_post_operation(
-        self,
-        geometry: ElementGeometry,
-        element_index: Tuple[int, int, int],
-        element_type: Union[FaceType, CellType],
-        src_vecs_accesses: List[List[Field.Access]],
-        dst_vecs_accesses: List[List[Field.Access]],
-    ) -> List[ps.astnodes.Node]:
-        tmp_symbols = sp.numbered_symbols(self.result_prefix)
-
-        # Add and store result to destination.
-        store_dst_vecs = [
-            SympyAssignment(a, a + s) for a, s in zip(dst_vecs_accesses[0], tmp_symbols)
-        ]
-        return store_dst_vecs
+    @property
+    def kernel_type(self) -> KernelType:
+        return AssembleDiagonal()
 
     def includes(self) -> Set[str]:
         return {"hyteg/solvers/Smoothables.hpp"} | self.dst.includes()
@@ -483,7 +676,7 @@ class AssembleDiagonal(KernelType):
         ]
 
 
-class Assemble(KernelType):
+class AssembleWrapper(KernelWrapperType):
     def __init__(
         self,
         src_space: FunctionSpace,
@@ -508,7 +701,6 @@ class Assemble(KernelType):
 
         self.type_descriptor = type_descriptor
         self.dims = dims
-        self.result_prefix = "elMat_"
 
         def macro_loop(dim: int) -> str:
             Macro = {2: "Face", 3: "Cell"}[dim]
@@ -563,116 +755,9 @@ class Assemble(KernelType):
             f'this->stopTiming( "{self.name}" );'
         )
 
-    def kernel_operation(
-        self,
-        src_vecs: List[sp.MatrixBase],
-        dst_vecs: List[sp.MatrixBase],
-        mat: sp.MatrixBase,
-        rows: int,
-    ) -> List[SympyAssignment]:
-        return [
-            SympyAssignment(sp.Symbol(f"{self.result_prefix}{r}_{c}"), mat[r, c])
-            for r in range(mat.shape[0])
-            for c in range(mat.shape[1])
-        ]
-
-    def kernel_post_operation(
-        self,
-        geometry: ElementGeometry,
-        element_index: Tuple[int, int, int],
-        element_type: Union[FaceType, CellType],
-        src_vecs_accesses: List[List[Field.Access]],
-        dst_vecs_accesses: List[List[Field.Access]],
-    ) -> List[ps.astnodes.Node]:
-        src, dst = dst_vecs_accesses
-        el_mat = sp.Matrix(
-            [
-                [sp.Symbol(f"{self.result_prefix}{r}_{c}") for c in range(len(src))]
-                for r in range(len(dst))
-            ]
-        )
-
-        # apply basis/dof transformations
-        transform_src_code, transform_src_mat = self.src.dof_transformation(
-            geometry, element_index, element_type
-        )
-        transform_dst_code, transform_dst_mat = self.dst.dof_transformation(
-            geometry, element_index, element_type
-        )
-        transformed_el_mat = transform_dst_mat.T * el_mat * transform_src_mat
-
-        nr = len(dst)
-        nc = len(src)
-        mat_size = nr * nc
-
-        rowIdx, colIdx = create_generic_fields(["rowIdx", "colIdx"], dtype=np.uint64)
-        # NOTE: The type is 'hyteg_type().pystencils_type', i.e. `real_t`
-        #       (not `self._type_descriptor.pystencils_type`)
-        #       because this function is implemented manually in HyTeG with
-        #       this signature. Passing `np.float64` is not ideal (if `real_t !=
-        #       double`) but it makes sure that casts are inserted if necessary
-        #       (though some might be superfluous).
-        mat = create_generic_fields(
-            ["mat"], dtype=hyteg_type(HOGPrecision.REAL_T).pystencils_type
-        )[0]
-        rowIdxSymb = FieldPointerSymbol(rowIdx.name, rowIdx.dtype, False)
-        colIdxSymb = FieldPointerSymbol(colIdx.name, colIdx.dtype, False)
-        matSymb = FieldPointerSymbol(mat.name, mat.dtype, False)
-
-        body: List[ps.astnodes.Node] = [
-            CustomCodeNode(
-                f"std::vector< uint_t > {rowIdxSymb.name}( {nr} );\n"
-                f"std::vector< uint_t > {colIdxSymb.name}( {nc} );\n"
-                f"std::vector< {mat.dtype} > {matSymb.name}( {mat_size} );",
-                [],
-                [rowIdxSymb, colIdxSymb, matSymb],
-            ),
-            ps.astnodes.EmptyLine(),
-        ]
-        body += [
-            SympyAssignment(
-                rowIdx.absolute_access((k,), (0,)), CastFunc(dst_access, np.uint64)
-            )
-            for k, dst_access in enumerate(dst)
-        ]
-        body += [
-            SympyAssignment(
-                colIdx.absolute_access((k,), (0,)), CastFunc(src_access, np.uint64)
-            )
-            for k, src_access in enumerate(src)
-        ]
-
-        body += [
-            ps.astnodes.EmptyLine(),
-            ps.astnodes.SourceCodeComment("Apply basis transformation"),
-            *sorted({transform_src_code, transform_dst_code}, key=lambda n: n._code),
-            ps.astnodes.EmptyLine(),
-        ]
-
-        body += [
-            SympyAssignment(
-                mat.absolute_access((r * nc + c,), (0,)),
-                CastFunc(transformed_el_mat[r, c], mat.dtype),
-            )
-            for r in range(nr)
-            for c in range(nc)
-        ]
-
-        body += [
-            ps.astnodes.EmptyLine(),
-            CustomCodeNode(
-                f"mat->addValues( {rowIdxSymb.name}, {colIdxSymb.name}, {matSymb.name} );",
-                [
-                    TypedSymbol("mat", "std::shared_ptr< SparseMatrixProxy >"),
-                    rowIdxSymb,
-                    colIdxSymb,
-                    matSymb,
-                ],
-                [],
-            ),
-        ]
-
-        return body
+    @property
+    def kernel_type(self) -> KernelType:
+        return Assemble(self.src.fe_space, self.dst.fe_space)
 
     def includes(self) -> Set[str]:
         return (

hog/operator_generation/optimizer.py

@@ -108,8 +108,8 @@ class Optimizer:
     def __getitem__(self, opt):
         return opt in self._opts
 
-    def check_opts_validity(self, loop_strategy: LoopStrategy) -> None:
-        """Checks if the desired optimizations are valid for the given loop strategy."""
+    def check_opts_validity(self) -> None:
+        """Checks if the desired optimizations are valid."""
 
         if Opts.VECTORIZE512 in self._opts and not Opts.VECTORIZE in self._opts:
             raise HOGException("Optimization VECTORIZE512 requires VECTORIZE.")
@@ -191,6 +191,7 @@ class Optimizer:
 
             with TimedLogger("simplifying conditionals", logging.DEBUG):
                 simplify_conditionals(loop, loop_counter_simplification=True)
+
         if self[Opts.MOVECONSTANTS]:
             with TimedLogger("moving constants out of loop", logging.DEBUG):
                 # This has to be done twice because sometimes constants are not moved completely to the surrounding block but

hog_tests/operator_generation/test_opgen_smoke.py

+# HyTeG Operator Generator
+# Copyright (C) 2024  HyTeG Team
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+from sympy.core.cache import clear_cache
+
+from hog.blending import IdentityMap
+from hog.operator_generation.loop_strategies import CUBES
+from hog.operator_generation.optimizer import Opts
+from hog.element_geometry import TriangleElement
+from hog.function_space import LagrangianFunctionSpace
+from hog.operator_generation.operators import (
+    HyTeGElementwiseOperator,
+    MacroIntegrationDomain,
+)
+from hog.symbolizer import Symbolizer
+from hog.quadrature import Quadrature, select_quadrule
+from hog.forms import div_k_grad, mass
+from hog.operator_generation.kernel_types import ApplyWrapper
+from hog.operator_generation.types import hyteg_type
+from hog.blending import AnnulusMap
+
+
+def test_opgen_smoke():
+    """
+    Just a simple smoke test to check that an operator can be generated.
+
+    If something is really broken, this will make the CI fail early.
+
+    We are generating a matvec method here for
+
+        ∫ k ∇u · ∇v dx + ∫ uv dx
+
+    with either integral being evaluated in their own kernel.
+
+    That may not be reasonable but tests some features.
+    """
+    clear_cache()
+
+    symbolizer = Symbolizer()
+    volume_geometry = TriangleElement()
+
+    name = f"P2DivKGradBlendingPlusMass"
+
+    trial = LagrangianFunctionSpace(2, symbolizer)
+    test = LagrangianFunctionSpace(2, symbolizer)
+    coeff = LagrangianFunctionSpace(2, symbolizer)
+    quad = Quadrature(select_quadrule(2, volume_geometry), volume_geometry)
+
+    divkgrad = div_k_grad(trial, test, volume_geometry, symbolizer, AnnulusMap(), coeff)
+    m = mass(trial, test, volume_geometry, symbolizer, AnnulusMap())
+
+    type_descriptor = hyteg_type()
+
+    kernel_types = [
+        ApplyWrapper(
+            test,
+            trial,
+            type_descriptor=type_descriptor,
+            dims=[2],
+        )
+    ]
+
+    opts = {Opts.MOVECONSTANTS, Opts.VECTORIZE, Opts.TABULATE, Opts.QUADLOOPS}
+
+    operator = HyTeGElementwiseOperator(
+        name,
+        symbolizer=symbolizer,
+        kernel_wrapper_types=kernel_types,
+        opts=opts,
+        type_descriptor=type_descriptor,
+    )
+
+    operator.add_integral(
+        name="div_k_grad",
+        dim=volume_geometry.dimensions,
+        geometry=volume_geometry,
+        integration_domain=MacroIntegrationDomain.VOLUME,
+        quad=quad,
+        blending=AnnulusMap(),
+        form=divkgrad,
+        loop_strategy=CUBES(),
+    )
+
+    operator.add_integral(
+        name="mass",
+        dim=volume_geometry.dimensions,
+        geometry=volume_geometry,
+        integration_domain=MacroIntegrationDomain.VOLUME,
+        quad=quad,
+        blending=AnnulusMap(),
+        form=m,
+        loop_strategy=CUBES(),
+    )
+
+    operator.generate_class_code(
+        ".",
+    )
+
+
+if __name__ == "__main__":
+    test_opgen_smoke()