feat(kernels): add AssembleNonInvertedDiagonal

refactor: rename stiffness operators

generate_all_operators.py

@@ -53,7 +53,7 @@ from hog.forms import (
     supg_advection,
     grad_rho_by_rho_dot_u,
     nonlinear_diffusion_stiffness,
-    nonlinear_diffusion_newton_galerkin_approximated_mass_lumped,
+    nonlinear_diffusion_newton_galerkin_approximated_mass_lumped
 )
 from hog.forms_boundary import (
     mass_boundary,
@@ -74,7 +74,7 @@ from hog.operator_generation.kernel_types import (
     ApplyWrapper,
     AssembleWrapper,
     AssembleDiagonalWrapper,
-    KernelWrapperType, LumpWrapper,
+    KernelWrapperType, LumpWrapper, AssembleNonInverseDiagonalWrapper
 )
 from hog.operator_generation.loop_strategies import (
     LoopStrategy,
@@ -574,13 +574,21 @@ class OperatorInfo:
                     )
                 )
 
-            self.kernel_types.append(
-                LumpWrapper(
-                    self.trial_space,
-                    type_descriptor=self.type_descriptor,
-                    dims=dims,
+                self.kernel_types.append(
+                    AssembleNonInverseDiagonalWrapper(
+                        self.trial_space,
+                        type_descriptor=self.type_descriptor,
+                        dims=dims,
+                    )
+                )
+
+                self.kernel_types.append(
+                    LumpWrapper(
+                        self.trial_space,
+                        type_descriptor=self.type_descriptor,
+                        dims=dims,
+                    )
                 )
-            )
 
 
 def all_operators(
@@ -635,12 +643,12 @@ def all_operators(
                             coefficient_function_space=P1),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
 
-    ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinStiffness", TrialSpace(P1),
+    ops.append(OperatorInfo("P1", "NonlinearDiffusionStiffness", TrialSpace(P1),
                             TestSpace(P1), form=partial(nonlinear_diffusion_stiffness,
                             coefficient_function_space=P1),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
 
-    ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinStiffnessDifferentiated", TrialSpace(P1),
+    ops.append(OperatorInfo("P1", "NonlinearDiffusionStiffnessDifferentiated", TrialSpace(P1),
                             TestSpace(P1), form=partial(nonlinear_diffusion_stiffness,
                             coefficient_function_space=P1, use_a_derivative=True),
                             type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))

hog/operator_generation/kernel_types.py

@@ -731,6 +731,130 @@ class AssembleDiagonalWrapper(KernelWrapperType):
             ),
         ]
 
+class AssembleNonInverseDiagonalWrapper(KernelWrapperType):
+    def __init__(
+        self,
+        fe_space: FunctionSpace,
+        type_descriptor: HOGType,
+        dst_field: str = "diag_",
+        dims: List[int] = [2, 3],
+    ):
+        self.name = "computeDiagonalOperatorValues"
+        self.dst: FunctionSpaceImpl = create_impl(
+            fe_space, dst_field, type_descriptor, is_pointer=True
+        )
+        self.src_fields = []
+        self.dst_fields = [self.dst]
+        self.dims = dims
+
+        def macro_loop(dim: int) -> str:
+            Macro = {2: "Face", 3: "Cell"}[dim]
+            macro = {2: "face", 3: "cell"}[dim]
+
+            if dim in dims:
+                return (
+                    f"for ( auto& it : storage_->get{Macro}s() )\n"
+                    f"{{\n"
+                    f"    {Macro}& {macro} = *it.second;\n"
+                    f"\n"
+                    f"    // get hold of the actual numerical data\n"
+                    f"{indent(self.dst.pointer_retrieval(dim), INDENT)}\n"
+                    f"    $pointer_retrieval_{dim}D\n"
+                    f"\n"
+                    f"    $scalar_parameter_setup_{dim}D\n"
+                    f"\n"
+                    f'    this->timingTree_->start( "kernel" );\n'
+                    f"    $kernel_function_call_{dim}D\n"
+                    f'    this->timingTree_->stop( "kernel" );\n'
+                    f"}}\n"
+                    f"\n"
+                    f"// Push result to lower-dimensional primitives\n"
+                    f"//\n"
+                    f'this->timingTree_->start( "post-communication" );\n'
+                    f"// Note: We could avoid communication here by implementing the apply() also for the respective\n"
+                    f"//       lower dimensional primitives!\n"
+                    f"{self.dst.post_communication(dim, 'level', transform_basis=False)}\n"
+                    f'this->timingTree_->stop( "post-communication" );'
+                )
+            else:
+                return 'WALBERLA_ABORT( "Not implemented." );'
+
+        fnc_name = "diagonal entries"
+
+        self._template = Template(
+            f'this->startTiming( "{self.name}" );\n'
+            f"\n"
+            f"if ( {self.dst.name} == nullptr )\n"
+            f"{{\n"
+            f"    {self.dst.name} =\n"
+            f'        std::make_shared< {self.dst.func_type_string()} >( "{fnc_name}", storage_, minLevel_, maxLevel_ );\n'
+            f"}}\n"
+            f"\n"
+            f"for ( uint_t level = minLevel_; level <= maxLevel_; level++ )\n"
+            f"{{\n"
+            f"    {self.dst.name}->setToZero( level );\n"
+            f"\n"
+            f"    if ( storage_->hasGlobalCells() )\n"
+            f"    {{\n"
+            f'        this->timingTree_->start( "pre-communication" );\n'
+            f"        $comm_fe_functions_3D\n"
+            f'        this->timingTree_->stop( "pre-communication" );\n'
+            f"\n"
+            f"{indent(macro_loop(3), 2 * INDENT)}\n"
+            f"    }}\n"
+            f"    else\n"
+            f"    {{\n"
+            f'        this->timingTree_->start( "pre-communication" );\n'
+            f"        $comm_fe_functions_2D\n"
+            f'        this->timingTree_->stop( "pre-communication" );\n'
+            f"\n"
+            f"{indent(macro_loop(2), 2 * INDENT)}\n"
+            f"    }}\n"
+            f"\n"
+            f"}}\n"
+            f"\n"
+            f'this->stopTiming( "{self.name}" );'
+        )
+
+    @property
+    def kernel_type(self) -> KernelType:
+        return AssembleDiagonal()
+
+    def includes(self) -> Set[str]:
+        return {"hyteg/solvers/Smoothables.hpp"} | self.dst.includes()
+
+    def base_classes(self) -> List[str]:
+        return [
+            f"public OperatorWithDiagonal< {self.dst.func_type_string()} >",
+        ]
+
+    def wrapper_methods(self) -> List[CppMethod]:
+        return [
+            CppMethod(
+                name=self.name,
+                arguments=[],
+                return_type="void",
+                content=self._template.template,
+            ),
+            CppMethod(
+                name="getDiagonalValues",
+                arguments=[],
+                return_type=f"std::shared_ptr< {self.dst.func_type_string()} >",
+                is_const=True,
+                content=f"return {self.dst.name};",
+            ),
+        ]
+
+    def member_variables(self) -> List[CppMemberVariable]:
+        return [
+            CppMemberVariable(
+                variable=CppVariable(
+                    name=self.dst.name,
+                    type=f"std::shared_ptr< {self.dst.func_type_string()} >",
+                )
+            ),
+        ]
+
 
 class AssembleWrapper(KernelWrapperType):
     def __init__(