diff --git a/hog/operator_generation/function_space_impls.py b/hog/operator_generation/function_space_impls.py
index 7b3a51fe8857c1a2c7b476ee7a9a143cbfda8e4b..fdad32bea2f71bb248b8c633c674dd4cba2c925e 100644
--- a/hog/operator_generation/function_space_impls.py
+++ b/hog/operator_generation/function_space_impls.py
@@ -36,8 +36,10 @@ from hog.function_space import (
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
@@ -98,6 +100,8 @@ class FunctionSpaceImpl(ABC):
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):
@@ -241,6 +245,113 @@ class FunctionSpaceImpl(ABC):
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 ""
+ # 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:
+ return ""
+ # 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:
+ return ""
+ # 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} = {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]:
+ vertex_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 vertex_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__(
diff --git a/hog/operator_generation/indexing.py b/hog/operator_generation/indexing.py
index e33704ecc6112ec192671029c69e88ed92e66095..91620ad4bd8f2a53f60e2831b3bda3d507bb7d8a 100644
--- a/hog/operator_generation/indexing.py
+++ b/hog/operator_generation/indexing.py
@@ -150,7 +150,7 @@ def generalized_macro_cell_index(
def num_microvertices_per_face_from_width(width: int) -> int:
"""Computes the number of microvertices in a refined macro triangle. width depends on the level and quantifies the amount of primitives in one direction of the refined triangle."""
- return (width * (width + 1)) >> 1
+ return int_div((width * (width + 1)), 2)
def num_microvertices_per_cell_from_width(width: int) -> int:
@@ -197,20 +197,20 @@ def num_faces_per_row_by_type(
def num_cells_per_row_by_type(
- level: int, cellType: Union[None, EdgeType, FaceType, CellType]
+ level: int, cellType: Union[None, EdgeType, FaceType, CellType], num_microedges_per_edge: sp.Symbol
) -> int:
if cellType == CellType.WHITE_UP:
- return num_microedges_per_edge(level)
+ return num_microedges_per_edge
elif cellType == CellType.GREEN_UP:
- return num_microedges_per_edge(level) - 1
+ return num_microedges_per_edge - 1
elif cellType == CellType.BLUE_UP:
- return num_microedges_per_edge(level) - 1
+ return num_microedges_per_edge - 1
elif cellType == CellType.WHITE_DOWN:
- return num_microedges_per_edge(level) - 2
+ return num_microedges_per_edge - 2
elif cellType == CellType.BLUE_DOWN:
- return num_microedges_per_edge(level) - 1
+ return num_microedges_per_edge - 1
elif cellType == CellType.GREEN_DOWN:
- return num_microedges_per_edge(level) - 1
+ return num_microedges_per_edge - 1
else:
raise HOGException(f"Unexpected cell type: {cellType}")
@@ -221,9 +221,9 @@ def num_micro_faces_per_macro_face(level: int, faceType: FaceType) -> int:
)
-def num_micro_cells_per_macro_cell(level: int, cellType: CellType) -> int:
+def num_micro_cells_per_macro_cell(level: int, cellType: CellType, num_microedges_per_edge: sp.Symbol) -> int:
return num_microvertices_per_cell_from_width(
- num_cells_per_row_by_type(level, cellType)
+ num_cells_per_row_by_type(level, cellType, num_microedges_per_edge)
)
@@ -254,16 +254,16 @@ def facedof_index(
level: int,
index: Tuple[int, int, int],
faceType: Union[None, EdgeType, FaceType, CellType],
+ num_microfaces_per_face: sp.Symbol,
+ num_microedges_per_edge: sp.Symbol
) -> int:
"""Indexes triangles/faces. Used to compute offsets in volume dof indexing in 2D and AoS layout."""
- x, y, _ = index
- width = num_faces_per_row_by_type(level, faceType)
+ x, y = index
+ # width = num_faces_per_row_by_type(level, faceType)
if faceType == FaceType.GRAY:
- return linear_macro_face_index(width, x, y)
+ return linear_macro_face_index(num_microedges_per_edge, x, y)
elif faceType == FaceType.BLUE:
- return num_micro_faces_per_macro_face(
- level, FaceType.GRAY
- ) + linear_macro_face_index(width, x, y)
+ return num_microvertices_per_face_from_width(num_microedges_per_edge) + linear_macro_face_index(num_microedges_per_edge - 1, x, y)
else:
raise HOGException(f"Unexpected face type: {faceType}")
@@ -272,44 +272,45 @@ def celldof_index(
level: int,
index: Tuple[int, int, int],
cellType: Union[None, EdgeType, FaceType, CellType],
+ num_microedges_per_edge: sp.Symbol,
) -> int:
"""Indexes cells/tetrahedra. Used to compute offsets in volume dof indexing in 3D and AoS layout."""
- x, y, z = index
- width = num_cells_per_row_by_type(level, cellType) # gives expr(level)
+ x, y, z = index
+ width = num_cells_per_row_by_type(level, cellType, num_microedges_per_edge) # gives expr(level)
if cellType == CellType.WHITE_UP:
return linear_macro_cell_index(width, x, y, z)
elif cellType == CellType.BLUE_UP:
return num_micro_cells_per_macro_cell( # gives expr(level)
- level, CellType.WHITE_UP
+ level, CellType.WHITE_UP, num_microedges_per_edge
) + linear_macro_cell_index(width, x, y, z)
elif cellType == CellType.GREEN_UP:
return (
- num_micro_cells_per_macro_cell(level, CellType.WHITE_UP)
- + num_micro_cells_per_macro_cell(level, CellType.BLUE_UP)
+ num_micro_cells_per_macro_cell(level, CellType.WHITE_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.BLUE_UP, num_microedges_per_edge)
+ linear_macro_cell_index(width, x, y, z)
)
elif cellType == CellType.WHITE_DOWN:
return (
- num_micro_cells_per_macro_cell(level, CellType.WHITE_UP)
- + num_micro_cells_per_macro_cell(level, CellType.BLUE_UP)
- + num_micro_cells_per_macro_cell(level, CellType.GREEN_UP)
+ num_micro_cells_per_macro_cell(level, CellType.WHITE_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.BLUE_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.GREEN_UP, num_microedges_per_edge)
+ linear_macro_cell_index(width, x, y, z)
)
elif cellType == CellType.BLUE_DOWN:
return (
- num_micro_cells_per_macro_cell(level, CellType.WHITE_UP)
- + num_micro_cells_per_macro_cell(level, CellType.BLUE_UP)
- + num_micro_cells_per_macro_cell(level, CellType.GREEN_UP)
- + num_micro_cells_per_macro_cell(level, CellType.WHITE_DOWN)
+ num_micro_cells_per_macro_cell(level, CellType.WHITE_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.BLUE_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.GREEN_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.WHITE_DOWN, num_microedges_per_edge)
+ linear_macro_cell_index(width, x, y, z)
)
elif cellType == CellType.GREEN_DOWN:
return (
- num_micro_cells_per_macro_cell(level, CellType.WHITE_UP)
- + num_micro_cells_per_macro_cell(level, CellType.BLUE_UP)
- + num_micro_cells_per_macro_cell(level, CellType.GREEN_UP)
- + num_micro_cells_per_macro_cell(level, CellType.WHITE_DOWN)
- + num_micro_cells_per_macro_cell(level, CellType.BLUE_DOWN)
+ num_micro_cells_per_macro_cell(level, CellType.WHITE_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.BLUE_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.GREEN_UP, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.WHITE_DOWN, num_microedges_per_edge)
+ + num_micro_cells_per_macro_cell(level, CellType.BLUE_DOWN, num_microedges_per_edge)
+ linear_macro_cell_index(width, x, y, z)
)
else:
@@ -326,6 +327,7 @@ class IndexingInfo:
self.micro_edges_per_macro_edge = TypedSymbol("micro_edges_per_macro_edge", int)
self.num_microfaces_per_face = TypedSymbol("num_microfaces_per_face", int)
self.num_microcells_per_cell = TypedSymbol("num_microcells_per_cell", int)
+ # self.num_microedges_per_edge = TypedSymbol("num_microedges_per_edge", int)
self.micro_edges_per_macro_edge_float = sp.Symbol(
"micro_edges_per_macro_edge_float"
)
@@ -445,7 +447,7 @@ class DoFIndex:
numMicroVolumes = indexing_info.num_microfaces_per_face
microVolume = facedof_index(
- indexing_info.level, self.primitive_index, self.dof_sub_type
+ indexing_info.level, self.primitive_index, self.dof_sub_type, indexing_info.num_microfaces_per_face, indexing_info.micro_edges_per_macro_edge
)
if self.mem_layout == VolumeDoFMemoryLayout.SoA:
@@ -461,7 +463,7 @@ class DoFIndex:
numMicroVolumes = indexing_info.num_microcells_per_cell
microVolume = celldof_index(
- indexing_info.level, self.primitive_index, self.dof_sub_type
+ indexing_info.level, self.primitive_index, self.dof_sub_type, indexing_info.micro_edges_per_macro_edge
)
if self.mem_layout == VolumeDoFMemoryLayout.SoA:
diff --git a/hog/operator_generation/operators.py b/hog/operator_generation/operators.py
index 18dbfede7ce05b0f9b57714e3c3203fd8f7a8918..b137074b28e9c55a3edc791822939857d48b306a 100644
--- a/hog/operator_generation/operators.py
+++ b/hog/operator_generation/operators.py
@@ -274,6 +274,7 @@ class HyTeGElementwiseOperator:
"hyteg/types/types.hpp",
}
+ VAR_NAME_MICRO_FACES_PER_MACRO_FACE = "num_microfaces_per_face"
VAR_NAME_MICRO_EDGES_PER_MACRO_EDGE = "micro_edges_per_macro_edge"
VAR_NAME_MICRO_EDGES_PER_MACRO_EDGE_FLOAT = "micro_edges_per_macro_edge_float"
@@ -1615,6 +1616,7 @@ class HyTeGElementwiseOperator:
# Setting up required scalar parameters.
scalar_parameter_setup = (
f"const auto {self.VAR_NAME_MICRO_EDGES_PER_MACRO_EDGE} = (int64_t) levelinfo::num_microedges_per_edge( level );\n"
+ f"const auto {self.VAR_NAME_MICRO_FACES_PER_MACRO_FACE} = (int64_t) levelinfo::num_microfaces_per_face( level );\n"
f"const auto {self.VAR_NAME_MICRO_EDGES_PER_MACRO_EDGE_FLOAT} = ({self._type_descriptor.pystencils_type}) levelinfo::num_microedges_per_edge( level );\n"
)
scalar_parameter_setup += "\n".join(