diff --git a/hog/operator_generation/__init__.py b/hog/operator_generation/__init__.py index 7de6b0e232ff9b89bccc8a525c7841951ce1c2ca..c4cab3c872f55c3312a8a9216b435f2daaf541dd 100644 --- a/hog/operator_generation/__init__.py +++ b/hog/operator_generation/__init__.py @@ -15,7 +15,7 @@ # along with this program. If not, see <https://www.gnu.org/licenses/>. __all__ = [ - "function_space_impls", + "function_space_implementations", "indexing", "kernel_types", "loop_strategies", diff --git a/hog/operator_generation/function_space_implementations/__init__.py b/hog/operator_generation/function_space_implementations/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..b133d9c3e2a9fb91d0a0adb2c637c1069dc457c8 --- /dev/null +++ b/hog/operator_generation/function_space_implementations/__init__.py @@ -0,0 +1,24 @@ +# 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/>. + +__all__ = [ + "function_space_impl_base", + "function_space_impl_factory", + "p0_space_impl", + "p1_space_impl", + "p2_space_impl", + "n1e1_space_impl", +] diff --git a/hog/operator_generation/function_space_implementations/function_space_impl_base.py b/hog/operator_generation/function_space_implementations/function_space_impl_base.py new file mode 100644 index 0000000000000000000000000000000000000000..d410fbc9de96af4b96d3e33ee855d78e8299068a --- /dev/null +++ b/hog/operator_generation/function_space_implementations/function_space_impl_base.py @@ -0,0 +1,193 @@ +# 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 abc import ABC, abstractmethod + +import itertools +import numpy as np +import sympy as sp +from typing import List, Tuple, Set, Union, Dict + +import pystencils as ps +from pystencils import Field, FieldType +from pystencils.backends.cbackend import CustomCodeNode + +from hog.element_geometry import ElementGeometry +from hog.exception import HOGException +from hog.function_space import ( + FunctionSpace, + LagrangianFunctionSpace, + TensorialVectorFunctionSpace, + N1E1Space, +) +from hog.operator_generation.indexing import ( + CellType, + FaceType, + VolumeDoFMemoryLayout, + micro_element_to_vertex_indices, + micro_vertex_to_edge_indices, + micro_element_to_volume_indices, + IndexingInfo, +) +from hog.operator_generation.types import HOGType + + +class FunctionSpaceImpl(ABC): + """A FunctionSpaceImpl is the counterpart of a Function in HyTeG. + + An instance of this class represents an instance of one of HyTeG's Function + classes. It is associated with a mathematical function space, an identifier + (variable name) and can optionally be a pointer. This class is intended + to abstract printing code for e.g. communication in an FE space and C++ + implementation agnostic way. + + It is impossible to create an instance of this abstract base class directly. + Preferrably, use `function_space_impl_factory.create_impl` which selects + the correct derived class for the `FunctionSpace`. + """ + + def __init__( + self, + fe_space: FunctionSpace, + name: str, + type_descriptor: HOGType, + is_pointer: bool, + ): + """Records the passed parameters in member variables. + + It is impossible to create an instance of this abstract base class + directly. This __init__ method is to be called by the derived classes. + Preferrably, use `function_space_impl_factory.create_impl` which + selects the correct derived class for the `FunctionSpace`. + """ + self.fe_space = fe_space + self.name = name + self.type_descriptor = type_descriptor + self.is_pointer = is_pointer + + def _create_field(self, name: str) -> Field: + """Creates a pystencils field with a given name and stride 1.""" + f = Field.create_generic( + name, + 1, + dtype=self.type_descriptor.pystencils_type, + field_type=FieldType.CUSTOM, + ) + f.strides = tuple([1 for _ in f.strides]) + return f + + @abstractmethod + def pre_communication(self, dim: int) -> str: + """C++ code for the function space communication prior to the macro-primitive loop.""" + ... + + @abstractmethod + def zero_halos(self, dim: int) -> str: + """C++ code to zero the halos on a macro in HyTeG.""" + ... + + @abstractmethod + def post_communication( + self, dim: int, params: str, transform_basis: bool = True + ) -> str: + """C++ code for the function space communication after the macro-primitive loop.""" + ... + + @abstractmethod + def pointer_retrieval(self, dim: int) -> str: + """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" + ... + + @abstractmethod + def invert_elementwise(self, dim: int) -> str: + """C++ code for inverting each DoF of the linalg vector.""" + ... + + @abstractmethod + def local_dofs( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + indexing_info: IndexingInfo, + element_vertex_ordering: List[int], + ) -> List[Field.Access]: + """ + Returns a list of local dof values on the current element. + + The element_vertex_ordering is a list that specifies the ordering of the reference vertices. + The ordering in which the DoFs are returned depends on this list. The "default" ordering is + [0, 1, ..., num_vertices - 1]. + """ + ... + + @abstractmethod + def func_type_string(self) -> str: + """Returns the C++ function class type as a string.""" + ... + + @abstractmethod + def includes(self) -> Set[str]: + """Returns the import location for the function space in HyTeG.""" + ... + + @abstractmethod + def dof_transformation( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + element_vertex_ordering: List[int], + ) -> Tuple[CustomCodeNode, sp.MatrixBase]: + """Returns HyTeG code that computes the basis/DoF transformation and a symbolic expression of the result. + + In general FEM spaces, neighboring elements must agree on the + orientation of shared mesh entities. For example, in N1E1 the + orientations of edges define the sign of the DoFs, while in P3 + reorienting an edge swaps the two DoFs located on this edge. Basis + transformations between these spaces allow us to view DoFs from + neighboring elements in our local orientation. + Recommended reading: + Scroggs et al., "Construction of Arbitrary Order Finite Element Degree- + of-Freedom Maps on Polygonal and Polyhedral Cell Meshes," 2022, doi: + https://doi.org/10.1145/3524456. + + This function returns a matrix that, when applied to a local DoF vector, + transforms the DoFs from the owning primitives to the basis of the + current macro primitive. + For example: + - P1/P2: The identity. + - P3: A permutation matrix. + - N1E1: A diagonal matrix with |aᵢᵢ| = 1. + + If the micro element is in the interior of the macro-cell, the + transformation is the identity. In matrix-free computations the + communication is responsible for the basis transformations. Only when + assembling operators into matrices, these transformations must be + "baked into" the matrix because vectors are assembled locally and our + communication routine is not performed during the operator application. + + The element_vertex_ordering is a list that specifies the ordering of the reference vertices. + The ordering in which the DoFs are returned depends on this list. The "default" ordering is + [0, 1, ..., num_vertices - 1]. + """ + ... + + def _deref(self) -> str: + if self.is_pointer: + return f"(*{self.name})" + else: + return self.name diff --git a/hog/operator_generation/function_space_implementations/function_space_impl_factory.py b/hog/operator_generation/function_space_implementations/function_space_impl_factory.py new file mode 100644 index 0000000000000000000000000000000000000000..c3762ee3c596c3279b34c93a2366528934cb81d5 --- /dev/null +++ b/hog/operator_generation/function_space_implementations/function_space_impl_factory.py @@ -0,0 +1,98 @@ +# 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 hog.exception import HOGException + +from hog.function_space import ( + FunctionSpace, + LagrangianFunctionSpace, + TensorialVectorFunctionSpace, + N1E1Space, +) + +from hog.operator_generation.function_space_implementations.function_space_impl_base import ( + FunctionSpaceImpl, +) + +from hog.operator_generation.function_space_implementations.p0_space_impl import ( + P0FunctionSpaceImpl, +) +from hog.operator_generation.function_space_implementations.p1_space_impl import ( + P1FunctionSpaceImpl, + P1VectorFunctionSpaceImpl, +) +from hog.operator_generation.function_space_implementations.p2_space_impl import ( + P2FunctionSpaceImpl, + P2VectorFunctionSpaceImpl, +) +from hog.operator_generation.function_space_implementations.n1e1_space_impl import ( + N1E1FunctionSpaceImpl, +) + +from hog.operator_generation.types import HOGType + + +def create_impl( + func_space: FunctionSpace, + name: str, + type_descriptor: HOGType, + is_pointer: bool = False, +) -> FunctionSpaceImpl: + """Takes a mathematical function space and produces the corresponding function space implementation. + + :param func_space: The mathematical function space. + :param name: The C++ variable name/identifier of this function. + :param type_descriptor: The value type of this function. + :param is_pointer: Whether the C++ variable is of a pointer type. Used + to print member accesses. + """ + impl_class: type + + if isinstance(func_space, LagrangianFunctionSpace): + if func_space.degree == 2: + impl_class = P2FunctionSpaceImpl + elif func_space.degree == 1: + impl_class = P1FunctionSpaceImpl + elif func_space.degree == 0: + impl_class = P0FunctionSpaceImpl + else: + raise HOGException("Lagrangian function space must be of order 1 or 2.") + elif isinstance(func_space, TensorialVectorFunctionSpace): + if isinstance(func_space.component_function_space, LagrangianFunctionSpace): + if func_space.component_function_space.degree == 1: + if func_space.single_component is None: + impl_class = P1VectorFunctionSpaceImpl + else: + impl_class = P1FunctionSpaceImpl + elif func_space.component_function_space.degree == 2: + if func_space.single_component is None: + impl_class = P2VectorFunctionSpaceImpl + else: + impl_class = P2FunctionSpaceImpl + else: + raise HOGException( + "TensorialVectorFunctionSpaces not supported for the chosen components." + ) + else: + raise HOGException( + "TensorialVectorFunctionSpaces are only supported with Lagrangian component spaces." + ) + elif isinstance(func_space, N1E1Space): + impl_class = N1E1FunctionSpaceImpl + else: + raise HOGException("Unexpected function space") + + return impl_class(func_space, name, type_descriptor, is_pointer) diff --git a/hog/operator_generation/function_space_implementations/n1e1_space_impl.py b/hog/operator_generation/function_space_implementations/n1e1_space_impl.py new file mode 100644 index 0000000000000000000000000000000000000000..16909cd489a0d04d10837b66c2f63b5faba0788b --- /dev/null +++ b/hog/operator_generation/function_space_implementations/n1e1_space_impl.py @@ -0,0 +1,169 @@ +# 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/>. + +import itertools +import numpy as np +import sympy as sp +from typing import List, Tuple, Set, Union, Dict + +import pystencils as ps +from pystencils import Field, FieldType +from pystencils.backends.cbackend import CustomCodeNode + +from hog.element_geometry import ElementGeometry +from hog.exception import HOGException +from hog.function_space import ( + FunctionSpace, + LagrangianFunctionSpace, + TensorialVectorFunctionSpace, + N1E1Space, +) +from hog.operator_generation.indexing import ( + CellType, + FaceType, + VolumeDoFMemoryLayout, + micro_element_to_vertex_indices, + micro_vertex_to_edge_indices, + micro_element_to_volume_indices, + IndexingInfo, +) +from hog.operator_generation.types import HOGType +from hog.operator_generation.function_space_implementations.function_space_impl_base import ( + FunctionSpaceImpl, +) + + +class N1E1FunctionSpaceImpl(FunctionSpaceImpl): + def __init__( + self, + fe_space: FunctionSpace, + name: str, + type_descriptor: HOGType, + is_pointer: bool = False, + ): + super().__init__(fe_space, name, type_descriptor, is_pointer) + self.field = self._create_field(name) + + def pre_communication(self, dim: int) -> str: + if dim == 2: + return "" + else: + return ( + f"{self._deref()}.communicate< Face, Cell >( level );\n" + f"{self._deref()}.communicate< Edge, Cell >( level );" + ) + + def zero_halos(self, dim: int) -> str: + if dim == 2: + return "" + else: + return f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}.getDoFs()->getCellDataID() );" + + def post_communication( + self, dim: int, params: str, transform_basis: bool = True + ) -> str: + if dim == 2: + return "" + else: + dofs = "" + if not transform_basis: + dofs = "getDoFs()->" + return ( + f"{self._deref()}.{dofs}communicateAdditively< Cell, Face >( {params} );\n" + f"{self._deref()}.{dofs}communicateAdditively< Cell, Edge >( {params} );" + ) + + def pointer_retrieval(self, dim: int) -> str: + """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" + Macro = {2: "Face", 3: "Cell"}[dim] + macro = {2: "face", 3: "cell"}[dim] + + return f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {macro}.getData( {self._deref()}.getDoFs()->get{Macro}DataID() )->getPointer( level );" + + def invert_elementwise(self, dim: int) -> str: + return f"{self._deref()}.getDoFs()->invertElementwise( level );" + + def local_dofs( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + indexing_info: IndexingInfo, + element_vertex_ordering: List[int], + ) -> List[Field.Access]: + vertex_dof_indices = micro_element_to_vertex_indices( + geometry, element_type, element_index + ) + + vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] + + edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices) + edge_array_indices = [ + dof_idx.array_index(geometry, indexing_info) for dof_idx in edge_dof_indices + ] + + return [self.field.absolute_access((idx,), (0,)) for idx in edge_array_indices] + + def func_type_string(self) -> str: + return f"n1e1::N1E1VectorFunction< {self.type_descriptor.pystencils_type} >" + + def includes(self) -> Set[str]: + return { + "hyteg/n1e1functionspace/N1E1VectorFunction.hpp", + "hyteg/n1e1functionspace/N1E1MacroCell.hpp", + } + + def dof_transformation( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + element_vertex_ordering: List[int], + ) -> Tuple[CustomCodeNode, sp.MatrixBase]: + if element_vertex_ordering != [0, 1, 2, 3]: + raise HOGException( + "Element vertex re-ordering not supported for Nédélec elements (yet)." + ) + + Macro = {2: "Face", 3: "Cell"}[geometry.dimensions] + macro = {2: "face", 3: "cell"}[geometry.dimensions] + name = "basisTransformation" + n_dofs = self.fe_space.num_dofs(geometry) + + # NOTE: Types are added manually because pystencils won't add types to accessed/ + # defined symbols of custom code nodes. As a result the symbols + # in the matrix will be typed but not their definition, leading to + # undefined symbols. + # NOTE: The type is `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). + symbols = [ + ps.TypedSymbol(f"{name}.diagonal()({i})", np.float64) for i in range(n_dofs) + ] + return ( + CustomCodeNode( + f"const Eigen::DiagonalMatrix< real_t, {n_dofs} > {name} = " + f"n1e1::macro{macro}::basisTransformation( level, {macro}, {{{element_index[0]}, {element_index[1]}, {element_index[2]}}}, {macro}dof::{Macro}Type::{element_type.name} );", + [ + ps.TypedSymbol("level", "const uint_t"), + ps.TypedSymbol(f"{macro}", f"const {Macro}&"), + ], + symbols, + ), + sp.diag(*symbols), + ) diff --git a/hog/operator_generation/function_space_implementations/p0_space_impl.py b/hog/operator_generation/function_space_implementations/p0_space_impl.py new file mode 100644 index 0000000000000000000000000000000000000000..48987bb8bddb0fdb2ea2b4589c98750b7ebb3b14 --- /dev/null +++ b/hog/operator_generation/function_space_implementations/p0_space_impl.py @@ -0,0 +1,118 @@ +# 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/>. + +import itertools +import numpy as np +import sympy as sp +from typing import List, Tuple, Set, Union, Dict + +import pystencils as ps +from pystencils import Field, FieldType +from pystencils.backends.cbackend import CustomCodeNode + +from hog.element_geometry import ElementGeometry +from hog.exception import HOGException +from hog.function_space import ( + FunctionSpace, + LagrangianFunctionSpace, + TensorialVectorFunctionSpace, + N1E1Space, +) +from hog.operator_generation.indexing import ( + CellType, + FaceType, + VolumeDoFMemoryLayout, + micro_element_to_vertex_indices, + micro_vertex_to_edge_indices, + micro_element_to_volume_indices, + IndexingInfo, +) +from hog.operator_generation.types import HOGType +from hog.operator_generation.function_space_implementations.function_space_impl_base import ( + FunctionSpaceImpl, +) + + +class P0FunctionSpaceImpl(FunctionSpaceImpl): + def __init__( + self, + fe_space: FunctionSpace, + name: str, + type_descriptor: HOGType, + is_pointer: bool = False, + ): + super().__init__(fe_space, name, type_descriptor, is_pointer) + self.field = self._create_field(name) + + def pre_communication(self, dim: int) -> str: + return "" + + def zero_halos(self, dim: int) -> str: + return "" + + def post_communication( + self, dim: int, params: str, transform_basis: bool = True + ) -> str: + return "" + + def pointer_retrieval(self, dim: int) -> str: + """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" + Macro = {2: "Face", 3: "Cell"}[dim] + macro = {2: "face", 3: "cell"}[dim] + + return f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {self._deref()}.getDGFunction()->volumeDoFFunction()->dofMemory( it.first, level );" + + def invert_elementwise(self, dim: int) -> str: + return f"{self._deref()}.invertElementwise( level );" + + def local_dofs( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + indexing_info: IndexingInfo, + element_vertex_ordering: List[int], + ) -> List[Field.Access]: + volume_dof_indices = micro_element_to_volume_indices( + element_type, element_index, 1, VolumeDoFMemoryLayout.SoA + ) + + vertex_array_indices = [ + dof_idx.array_index(geometry, indexing_info) + for dof_idx in volume_dof_indices + ] + + return [ + self.field.absolute_access((idx,), (0,)) for idx in vertex_array_indices + ] + + def func_type_string(self) -> str: + return f"P0Function< {self.type_descriptor.pystencils_type} >" + + def includes(self) -> Set[str]: + return {f"hyteg/p0functionspace/P0Function.hpp"} + + def dof_transformation( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + element_vertex_ordering: List[int], + ) -> Tuple[CustomCodeNode, sp.MatrixBase]: + return ( + CustomCodeNode("", [], []), + sp.Identity(self.fe_space.num_dofs(geometry)), + ) diff --git a/hog/operator_generation/function_space_implementations/p1_space_impl.py b/hog/operator_generation/function_space_implementations/p1_space_impl.py new file mode 100644 index 0000000000000000000000000000000000000000..0f993dce5f4820d6ab40de9f9bce4ea9d0d56ea9 --- /dev/null +++ b/hog/operator_generation/function_space_implementations/p1_space_impl.py @@ -0,0 +1,333 @@ +# 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/>. + +import itertools +import numpy as np +import sympy as sp +from typing import List, Tuple, Set, Union, Dict + +import pystencils as ps +from pystencils import Field, FieldType +from pystencils.backends.cbackend import CustomCodeNode + +from hog.element_geometry import ElementGeometry +from hog.exception import HOGException +from hog.function_space import ( + FunctionSpace, + LagrangianFunctionSpace, + TensorialVectorFunctionSpace, + N1E1Space, +) +from hog.operator_generation.indexing import ( + CellType, + FaceType, + VolumeDoFMemoryLayout, + micro_element_to_vertex_indices, + micro_vertex_to_edge_indices, + micro_element_to_volume_indices, + IndexingInfo, +) +from hog.operator_generation.types import HOGType +from hog.operator_generation.function_space_implementations.function_space_impl_base import ( + FunctionSpaceImpl, +) + + +class P1FunctionSpaceImpl(FunctionSpaceImpl): + def __init__( + self, + fe_space: FunctionSpace, + name: str, + type_descriptor: HOGType, + is_pointer: bool = False, + ): + super().__init__(fe_space, name, type_descriptor, is_pointer) + self.field = self._create_field(name) + + def pre_communication(self, dim: int) -> str: + if dim == 2: + return f"communication::syncFunctionBetweenPrimitives( {self.name}, level, communication::syncDirection_t::LOW2HIGH );" + else: + return ( + f"{self._deref()}.communicate< Face, Cell >( level );\n" + f"{self._deref()}.communicate< Edge, Cell >( level );\n" + f"{self._deref()}.communicate< Vertex, Cell >( level );" + ) + + def zero_halos(self, dim: int) -> str: + if dim == 2: + return ( + f"for ( const auto& idx : vertexdof::macroface::Iterator( level ) ) {{\n" + f" if ( vertexdof::macroface::isVertexOnBoundary( level, idx ) ) {{\n" + f" auto arrayIdx = vertexdof::macroface::index( level, idx.x(), idx.y() );\n" + f" _data_{self.name}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" }}\n" + f"}}" + ) + else: + return ( + f"for ( const auto& idx : vertexdof::macrocell::Iterator( level ) ) {{\n" + f" if ( !vertexdof::macrocell::isOnCellFace( idx, level ).empty() ) {{\n" + f" auto arrayIdx = vertexdof::macrocell::index( level, idx.x(), idx.y(), idx.z() );\n" + f" _data_{self.name}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" }}\n" + f"}}" + ) + + def post_communication( + self, dim: int, params: str, transform_basis: bool = True + ) -> str: + if dim == 2: + return ( + f"{self._deref()}.communicateAdditively < Face, Edge > ( {params} );\n" + f"{self._deref()}.communicateAdditively < Face, Vertex > ( {params} );" + ) + else: + return ( + f"{self._deref()}.communicateAdditively< Cell, Face >( {params} );\n" + f"{self._deref()}.communicateAdditively< Cell, Edge >( {params} );\n" + f"{self._deref()}.communicateAdditively< Cell, Vertex >( {params} );" + ) + + def pointer_retrieval(self, dim: int) -> str: + """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" + Macro = {2: "Face", 3: "Cell"}[dim] + macro = {2: "face", 3: "cell"}[dim] + + return f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {macro}.getData( {self._deref()}.get{Macro}DataID() )->getPointer( level );" + + def invert_elementwise(self, dim: int) -> str: + return f"{self._deref()}.invertElementwise( level );" + + def local_dofs( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + indexing_info: IndexingInfo, + element_vertex_ordering: List[int], + ) -> List[Field.Access]: + vertex_dof_indices = micro_element_to_vertex_indices( + geometry, element_type, element_index + ) + + vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] + + vertex_array_indices = [ + dof_idx.array_index(geometry, indexing_info) + for dof_idx in vertex_dof_indices + ] + + return [ + self.field.absolute_access((idx,), (0,)) for idx in vertex_array_indices + ] + + def func_type_string(self) -> str: + return f"P1Function< {self.type_descriptor.pystencils_type} >" + + def includes(self) -> Set[str]: + return {f"hyteg/p1functionspace/P1Function.hpp"} + + def dof_transformation( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + element_vertex_ordering: List[int], + ) -> Tuple[CustomCodeNode, sp.MatrixBase]: + return ( + CustomCodeNode("", [], []), + sp.Identity(self.fe_space.num_dofs(geometry)), + ) + + +class P1VectorFunctionSpaceImpl(FunctionSpaceImpl): + def __init__( + self, + fe_space: FunctionSpace, + name: str, + type_descriptor: HOGType, + is_pointer: bool = False, + ): + super().__init__(fe_space, name, type_descriptor, is_pointer) + self.fields: Dict[int, Field] = {} + + def _field_name(self, component: int) -> str: + return self.name + f"_{component}" + + def _raw_pointer_name(self, component: int) -> str: + return f"_data_" + self._field_name(component) + + def pre_communication(self, dim: int) -> str: + if dim == 2: + return f"communication::syncVectorFunctionBetweenPrimitives( {self.name}, level, communication::syncDirection_t::LOW2HIGH );" + else: + ret_str = "" + for i in range(dim): + ret_str += ( + f"{self._deref()}[{i}].communicate< Face, Cell >( level );\n" + f"{self._deref()}[{i}].communicate< Edge, Cell >( level );\n" + f"{self._deref()}[{i}].communicate< Vertex, Cell >( level );" + ) + return ret_str + + def zero_halos(self, dim: int) -> str: + if dim == 2: + return ( + f"for ( const auto& idx : vertexdof::macroface::Iterator( level ) ) {{\n" + f" if ( vertexdof::macroface::isVertexOnBoundary( level, idx ) ) {{\n" + f" auto arrayIdx = vertexdof::macroface::index( level, idx.x(), idx.y() );\n" + f" {self._raw_pointer_name(0)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" {self._raw_pointer_name(1)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" }}\n" + f"}}" + ) + else: + return ( + f"for ( const auto& idx : vertexdof::macrocell::Iterator( level ) ) {{\n" + f" if ( !vertexdof::macrocell::isOnCellFace( idx, level ).empty() ) {{\n" + f" auto arrayIdx = vertexdof::macrocell::index( level, idx.x(), idx.y(), idx.z() );\n" + f" {self._raw_pointer_name(0)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" {self._raw_pointer_name(1)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" {self._raw_pointer_name(2)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" }}\n" + f"}}" + ) + + def post_communication( + self, dim: int, params: str, transform_basis: bool = True + ) -> str: + if dim == 2: + ret_str = "" + for i in range(dim): + ret_str += ( + f"{self._deref()}[{i}].communicateAdditively < Face, Edge > ( {params} );\n" + f"{self._deref()}[{i}].communicateAdditively < Face, Vertex > ( {params} );\n" + ) + return ret_str + else: + ret_str = "" + for i in range(dim): + ret_str += ( + f"{self._deref()}[{i}].communicateAdditively< Cell, Face >( {params} );\n" + f"{self._deref()}[{i}].communicateAdditively< Cell, Edge >( {params} );\n" + f"{self._deref()}[{i}].communicateAdditively< Cell, Vertex >( {params} );" + ) + return ret_str + + def pointer_retrieval(self, dim: int) -> str: + """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" + Macro = {2: "Face", 3: "Cell"}[dim] + macro = {2: "face", 3: "cell"}[dim] + + ret_str = "" + for i in range(dim): + ret_str += f"{self.type_descriptor.pystencils_type}* {self._raw_pointer_name(i)} = {macro}.getData( {self._deref()}[{i}].get{Macro}DataID() )->getPointer( level );\n" + return ret_str + + def invert_elementwise(self, dim: int) -> str: + ret_str = "" + for i in range(dim): + ret_str += f"{self._deref()}[{i}].invertElementwise( level );\n" + return ret_str + + def local_dofs( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + indexing_info: IndexingInfo, + element_vertex_ordering: List[int], + ) -> List[Field.Access]: + """ + Returns the element-local DoFs (field accesses) in a list (i.e., linearized). + + The order here has to match that of the function space implementation. + TODO: This is a little concerning since that order has to match in two very different parts of this + implementation. Here, and in the TensorialVectorFunctionSpace. Maybe we can define this order in a single + location. + We return them in "SoA" order: first all DoFs corresponding to the function with the first component != 0, + then all DoFs corresponding to the function with the second component != 0 etc. + + For instance, in 2D we get the DoFs corresponding to the 6 shape functions in the following order: + + First component != 0 + + list[0]: ⎡-x_ref_0 - x_ref_1 + 1⎤ + ⎢ ⎥ + ⎣ 0 ⎦ + + list[1]: ⎡x_ref_0⎤ + ⎢ ⎥ + ⎣ 0 ⎦ + + list[2]: ⎡x_ref_1⎤ + ⎢ ⎥ + ⎣ 0 ⎦ + + Second component != 0 + + list[3]: ⎡ 0 ⎤ + ⎢ ⎥ + ⎣-x_ref_0 - x_ref_1 + 1⎦ + + list[4]: ⎡ 0 ⎤ + ⎢ ⎥ + ⎣x_ref_0⎦ + + list[5]: ⎡ 0 ⎤ + ⎢ ⎥ + ⎣x_ref_1⎦ + """ + + for c in range(geometry.dimensions): + if c not in self.fields: + self.fields[c] = self._create_field(self._field_name(c)) + + vertex_dof_indices = micro_element_to_vertex_indices( + geometry, element_type, element_index + ) + + vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] + + vertex_array_indices = [ + dof_idx.array_index(geometry, indexing_info) + for dof_idx in vertex_dof_indices + ] + + return [ + self.fields[c].absolute_access((idx,), (0,)) + for c in range(geometry.dimensions) + for idx in vertex_array_indices + ] + + def func_type_string(self) -> str: + return f"P1VectorFunction< {self.type_descriptor.pystencils_type} >" + + def includes(self) -> Set[str]: + return {f"hyteg/p1functionspace/P1VectorFunction.hpp"} + + def dof_transformation( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + element_vertex_ordering: List[int], + ) -> Tuple[CustomCodeNode, sp.MatrixBase]: + return ( + CustomCodeNode("", [], []), + sp.Identity(self.fe_space.num_dofs(geometry)), + ) diff --git a/hog/operator_generation/function_space_implementations/p2_space_impl.py b/hog/operator_generation/function_space_implementations/p2_space_impl.py new file mode 100644 index 0000000000000000000000000000000000000000..01eec61b8cca34a84bb54bf04f87cba9f34d5219 --- /dev/null +++ b/hog/operator_generation/function_space_implementations/p2_space_impl.py @@ -0,0 +1,372 @@ +# 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/>. + +import itertools +import numpy as np +import sympy as sp +from typing import List, Tuple, Set, Union, Dict + +import pystencils as ps +from pystencils import Field, FieldType +from pystencils.backends.cbackend import CustomCodeNode + +from hog.element_geometry import ElementGeometry +from hog.exception import HOGException +from hog.function_space import ( + FunctionSpace, + LagrangianFunctionSpace, + TensorialVectorFunctionSpace, + N1E1Space, +) +from hog.operator_generation.indexing import ( + CellType, + FaceType, + VolumeDoFMemoryLayout, + micro_element_to_vertex_indices, + micro_vertex_to_edge_indices, + micro_element_to_volume_indices, + IndexingInfo, +) +from hog.operator_generation.types import HOGType +from hog.operator_generation.function_space_implementations.function_space_impl_base import ( + FunctionSpaceImpl, +) + + +class P2FunctionSpaceImpl(FunctionSpaceImpl): + def __init__( + self, + fe_space: FunctionSpace, + name: str, + type_descriptor: HOGType, + is_pointer: bool = False, + ): + super().__init__(fe_space, name, type_descriptor, is_pointer) + self.f_vertex = self._create_field(name + "Vertex") + self.f_edge = self._create_field(name + "Edge") + + def pre_communication(self, dim: int) -> str: + if dim == 2: + return f"communication::syncFunctionBetweenPrimitives( {self.name}, level, communication::syncDirection_t::LOW2HIGH );" + else: + return ( + f"{self._deref()}.communicate< Face, Cell >( level );\n" + f"{self._deref()}.communicate< Edge, Cell >( level );\n" + f"{self._deref()}.communicate< Vertex, Cell >( level );" + ) + + def zero_halos(self, dim: int) -> str: + if dim == 2: + return ( + f"for ( const auto& idx : vertexdof::macroface::Iterator( level ) )\n" + f"{{\n" + f" if ( vertexdof::macroface::isVertexOnBoundary( level, idx ) )\n" + f" {{\n" + f" auto arrayIdx = vertexdof::macroface::index( level, idx.x(), idx.y() );\n" + f" _data_{self.name }Vertex[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" + f" }}\n" + f"}}\n" + f"for ( const auto& idx : edgedof::macroface::Iterator( level ) )\n" + f"{{\n" + f" for ( const auto& orientation : edgedof::faceLocalEdgeDoFOrientations )\n" + f" {{\n" + f" if ( !edgedof::macroface::isInnerEdgeDoF( level, idx, orientation ) )\n" + f" {{\n" + f" auto arrayIdx = edgedof::macroface::index( level, idx.x(), idx.y(), orientation );\n" + f" _data_{self.name }Edge[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" + f" }}\n" + f" }}\n" + f"}}" + ) + else: + return ( + f"for ( const auto& idx : vertexdof::macrocell::Iterator( level ) )\n" + f"{{\n" + f"if ( !vertexdof::macrocell::isOnCellFace( idx, level ).empty() )\n" + f" {{\n" + f" auto arrayIdx = vertexdof::macrocell::index( level, idx.x(), idx.y(), idx.z() );\n" + f" _data_{self.name}Vertex[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" }}\n" + f"}}\n" + f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}.getEdgeDoFFunction().getCellDataID() );" + ) + + def post_communication( + self, dim: int, params: str, transform_basis: bool = True + ) -> str: + if dim == 2: + return ( + f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Face, Edge >( {params} );\n" + f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Face, Vertex >( {params} );\n" + f"{self._deref()}.getEdgeDoFFunction().communicateAdditively< Face, Edge >( {params} );" + ) + else: + return ( + f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Cell, Face >( {params} );\n" + f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Cell, Edge >( {params} );\n" + f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Cell, Vertex >( {params} );\n" + f"{self._deref()}.getEdgeDoFFunction().communicateAdditively< Cell, Face >( {params} );\n" + f"{self._deref()}.getEdgeDoFFunction().communicateAdditively< Cell, Edge >( {params} );" + ) + + def pointer_retrieval(self, dim: int) -> str: + """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" + Macro = {2: "Face", 3: "Cell"}[dim] + macro = {2: "face", 3: "cell"}[dim] + + return ( + f"{self.type_descriptor.pystencils_type}* _data_{self.name}Vertex = {macro}.getData( {self._deref()}.getVertexDoFFunction().get{Macro}DataID() )->getPointer( level );\n" + f"{self.type_descriptor.pystencils_type}* _data_{self.name}Edge = {macro}.getData( {self._deref()}.getEdgeDoFFunction().get{Macro}DataID() )->getPointer( level );" + ) + + def invert_elementwise(self, dim: int) -> str: + return f"{self._deref()}.invertElementwise( level );" + + def local_dofs( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + indexing_info: IndexingInfo, + element_vertex_ordering: List[int], + ) -> List[Field.Access]: + vertex_dof_indices = micro_element_to_vertex_indices( + geometry, element_type, element_index + ) + + vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] + + edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices) + + vrtx_array_idcs = [ + dof_idx.array_index(geometry, indexing_info) + for dof_idx in vertex_dof_indices + ] + edge_array_idcs = [ + dof_idx.array_index(geometry, indexing_info) for dof_idx in edge_dof_indices + ] + + # NOTE: The order of DoFs must match the order of shape functions + # defined in the FunctionSpace. Holds within the individual lists + # and overall. + return [ + self.f_vertex.absolute_access((idx,), (0,)) for idx in vrtx_array_idcs + ] + [self.f_edge.absolute_access((idx,), (0,)) for idx in edge_array_idcs] + + def func_type_string(self) -> str: + return f"P2Function< {self.type_descriptor.pystencils_type} >" + + def includes(self) -> Set[str]: + return {f"hyteg/p2functionspace/P2Function.hpp"} + + def dof_transformation( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + element_vertex_ordering: List[int], + ) -> Tuple[CustomCodeNode, sp.MatrixBase]: + return ( + CustomCodeNode("", [], []), + sp.Identity(self.fe_space.num_dofs(geometry)), + ) + + +class P2VectorFunctionSpaceImpl(FunctionSpaceImpl): + def __init__( + self, + fe_space: FunctionSpace, + name: str, + type_descriptor: HOGType, + is_pointer: bool = False, + ): + super().__init__(fe_space, name, type_descriptor, is_pointer) + self.fields: Dict[Tuple[str, int], Field] = {} + + def _field_name(self, component: int, dof_type: str) -> str: + """dof_type should either be 'vertex' or 'edge'""" + return self.name + f"_{dof_type}_{component}" + + def _raw_pointer_name(self, component: int, dof_type: str) -> str: + """dof_type should either be 'vertex' or 'edge'""" + return f"_data_" + self._field_name(component, dof_type) + + def pre_communication(self, dim: int) -> str: + if dim == 2: + return f"communication::syncVectorFunctionBetweenPrimitives( {self.name}, level, communication::syncDirection_t::LOW2HIGH );" + else: + ret_str = "" + for i in range(dim): + ret_str += ( + f"{self._deref()}[{i}].communicate< Face, Cell >( level );\n" + f"{self._deref()}[{i}].communicate< Edge, Cell >( level );\n" + f"{self._deref()}[{i}].communicate< Vertex, Cell >( level );" + ) + return ret_str + + def zero_halos(self, dim: int) -> str: + if dim == 2: + return ( + f"for ( const auto& idx : vertexdof::macroface::Iterator( level ) )\n" + f"{{\n" + f" if ( vertexdof::macroface::isVertexOnBoundary( level, idx ) )\n" + f" {{\n" + f" auto arrayIdx = vertexdof::macroface::index( level, idx.x(), idx.y() );\n" + f" {self._raw_pointer_name(0, 'vertex')}[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" + f" {self._raw_pointer_name(1, 'vertex')}[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" + f" }}\n" + f"}}\n" + f"for ( const auto& idx : edgedof::macroface::Iterator( level ) )\n" + f"{{\n" + f" for ( const auto& orientation : edgedof::faceLocalEdgeDoFOrientations )\n" + f" {{\n" + f" if ( !edgedof::macroface::isInnerEdgeDoF( level, idx, orientation ) )\n" + f" {{\n" + f" auto arrayIdx = edgedof::macroface::index( level, idx.x(), idx.y(), orientation );\n" + f" {self._raw_pointer_name(0, 'edge')}[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" + f" {self._raw_pointer_name(1, 'edge')}[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" + f" }}\n" + f" }}\n" + f"}}" + ) + else: + return ( + f"for ( const auto& idx : vertexdof::macrocell::Iterator( level ) )\n" + f"{{\n" + f"if ( !vertexdof::macrocell::isOnCellFace( idx, level ).empty() )\n" + f" {{\n" + f" auto arrayIdx = vertexdof::macrocell::index( level, idx.x(), idx.y(), idx.z() );\n" + f" {self._raw_pointer_name(0, 'vertex')}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" {self._raw_pointer_name(1, 'vertex')}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" {self._raw_pointer_name(2, 'vertex')}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" + f" }}\n" + f"}}\n" + f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}[0].getEdgeDoFFunction().getCellDataID() );\n" + f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}[1].getEdgeDoFFunction().getCellDataID() );\n" + f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}[2].getEdgeDoFFunction().getCellDataID() );\n" + ) + + def post_communication( + self, dim: int, params: str, transform_basis: bool = True + ) -> str: + if dim == 2: + ret_str = "" + for i in range(dim): + ret_str += ( + f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Face, Edge >( {params} );\n" + f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Face, Vertex >( {params} );\n" + f"{self._deref()}[{i}].getEdgeDoFFunction().communicateAdditively< Face, Edge >( {params} );" + ) + return ret_str + else: + ret_str = "" + for i in range(dim): + ret_str += ( + f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Cell, Face >( {params} );\n" + f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Cell, Edge >( {params} );\n" + f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Cell, Vertex >( {params} );\n" + f"{self._deref()}[{i}].getEdgeDoFFunction().communicateAdditively< Cell, Face >( {params} );\n" + f"{self._deref()}[{i}].getEdgeDoFFunction().communicateAdditively< Cell, Edge >( {params} );" + ) + return ret_str + + def pointer_retrieval(self, dim: int) -> str: + """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" + Macro = {2: "Face", 3: "Cell"}[dim] + macro = {2: "face", 3: "cell"}[dim] + + ret_str = "" + for i in range(dim): + ret_str += f"{self.type_descriptor.pystencils_type}* {self._raw_pointer_name(i, 'vertex')} = {macro}.getData( {self._deref()}[{i}].getVertexDoFFunction().get{Macro}DataID() )->getPointer( level );\n" + ret_str += f"{self.type_descriptor.pystencils_type}* {self._raw_pointer_name(i, 'edge')} = {macro}.getData( {self._deref()}[{i}].getEdgeDoFFunction().get{Macro}DataID() )->getPointer( level );\n" + return ret_str + + def invert_elementwise(self, dim: int) -> str: + ret_str = "" + for i in range(dim): + ret_str += f"{self._deref()}[{i}].invertElementwise( level );\n" + return ret_str + + def local_dofs( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + indexing_info: IndexingInfo, + element_vertex_ordering: List[int], + ) -> List[Field.Access]: + """ + Returns the element-local DoFs (field accesses) in a list (i.e., linearized). + + See P1VectorFunctionSpaceImpl::local_dofs() for details. + """ + + for dt in ["vertex", "edge"]: + for c in range(geometry.dimensions): + if (dt, c) not in self.fields: + self.fields[(dt, c)] = self._create_field(self._field_name(c, dt)) + + vertex_dof_indices = micro_element_to_vertex_indices( + geometry, element_type, element_index + ) + + vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] + + edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices) + + vertex_array_indices = [ + dof_idx.array_index(geometry, indexing_info) + for dof_idx in vertex_dof_indices + ] + + edge_array_indices = [ + dof_idx.array_index(geometry, indexing_info) for dof_idx in edge_dof_indices + ] + + loc_dofs = [] + + for c in range(geometry.dimensions): + loc_dofs += [ + self.fields[("vertex", c)].absolute_access((idx,), (0,)) + for idx in vertex_array_indices + ] + loc_dofs += [ + self.fields[("edge", c)].absolute_access((idx,), (0,)) + for idx in edge_array_indices + ] + + return loc_dofs + + def func_type_string(self) -> str: + return f"P2VectorFunction< {self.type_descriptor.pystencils_type} >" + + def includes(self) -> Set[str]: + return { + f"hyteg/p2functionspace/P2VectorFunction.hpp", + f"hyteg/p2functionspace/P2Function.hpp", + } + + def dof_transformation( + self, + geometry: ElementGeometry, + element_index: Tuple[int, int, int], + element_type: Union[FaceType, CellType], + element_vertex_ordering: List[int], + ) -> Tuple[CustomCodeNode, sp.MatrixBase]: + return ( + CustomCodeNode("", [], []), + sp.Identity(self.fe_space.num_dofs(geometry)), + ) diff --git a/hog/operator_generation/function_space_impls.py b/hog/operator_generation/function_space_impls.py deleted file mode 100644 index 956001b75dd2c2f62b33f1695891047f03c923ec..0000000000000000000000000000000000000000 --- a/hog/operator_generation/function_space_impls.py +++ /dev/null @@ -1,1053 +0,0 @@ -# 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 abc import ABC, abstractmethod - -import itertools -import numpy as np -import sympy as sp -from typing import List, Tuple, Set, Union, Dict - -import pystencils as ps -from pystencils import Field, FieldType -from pystencils.backends.cbackend import CustomCodeNode - -from hog.element_geometry import ElementGeometry -from hog.exception import HOGException -from hog.function_space import ( - FunctionSpace, - LagrangianFunctionSpace, - TensorialVectorFunctionSpace, - N1E1Space, -) -from hog.operator_generation.indexing import ( - CellType, - FaceType, - VolumeDoFMemoryLayout, - micro_element_to_vertex_indices, - micro_vertex_to_edge_indices, - micro_element_to_volume_indices, - IndexingInfo, -) -from hog.operator_generation.types import HOGType - - -class FunctionSpaceImpl(ABC): - """A FunctionSpaceImpl is the counterpart of a Function in HyTeG. - - An instance of this class represents an instance of one of HyTeG's Function - classes. It is associated with a mathematical function space, an identifier - (variable name) and can optionally be a pointer. This class is intended - to abstract printing code for e.g. communication in an FE space and C++ - implementation agnostic way. - - It is impossible to create an instance of this abstract base class directly. - Preferrably, use the static method `create_impl` which selects the correct - derived class for the `FunctionSpace`. - """ - - def __init__( - self, - fe_space: FunctionSpace, - name: str, - type_descriptor: HOGType, - is_pointer: bool, - ): - """Records the passed parameters in member variables. - - It is impossible to create an instance of this abstract base class - directly. This __init__ method is to be called by the derived classes. - Preferrably, use the static method `create_impl` which selects the - correct derived class for the `FunctionSpace`. - """ - self.fe_space = fe_space - self.name = name - self.type_descriptor = type_descriptor - self.is_pointer = is_pointer - - @staticmethod - def create_impl( # type: ignore[no-untyped-def] # returns Self but Self is kind of new - func_space: FunctionSpace, - name: str, - type_descriptor: HOGType, - is_pointer: bool = False, - ): - """Takes a mathematical function space and produces the corresponding function space implementation. - - :param func_space: The mathematical function space. - :param name: The C++ variable name/identifier of this function. - :param type_descriptor: The value type of this function. - :param is_pointer: Whether the C++ variable is of a pointer type. Used - to print member accesses. - """ - impl_class: type - - if isinstance(func_space, LagrangianFunctionSpace): - if func_space.degree == 2: - impl_class = P2FunctionSpaceImpl - elif func_space.degree == 1: - impl_class = P1FunctionSpaceImpl - elif func_space.degree == 0: - impl_class = P0FunctionSpaceImpl - else: - raise HOGException("Lagrangian function space must be of order 1 or 2.") - elif isinstance(func_space, TensorialVectorFunctionSpace): - if isinstance(func_space.component_function_space, LagrangianFunctionSpace): - if func_space.component_function_space.degree == 1: - if func_space.single_component is None: - impl_class = P1VectorFunctionSpaceImpl - else: - impl_class = P1FunctionSpaceImpl - elif func_space.component_function_space.degree == 2: - if func_space.single_component is None: - impl_class = P2VectorFunctionSpaceImpl - else: - impl_class = P2FunctionSpaceImpl - else: - raise HOGException( - "TensorialVectorFunctionSpaces not supported for the chosen components." - ) - else: - raise HOGException( - "TensorialVectorFunctionSpaces are only supported with Lagrangian component spaces." - ) - elif isinstance(func_space, N1E1Space): - impl_class = N1E1FunctionSpaceImpl - else: - raise HOGException("Unexpected function space") - - return impl_class(func_space, name, type_descriptor, is_pointer) - - def _create_field(self, name: str) -> Field: - """Creates a pystencils field with a given name and stride 1.""" - f = Field.create_generic( - name, - 1, - dtype=self.type_descriptor.pystencils_type, - field_type=FieldType.CUSTOM, - ) - f.strides = tuple([1 for _ in f.strides]) - return f - - @abstractmethod - def pre_communication(self, dim: int) -> str: - """C++ code for the function space communication prior to the macro-primitive loop.""" - ... - - @abstractmethod - def zero_halos(self, dim: int) -> str: - """C++ code to zero the halos on a macro in HyTeG.""" - ... - - @abstractmethod - def post_communication( - self, dim: int, params: str, transform_basis: bool = True - ) -> str: - """C++ code for the function space communication after the macro-primitive loop.""" - ... - - @abstractmethod - def pointer_retrieval(self, dim: int) -> str: - """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" - ... - - @abstractmethod - def invert_elementwise(self, dim: int) -> str: - """C++ code for inverting each DoF of the linalg vector.""" - ... - - @abstractmethod - def local_dofs( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - indexing_info: IndexingInfo, - element_vertex_ordering: List[int], - ) -> List[Field.Access]: - """ - Returns a list of local dof values on the current element. - - The element_vertex_ordering is a list that specifies the ordering of the reference vertices. - The ordering in which the DoFs are returned depends on this list. The "default" ordering is - [0, 1, ..., num_vertices - 1]. - """ - ... - - @abstractmethod - def func_type_string(self) -> str: - """Returns the C++ function class type as a string.""" - ... - - @abstractmethod - def includes(self) -> Set[str]: - """Returns the import location for the function space in HyTeG.""" - ... - - @abstractmethod - def dof_transformation( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - element_vertex_ordering: List[int], - ) -> Tuple[CustomCodeNode, sp.MatrixBase]: - """Returns HyTeG code that computes the basis/DoF transformation and a symbolic expression of the result. - - In general FEM spaces, neighboring elements must agree on the - orientation of shared mesh entities. For example, in N1E1 the - orientations of edges define the sign of the DoFs, while in P3 - reorienting an edge swaps the two DoFs located on this edge. Basis - transformations between these spaces allow us to view DoFs from - neighboring elements in our local orientation. - Recommended reading: - Scroggs et al., "Construction of Arbitrary Order Finite Element Degree- - of-Freedom Maps on Polygonal and Polyhedral Cell Meshes," 2022, doi: - https://doi.org/10.1145/3524456. - - This function returns a matrix that, when applied to a local DoF vector, - transforms the DoFs from the owning primitives to the basis of the - current macro primitive. - For example: - - P1/P2: The identity. - - P3: A permutation matrix. - - N1E1: A diagonal matrix with |aᵢᵢ| = 1. - - If the micro element is in the interior of the macro-cell, the - transformation is the identity. In matrix-free computations the - communication is responsible for the basis transformations. Only when - assembling operators into matrices, these transformations must be - "baked into" the matrix because vectors are assembled locally and our - communication routine is not performed during the operator application. - - The element_vertex_ordering is a list that specifies the ordering of the reference vertices. - The ordering in which the DoFs are returned depends on this list. The "default" ordering is - [0, 1, ..., num_vertices - 1]. - """ - ... - - def _deref(self) -> str: - if self.is_pointer: - return f"(*{self.name})" - else: - return self.name - -class P0FunctionSpaceImpl(FunctionSpaceImpl): - def __init__( - self, - fe_space: FunctionSpace, - name: str, - type_descriptor: HOGType, - is_pointer: bool = False, - ): - super().__init__(fe_space, name, type_descriptor, is_pointer) - self.field = self._create_field(name) - - def pre_communication(self, dim: int) -> str: - return "" - - def zero_halos(self, dim: int) -> str: - return "" - - def post_communication( - self, dim: int, params: str, transform_basis: bool = True - ) -> str: - return "" - - def pointer_retrieval(self, dim: int) -> str: - """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" - Macro = {2: "Face", 3: "Cell"}[dim] - macro = {2: "face", 3: "cell"}[dim] - - return f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {self._deref()}.getDGFunction()->volumeDoFFunction()->dofMemory( it.first, level );" - - def invert_elementwise(self, dim: int) -> str: - return f"{self._deref()}.invertElementwise( level );" - - def local_dofs( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - indexing_info: IndexingInfo, - element_vertex_ordering: List[int], - ) -> List[Field.Access]: - volume_dof_indices = micro_element_to_volume_indices( - element_type, element_index, 1, VolumeDoFMemoryLayout.SoA - ) - - vertex_array_indices = [ - dof_idx.array_index(geometry, indexing_info) - for dof_idx in volume_dof_indices - ] - - return [ - self.field.absolute_access((idx,), (0,)) for idx in vertex_array_indices - ] - - def func_type_string(self) -> str: - return f"P0Function< {self.type_descriptor.pystencils_type} >" - - def includes(self) -> Set[str]: - return {f"hyteg/p0functionspace/P0Function.hpp"} - - def dof_transformation( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - element_vertex_ordering: List[int], - ) -> Tuple[CustomCodeNode, sp.MatrixBase]: - return ( - CustomCodeNode("", [], []), - sp.Identity(self.fe_space.num_dofs(geometry)), - ) - -class P1FunctionSpaceImpl(FunctionSpaceImpl): - def __init__( - self, - fe_space: FunctionSpace, - name: str, - type_descriptor: HOGType, - is_pointer: bool = False, - ): - super().__init__(fe_space, name, type_descriptor, is_pointer) - self.field = self._create_field(name) - - def pre_communication(self, dim: int) -> str: - if dim == 2: - return f"communication::syncFunctionBetweenPrimitives( {self.name}, level, communication::syncDirection_t::LOW2HIGH );" - else: - return ( - f"{self._deref()}.communicate< Face, Cell >( level );\n" - f"{self._deref()}.communicate< Edge, Cell >( level );\n" - f"{self._deref()}.communicate< Vertex, Cell >( level );" - ) - - def zero_halos(self, dim: int) -> str: - if dim == 2: - return ( - f"for ( const auto& idx : vertexdof::macroface::Iterator( level ) ) {{\n" - f" if ( vertexdof::macroface::isVertexOnBoundary( level, idx ) ) {{\n" - f" auto arrayIdx = vertexdof::macroface::index( level, idx.x(), idx.y() );\n" - f" _data_{self.name}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" }}\n" - f"}}" - ) - else: - return ( - f"for ( const auto& idx : vertexdof::macrocell::Iterator( level ) ) {{\n" - f" if ( !vertexdof::macrocell::isOnCellFace( idx, level ).empty() ) {{\n" - f" auto arrayIdx = vertexdof::macrocell::index( level, idx.x(), idx.y(), idx.z() );\n" - f" _data_{self.name}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" }}\n" - f"}}" - ) - - def post_communication( - self, dim: int, params: str, transform_basis: bool = True - ) -> str: - if dim == 2: - return ( - f"{self._deref()}.communicateAdditively < Face, Edge > ( {params} );\n" - f"{self._deref()}.communicateAdditively < Face, Vertex > ( {params} );" - ) - else: - return ( - f"{self._deref()}.communicateAdditively< Cell, Face >( {params} );\n" - f"{self._deref()}.communicateAdditively< Cell, Edge >( {params} );\n" - f"{self._deref()}.communicateAdditively< Cell, Vertex >( {params} );" - ) - - def pointer_retrieval(self, dim: int) -> str: - """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" - Macro = {2: "Face", 3: "Cell"}[dim] - macro = {2: "face", 3: "cell"}[dim] - - return f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {macro}.getData( {self._deref()}.get{Macro}DataID() )->getPointer( level );" - - def invert_elementwise(self, dim: int) -> str: - return f"{self._deref()}.invertElementwise( level );" - - def local_dofs( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - indexing_info: IndexingInfo, - element_vertex_ordering: List[int], - ) -> List[Field.Access]: - vertex_dof_indices = micro_element_to_vertex_indices( - geometry, element_type, element_index - ) - - vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] - - vertex_array_indices = [ - dof_idx.array_index(geometry, indexing_info) - for dof_idx in vertex_dof_indices - ] - - return [ - self.field.absolute_access((idx,), (0,)) for idx in vertex_array_indices - ] - - def func_type_string(self) -> str: - return f"P1Function< {self.type_descriptor.pystencils_type} >" - - def includes(self) -> Set[str]: - return {f"hyteg/p1functionspace/P1Function.hpp"} - - def dof_transformation( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - element_vertex_ordering: List[int], - ) -> Tuple[CustomCodeNode, sp.MatrixBase]: - return ( - CustomCodeNode("", [], []), - sp.Identity(self.fe_space.num_dofs(geometry)), - ) - - -class P2FunctionSpaceImpl(FunctionSpaceImpl): - def __init__( - self, - fe_space: FunctionSpace, - name: str, - type_descriptor: HOGType, - is_pointer: bool = False, - ): - super().__init__(fe_space, name, type_descriptor, is_pointer) - self.f_vertex = self._create_field(name + "Vertex") - self.f_edge = self._create_field(name + "Edge") - - def pre_communication(self, dim: int) -> str: - if dim == 2: - return f"communication::syncFunctionBetweenPrimitives( {self.name}, level, communication::syncDirection_t::LOW2HIGH );" - else: - return ( - f"{self._deref()}.communicate< Face, Cell >( level );\n" - f"{self._deref()}.communicate< Edge, Cell >( level );\n" - f"{self._deref()}.communicate< Vertex, Cell >( level );" - ) - - def zero_halos(self, dim: int) -> str: - if dim == 2: - return ( - f"for ( const auto& idx : vertexdof::macroface::Iterator( level ) )\n" - f"{{\n" - f" if ( vertexdof::macroface::isVertexOnBoundary( level, idx ) )\n" - f" {{\n" - f" auto arrayIdx = vertexdof::macroface::index( level, idx.x(), idx.y() );\n" - f" _data_{self.name }Vertex[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" - f" }}\n" - f"}}\n" - f"for ( const auto& idx : edgedof::macroface::Iterator( level ) )\n" - f"{{\n" - f" for ( const auto& orientation : edgedof::faceLocalEdgeDoFOrientations )\n" - f" {{\n" - f" if ( !edgedof::macroface::isInnerEdgeDoF( level, idx, orientation ) )\n" - f" {{\n" - f" auto arrayIdx = edgedof::macroface::index( level, idx.x(), idx.y(), orientation );\n" - f" _data_{self.name }Edge[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" - f" }}\n" - f" }}\n" - f"}}" - ) - else: - return ( - f"for ( const auto& idx : vertexdof::macrocell::Iterator( level ) )\n" - f"{{\n" - f"if ( !vertexdof::macrocell::isOnCellFace( idx, level ).empty() )\n" - f" {{\n" - f" auto arrayIdx = vertexdof::macrocell::index( level, idx.x(), idx.y(), idx.z() );\n" - f" _data_{self.name}Vertex[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" }}\n" - f"}}\n" - f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}.getEdgeDoFFunction().getCellDataID() );" - ) - - def post_communication( - self, dim: int, params: str, transform_basis: bool = True - ) -> str: - if dim == 2: - return ( - f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Face, Edge >( {params} );\n" - f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Face, Vertex >( {params} );\n" - f"{self._deref()}.getEdgeDoFFunction().communicateAdditively< Face, Edge >( {params} );" - ) - else: - return ( - f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Cell, Face >( {params} );\n" - f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Cell, Edge >( {params} );\n" - f"{self._deref()}.getVertexDoFFunction().communicateAdditively< Cell, Vertex >( {params} );\n" - f"{self._deref()}.getEdgeDoFFunction().communicateAdditively< Cell, Face >( {params} );\n" - f"{self._deref()}.getEdgeDoFFunction().communicateAdditively< Cell, Edge >( {params} );" - ) - - def pointer_retrieval(self, dim: int) -> str: - """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" - Macro = {2: "Face", 3: "Cell"}[dim] - macro = {2: "face", 3: "cell"}[dim] - - return ( - f"{self.type_descriptor.pystencils_type}* _data_{self.name}Vertex = {macro}.getData( {self._deref()}.getVertexDoFFunction().get{Macro}DataID() )->getPointer( level );\n" - f"{self.type_descriptor.pystencils_type}* _data_{self.name}Edge = {macro}.getData( {self._deref()}.getEdgeDoFFunction().get{Macro}DataID() )->getPointer( level );" - ) - - def invert_elementwise(self, dim: int) -> str: - return f"{self._deref()}.invertElementwise( level );" - - def local_dofs( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - indexing_info: IndexingInfo, - element_vertex_ordering: List[int], - ) -> List[Field.Access]: - vertex_dof_indices = micro_element_to_vertex_indices( - geometry, element_type, element_index - ) - - vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] - - edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices) - - vrtx_array_idcs = [ - dof_idx.array_index(geometry, indexing_info) - for dof_idx in vertex_dof_indices - ] - edge_array_idcs = [ - dof_idx.array_index(geometry, indexing_info) for dof_idx in edge_dof_indices - ] - - # NOTE: The order of DoFs must match the order of shape functions - # defined in the FunctionSpace. Holds within the individual lists - # and overall. - return [ - self.f_vertex.absolute_access((idx,), (0,)) for idx in vrtx_array_idcs - ] + [self.f_edge.absolute_access((idx,), (0,)) for idx in edge_array_idcs] - - def func_type_string(self) -> str: - return f"P2Function< {self.type_descriptor.pystencils_type} >" - - def includes(self) -> Set[str]: - return {f"hyteg/p2functionspace/P2Function.hpp"} - - def dof_transformation( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - element_vertex_ordering: List[int], - ) -> Tuple[CustomCodeNode, sp.MatrixBase]: - return ( - CustomCodeNode("", [], []), - sp.Identity(self.fe_space.num_dofs(geometry)), - ) - - -class P1VectorFunctionSpaceImpl(FunctionSpaceImpl): - def __init__( - self, - fe_space: FunctionSpace, - name: str, - type_descriptor: HOGType, - is_pointer: bool = False, - ): - super().__init__(fe_space, name, type_descriptor, is_pointer) - self.fields: Dict[int, Field] = {} - - def _field_name(self, component: int) -> str: - return self.name + f"_{component}" - - def _raw_pointer_name(self, component: int) -> str: - return f"_data_" + self._field_name(component) - - def pre_communication(self, dim: int) -> str: - if dim == 2: - return f"communication::syncVectorFunctionBetweenPrimitives( {self.name}, level, communication::syncDirection_t::LOW2HIGH );" - else: - ret_str = "" - for i in range(dim): - ret_str += ( - f"{self._deref()}[{i}].communicate< Face, Cell >( level );\n" - f"{self._deref()}[{i}].communicate< Edge, Cell >( level );\n" - f"{self._deref()}[{i}].communicate< Vertex, Cell >( level );" - ) - return ret_str - - def zero_halos(self, dim: int) -> str: - if dim == 2: - return ( - f"for ( const auto& idx : vertexdof::macroface::Iterator( level ) ) {{\n" - f" if ( vertexdof::macroface::isVertexOnBoundary( level, idx ) ) {{\n" - f" auto arrayIdx = vertexdof::macroface::index( level, idx.x(), idx.y() );\n" - f" {self._raw_pointer_name(0)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" {self._raw_pointer_name(1)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" }}\n" - f"}}" - ) - else: - return ( - f"for ( const auto& idx : vertexdof::macrocell::Iterator( level ) ) {{\n" - f" if ( !vertexdof::macrocell::isOnCellFace( idx, level ).empty() ) {{\n" - f" auto arrayIdx = vertexdof::macrocell::index( level, idx.x(), idx.y(), idx.z() );\n" - f" {self._raw_pointer_name(0)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" {self._raw_pointer_name(1)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" {self._raw_pointer_name(2)}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" }}\n" - f"}}" - ) - - def post_communication( - self, dim: int, params: str, transform_basis: bool = True - ) -> str: - if dim == 2: - ret_str = "" - for i in range(dim): - ret_str += ( - f"{self._deref()}[{i}].communicateAdditively < Face, Edge > ( {params} );\n" - f"{self._deref()}[{i}].communicateAdditively < Face, Vertex > ( {params} );\n" - ) - return ret_str - else: - ret_str = "" - for i in range(dim): - ret_str += ( - f"{self._deref()}[{i}].communicateAdditively< Cell, Face >( {params} );\n" - f"{self._deref()}[{i}].communicateAdditively< Cell, Edge >( {params} );\n" - f"{self._deref()}[{i}].communicateAdditively< Cell, Vertex >( {params} );" - ) - return ret_str - - def pointer_retrieval(self, dim: int) -> str: - """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" - Macro = {2: "Face", 3: "Cell"}[dim] - macro = {2: "face", 3: "cell"}[dim] - - ret_str = "" - for i in range(dim): - ret_str += f"{self.type_descriptor.pystencils_type}* {self._raw_pointer_name(i)} = {macro}.getData( {self._deref()}[{i}].get{Macro}DataID() )->getPointer( level );\n" - return ret_str - - def invert_elementwise(self, dim: int) -> str: - ret_str = "" - for i in range(dim): - ret_str += f"{self._deref()}[{i}].invertElementwise( level );\n" - return ret_str - - def local_dofs( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - indexing_info: IndexingInfo, - element_vertex_ordering: List[int], - ) -> List[Field.Access]: - """ - Returns the element-local DoFs (field accesses) in a list (i.e., linearized). - - The order here has to match that of the function space implementation. - TODO: This is a little concerning since that order has to match in two very different parts of this - implementation. Here, and in the TensorialVectorFunctionSpace. Maybe we can define this order in a single - location. - We return them in "SoA" order: first all DoFs corresponding to the function with the first component != 0, - then all DoFs corresponding to the function with the second component != 0 etc. - - For instance, in 2D we get the DoFs corresponding to the 6 shape functions in the following order: - - First component != 0 - - list[0]: ⎡-x_ref_0 - x_ref_1 + 1⎤ - ⎢ ⎥ - ⎣ 0 ⎦ - - list[1]: ⎡x_ref_0⎤ - ⎢ ⎥ - ⎣ 0 ⎦ - - list[2]: ⎡x_ref_1⎤ - ⎢ ⎥ - ⎣ 0 ⎦ - - Second component != 0 - - list[3]: ⎡ 0 ⎤ - ⎢ ⎥ - ⎣-x_ref_0 - x_ref_1 + 1⎦ - - list[4]: ⎡ 0 ⎤ - ⎢ ⎥ - ⎣x_ref_0⎦ - - list[5]: ⎡ 0 ⎤ - ⎢ ⎥ - ⎣x_ref_1⎦ - """ - - for c in range(geometry.dimensions): - if c not in self.fields: - self.fields[c] = self._create_field(self._field_name(c)) - - vertex_dof_indices = micro_element_to_vertex_indices( - geometry, element_type, element_index - ) - - vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] - - vertex_array_indices = [ - dof_idx.array_index(geometry, indexing_info) - for dof_idx in vertex_dof_indices - ] - - return [ - self.fields[c].absolute_access((idx,), (0,)) - for c in range(geometry.dimensions) - for idx in vertex_array_indices - ] - - def func_type_string(self) -> str: - return f"P1VectorFunction< {self.type_descriptor.pystencils_type} >" - - def includes(self) -> Set[str]: - return {f"hyteg/p1functionspace/P1VectorFunction.hpp"} - - def dof_transformation( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - element_vertex_ordering: List[int], - ) -> Tuple[CustomCodeNode, sp.MatrixBase]: - return ( - CustomCodeNode("", [], []), - sp.Identity(self.fe_space.num_dofs(geometry)), - ) - - -class P2VectorFunctionSpaceImpl(FunctionSpaceImpl): - def __init__( - self, - fe_space: FunctionSpace, - name: str, - type_descriptor: HOGType, - is_pointer: bool = False, - ): - super().__init__(fe_space, name, type_descriptor, is_pointer) - self.fields: Dict[Tuple[str, int], Field] = {} - - def _field_name(self, component: int, dof_type: str) -> str: - """dof_type should either be 'vertex' or 'edge'""" - return self.name + f"_{dof_type}_{component}" - - def _raw_pointer_name(self, component: int, dof_type: str) -> str: - """dof_type should either be 'vertex' or 'edge'""" - return f"_data_" + self._field_name(component, dof_type) - - def pre_communication(self, dim: int) -> str: - if dim == 2: - return f"communication::syncVectorFunctionBetweenPrimitives( {self.name}, level, communication::syncDirection_t::LOW2HIGH );" - else: - ret_str = "" - for i in range(dim): - ret_str += ( - f"{self._deref()}[{i}].communicate< Face, Cell >( level );\n" - f"{self._deref()}[{i}].communicate< Edge, Cell >( level );\n" - f"{self._deref()}[{i}].communicate< Vertex, Cell >( level );" - ) - return ret_str - - def zero_halos(self, dim: int) -> str: - if dim == 2: - return ( - f"for ( const auto& idx : vertexdof::macroface::Iterator( level ) )\n" - f"{{\n" - f" if ( vertexdof::macroface::isVertexOnBoundary( level, idx ) )\n" - f" {{\n" - f" auto arrayIdx = vertexdof::macroface::index( level, idx.x(), idx.y() );\n" - f" {self._raw_pointer_name(0, 'vertex')}[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" - f" {self._raw_pointer_name(1, 'vertex')}[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" - f" }}\n" - f"}}\n" - f"for ( const auto& idx : edgedof::macroface::Iterator( level ) )\n" - f"{{\n" - f" for ( const auto& orientation : edgedof::faceLocalEdgeDoFOrientations )\n" - f" {{\n" - f" if ( !edgedof::macroface::isInnerEdgeDoF( level, idx, orientation ) )\n" - f" {{\n" - f" auto arrayIdx = edgedof::macroface::index( level, idx.x(), idx.y(), orientation );\n" - f" {self._raw_pointer_name(0, 'edge')}[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" - f" {self._raw_pointer_name(1, 'edge')}[arrayIdx] = walberla::numeric_cast< {self.type_descriptor.pystencils_type} >( 0 );\n" - f" }}\n" - f" }}\n" - f"}}" - ) - else: - return ( - f"for ( const auto& idx : vertexdof::macrocell::Iterator( level ) )\n" - f"{{\n" - f"if ( !vertexdof::macrocell::isOnCellFace( idx, level ).empty() )\n" - f" {{\n" - f" auto arrayIdx = vertexdof::macrocell::index( level, idx.x(), idx.y(), idx.z() );\n" - f" {self._raw_pointer_name(0, 'vertex')}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" {self._raw_pointer_name(1, 'vertex')}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" {self._raw_pointer_name(2, 'vertex')}[arrayIdx] = {self.type_descriptor.pystencils_type}( 0 );\n" - f" }}\n" - f"}}\n" - f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}[0].getEdgeDoFFunction().getCellDataID() );\n" - f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}[1].getEdgeDoFFunction().getCellDataID() );\n" - f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}[2].getEdgeDoFFunction().getCellDataID() );\n" - ) - - def post_communication( - self, dim: int, params: str, transform_basis: bool = True - ) -> str: - if dim == 2: - ret_str = "" - for i in range(dim): - ret_str += ( - f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Face, Edge >( {params} );\n" - f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Face, Vertex >( {params} );\n" - f"{self._deref()}[{i}].getEdgeDoFFunction().communicateAdditively< Face, Edge >( {params} );" - ) - return ret_str - else: - ret_str = "" - for i in range(dim): - ret_str += ( - f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Cell, Face >( {params} );\n" - f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Cell, Edge >( {params} );\n" - f"{self._deref()}[{i}].getVertexDoFFunction().communicateAdditively< Cell, Vertex >( {params} );\n" - f"{self._deref()}[{i}].getEdgeDoFFunction().communicateAdditively< Cell, Face >( {params} );\n" - f"{self._deref()}[{i}].getEdgeDoFFunction().communicateAdditively< Cell, Edge >( {params} );" - ) - return ret_str - - def pointer_retrieval(self, dim: int) -> str: - """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" - Macro = {2: "Face", 3: "Cell"}[dim] - macro = {2: "face", 3: "cell"}[dim] - - ret_str = "" - for i in range(dim): - ret_str += f"{self.type_descriptor.pystencils_type}* {self._raw_pointer_name(i, 'vertex')} = {macro}.getData( {self._deref()}[{i}].getVertexDoFFunction().get{Macro}DataID() )->getPointer( level );\n" - ret_str += f"{self.type_descriptor.pystencils_type}* {self._raw_pointer_name(i, 'edge')} = {macro}.getData( {self._deref()}[{i}].getEdgeDoFFunction().get{Macro}DataID() )->getPointer( level );\n" - return ret_str - - def invert_elementwise(self, dim: int) -> str: - ret_str = "" - for i in range(dim): - ret_str += f"{self._deref()}[{i}].invertElementwise( level );\n" - return ret_str - - def local_dofs( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - indexing_info: IndexingInfo, - element_vertex_ordering: List[int], - ) -> List[Field.Access]: - """ - Returns the element-local DoFs (field accesses) in a list (i.e., linearized). - - See P1VectorFunctionSpaceImpl::local_dofs() for details. - """ - - for dt in ["vertex", "edge"]: - for c in range(geometry.dimensions): - if (dt, c) not in self.fields: - self.fields[(dt, c)] = self._create_field(self._field_name(c, dt)) - - vertex_dof_indices = micro_element_to_vertex_indices( - geometry, element_type, element_index - ) - - vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] - - edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices) - - vertex_array_indices = [ - dof_idx.array_index(geometry, indexing_info) - for dof_idx in vertex_dof_indices - ] - - edge_array_indices = [ - dof_idx.array_index(geometry, indexing_info) for dof_idx in edge_dof_indices - ] - - loc_dofs = [] - - for c in range(geometry.dimensions): - loc_dofs += [ - self.fields[("vertex", c)].absolute_access((idx,), (0,)) - for idx in vertex_array_indices - ] - loc_dofs += [ - self.fields[("edge", c)].absolute_access((idx,), (0,)) - for idx in edge_array_indices - ] - - return loc_dofs - - def func_type_string(self) -> str: - return f"P2VectorFunction< {self.type_descriptor.pystencils_type} >" - - def includes(self) -> Set[str]: - return { - f"hyteg/p2functionspace/P2VectorFunction.hpp", - f"hyteg/p2functionspace/P2Function.hpp", - } - - def dof_transformation( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - element_vertex_ordering: List[int], - ) -> Tuple[CustomCodeNode, sp.MatrixBase]: - return ( - CustomCodeNode("", [], []), - sp.Identity(self.fe_space.num_dofs(geometry)), - ) - - -class N1E1FunctionSpaceImpl(FunctionSpaceImpl): - def __init__( - self, - fe_space: FunctionSpace, - name: str, - type_descriptor: HOGType, - is_pointer: bool = False, - ): - super().__init__(fe_space, name, type_descriptor, is_pointer) - self.field = self._create_field(name) - - def pre_communication(self, dim: int) -> str: - if dim == 2: - return "" - else: - return ( - f"{self._deref()}.communicate< Face, Cell >( level );\n" - f"{self._deref()}.communicate< Edge, Cell >( level );" - ) - - def zero_halos(self, dim: int) -> str: - if dim == 2: - return "" - else: - return f"edgedof::macrocell::setBoundaryToZero( level, cell, {self._deref()}.getDoFs()->getCellDataID() );" - - def post_communication( - self, dim: int, params: str, transform_basis: bool = True - ) -> str: - if dim == 2: - return "" - else: - dofs = "" - if not transform_basis: - dofs = "getDoFs()->" - return ( - f"{self._deref()}.{dofs}communicateAdditively< Cell, Face >( {params} );\n" - f"{self._deref()}.{dofs}communicateAdditively< Cell, Edge >( {params} );" - ) - - def pointer_retrieval(self, dim: int) -> str: - """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d).""" - Macro = {2: "Face", 3: "Cell"}[dim] - macro = {2: "face", 3: "cell"}[dim] - - return f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {macro}.getData( {self._deref()}.getDoFs()->get{Macro}DataID() )->getPointer( level );" - - def invert_elementwise(self, dim: int) -> str: - return f"{self._deref()}.getDoFs()->invertElementwise( level );" - - def local_dofs( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - indexing_info: IndexingInfo, - element_vertex_ordering: List[int], - ) -> List[Field.Access]: - vertex_dof_indices = micro_element_to_vertex_indices( - geometry, element_type, element_index - ) - - vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering] - - edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices) - edge_array_indices = [ - dof_idx.array_index(geometry, indexing_info) for dof_idx in edge_dof_indices - ] - - return [self.field.absolute_access((idx,), (0,)) for idx in edge_array_indices] - - def func_type_string(self) -> str: - return f"n1e1::N1E1VectorFunction< {self.type_descriptor.pystencils_type} >" - - def includes(self) -> Set[str]: - return { - "hyteg/n1e1functionspace/N1E1VectorFunction.hpp", - "hyteg/n1e1functionspace/N1E1MacroCell.hpp", - } - - def dof_transformation( - self, - geometry: ElementGeometry, - element_index: Tuple[int, int, int], - element_type: Union[FaceType, CellType], - element_vertex_ordering: List[int], - ) -> Tuple[CustomCodeNode, sp.MatrixBase]: - - if element_vertex_ordering != [0, 1, 2, 3]: - raise HOGException( - "Element vertex re-ordering not supported for Nédélec elements (yet)." - ) - - Macro = {2: "Face", 3: "Cell"}[geometry.dimensions] - macro = {2: "face", 3: "cell"}[geometry.dimensions] - name = "basisTransformation" - n_dofs = self.fe_space.num_dofs(geometry) - - # NOTE: Types are added manually because pystencils won't add types to accessed/ - # defined symbols of custom code nodes. As a result the symbols - # in the matrix will be typed but not their definition, leading to - # undefined symbols. - # NOTE: The type is `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). - symbols = [ - ps.TypedSymbol(f"{name}.diagonal()({i})", np.float64) for i in range(n_dofs) - ] - return ( - CustomCodeNode( - f"const Eigen::DiagonalMatrix< real_t, {n_dofs} > {name} = " - f"n1e1::macro{macro}::basisTransformation( level, {macro}, {{{element_index[0]}, {element_index[1]}, {element_index[2]}}}, {macro}dof::{Macro}Type::{element_type.name} );", - [ - ps.TypedSymbol("level", "const uint_t"), - ps.TypedSymbol(f"{macro}", f"const {Macro}&"), - ], - symbols, - ), - sp.diag(*symbols), - ) diff --git a/hog/operator_generation/kernel_types.py b/hog/operator_generation/kernel_types.py index f0ba9f13523100cb3c2862dd1bcce32faea0644f..712dc0365b750c6bc6b95088e39133bcf71fde6b 100644 --- a/hog/operator_generation/kernel_types.py +++ b/hog/operator_generation/kernel_types.py @@ -38,7 +38,12 @@ from hog.cpp_printing import ( from hog.element_geometry import ElementGeometry from hog.exception import HOGException from hog.function_space import FunctionSpace, TrialSpace, TestSpace -from hog.operator_generation.function_space_impls import FunctionSpaceImpl +from hog.operator_generation.function_space_implementations.function_space_impl_base import ( + FunctionSpaceImpl, +) +from hog.operator_generation.function_space_implementations.function_space_impl_factory import ( + create_impl, +) 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 @@ -377,12 +382,8 @@ class ApplyWrapper(KernelWrapperType): dims: List[int] = [2, 3], ): self.name = "apply" - self.src: FunctionSpaceImpl = FunctionSpaceImpl.create_impl( - src_space, "src", type_descriptor - ) - self.dst: FunctionSpaceImpl = FunctionSpaceImpl.create_impl( - dst_space, "dst", type_descriptor - ) + self.src: FunctionSpaceImpl = create_impl(src_space, "src", type_descriptor) + self.dst: FunctionSpaceImpl = create_impl(dst_space, "dst", type_descriptor) self.src_fields = [self.src] self.dst_fields = [self.dst] self.dims = dims @@ -562,7 +563,7 @@ class AssembleDiagonalWrapper(KernelWrapperType): dims: List[int] = [2, 3], ): self.name = "computeInverseDiagonalOperatorValues" - self.dst: FunctionSpaceImpl = FunctionSpaceImpl.create_impl( + self.dst: FunctionSpaceImpl = create_impl( fe_space, dst_field, type_descriptor, is_pointer=True ) self.src_fields = [] @@ -690,12 +691,8 @@ class AssembleWrapper(KernelWrapperType): ): idx_t = HOGType("idx_t", np.int64) self.name = "toMatrix" - self.src: FunctionSpaceImpl = FunctionSpaceImpl.create_impl( - src_space, "src", idx_t - ) - self.dst: FunctionSpaceImpl = FunctionSpaceImpl.create_impl( - dst_space, "dst", idx_t - ) + self.src: FunctionSpaceImpl = create_impl(src_space, "src", idx_t) + self.dst: FunctionSpaceImpl = create_impl(dst_space, "dst", idx_t) # Treating both src and dst as dst_fields because src_fields are loaded # explicitly from memory prior to the kernel operation but we do not diff --git a/hog/operator_generation/operators.py b/hog/operator_generation/operators.py index b137074b28e9c55a3edc791822939857d48b306a..ff104c33a346c2c2f48e5a831f0d43c0d03106a2 100644 --- a/hog/operator_generation/operators.py +++ b/hog/operator_generation/operators.py @@ -45,7 +45,12 @@ from hog.hyteg_code_generation import ( GCC_WARNING_WORKAROUND, ) from hog.logger import TimedLogger -from hog.operator_generation.function_space_impls import FunctionSpaceImpl +from hog.operator_generation.function_space_implementations.function_space_impl_base import ( + FunctionSpaceImpl, +) +from hog.operator_generation.function_space_implementations.function_space_impl_factory import ( + create_impl, +) from hog.operator_generation.pystencils_extensions import create_generic_fields import pystencils as ps @@ -219,7 +224,6 @@ def micro_vertex_permutation_for_facet( """ if volume_geometry == TriangleElement(): - if element_type == FaceType.BLUE: return [0, 1, 2] @@ -232,7 +236,6 @@ def micro_vertex_permutation_for_facet( return shuffle_order_gray[facet_id] elif volume_geometry == TetrahedronElement(): - if element_type == CellType.WHITE_DOWN: return [0, 1, 2, 3] @@ -557,7 +560,6 @@ class HyTeGElementwiseOperator: with TimedLogger( f"Generating kernels for operator {self.name}", level=logging.INFO ): - # Generate each kernel type (apply, gemv, ...). self.generate_kernels() @@ -602,7 +604,6 @@ class HyTeGElementwiseOperator: with TimedLogger("Generating C code from kernel AST(s)"): # Add all kernels to the class. for operator_method in self.operator_methods: - num_integrals = len(operator_method.integration_infos) if num_integrals != len( @@ -798,7 +799,6 @@ class HyTeGElementwiseOperator: ) blending_includes = set() for dim, integration_infos in self.integration_infos.items(): - if not all( [ integration_infos[0].blending.coupling_includes() @@ -1194,7 +1194,6 @@ class HyTeGElementwiseOperator: element_types = list(integration_info.loop_strategy.element_loops.keys()) for element_type in element_types: - # Re-ordering micro-element vertices for the handling of domain boundary integrals. # # Boundary integrals are handled by looping over all (volume-)elements that have a facet at one of the @@ -1283,7 +1282,7 @@ class HyTeGElementwiseOperator: coeffs = dict( ( dof_symbol.function_id, - FunctionSpaceImpl.create_impl( + create_impl( dof_symbol.function_space, dof_symbol.function_id, self._type_descriptor, @@ -1299,7 +1298,7 @@ class HyTeGElementwiseOperator: sp.Symbol(dof_symbol.name), coeffs[dof_symbol.function_id].local_dofs( geometry, - element_index, + element_index, # type: ignore[arg-type] # list of sympy expressions also works element_type, indexing_info, element_vertex_order, @@ -1530,7 +1529,6 @@ class HyTeGElementwiseOperator: for kernel_wrapper_type in self.kernel_wrapper_types: for dim, integration_infos in self.integration_infos.items(): - kernel_functions = [] kernel_op_counts = [] platform_dep_kernels = [] @@ -1546,13 +1544,11 @@ class HyTeGElementwiseOperator: ) for integration_info in integration_infos: - # generate AST of kernel loop with TimedLogger( f"Generating kernel {integration_info.name} ({kernel_wrapper_type.name}, {dim}D)", logging.INFO, ): - ( function_body, kernel_op_count, @@ -1655,7 +1651,6 @@ class HyTeGElementwiseOperator: for kernel_function, integration_info in zip( kernel_functions, integration_infos ): - pre_call_code = "" post_call_code = "" @@ -1663,7 +1658,6 @@ class HyTeGElementwiseOperator: integration_info.integration_domain == MacroIntegrationDomain.DOMAIN_BOUNDARY ): - if not isinstance(integration_info.loop_strategy, BOUNDARY): raise HOGException( "The loop strategy should be BOUNDARY for boundary integrals." diff --git a/hyteg_integration_tests/src/OperatorGenerationTest.hpp b/hyteg_integration_tests/src/OperatorGenerationTest.hpp index 8719f6d5ef62fbed630bedb3fc89c39deb71f1d7..0e66c402acf55b0bfc0a6f673519981e03265c75 100644 --- a/hyteg_integration_tests/src/OperatorGenerationTest.hpp +++ b/hyteg_integration_tests/src/OperatorGenerationTest.hpp @@ -220,7 +220,7 @@ void compareApply( OperatorFactory< RefOpType > refOpFactory, real_t maxAbs = 1.0; if constexpr ( std::is_same< RefDstFncType, n1e1::N1E1VectorFunction< RefType > >::value ) { - maxAbs = err.getDoFs()->getMaxMagnitude( level ); + maxAbs = err.getDoFs()->getMaxDoFMagnitude( level ); } else if constexpr ( std::is_same_v< RefDstFncType, P1VectorFunction< RefType > > || std::is_same_v< RefDstFncType, P2VectorFunction< RefType > > ) { @@ -228,7 +228,7 @@ void compareApply( OperatorFactory< RefOpType > refOpFactory, } else { - maxAbs = err.getMaxMagnitude( level ); + maxAbs = err.getMaxDoFMagnitude( level ); } const real_t threshold = thresholdOverMachineEps * precisionDict::epsilon< TestType >(); @@ -306,7 +306,7 @@ void compareGEMV( const uint_t level, real_t maxAbs = 1.0; if constexpr ( std::is_same< DstFncType, n1e1::N1E1VectorFunction< real_t > >::value ) { - maxAbs = err.getDoFs()->getMaxMagnitude( level ); + maxAbs = err.getDoFs()->getMaxDoFMagnitude( level ); } else if constexpr ( std::is_same_v< DstFncType, P1VectorFunction< real_t > > || std::is_same_v< DstFncType, P2VectorFunction< real_t > > ) { @@ -314,7 +314,7 @@ void compareGEMV( const uint_t level, } else { - maxAbs = err.getMaxMagnitude( level ); + maxAbs = err.getMaxDoFMagnitude( level ); } const real_t threshold = thresholdOverMachineEps * precisionDict::epsilon< TestType >(); @@ -393,7 +393,7 @@ void compareInvDiag( OperatorFactory< RefOpType > refOpFactory, real_t maxAbs = 1.0; if constexpr ( std::is_same< RefDiagType, n1e1::N1E1VectorFunction< RefType > >::value ) { - maxAbs = err.getDoFs()->getMaxMagnitude( level ); + maxAbs = err.getDoFs()->getMaxDoFMagnitude( level ); } else if constexpr ( std::is_same_v< RefDiagType, P1VectorFunction< RefType > > || std::is_same_v< RefDiagType, P2VectorFunction< RefType > > ) { @@ -401,7 +401,7 @@ void compareInvDiag( OperatorFactory< RefOpType > refOpFactory, } else { - maxAbs = err.getMaxMagnitude( level ); + maxAbs = err.getMaxDoFMagnitude( level ); } const real_t threshold = thresholdOverMachineEps * precisionDict::epsilon< TestType >();