From 236cd1c6715ed55039cb92e7d90d4d5e0c69d961 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Wed, 13 Mar 2024 16:16:54 +0100
Subject: [PATCH 01/33] Cherry-picked be7e45e5bb08baf12bb18edfafbfc6052062c5c0
from old version: Draft for an updated version of the TensorialVectorSpace as
a first step towards vector-vector operators.
---
hog/function_space.py | 94 ++++++++++++++++++++-----------
hog/math_helpers.py | 10 +++-
hog_tests/test_function_spaces.py | 66 ++++++++++++++++++++++
3 files changed, 133 insertions(+), 37 deletions(-)
create mode 100644 hog_tests/test_function_spaces.py
diff --git a/hog/function_space.py b/hog/function_space.py
index 08fe799..aa70ef3 100644
--- a/hog/function_space.py
+++ b/hog/function_space.py
@@ -57,7 +57,7 @@ class FunctionSpace(Protocol):
geometry: ElementGeometry,
domain: str = "reference",
dof_map: Optional[List[int]] = None,
- ) -> List[sp.Expr]:
+ ) -> List[sp.MatrixBase]:
"""The basis functions of this FEM space."""
...
@@ -69,6 +69,17 @@ class FunctionSpace(Protocol):
) -> List[sp.MatrixBase]:
"""Returns a list containing the gradients of the shape functions on the element.
+ Particularly, for each (vector- or scalar-valued) shape function N = (N_1, ..., N_n) returns the _transposed_
+ Jacobian
+
+ ⎡∂Nâ‚/∂x₠··· ∂Nâ‚™/∂xâ‚⎤
+ ⎢ · · ⎥
+ grad(N) = J(N)ᵀ = ⎢ · · ⎥
+ ⎢ · · ⎥
+ ⎣∂Nâ‚/∂xâ‚– ··· ∂Nâ‚™/∂xâ‚– ⎦
+
+ i.e., the returned gradient is a column vector for scalar shape functions.
+
:param dof_map: this list can be used to specify (remap) the DoF ordering of the element
"""
if domain in ["ref", "reference"]:
@@ -131,7 +142,7 @@ class LagrangianFunctionSpace(FunctionSpace):
geometry: ElementGeometry,
domain: str = "reference",
dof_map: Optional[List[int]] = None,
- ) -> List[sp.Expr]:
+ ) -> List[sp.MatrixBase]:
"""Returns a list containing the shape functions on the element.
:param dof_map: this list can be used to specify (remap) the DoF ordering of the element
@@ -249,6 +260,7 @@ class LagrangianFunctionSpace(FunctionSpace):
if dof_map:
raise HOGException("DoF reordering not implemented.")
+ basis_functions = [sp.Matrix([b]) for b in basis_functions]
return basis_functions
raise HOGException(f"Shape function not available for domain type {domain}")
@@ -266,9 +278,35 @@ class LagrangianFunctionSpace(FunctionSpace):
class TensorialVectorFunctionSpace(FunctionSpace):
- def __init__(self, component: int, function_space: FunctionSpace):
+ """
+ Given a scalar function space, this class can be used to construct a tensorial vector-function space.
+
+ Each shape function of the resulting tensorial space has only one non-zero component that is set to one of the
+ shape functions of the given scalar function space.
+
+ For instance, a tensorial function space of dim d = 3 that is based on a P1 (Lagrangian) scalar function space with
+ n shape functions N_1, ..., N_n has the shape functions
+
+ (N_1, 0, 0),
+ ...
+ (N_n, 0, 0),
+ (0, N_1, 0),
+ ...
+ (0, N_n, 0),
+ (0, 0, N_1),
+ ...
+ (0, 0, N_n).
+ """
+
+ def __init__(self, function_space: FunctionSpace, dim: int = None):
+ """
+ Initializes a tensorial vector function space from a scalar function space.
+
+ :param function_space: the (scalar) component function space
+ :param dim: optional parameter that is set to the dimensionality of the element later if left None here
+ """
self._component_function_space = function_space
- self._component = component
+ self._dim = None
@property
def is_vectorial(self) -> bool:
@@ -286,24 +324,22 @@ class TensorialVectorFunctionSpace(FunctionSpace):
def family(self) -> str:
return self._component_function_space.family
- @property
- def component(self) -> int:
- return self._component
-
@property
def is_continuous(self) -> bool:
return self._component_function_space.is_continuous
- def _to_vector(self, phi: sp.Expr, geometry: ElementGeometry) -> sp.MatrixBase:
- dimensions = geometry.dimensions
+ def _to_vector(
+ self, phi: sp.MatrixBase, component: int, dimensions: int
+ ) -> sp.MatrixBase:
+ if phi.shape != (1, 1):
+ raise HOGException("Component of tensorial space must be scalar.")
return sp.Matrix(
- [
- [phi if c == self._component else sp.sympify(0)]
- for c in range(dimensions)
- ]
+ [[phi if c == component else sp.sympify(0)] for c in range(dimensions)]
)
- def _to_matrix(self, grad_phi: sp.Expr, geometry: ElementGeometry) -> sp.MatrixBase:
+ def _to_matrix(
+ self, grad_phi: sp.MatrixBase, geometry: ElementGeometry
+ ) -> sp.MatrixBase:
grad_phi = grad_phi.transpose().tolist()[0]
dimensions = geometry.dimensions
zero = [sp.sympify(0) for i in range(dimensions)]
@@ -316,40 +352,30 @@ class TensorialVectorFunctionSpace(FunctionSpace):
geometry: ElementGeometry,
domain: str = "reference",
dof_map: Optional[List[int]] = None,
- ) -> List[sp.Expr]:
- shape_functions = self._component_function_space.shape(
- geometry, domain, dof_map
- )
- return [self._to_vector(phi, geometry) for phi in shape_functions]
-
- def grad_shape(
- self,
- geometry: ElementGeometry,
- domain: str = "reference",
- dof_map: Optional[List[int]] = None,
) -> List[sp.MatrixBase]:
- grad_shape_functions = self._component_function_space.grad_shape(
+ # Defaulting to the dimensionality of the geometry if not specified otherwise.
+ dim = geometry.dimensions
+ if self._dim:
+ dim = self._dim
+
+ shape_functions = self._component_function_space.shape(
geometry, domain, dof_map
)
return [
- self._to_matrix(grad_shape, geometry) for grad_shape in grad_shape_functions
+ self._to_vector(phi, c, dim) for c in range(dim) for phi in shape_functions
]
- def num_dofs(self, geometry: ElementGeometry) -> int:
- """Returns the number of DoFs per element."""
- return len(self.shape(geometry))
-
def __eq__(self, other: Any) -> bool:
if type(self) != type(other):
return False
if not hasattr(other, "_component_function_space"):
return False
- if self._component != other._component:
+ if self._dim != other._dim:
return False
return self._component_function_space == other._component_function_space
def __str__(self):
- return f"Vectorial({self._component_function_space})"
+ return f"TensorialVectorSpace({self._component_function_space})"
def __repr__(self):
return str(self)
diff --git a/hog/math_helpers.py b/hog/math_helpers.py
index 6274055..eca3b78 100644
--- a/hog/math_helpers.py
+++ b/hog/math_helpers.py
@@ -14,7 +14,7 @@
# 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 typing import List
+from typing import List, Union
import sympy as sp
from hog.exception import HOGException
@@ -46,9 +46,13 @@ def abs(expr: sp.Expr) -> sp.Expr:
return sp.Abs(expr)
-def grad(f: sp.Expr, symbols: List[sp.Symbol]) -> sp.MatrixBase:
+def grad(f: Union[sp.Expr, sp.MatrixBase], symbols: List[sp.Symbol]) -> sp.MatrixBase:
"""Returns the gradient of the passed sympy expression with respect to the passed symbols."""
- return sp.simplify(sp.Matrix([[sp.diff(f, s)] for s in symbols]))
+ if isinstance(f, sp.MatrixBase):
+ return sp.simplify(f.jacobian(symbols).T)
+ elif isinstance(f, sp.Expr):
+ return sp.simplify(sp.Matrix([[sp.diff(f, s)] for s in symbols]))
+ raise HOGException("Invalid data type in grad().")
def curl(u: sp.Matrix, symbols: List[sp.Symbol]) -> sp.Expr:
diff --git a/hog_tests/test_function_spaces.py b/hog_tests/test_function_spaces.py
new file mode 100644
index 0000000..1ffa604
--- /dev/null
+++ b/hog_tests/test_function_spaces.py
@@ -0,0 +1,66 @@
+# 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 sympy as sp
+from hog.element_geometry import TriangleElement
+from hog.function_space import LagrangianFunctionSpace, TensorialVectorFunctionSpace
+from hog.symbolizer import Symbolizer
+from hog.exception import HOGException
+
+
+def test_function_spaces():
+ symbolizer = Symbolizer()
+ geometry = TriangleElement()
+
+ print()
+
+ f = LagrangianFunctionSpace(1, symbolizer)
+ f_shape = f.shape(geometry)
+ f_grad_shape = f.grad_shape(geometry)
+ print(f)
+ print()
+ print("- shape")
+ for s in f_shape:
+ if s.shape != (1, 1):
+ raise HOGException("Shape function has wrong shape (sic).")
+ sp.pprint(s)
+ print()
+ print("- grad shape")
+ for g in f_grad_shape:
+ if g.shape != (geometry.dimensions, 1):
+ raise HOGException("Gradient has wrong shape.")
+ sp.pprint(g)
+
+ print()
+ print()
+
+ f_scalar = LagrangianFunctionSpace(1, symbolizer)
+ f = TensorialVectorFunctionSpace(f_scalar)
+ f_shape = f.shape(geometry)
+ f_grad_shape = f.grad_shape(geometry)
+ print(f)
+ print()
+ print("- shape")
+ for s in f_shape:
+ if s.shape != (geometry.dimensions, 1):
+ raise HOGException("Shape function has wrong shape (sic).")
+ sp.pprint(s)
+ print()
+ print("- grad shape")
+ for g in f_grad_shape:
+ if g.shape != (geometry.dimensions, geometry.dimensions):
+ raise HOGException("Gradient has wrong shape.")
+ sp.pprint(g)
--
GitLab
From c72965a51de26ca0e8a6229f573400eae8392471 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Fri, 26 Apr 2024 16:02:47 +0200
Subject: [PATCH 02/33] First apparently working version of tensorial vector
function operators.
Instead of having dedicated forms, the diffusion form has been modified slightly by replacing the dot product with a double contraction. This should ensure compatibility with the scalar version.
* only P1
* only 2D (3D should be straightforward)
* only briefly tested vector Laplace
* some cleanup pending
---
generate_all_operators.py | 12 +-
hog/forms.py | 6 +-
hog/function_space.py | 25 +--
.../function_space_impls.py | 195 +++++++++++++++++-
hog/operator_generation/kernel_types.py | 4 +-
5 files changed, 213 insertions(+), 29 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 47c7c9e..c9c1061 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -43,7 +43,12 @@ from hog.forms import (
nonlinear_diffusion_newton_galerkin,
)
from hog.forms_vectorial import curl_curl, curl_curl_plus_mass, mass_n1e1
-from hog.function_space import FunctionSpace, LagrangianFunctionSpace, N1E1Space
+from hog.function_space import (
+ FunctionSpace,
+ LagrangianFunctionSpace,
+ N1E1Space,
+ TensorialVectorFunctionSpace,
+)
from hog.logger import get_logger, TimedLogger
from hog.operator_generation.kernel_types import (
Apply,
@@ -516,6 +521,7 @@ def all_operators(
blending: GeometryMap,
) -> List[OperatorInfo]:
P1 = LagrangianFunctionSpace(1, symbolizer)
+ P1Vector = TensorialVectorFunctionSpace(P1)
P2 = LagrangianFunctionSpace(2, symbolizer)
N1E1 = N1E1Space(symbolizer)
@@ -556,6 +562,10 @@ def all_operators(
coefficient_function_space=P1, onlyNewtonGalerkinPartOfForm=False),
type_descriptor=type_descriptor, opts=opts, blending=blending))
+ ops.append(OperatorInfo(mapping="P1Vector", name="Diffusion", trial_space=P1Vector, test_space=P1Vector,
+ form=diffusion, type_descriptor=type_descriptor, opts=opts, blending=blending,
+ geometries=[TriangleElement()]))
+
# fmt: on
for c in [0, 1, 2]:
diff --git a/hog/forms.py b/hog/forms.py
index dc98e44..0e2378c 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -79,7 +79,7 @@ Weak formulation
u: trial function (space: {trial})
v: test function (space: {test})
- ∫ ∇u · ∇v
+ ∫ ∇u : ∇v
"""
if trial != test:
@@ -124,7 +124,7 @@ Weak formulation
jac_affine_inv.T * grad_psi,
)
form = (
- dot(
+ double_contraction(
jac_blending_inv.T
* sp.Matrix(jac_affine_inv_T_grad_phi_symbols),
jac_blending_inv.T
@@ -137,7 +137,7 @@ Weak formulation
jac_affine_inv_grad_phi_jac_affine_inv_grad_psi_det_symbol = (
tabulation.register_factor(
"jac_affine_inv_grad_phi_jac_affine_inv_grad_psi_det",
- dot(
+ double_contraction(
jac_affine_inv.T * grad_phi,
jac_affine_inv.T * grad_psi,
)
diff --git a/hog/function_space.py b/hog/function_space.py
index aa70ef3..67dfe74 100644
--- a/hog/function_space.py
+++ b/hog/function_space.py
@@ -298,15 +298,13 @@ class TensorialVectorFunctionSpace(FunctionSpace):
(0, 0, N_n).
"""
- def __init__(self, function_space: FunctionSpace, dim: int = None):
+ def __init__(self, function_space: FunctionSpace):
"""
Initializes a tensorial vector function space from a scalar function space.
:param function_space: the (scalar) component function space
- :param dim: optional parameter that is set to the dimensionality of the element later if left None here
"""
self._component_function_space = function_space
- self._dim = None
@property
def is_vectorial(self) -> bool:
@@ -328,6 +326,10 @@ class TensorialVectorFunctionSpace(FunctionSpace):
def is_continuous(self) -> bool:
return self._component_function_space.is_continuous
+ @property
+ def component_function_space(self) -> FunctionSpace:
+ return self._component_function_space
+
def _to_vector(
self, phi: sp.MatrixBase, component: int, dimensions: int
) -> sp.MatrixBase:
@@ -337,30 +339,19 @@ class TensorialVectorFunctionSpace(FunctionSpace):
[[phi if c == component else sp.sympify(0)] for c in range(dimensions)]
)
- def _to_matrix(
- self, grad_phi: sp.MatrixBase, geometry: ElementGeometry
- ) -> sp.MatrixBase:
- grad_phi = grad_phi.transpose().tolist()[0]
- dimensions = geometry.dimensions
- zero = [sp.sympify(0) for i in range(dimensions)]
- return sp.Matrix(
- [(grad_phi if c == self._component else zero) for c in range(dimensions)]
- )
-
def shape(
self,
geometry: ElementGeometry,
domain: str = "reference",
dof_map: Optional[List[int]] = None,
) -> List[sp.MatrixBase]:
- # Defaulting to the dimensionality of the geometry if not specified otherwise.
+
dim = geometry.dimensions
- if self._dim:
- dim = self._dim
shape_functions = self._component_function_space.shape(
geometry, domain, dof_map
)
+
return [
self._to_vector(phi, c, dim) for c in range(dim) for phi in shape_functions
]
@@ -370,8 +361,6 @@ class TensorialVectorFunctionSpace(FunctionSpace):
return False
if not hasattr(other, "_component_function_space"):
return False
- if self._dim != other._dim:
- return False
return self._component_function_space == other._component_function_space
def __str__(self):
diff --git a/hog/operator_generation/function_space_impls.py b/hog/operator_generation/function_space_impls.py
index ac8f8a9..fc92240 100644
--- a/hog/operator_generation/function_space_impls.py
+++ b/hog/operator_generation/function_space_impls.py
@@ -27,7 +27,12 @@ from pystencils.backends.cbackend import CustomCodeNode
from hog.element_geometry import ElementGeometry
from hog.exception import HOGException
-from hog.function_space import FunctionSpace, LagrangianFunctionSpace, N1E1Space
+from hog.function_space import (
+ FunctionSpace,
+ LagrangianFunctionSpace,
+ TensorialVectorFunctionSpace,
+ N1E1Space,
+)
from hog.operator_generation.indexing import (
CellType,
FaceType,
@@ -95,6 +100,18 @@ class FunctionSpaceImpl(ABC):
impl_class = P1FunctionSpaceImpl
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:
+ impl_class = P1VectorFunctionSpaceImpl
+ else:
+ raise HOGException(
+ "TensorialVectorFunctionSpaces are only supported for P1 components."
+ )
+ else:
+ raise HOGException(
+ "TensorialVectorFunctionSpaces are only supported with Lagrangian component spaces."
+ )
elif isinstance(func_space, N1E1Space):
impl_class = N1E1FunctionSpaceImpl
else:
@@ -136,7 +153,7 @@ class FunctionSpaceImpl(ABC):
...
@abstractmethod
- def invert_elementwise(self) -> str:
+ def invert_elementwise(self, dim: int) -> str:
"""C++ code for inverting each DoF of the linalg vector."""
...
@@ -268,7 +285,7 @@ class P1FunctionSpaceImpl(FunctionSpaceImpl):
return f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {macro}.getData( {self._deref()}.get{Macro}DataID() )->getPointer( level );"
- def invert_elementwise(self) -> str:
+ def invert_elementwise(self, dim: int) -> str:
return f"{self._deref()}.invertElementwise( level );"
def local_dofs(
@@ -394,7 +411,7 @@ class P2FunctionSpaceImpl(FunctionSpaceImpl):
f"{self.type_descriptor.pystencils_type}* _data_{self.name}Edge = {macro}.getData( {self._deref()}.getEdgeDoFFunction().get{Macro}DataID() )->getPointer( level );"
)
- def invert_elementwise(self) -> str:
+ def invert_elementwise(self, dim: int) -> str:
return f"{self._deref()}.invertElementwise( level );"
def local_dofs(
@@ -442,6 +459,174 @@ class P2FunctionSpaceImpl(FunctionSpaceImpl):
)
+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 = {}
+
+ 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:
+ 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" {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:
+ 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]
+
+ 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,
+ ) -> 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_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],
+ ) -> Tuple[CustomCodeNode, sp.MatrixBase]:
+ return (
+ CustomCodeNode("", [], []),
+ sp.Identity(self.fe_space.num_dofs(geometry)),
+ )
+
+
class N1E1FunctionSpaceImpl(FunctionSpaceImpl):
def __init__(
self,
@@ -489,7 +674,7 @@ class N1E1FunctionSpaceImpl(FunctionSpaceImpl):
return f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {macro}.getData( {self._deref()}.getDoFs()->get{Macro}DataID() )->getPointer( level );"
- def invert_elementwise(self) -> str:
+ def invert_elementwise(self, dim: int) -> str:
return f"{self._deref()}.getDoFs()->invertElementwise( level );"
def local_dofs(
diff --git a/hog/operator_generation/kernel_types.py b/hog/operator_generation/kernel_types.py
index a474d9c..c2c781a 100644
--- a/hog/operator_generation/kernel_types.py
+++ b/hog/operator_generation/kernel_types.py
@@ -387,7 +387,7 @@ class AssembleDiagonal(KernelType):
f"\n"
f"for ( uint_t level = minLevel_; level <= maxLevel_; level++ )\n"
f"{{\n"
- f" {self.dst.name}->setToZero( level );\n"
+ f" {self.dst.name}->interpolate( 0, level );\n"
f"\n"
f" if ( storage_->hasGlobalCells() )\n"
f" {{\n"
@@ -406,7 +406,7 @@ class AssembleDiagonal(KernelType):
f"{indent(macro_loop(2), 2 * INDENT)}\n"
f" }}\n"
f"\n"
- f"{indent(self.dst.invert_elementwise(), INDENT)}\n"
+ f"{indent(self.dst.invert_elementwise(dim=0), INDENT)}\n"
f"}}\n"
f"\n"
f'this->stopTiming( "{self.name}" );'
--
GitLab
From 1b4059c8e08adf2a6aaf8a02815f6eee3d660f0f Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Tue, 30 Apr 2024 08:54:21 +0200
Subject: [PATCH 03/33] Extended vector ops to 3D (only P1 for now).
---
generate_all_operators.py | 3 +-
.../function_space_impls.py | 36 +++++++++++--------
2 files changed, 22 insertions(+), 17 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index c9c1061..4dfe899 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -563,8 +563,7 @@ def all_operators(
type_descriptor=type_descriptor, opts=opts, blending=blending))
ops.append(OperatorInfo(mapping="P1Vector", name="Diffusion", trial_space=P1Vector, test_space=P1Vector,
- form=diffusion, type_descriptor=type_descriptor, opts=opts, blending=blending,
- geometries=[TriangleElement()]))
+ form=diffusion, type_descriptor=type_descriptor, opts=opts, blending=blending))
# fmt: on
diff --git a/hog/operator_generation/function_space_impls.py b/hog/operator_generation/function_space_impls.py
index fc92240..276b362 100644
--- a/hog/operator_generation/function_space_impls.py
+++ b/hog/operator_generation/function_space_impls.py
@@ -19,7 +19,7 @@ from abc import ABC, abstractmethod
import itertools
import numpy as np
import sympy as sp
-from typing import List, Tuple, Set, Union
+from typing import List, Tuple, Set, Union, Dict
import pystencils as ps
from pystencils import Field, FieldType
@@ -468,7 +468,7 @@ class P1VectorFunctionSpaceImpl(FunctionSpaceImpl):
is_pointer: bool = False,
):
super().__init__(fe_space, name, type_descriptor, is_pointer)
- self.fields = {}
+ self.fields: Dict = {}
def _field_name(self, component: int) -> str:
return self.name + f"_{component}"
@@ -480,11 +480,14 @@ class P1VectorFunctionSpaceImpl(FunctionSpaceImpl):
if dim == 2:
return f"communication::syncVectorFunctionBetweenPrimitives( {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 );"
- )
+ 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:
@@ -516,16 +519,19 @@ class P1VectorFunctionSpaceImpl(FunctionSpaceImpl):
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"
- )
+ f"{self._deref()}[{i}].communicateAdditively < Face, Edge > ( {params} );\n"
+ f"{self._deref()}[{i}].communicateAdditively < Face, Vertex > ( {params} );\n"
+ )
return ret_str
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} );"
- )
+ 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)."""
--
GitLab
From 55ae07cc89c29914f635ba813fd2610b59a847e1 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Mon, 6 May 2024 14:17:24 +0200
Subject: [PATCH 04/33] Extended vector ops to P2.
---
hog/forms.py | 5 +-
.../function_space_impls.py | 187 +++++++++++++++++-
2 files changed, 189 insertions(+), 3 deletions(-)
diff --git a/hog/forms.py b/hog/forms.py
index 0e2378c..3b6c72a 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -80,6 +80,10 @@ Weak formulation
v: test function (space: {test})
∫ ∇u : ∇v
+
+ Note that the double contraction (:) reduces to the dot product for scalar function spaces, i.e. the form becomes
+
+ ∫ ∇u · ∇v
"""
if trial != test:
@@ -90,7 +94,6 @@ Weak formulation
with TimedLogger("assembling diffusion matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
diff --git a/hog/operator_generation/function_space_impls.py b/hog/operator_generation/function_space_impls.py
index 276b362..130539f 100644
--- a/hog/operator_generation/function_space_impls.py
+++ b/hog/operator_generation/function_space_impls.py
@@ -104,9 +104,11 @@ class FunctionSpaceImpl(ABC):
if isinstance(func_space.component_function_space, LagrangianFunctionSpace):
if func_space.component_function_space.degree == 1:
impl_class = P1VectorFunctionSpaceImpl
+ elif func_space.component_function_space.degree == 2:
+ impl_class = P2VectorFunctionSpaceImpl
else:
raise HOGException(
- "TensorialVectorFunctionSpaces are only supported for P1 components."
+ "TensorialVectorFunctionSpaces not supported for the chosen components."
)
else:
raise HOGException(
@@ -468,7 +470,7 @@ class P1VectorFunctionSpaceImpl(FunctionSpaceImpl):
is_pointer: bool = False,
):
super().__init__(fe_space, name, type_descriptor, is_pointer)
- self.fields: Dict = {}
+ self.fields: Dict[int, Field] = {}
def _field_name(self, component: int) -> str:
return self.name + f"_{component}"
@@ -633,6 +635,187 @@ class P1VectorFunctionSpaceImpl(FunctionSpaceImpl):
)
+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,
+ ) -> 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
+ )
+
+ 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],
+ ) -> Tuple[CustomCodeNode, sp.MatrixBase]:
+ return (
+ CustomCodeNode("", [], []),
+ sp.Identity(self.fe_space.num_dofs(geometry)),
+ )
+
+
class N1E1FunctionSpaceImpl(FunctionSpaceImpl):
def __init__(
self,
--
GitLab
From 30200e188d840f3858f874b42070ecc1eaab2b61 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Mon, 13 May 2024 13:47:32 +0200
Subject: [PATCH 05/33] Some debugging output in generate_all_operators.py
(mainly for the CI).
---
generate_all_operators.py | 3 +++
1 file changed, 3 insertions(+)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 4dfe899..d0d76fb 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -15,6 +15,7 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>.
import argparse
+import sys
from argparse import ArgumentParser
from dataclasses import dataclass, field
from functools import partial
@@ -280,6 +281,8 @@ def main():
logger = get_logger(log_level)
TimedLogger.set_log_level(log_level)
+ logger.info(f"Running {__file__} with arguments: {' '.join(sys.argv[1:])}")
+
symbolizer = Symbolizer()
if args.list_quadratures:
--
GitLab
From ac8f91ff301632c36cc0713c7def643289583e1e Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Wed, 15 May 2024 16:02:26 +0200
Subject: [PATCH 06/33] Using re.fullmatch also for the FE spaces in the
generate_all_operators.py script. Otherwise, e.g., 'P1Vector' is also chosen
if only 'P1' is specified.
---
generate_all_operators.py | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index d0d76fb..b9b7d1c 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -357,7 +357,8 @@ def main():
filtered_operators = list(
filter(
- lambda o: re_mapping.search(o.mapping) and re.fullmatch(args.form, o.name),
+ lambda o: re.fullmatch(args.space_mapping, o.mapping)
+ and re.fullmatch(args.form, o.name),
operators,
)
)
--
GitLab
From 98f9c4396d97df9624d372276fcd1d71ec34ce66 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Wed, 15 May 2024 16:03:28 +0200
Subject: [PATCH 07/33] Updating the epsilon form to also work with (tensorial)
vector spaces to test offdiagonal blocks.
---
hog/forms.py | 34 +++++++++++++++++++++++++---------
1 file changed, 25 insertions(+), 9 deletions(-)
diff --git a/hog/forms.py b/hog/forms.py
index 3b6c72a..d9ba81b 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -525,11 +525,22 @@ def epsilon(
variable_viscosity: bool = True,
coefficient_function_space: Optional[FunctionSpace] = None,
) -> Form:
+
+ if trial.is_vectorial ^ test.is_vectorial:
+ raise HOGException(
+ "Either both (trial and test) spaces or none should be vectorial."
+ )
+
+ vectorial_spaces = trial.is_vectorial
+ docstring_components = (
+ ""
+ if vectorial_spaces
+ else f"\nComponent trial: {component_trial}\nComponent test: {component_test}"
+ )
+
docstring = f"""
"Epsilon" operator.
-
-Component trial: {component_trial}
-Component test: {component_test}
+{docstring_components}
Geometry map: {blending}
Weak formulation
@@ -544,12 +555,17 @@ where
ε(w) := (1/2) (∇w + (∇w)ᵀ)
"""
- if variable_viscosity == False:
+ if not variable_viscosity:
raise HOGException("Constant viscosity currently not supported.")
# TODO fix issue with undeclared p_affines
- if geometry.dimensions < 3 and (component_trial > 1 or component_test > 1):
+ if (
+ not vectorial_spaces
+ and geometry.dimensions < 3
+ and (component_trial > 1 or component_test > 1)
+ ):
return create_empty_element_matrix(trial, test, geometry)
+
with TimedLogger("assembling epsilon matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
@@ -602,7 +618,7 @@ where
if isinstance(trial, EnrichedGalerkinFunctionSpace):
# for EDG, the shape function is already vectorial and does not have to be multiplied by e_vec
grad_phi_vec = jac_affine * grad_phi_vec
- else:
+ elif not vectorial_spaces:
grad_phi_vec = (
(e_vec(geometry.dimensions, component_trial) * phi)
.jacobian(ref_symbols_list)
@@ -612,7 +628,7 @@ where
if isinstance(test, EnrichedGalerkinFunctionSpace):
# for EDG, the shape function is already vectorial and does not have to be multiplied by e_vec
grad_psi_vec = jac_affine * grad_psi_vec
- else:
+ elif not vectorial_spaces:
grad_psi_vec = (
(e_vec(geometry.dimensions, component_test) * psi)
.jacobian(ref_symbols_list)
@@ -697,7 +713,7 @@ where
return Form(
mat,
tabulation,
- symmetric=component_trial == component_test,
+ symmetric=vectorial_spaces or (component_trial == component_test),
docstring=docstring,
)
@@ -727,7 +743,7 @@ Weak formulation
if trial != test:
TimedLogger(
"Trial and test space can be different, but please make sure this is intensional!",
- level=logging.INFO
+ level=logging.INFO,
).log()
with TimedLogger("assembling k-mass matrix", level=logging.DEBUG):
--
GitLab
From 01cb0a1ee9d74a8ddcbb07aae9d350b4d67e5cbc Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Wed, 15 May 2024 16:15:39 +0200
Subject: [PATCH 08/33] Attempt to add a P1VectorLaplace HyTeG-integration
test.
---
hyteg_integration_tests/src/CMakeLists.txt | 9 +--
.../src/DiffusionP1Vector.cpp | 61 +++++++++++++++++++
2 files changed, 66 insertions(+), 4 deletions(-)
create mode 100644 hyteg_integration_tests/src/DiffusionP1Vector.cpp
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index 73be793..b84602b 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -44,10 +44,11 @@ add_operator_test(FILE CurlCurl.cpp FORM CurlCurl ABBR CVIP
add_operator_test(FILE CurlCurlPlusMass.cpp DEF TEST_ASSEMBLE FORM CurlCurlPlusMass ABBR noopts GEN_ARGS --quad-degree 2 2)
add_operator_test(FILE CurlCurlPlusMass.cpp FORM CurlCurlPlusMass ABBR VI GEN_ARGS --quad-degree 2 2 -o MOVECONSTANTS VECTORIZE)
-add_operator_test(FILE DiffusionP1.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1 GEN_ARGS -s P1 --quad-degree 0 0)
-add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1VIP GEN_ARGS -s P1 --quad-degree 0 0 -o MOVECONSTANTS VECTORIZE POLYCSE)
-add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1CVI GEN_ARGS -s P1 --quad-degree 0 0 --loop-strategy CUBES -o MOVECONSTANTS VECTORIZE)
-add_operator_test(FILE DiffusionP2.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 2 GEN_ARGS -s P2 --quad-degree 2 2)
+add_operator_test(FILE DiffusionP1.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1 GEN_ARGS -s P1 --quad-degree 0 0)
+add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1VIP GEN_ARGS -s P1 --quad-degree 0 0 -o MOVECONSTANTS VECTORIZE POLYCSE)
+add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1CVI GEN_ARGS -s P1 --quad-degree 0 0 --loop-strategy CUBES -o MOVECONSTANTS VECTORIZE)
+add_operator_test(FILE DiffusionP2.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 2 GEN_ARGS -s P2 --quad-degree 2 2)
+add_operator_test(FILE DiffusionP1Vector.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1Vec GEN_ARGS -s P1Vector --quad-degree 0 0)
add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR noopts GEN_ARGS -s P2 --quad-degree 4 4)
add_operator_test(FILE DivKGrad.cpp FORM DivKGrad ABBR VIQ GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS VECTORIZE QUADLOOPS )
diff --git a/hyteg_integration_tests/src/DiffusionP1Vector.cpp b/hyteg_integration_tests/src/DiffusionP1Vector.cpp
new file mode 100644
index 0000000..5c3424b
--- /dev/null
+++ b/hyteg_integration_tests/src/DiffusionP1Vector.cpp
@@ -0,0 +1,61 @@
+/*
+ * HyTeG Operator Generator
+ * Copyright (C) 2017-2024 Nils Kohl, Fabian Böhm, Daniel Bauer
+ *
+ * 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/>.
+ */
+
+#include "core/DataTypes.h"
+
+#include "hyteg/elementwiseoperators/P1ElementwiseOperator.hpp"
+#include "mixed_operator/VectorLaplaceOperator.hpp"
+
+#include "Diffusion/TestOpDiffusion.hpp"
+#include "OperatorGenerationTest.hpp"
+
+using namespace hyteg;
+using walberla::real_t;
+
+int main( int argc, char* argv[] )
+{
+ walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+ walberla::MPIManager::instance()->useWorldComm();
+
+ const uint_t level = 3;
+
+ for ( uint_t d = 2; d <= 3; ++d )
+ {
+ StorageSetup storageSetup;
+ if ( d == 2 )
+ {
+ storageSetup = StorageSetup(
+ "quad_4el", MeshInfo::fromGmshFile( "../hyteg/data/meshes/quad_4el.msh" ), GeometryMap::Type::IDENTITY );
+ }
+ else
+ {
+ storageSetup = StorageSetup(
+ "cube_6el", MeshInfo::fromGmshFile( "../hyteg/data/meshes/3D/cube_6el.msh" ), GeometryMap::Type::IDENTITY );
+ }
+
+ testOperators< P1VectorFunction< real_t >, P1ElementwiseVectorLaplaceOperator, operatorgeneration::TestOpDiffusion >(
+ level, storageSetup, 225, 9.0e6 );
+
+#ifdef TEST_ASSEMBLE
+ compareAssembledMatrix< P1ElementwiseLaplaceOperator, operatorgeneration::TestOpDiffusion >(
+ level, storageSetup, storageSetup.description() + " Assembly", 360 );
+#endif
+ }
+
+ return EXIT_SUCCESS;
+}
--
GitLab
From 05b9f9384b03a8873d487b6e3a312ff279a24e8f Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 09:51:36 +0200
Subject: [PATCH 09/33] Fixed regex default argument for --space-mapping
option.
---
generate_all_operators.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index b9b7d1c..32d9671 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -120,7 +120,7 @@ def parse_arguments():
parser.add_argument(
"-s",
"--space-mapping",
- default="",
+ default=".*?",
help="Filter by function space, e.g. 'P1' or 'P1toP2' (regex).",
# TODO add choices list and type=string.lower
)
--
GitLab
From c24f1a3b545d227fb7e0442b6b2c37a41b4b2261 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 11:53:31 +0200
Subject: [PATCH 10/33] Ignoring case in re.fullmatch.
---
generate_all_operators.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 32d9671..ee34f56 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -357,8 +357,8 @@ def main():
filtered_operators = list(
filter(
- lambda o: re.fullmatch(args.space_mapping, o.mapping)
- and re.fullmatch(args.form, o.name),
+ lambda o: re.fullmatch(args.space_mapping, o.mapping, re.IGNORECASE)
+ and re.fullmatch(args.form, o.name, re.IGNORECASE),
operators,
)
)
--
GitLab
From ad2923f5c83775c10a26792310cf954628c98750 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 11:54:10 +0200
Subject: [PATCH 11/33] HyTeG integration tests: Fixing space mappings. Added
P2Vector epsilon test.
---
hyteg_integration_tests/src/CMakeLists.txt | 10 +-
hyteg_integration_tests/src/EpsilonVector.cpp | 111 ++++++++++++++++++
2 files changed, 117 insertions(+), 4 deletions(-)
create mode 100644 hyteg_integration_tests/src/EpsilonVector.cpp
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index b84602b..40a39cc 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -54,10 +54,12 @@ add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR noopts
add_operator_test(FILE DivKGrad.cpp FORM DivKGrad ABBR VIQ GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS VECTORIZE QUADLOOPS )
add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR IQT GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS QUADLOOPS TABULATE )
-add_operator_test(FILE Div.cpp DEF REF_OP=P2ToP1ElementwiseDivxOperator FORM Div_0 ABBR VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
-add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM DivT_2 ABBR VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
-add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
-add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
+add_operator_test(FILE Div.cpp DEF REF_OP=P2ToP1ElementwiseDivxOperator FORM Div_0 ABBR VIQT GEN_ARGS -s P2ToP1 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
+add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM DivT_2 ABBR VIQT GEN_ARGS -s P1ToP2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
+
+add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
+add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
+add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
# tests with blending
diff --git a/hyteg_integration_tests/src/EpsilonVector.cpp b/hyteg_integration_tests/src/EpsilonVector.cpp
new file mode 100644
index 0000000..394e9d7
--- /dev/null
+++ b/hyteg_integration_tests/src/EpsilonVector.cpp
@@ -0,0 +1,111 @@
+/*
+ * HyTeG Operator Generator
+ * Copyright (C) 2017-2024 Nils Kohl, Fabian Böhm, Daniel Bauer
+ *
+ * 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/>.
+ */
+
+#include <type_traits>
+
+#include "core/DataTypes.h"
+
+#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
+#include "hyteg/forms/form_hyteg_generated/p2/p2_epsilonvar_affine_q4.hpp"
+#include "hyteg/p2functionspace/P2Function.hpp"
+
+#include "Epsilon/TestOpEpsilon.hpp"
+#include "OperatorGenerationTest.hpp"
+
+using namespace hyteg;
+using walberla::real_t;
+
+constexpr bool only3d = std::is_constructible< FORM, std::function< real_t( const Point3D& ) > >::value;
+
+real_t k( const hyteg::Point3D& x )
+{
+ return x[0] * x[0] - x[2] * x[2] - x[1] * x[1] + 3 * x[0] + x[1] - 5 * x[2] + 2;
+}
+
+P2ElementwiseAffineEpsilonOperator makeRefOp( std::shared_ptr< PrimitiveStorage > storage, uint_t minLevel, uint_t maxLevel )
+{
+ return P2ElementwiseAffineEpsilonOperator( storage, minLevel, maxLevel, k );
+};
+
+template < class Op >
+Op makeTestOp( std::shared_ptr< PrimitiveStorage > storage, uint_t minLevel, uint_t maxLevel )
+{
+ P2Function< real_t > k( "k", storage, minLevel, maxLevel );
+ for ( size_t lvl = minLevel; lvl <= maxLevel; ++lvl )
+ {
+ k.interpolate( ::k, lvl );
+ }
+ return Op( storage, minLevel, maxLevel, k );
+};
+
+int main( int argc, char* argv[] )
+{
+ walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+ walberla::MPIManager::instance()->useWorldComm();
+
+ const uint_t level = 3;
+
+ real_t thresholdOverMachineEpsApply = 225;
+ real_t thresholdOverMachineEpsInvDiag = 9.0e6;
+ real_t thresholdOverMachineEpsAssembly = 360;
+
+ for ( uint_t d = only3d ? 3 : 2; d <= 3; ++d )
+ {
+ StorageSetup storageSetup;
+ if ( d == 2 )
+ {
+ storageSetup = StorageSetup(
+ "quad_4el", MeshInfo::fromGmshFile( "../hyteg/data/meshes/quad_4el.msh" ), GeometryMap::Type::IDENTITY );
+ }
+ else
+ {
+ storageSetup = StorageSetup(
+ "cube_6el", MeshInfo::fromGmshFile( "../hyteg/data/meshes/3D/cube_6el.msh" ), GeometryMap::Type::IDENTITY );
+ }
+
+ compareApply< P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Apply",
+ thresholdOverMachineEpsApply );
+
+#ifdef TEST_DIAG
+ compareInvDiag< P2Function< real_t >, P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Inverse Diagonal",
+ thresholdOverMachineEpsInvDiag );
+#endif
+
+#ifdef TEST_ASSEMBLE
+ compareAssembledMatrix< P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Assembly",
+ thresholdOverMachineEpsAssembly );
+#endif
+ }
+
+ return EXIT_SUCCESS;
+}
--
GitLab
From bafe845030e8d8a6d9c0ce936665fee158db210b Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 12:11:28 +0200
Subject: [PATCH 12/33] HyTeG integration tests: Fixing DiffusionP1Vector.cpp.
---
hyteg_integration_tests/src/DiffusionP1Vector.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/hyteg_integration_tests/src/DiffusionP1Vector.cpp b/hyteg_integration_tests/src/DiffusionP1Vector.cpp
index 5c3424b..2cbd2c9 100644
--- a/hyteg_integration_tests/src/DiffusionP1Vector.cpp
+++ b/hyteg_integration_tests/src/DiffusionP1Vector.cpp
@@ -52,7 +52,7 @@ int main( int argc, char* argv[] )
level, storageSetup, 225, 9.0e6 );
#ifdef TEST_ASSEMBLE
- compareAssembledMatrix< P1ElementwiseLaplaceOperator, operatorgeneration::TestOpDiffusion >(
+ compareAssembledMatrix< P1ElementwiseVectorLaplaceOperator, operatorgeneration::TestOpDiffusion >(
level, storageSetup, storageSetup.description() + " Assembly", 360 );
#endif
}
--
GitLab
From 175ab6ea7edac101876b153ad1c966785899d5bf Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 12:58:13 +0200
Subject: [PATCH 13/33] HyTeG integration tests: Fixing max magnitude check for
vector functions.
---
.../src/OperatorGenerationTest.hpp | 30 +++++++++++++------
1 file changed, 21 insertions(+), 9 deletions(-)
diff --git a/hyteg_integration_tests/src/OperatorGenerationTest.hpp b/hyteg_integration_tests/src/OperatorGenerationTest.hpp
index 8e4e6ca..abae369 100644
--- a/hyteg_integration_tests/src/OperatorGenerationTest.hpp
+++ b/hyteg_integration_tests/src/OperatorGenerationTest.hpp
@@ -218,13 +218,17 @@ void compareApply( OperatorFactory< RefOpType > refOpFactory,
}
real_t maxAbs = 1.0;
- if constexpr ( !std::is_same< RefDstFncType, n1e1::N1E1VectorFunction< RefType > >::value )
+ if constexpr ( std::is_same< RefDstFncType, n1e1::N1E1VectorFunction< RefType > >::value )
{
- maxAbs = err.getMaxMagnitude( level );
+ maxAbs = err.getDoFs()->getMaxMagnitude( level );
+ }
+ else if constexpr ( std::is_same_v< RefDstFncType, P1VectorFunction< RefType > > || std::is_same_v< RefDstFncType, P2VectorFunction< RefType > > )
+ {
+ maxAbs = err.getMaxComponentMagnitude( level, All );
}
else
{
- maxAbs = err.getDoFs()->getMaxMagnitude( level );
+ maxAbs = err.getMaxMagnitude( level );
}
const real_t threshold = thresholdOverMachineEps * precisionDict::epsilon< TestType >();
@@ -300,13 +304,17 @@ void compareGEMV( const uint_t level,
// compare
err.assign( { 1.0, -1.0 }, { dst, refDst }, level, All );
real_t maxAbs = 1.0;
- if constexpr ( !std::is_same< DstFncType, n1e1::N1E1VectorFunction< real_t > >::value )
+ if constexpr ( std::is_same< DstFncType, n1e1::N1E1VectorFunction< real_t > >::value )
{
- maxAbs = err.getMaxMagnitude( level );
+ maxAbs = err.getDoFs()->getMaxMagnitude( level );
+ }
+ else if constexpr ( std::is_same_v< RefDstFncType, P1VectorFunction< real_t > > || std::is_same_v< RefDstFncType, P2VectorFunction< real_t > > )
+ {
+ maxAbs = err.getMaxComponentMagnitude( level, All );
}
else
{
- maxAbs = err.getDoFs()->getMaxMagnitude( level );
+ maxAbs = err.getMaxMagnitude( level );
}
const real_t threshold = thresholdOverMachineEps * precisionDict::epsilon< TestType >();
@@ -381,13 +389,17 @@ void compareInvDiag( OperatorFactory< RefOpType > refOpFactory,
}
real_t maxAbs = 1.0;
- if constexpr ( !std::is_same< RefDiagType, n1e1::N1E1VectorFunction< RefType > >::value )
+ if constexpr ( std::is_same< RefDiagType, n1e1::N1E1VectorFunction< RefType > >::value )
{
- maxAbs = err.getMaxMagnitude( level );
+ maxAbs = err.getDoFs()->getMaxMagnitude( level );
+ }
+ else if constexpr ( std::is_same_v< RefDiagType, P1VectorFunction< RefType > > || std::is_same_v< RefDiagType, P2VectorFunction< RefType > > )
+ {
+ maxAbs = err.getMaxComponentMagnitude( level, All );
}
else
{
- maxAbs = err.getDoFs()->getMaxMagnitude( level );
+ maxAbs = err.getMaxMagnitude( level );
}
const real_t threshold = thresholdOverMachineEps * precisionDict::epsilon< TestType >();
--
GitLab
From 752a12f651cf2f155fefc1a80be7e6d4eba2f02d Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 13:22:00 +0200
Subject: [PATCH 14/33] HyTeG integration tests: Minor test fix.
---
hyteg_integration_tests/src/OperatorGenerationTest.hpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/hyteg_integration_tests/src/OperatorGenerationTest.hpp b/hyteg_integration_tests/src/OperatorGenerationTest.hpp
index abae369..3eaec8c 100644
--- a/hyteg_integration_tests/src/OperatorGenerationTest.hpp
+++ b/hyteg_integration_tests/src/OperatorGenerationTest.hpp
@@ -308,7 +308,7 @@ void compareGEMV( const uint_t level,
{
maxAbs = err.getDoFs()->getMaxMagnitude( level );
}
- else if constexpr ( std::is_same_v< RefDstFncType, P1VectorFunction< real_t > > || std::is_same_v< RefDstFncType, P2VectorFunction< real_t > > )
+ else if constexpr ( std::is_same_v< DstFncType, P1VectorFunction< real_t > > || std::is_same_v< DstFncType, P2VectorFunction< real_t > > )
{
maxAbs = err.getMaxComponentMagnitude( level, All );
}
--
GitLab
From 0d385ac28db9f74cbb618adf0b68a7eec85f9a96 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 13:50:28 +0200
Subject: [PATCH 15/33] HyTeG integration tests: Trying to fix test builds by
linking mixed_operator lib.
---
hyteg_integration_tests/src/CMakeLists.txt | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index 40a39cc..8ab65f8 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -48,7 +48,7 @@ add_operator_test(FILE DiffusionP1.cpp DEF TEST_ASSEMBLE FORM Diffusion AB
add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1VIP GEN_ARGS -s P1 --quad-degree 0 0 -o MOVECONSTANTS VECTORIZE POLYCSE)
add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1CVI GEN_ARGS -s P1 --quad-degree 0 0 --loop-strategy CUBES -o MOVECONSTANTS VECTORIZE)
add_operator_test(FILE DiffusionP2.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 2 GEN_ARGS -s P2 --quad-degree 2 2)
-add_operator_test(FILE DiffusionP1Vector.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1Vec GEN_ARGS -s P1Vector --quad-degree 0 0)
+add_operator_test(FILE DiffusionP1Vector.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1Vec GEN_ARGS -s P1Vector --quad-degree 0 0 ARG_LIBS mixed_operator)
add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR noopts GEN_ARGS -s P2 --quad-degree 4 4)
add_operator_test(FILE DivKGrad.cpp FORM DivKGrad ABBR VIQ GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS VECTORIZE QUADLOOPS )
--
GitLab
From a7e0df987e6958cf162f00f35d01dc90d1d16d90 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 13:54:10 +0200
Subject: [PATCH 16/33] HyTeG integration tests: Trying to fix test builds by
linking mixed_operator lib part II
---
hyteg_integration_tests/src/CMakeLists.txt | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index 8ab65f8..fb37499 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -48,7 +48,7 @@ add_operator_test(FILE DiffusionP1.cpp DEF TEST_ASSEMBLE FORM Diffusion AB
add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1VIP GEN_ARGS -s P1 --quad-degree 0 0 -o MOVECONSTANTS VECTORIZE POLYCSE)
add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1CVI GEN_ARGS -s P1 --quad-degree 0 0 --loop-strategy CUBES -o MOVECONSTANTS VECTORIZE)
add_operator_test(FILE DiffusionP2.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 2 GEN_ARGS -s P2 --quad-degree 2 2)
-add_operator_test(FILE DiffusionP1Vector.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1Vec GEN_ARGS -s P1Vector --quad-degree 0 0 ARG_LIBS mixed_operator)
+add_operator_test(FILE DiffusionP1Vector.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1Vec GEN_ARGS -s P1Vector --quad-degree 0 0 LIBS mixed_operator)
add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR noopts GEN_ARGS -s P2 --quad-degree 4 4)
add_operator_test(FILE DivKGrad.cpp FORM DivKGrad ABBR VIQ GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS VECTORIZE QUADLOOPS )
--
GitLab
From a270737fe11ca32d825341fe99d6677aa8168ada Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 16 May 2024 17:11:44 +0200
Subject: [PATCH 17/33] Adding P2Vector epsilon op to available operators.
---
generate_all_operators.py | 19 +++++++++++++++++++
1 file changed, 19 insertions(+)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index ee34f56..280a298 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -527,6 +527,7 @@ def all_operators(
P1 = LagrangianFunctionSpace(1, symbolizer)
P1Vector = TensorialVectorFunctionSpace(P1)
P2 = LagrangianFunctionSpace(2, symbolizer)
+ P2Vector = TensorialVectorFunctionSpace(P2)
N1E1 = N1E1Space(symbolizer)
two_d = [TriangleElement()]
@@ -571,6 +572,24 @@ def all_operators(
# fmt: on
+ p2vec_epsilon = partial(
+ epsilon,
+ variable_viscosity=True,
+ coefficient_function_space=P2,
+ )
+ ops.append(
+ OperatorInfo(
+ mapping=f"P2Vector",
+ name=f"Epsilon",
+ trial_space=P2Vector,
+ test_space=P2Vector,
+ form=p2vec_epsilon,
+ type_descriptor=type_descriptor,
+ opts=opts,
+ blending=blending,
+ )
+ )
+
for c in [0, 1, 2]:
# fmt: off
div_geometries = [TriangleElement(), TetrahedronElement()]
--
GitLab
From 7126768dbbd310b13a63957192d6e11ca0650e2c Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Mon, 27 May 2024 14:13:50 +0200
Subject: [PATCH 18/33] Adding option to only generate 2D or 3D part in script.
Then removed VectorEpsilon test in 3D since it generates very long in CI.
---
generate_all_operators.py | 19 ++++++
hyteg_integration_tests/src/CMakeLists.txt | 6 +-
hyteg_integration_tests/src/EpsilonVector.cpp | 66 ++++++++-----------
3 files changed, 50 insertions(+), 41 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 280a298..83e1791 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -198,6 +198,14 @@ def parse_arguments():
help=f"Specify the name of the generated C++ class.",
)
+ parser.add_argument(
+ "--dimensions",
+ type=int,
+ nargs="+",
+ default=[2, 3],
+ help=f"Domain dimensions for which operators shall be generated.",
+ )
+
logging_group = parser.add_mutually_exclusive_group(required=False)
logging_group.add_argument(
"-v",
@@ -363,6 +371,17 @@ def main():
)
)
+ enabled_geometries = set()
+ if 2 in args.dimensions:
+ enabled_geometries.add(TriangleElement())
+ if 3 in args.dimensions:
+ enabled_geometries.add(TetrahedronElement())
+
+ for f in filtered_operators:
+ f.geometries = list(set(f.geometries) & enabled_geometries)
+
+ filtered_operators = [f for f in filtered_operators if f.geometries]
+
if len(filtered_operators) == 0:
raise HOGException(
f"Form {args.form} does not match any available forms. Run --list for a list."
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index fb37499..76a9e24 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -57,9 +57,9 @@ add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR IQT
add_operator_test(FILE Div.cpp DEF REF_OP=P2ToP1ElementwiseDivxOperator FORM Div_0 ABBR VIQT GEN_ARGS -s P2ToP1 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM DivT_2 ABBR VIQT GEN_ARGS -s P1ToP2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
-add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
-add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
-add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
+add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 )
+add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 )
+add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 --dimensions 2)
# tests with blending
diff --git a/hyteg_integration_tests/src/EpsilonVector.cpp b/hyteg_integration_tests/src/EpsilonVector.cpp
index 394e9d7..a6c99bd 100644
--- a/hyteg_integration_tests/src/EpsilonVector.cpp
+++ b/hyteg_integration_tests/src/EpsilonVector.cpp
@@ -30,8 +30,6 @@
using namespace hyteg;
using walberla::real_t;
-constexpr bool only3d = std::is_constructible< FORM, std::function< real_t( const Point3D& ) > >::value;
-
real_t k( const hyteg::Point3D& x )
{
return x[0] * x[0] - x[2] * x[2] - x[1] * x[1] + 3 * x[0] + x[1] - 5 * x[2] + 2;
@@ -64,48 +62,40 @@ int main( int argc, char* argv[] )
real_t thresholdOverMachineEpsInvDiag = 9.0e6;
real_t thresholdOverMachineEpsAssembly = 360;
- for ( uint_t d = only3d ? 3 : 2; d <= 3; ++d )
- {
- StorageSetup storageSetup;
- if ( d == 2 )
- {
- storageSetup = StorageSetup(
- "quad_4el", MeshInfo::fromGmshFile( "../hyteg/data/meshes/quad_4el.msh" ), GeometryMap::Type::IDENTITY );
- }
- else
- {
- storageSetup = StorageSetup(
- "cube_6el", MeshInfo::fromGmshFile( "../hyteg/data/meshes/3D/cube_6el.msh" ), GeometryMap::Type::IDENTITY );
- }
-
- compareApply< P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
- makeRefOp,
- makeTestOp< operatorgeneration::TestOpEpsilon >,
- level,
- storageSetup,
- storageSetup.description() + " Apply",
- thresholdOverMachineEpsApply );
+ // Testing only 2D since generation of 3D operator just takes too long.
+
+ StorageSetup storageSetup;
+ storageSetup = StorageSetup(
+ "quad_4el", MeshInfo::fromGmshFile( "../hyteg/data/meshes/quad_4el.msh" ), GeometryMap::Type::IDENTITY );
+ }
+
+ compareApply< P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Apply",
+ thresholdOverMachineEpsApply );
#ifdef TEST_DIAG
- compareInvDiag< P2Function< real_t >, P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
- makeRefOp,
- makeTestOp< operatorgeneration::TestOpEpsilon >,
- level,
- storageSetup,
- storageSetup.description() + " Inverse Diagonal",
- thresholdOverMachineEpsInvDiag );
+ compareInvDiag< P2Function< real_t >, P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Inverse Diagonal",
+ thresholdOverMachineEpsInvDiag );
#endif
#ifdef TEST_ASSEMBLE
- compareAssembledMatrix< P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
- makeRefOp,
- makeTestOp< operatorgeneration::TestOpEpsilon >,
- level,
- storageSetup,
- storageSetup.description() + " Assembly",
- thresholdOverMachineEpsAssembly );
+ compareAssembledMatrix< P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Assembly",
+ thresholdOverMachineEpsAssembly );
#endif
- }
return EXIT_SUCCESS;
}
--
GitLab
From e1d588f155469a0b66ed2292af186aca0aab634b Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Mon, 27 May 2024 14:22:38 +0200
Subject: [PATCH 19/33] Trying to satisfy the mypy type checker for the
dimensions parameter in the operator generation script.
---
generate_all_operators.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 83e1791..ae1c2f2 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -371,7 +371,7 @@ def main():
)
)
- enabled_geometries = set()
+ enabled_geometries: Set[TriangleElement | TetrahedronElement] = set()
if 2 in args.dimensions:
enabled_geometries.add(TriangleElement())
if 3 in args.dimensions:
--
GitLab
From 9d6bf5c0148e05a62d595317de98d40e4870ef73 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Mon, 27 May 2024 15:19:32 +0200
Subject: [PATCH 20/33] Fixing epsilon vector test.
---
hyteg_integration_tests/src/CMakeLists.txt | 2 +-
hyteg_integration_tests/src/EpsilonVector.cpp | 2 +-
2 files changed, 2 insertions(+), 2 deletions(-)
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index 76a9e24..c84aea7 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -59,7 +59,7 @@ add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM D
add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 )
add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 )
-add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 --dimensions 2)
+add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 --dimensions 2 LIBS constant_stencil_operator)
# tests with blending
diff --git a/hyteg_integration_tests/src/EpsilonVector.cpp b/hyteg_integration_tests/src/EpsilonVector.cpp
index a6c99bd..581de51 100644
--- a/hyteg_integration_tests/src/EpsilonVector.cpp
+++ b/hyteg_integration_tests/src/EpsilonVector.cpp
@@ -21,6 +21,7 @@
#include "core/DataTypes.h"
#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
+#include "hyteg/src/constant_stencil_operator/P2ConstantEpsilonOperator.hpp"
#include "hyteg/forms/form_hyteg_generated/p2/p2_epsilonvar_affine_q4.hpp"
#include "hyteg/p2functionspace/P2Function.hpp"
@@ -67,7 +68,6 @@ int main( int argc, char* argv[] )
StorageSetup storageSetup;
storageSetup = StorageSetup(
"quad_4el", MeshInfo::fromGmshFile( "../hyteg/data/meshes/quad_4el.msh" ), GeometryMap::Type::IDENTITY );
- }
compareApply< P2ElementwiseAffineEpsilonOperator, operatorgeneration::TestOpEpsilon >(
makeRefOp,
--
GitLab
From 53b7711f646494098caf2b3cd90de08b46235c46 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Mon, 27 May 2024 17:24:48 +0200
Subject: [PATCH 21/33] Fixing epsilon vector test pt. II (wrong include path).
---
hyteg_integration_tests/src/EpsilonVector.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/hyteg_integration_tests/src/EpsilonVector.cpp b/hyteg_integration_tests/src/EpsilonVector.cpp
index 581de51..5ad41ed 100644
--- a/hyteg_integration_tests/src/EpsilonVector.cpp
+++ b/hyteg_integration_tests/src/EpsilonVector.cpp
@@ -21,7 +21,7 @@
#include "core/DataTypes.h"
#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
-#include "hyteg/src/constant_stencil_operator/P2ConstantEpsilonOperator.hpp"
+#include "constant_stencil_operator/P2ConstantEpsilonOperator.hpp"
#include "hyteg/forms/form_hyteg_generated/p2/p2_epsilonvar_affine_q4.hpp"
#include "hyteg/p2functionspace/P2Function.hpp"
--
GitLab
From e84da78aaf5459865e132ec43fc20d2c95788dd8 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Tue, 28 May 2024 13:55:11 +0200
Subject: [PATCH 22/33] Passing correct geometries already in the operator info
in generate_all_operators.py.
---
generate_all_operators.py | 51 +++++++++++++++++++--------------------
1 file changed, 25 insertions(+), 26 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index ae1c2f2..36a1b66 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -353,16 +353,21 @@ def main():
args.quad_degree if args.quad_rule is None else args.quad_rule,
)
}
+
+ enabled_geometries: Set[TriangleElement | TetrahedronElement] = set()
+ if 2 in args.dimensions:
+ enabled_geometries.add(TriangleElement())
+ if 3 in args.dimensions:
+ enabled_geometries.add(TetrahedronElement())
+
operators = all_operators(
symbolizer,
[(opts, loop_strategy, ordered_opts_suffix(loop_strategy, opts))],
type_descriptor,
blending,
+ geometries=enabled_geometries,
)
- re_mapping = re.compile(args.space_mapping)
- re_form = re.compile(args.form)
-
filtered_operators = list(
filter(
lambda o: re.fullmatch(args.space_mapping, o.mapping, re.IGNORECASE)
@@ -371,15 +376,6 @@ def main():
)
)
- enabled_geometries: Set[TriangleElement | TetrahedronElement] = set()
- if 2 in args.dimensions:
- enabled_geometries.add(TriangleElement())
- if 3 in args.dimensions:
- enabled_geometries.add(TetrahedronElement())
-
- for f in filtered_operators:
- f.geometries = list(set(f.geometries) & enabled_geometries)
-
filtered_operators = [f for f in filtered_operators if f.geometries]
if len(filtered_operators) == 0:
@@ -542,6 +538,7 @@ def all_operators(
opts: List[Tuple[Set[Opts], LoopStrategy, str]],
type_descriptor: HOGType,
blending: GeometryMap,
+ geometries: Set[ElementGeometry],
) -> List[OperatorInfo]:
P1 = LagrangianFunctionSpace(1, symbolizer)
P1Vector = TensorialVectorFunctionSpace(P1)
@@ -549,8 +546,8 @@ def all_operators(
P2Vector = TensorialVectorFunctionSpace(P2)
N1E1 = N1E1Space(symbolizer)
- two_d = [TriangleElement()]
- three_d = [TetrahedronElement()]
+ two_d = list({TriangleElement()} | geometries)
+ three_d = list({TetrahedronElement()} | geometries)
ops: List[OperatorInfo] = []
@@ -564,30 +561,31 @@ def all_operators(
form=partial(curl_curl_plus_mass, alpha_fem_space=P1, beta_fem_space=P1),
type_descriptor=type_descriptor, geometries=three_d, opts=opts, blending=blending))
ops.append(OperatorInfo(mapping="P1", name="Diffusion", trial_space=P1, test_space=P1, form=diffusion,
- type_descriptor=type_descriptor, opts=opts, blending=blending))
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo(mapping="P1", name="DivKGrad", trial_space=P1, test_space=P1,
form=partial(div_k_grad, coefficient_function_space=P1),
- type_descriptor=type_descriptor, opts=opts, blending=blending))
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo(mapping="P2", name="Diffusion", trial_space=P2, test_space=P2, form=diffusion,
- type_descriptor=type_descriptor, opts=opts, blending=blending))
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo(mapping="P2", name="DivKGrad", trial_space=P2, test_space=P2,
form=partial(div_k_grad, coefficient_function_space=P2),
- type_descriptor=type_descriptor, opts=opts, blending=blending))
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+
ops.append(OperatorInfo(mapping="P2", name="ShearHeating", trial_space=P2, test_space=P2,
form=partial(shear_heating, viscosity_function_space=P2, velocity_function_space=P2),
- type_descriptor=type_descriptor, opts=opts, blending=blending))
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo(mapping="P1", name="NonlinearDiffusion", trial_space=P1, test_space=P1,
form=partial(nonlinear_diffusion, coefficient_function_space=P1),
- type_descriptor=type_descriptor, opts=opts, blending=blending))
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo(mapping="P1", name="NonlinearDiffusionNewtonGalerkin", trial_space=P1,
test_space=P1, form=partial(nonlinear_diffusion_newton_galerkin,
coefficient_function_space=P1, onlyNewtonGalerkinPartOfForm=False),
- type_descriptor=type_descriptor, opts=opts, blending=blending))
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo(mapping="P1Vector", name="Diffusion", trial_space=P1Vector, test_space=P1Vector,
- form=diffusion, type_descriptor=type_descriptor, opts=opts, blending=blending))
+ form=diffusion, type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
# fmt: on
@@ -604,6 +602,7 @@ def all_operators(
test_space=P2Vector,
form=p2vec_epsilon,
type_descriptor=type_descriptor,
+ geometries=list(geometries),
opts=opts,
blending=blending,
)
@@ -611,9 +610,9 @@ def all_operators(
for c in [0, 1, 2]:
# fmt: off
- div_geometries = [TriangleElement(), TetrahedronElement()]
+ div_geometries = geometries
if c == 2:
- div_geometries = [TetrahedronElement()]
+ div_geometries = three_d
ops.append(OperatorInfo(mapping=f"P2ToP1", name=f"Div_{c}", trial_space=P1, test_space=P2,
form=partial(divergence, transpose=False, component_index=c),
type_descriptor=type_descriptor, opts=opts, geometries=div_geometries,
@@ -644,9 +643,9 @@ def all_operators(
# fmt: off
ops.append(
OperatorInfo(mapping=f"P2", name=f"Epsilon_{r}_{c}", trial_space=P2, test_space=P2, form=p2_epsilon,
- type_descriptor=type_descriptor, opts=opts, blending=blending))
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo(mapping=f"P2", name=f"FullStokes_{r}_{c}", trial_space=P2, test_space=P2,
- form=p2_full_stokes, type_descriptor=type_descriptor, opts=opts,
+ form=p2_full_stokes, type_descriptor=type_descriptor, geometries=list(geometries), opts=opts,
blending=blending))
# fmt: on
for c, r in [(0, 2), (1, 2), (2, 2), (2, 1), (2, 0)]:
--
GitLab
From 26bb27787180b70e7e530d722f6565c2be1c8e08 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Tue, 28 May 2024 14:03:00 +0200
Subject: [PATCH 23/33] Typing issue in generate_all_operators.py.
---
generate_all_operators.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 36a1b66..a2539bc 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -610,7 +610,7 @@ def all_operators(
for c in [0, 1, 2]:
# fmt: off
- div_geometries = geometries
+ div_geometries = list(geometries)
if c == 2:
div_geometries = three_d
ops.append(OperatorInfo(mapping=f"P2ToP1", name=f"Div_{c}", trial_space=P1, test_space=P2,
--
GitLab
From 141eca525da85e88a4463be718a0e8bed78ebcfa Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Tue, 28 May 2024 14:09:45 +0200
Subject: [PATCH 24/33] Typing issue II in generate_all_operators.py.
---
generate_all_operators.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index a2539bc..6fda1e5 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -22,7 +22,7 @@ from functools import partial
import logging
import os
import re
-from typing import Callable, List, Optional, Sequence, Set, Tuple
+from typing import Callable, List, Optional, Sequence, Set, Tuple, Union
import numpy as np
import sympy as sp
@@ -538,7 +538,7 @@ def all_operators(
opts: List[Tuple[Set[Opts], LoopStrategy, str]],
type_descriptor: HOGType,
blending: GeometryMap,
- geometries: Set[ElementGeometry],
+ geometries: Set[Union[TriangleElement, TetrahedronElement]],
) -> List[OperatorInfo]:
P1 = LagrangianFunctionSpace(1, symbolizer)
P1Vector = TensorialVectorFunctionSpace(P1)
--
GitLab
From ba283efac2f848e82d4c22c68b5c9a1772c685f2 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Tue, 28 May 2024 15:31:21 +0200
Subject: [PATCH 25/33] Fixed set intersection in generate_all_operators.py.
---
generate_all_operators.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 6fda1e5..eb5fd23 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -546,8 +546,8 @@ def all_operators(
P2Vector = TensorialVectorFunctionSpace(P2)
N1E1 = N1E1Space(symbolizer)
- two_d = list({TriangleElement()} | geometries)
- three_d = list({TetrahedronElement()} | geometries)
+ two_d = list({TriangleElement()} & geometries)
+ three_d = list({TetrahedronElement()} & geometries)
ops: List[OperatorInfo] = []
--
GitLab
From ad544078e4ce80f3f00639c14df60cd260a06901 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Tue, 28 May 2024 15:48:22 +0200
Subject: [PATCH 26/33] Another type fix.
mypy is driving me insane.
---
generate_all_operators.py | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index eb5fd23..9113125 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -610,9 +610,10 @@ def all_operators(
for c in [0, 1, 2]:
# fmt: off
- div_geometries = list(geometries)
if c == 2:
div_geometries = three_d
+ else:
+ div_geometries = list(geometries)
ops.append(OperatorInfo(mapping=f"P2ToP1", name=f"Div_{c}", trial_space=P1, test_space=P2,
form=partial(divergence, transpose=False, component_index=c),
type_descriptor=type_descriptor, opts=opts, geometries=div_geometries,
--
GitLab
From 6133ddb2dd42738150a4dccc64b013d2a93b8815 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Tue, 28 May 2024 16:11:28 +0200
Subject: [PATCH 27/33] Another type fix (II).
mypy is driving me insane.
---
generate_all_operators.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 9113125..ea28b40 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -538,7 +538,7 @@ def all_operators(
opts: List[Tuple[Set[Opts], LoopStrategy, str]],
type_descriptor: HOGType,
blending: GeometryMap,
- geometries: Set[Union[TriangleElement, TetrahedronElement]],
+ geometries: Set,
) -> List[OperatorInfo]:
P1 = LagrangianFunctionSpace(1, symbolizer)
P1Vector = TensorialVectorFunctionSpace(P1)
--
GitLab
From 4e303f6e24d9f9c785ec28f8e5bd161bc9356008 Mon Sep 17 00:00:00 2001
From: Daniel Bauer <daniel.j.bauer@fau.de>
Date: Wed, 29 May 2024 11:01:49 +0200
Subject: [PATCH 28/33] fix typing
---
generate_all_operators.py | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index ea28b40..4159e66 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -538,7 +538,7 @@ def all_operators(
opts: List[Tuple[Set[Opts], LoopStrategy, str]],
type_descriptor: HOGType,
blending: GeometryMap,
- geometries: Set,
+ geometries: Set[Union[TriangleElement, TetrahedronElement]],
) -> List[OperatorInfo]:
P1 = LagrangianFunctionSpace(1, symbolizer)
P1Vector = TensorialVectorFunctionSpace(P1)
@@ -546,8 +546,8 @@ def all_operators(
P2Vector = TensorialVectorFunctionSpace(P2)
N1E1 = N1E1Space(symbolizer)
- two_d = list({TriangleElement()} & geometries)
- three_d = list({TetrahedronElement()} & geometries)
+ two_d = list(geometries & {TriangleElement()})
+ three_d = list(geometries & {TetrahedronElement()})
ops: List[OperatorInfo] = []
--
GitLab
From 69b2f873901beb3fb288263884c558790aa06217 Mon Sep 17 00:00:00 2001
From: Daniel Bauer <daniel.j.bauer@fau.de>
Date: Wed, 29 May 2024 11:14:23 +0200
Subject: [PATCH 29/33] fix inverse diagonal kernel
---
hog/operator_generation/kernel_types.py | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/hog/operator_generation/kernel_types.py b/hog/operator_generation/kernel_types.py
index c2c781a..778077c 100644
--- a/hog/operator_generation/kernel_types.py
+++ b/hog/operator_generation/kernel_types.py
@@ -396,6 +396,7 @@ class AssembleDiagonal(KernelType):
f' this->timingTree_->stop( "pre-communication" );\n'
f"\n"
f"{indent(macro_loop(3), 2 * INDENT)}\n"
+ f"{indent(self.dst.invert_elementwise(3), 2 * INDENT)}\n"
f" }}\n"
f" else\n"
f" {{\n"
@@ -404,9 +405,9 @@ class AssembleDiagonal(KernelType):
f' this->timingTree_->stop( "pre-communication" );\n'
f"\n"
f"{indent(macro_loop(2), 2 * INDENT)}\n"
+ f"{indent(self.dst.invert_elementwise(2), 2 * INDENT)}\n"
f" }}\n"
f"\n"
- f"{indent(self.dst.invert_elementwise(dim=0), INDENT)}\n"
f"}}\n"
f"\n"
f'this->stopTiming( "{self.name}" );'
--
GitLab
From e63f998a404c134f4908419ba34731b20d874edc Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 30 May 2024 13:19:35 +0200
Subject: [PATCH 30/33] Lowering quadrature degree of vector epsilon operator
in HyTeG integration test to match that of the reference operator. Otherwise,
the comparison tests do not pass.
---
hyteg_integration_tests/src/CMakeLists.txt | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index c84aea7..2f35d52 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -59,7 +59,7 @@ add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM D
add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 )
add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 )
-add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 --dimensions 2 LIBS constant_stencil_operator)
+add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2 --dimensions 2 LIBS constant_stencil_operator)
# tests with blending
--
GitLab
From 81c8b53c32a3f7ecc241ae4bb9e28b8dc2036c2d Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Thu, 30 May 2024 16:34:13 +0200
Subject: [PATCH 31/33] Added a dedicated P2Vector epsilon op blending test in
2D (testing combination vectorspace + blending).
---
hyteg_integration_tests/src/CMakeLists.txt | 33 +-----
.../src/EpsilonVectorAnnulus.cpp | 106 ++++++++++++++++++
2 files changed, 108 insertions(+), 31 deletions(-)
create mode 100644 hyteg_integration_tests/src/EpsilonVectorAnnulus.cpp
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index 2f35d52..03c3b8f 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -37,35 +37,6 @@ endfunction()
# tests without blending
-add_operator_test(FILE CurlCurl.cpp DEF TEST_ASSEMBLE FORM CurlCurl ABBR noopts GEN_ARGS --quad-degree 0 0)
-add_operator_test(FILE CurlCurl.cpp FORM CurlCurl ABBR VIQ GEN_ARGS --quad-degree 0 0 -o MOVECONSTANTS VECTORIZE QUADLOOPS)
-add_operator_test(FILE CurlCurl.cpp FORM CurlCurl ABBR CVIP GEN_ARGS --quad-degree 0 0 --loop-strategy CUBES -o MOVECONSTANTS VECTORIZE POLYCSE)
-add_operator_test(FILE CurlCurlPlusMass.cpp DEF TEST_ASSEMBLE FORM CurlCurlPlusMass ABBR noopts GEN_ARGS --quad-degree 2 2)
-add_operator_test(FILE CurlCurlPlusMass.cpp FORM CurlCurlPlusMass ABBR VI GEN_ARGS --quad-degree 2 2 -o MOVECONSTANTS VECTORIZE)
-
-add_operator_test(FILE DiffusionP1.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1 GEN_ARGS -s P1 --quad-degree 0 0)
-add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1VIP GEN_ARGS -s P1 --quad-degree 0 0 -o MOVECONSTANTS VECTORIZE POLYCSE)
-add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1CVI GEN_ARGS -s P1 --quad-degree 0 0 --loop-strategy CUBES -o MOVECONSTANTS VECTORIZE)
-add_operator_test(FILE DiffusionP2.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 2 GEN_ARGS -s P2 --quad-degree 2 2)
-add_operator_test(FILE DiffusionP1Vector.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1Vec GEN_ARGS -s P1Vector --quad-degree 0 0 LIBS mixed_operator)
-
-add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR noopts GEN_ARGS -s P2 --quad-degree 4 4)
-add_operator_test(FILE DivKGrad.cpp FORM DivKGrad ABBR VIQ GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS VECTORIZE QUADLOOPS )
-add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR IQT GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS QUADLOOPS TABULATE )
-
-add_operator_test(FILE Div.cpp DEF REF_OP=P2ToP1ElementwiseDivxOperator FORM Div_0 ABBR VIQT GEN_ARGS -s P2ToP1 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
-add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM DivT_2 ABBR VIQT GEN_ARGS -s P1ToP2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
-
-add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 )
-add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4 )
-add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2 --dimensions 2 LIBS constant_stencil_operator)
-
-# tests with blending
-
-add_operator_test(FILE DiffusionP1Annulus.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR b2 GEN_ARGS -s P1 -b AnnulusMap --quad-rule hillion_07 yu_3)
-add_operator_test(FILE DiffusionP1Annulus.cpp FORM Diffusion ABBR b2VIQT GEN_ARGS -s P1 -b AnnulusMap -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-rule hillion_07 yu_3)
-add_operator_test(FILE DiffusionP1IcosahedralShell.cpp FORM Diffusion ABBR b3VIQT GEN_ARGS -s P1 -b IcosahedralShellMap -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-rule hillion_07 yu_3)
-
-add_operator_test(FILE FullStokes.cpp DEF TEST_DIAG FORM=forms::p2_full_stokesvar_0_0_blending_q3 FORM FullStokes_0_0 ABBR 00b3VIQT GEN_ARGS -s P2 -b IcosahedralShellMap -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-rule yu_3 yu_3)
-add_operator_test(FILE FullStokes.cpp DEF FORM=forms::p2_full_stokesvar_2_1_blending_q3 FORM FullStokes_2_1 ABBR 21b3VIQT GEN_ARGS -s P2 -b IcosahedralShellMap -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-rule yu_3 yu_3)
+add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE --quad-degree 2 2 --dimensions 2 LIBS constant_stencil_operator)
+add_operator_test(FILE EpsilonVectorAnnulus.cpp DEF FORM Epsilon ABBR VecbVIQT GEN_ARGS -s P2Vector -b AnnulusMap -o MOVECONSTANTS QUADLOOPS TABULATE --quad-degree 2 2 --dimensions 2 LIBS constant_stencil_operator)
diff --git a/hyteg_integration_tests/src/EpsilonVectorAnnulus.cpp b/hyteg_integration_tests/src/EpsilonVectorAnnulus.cpp
new file mode 100644
index 0000000..89a053a
--- /dev/null
+++ b/hyteg_integration_tests/src/EpsilonVectorAnnulus.cpp
@@ -0,0 +1,106 @@
+/*
+ * HyTeG Operator Generator
+ * Copyright (C) 2017-2024 Nils Kohl, Fabian Böhm, Daniel Bauer
+ *
+ * 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/>.
+ */
+
+#include <type_traits>
+
+#include "core/DataTypes.h"
+
+#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
+#include "constant_stencil_operator/P2ConstantEpsilonOperator.hpp"
+#include "hyteg/forms/form_hyteg_generated/p2/p2_epsilonvar_affine_q4.hpp"
+#include "hyteg/p2functionspace/P2Function.hpp"
+
+#include "Epsilon/TestOpEpsilon.hpp"
+#include "OperatorGenerationTest.hpp"
+
+using namespace hyteg;
+using walberla::real_t;
+
+real_t k( const hyteg::Point3D& x )
+{
+ // The operator also works with a non-constant viscosity.
+ // However, the reference operator seems to have such a low quadrature order that the difference between the
+ // reference and generated operator is relatively large.
+ // Under refinement that error vanishes, but not quickly enough that we can afford refinement in this (quick) test.
+ // For simplicity, the viscosity is therefore set to 1 here, such that the comparison succeeds (just blending, no var
+ // viscosity). Note that there is a dedicated test with a variable viscosity, but without blending.
+ return 1.0;
+}
+
+P2ElementwiseBlendingEpsilonOperator makeRefOp( std::shared_ptr< PrimitiveStorage > storage, uint_t minLevel, uint_t maxLevel )
+{
+ return P2ElementwiseBlendingEpsilonOperator( storage, minLevel, maxLevel, k );
+};
+
+template < class Op >
+Op makeTestOp( std::shared_ptr< PrimitiveStorage > storage, uint_t minLevel, uint_t maxLevel )
+{
+ P2Function< real_t > k( "k", storage, minLevel, maxLevel );
+ for ( size_t lvl = minLevel; lvl <= maxLevel; ++lvl )
+ {
+ k.interpolate( ::k, lvl );
+ }
+ return Op( storage, minLevel, maxLevel, k );
+};
+
+int main( int argc, char* argv[] )
+{
+ walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+ walberla::MPIManager::instance()->useWorldComm();
+
+ const uint_t level = 3;
+
+ real_t thresholdOverMachineEpsApply = 2e3;
+ real_t thresholdOverMachineEpsInvDiag = 9.0e6;
+ real_t thresholdOverMachineEpsAssembly = 360;
+
+ // Testing only 2D since generation of 3D operator just takes too long.
+
+ StorageSetup storageSetup(
+ "Annulus", MeshInfo::meshAnnulus( 1.0, 2.0, MeshInfo::CRISS, 12, 2 ), GeometryMap::Type::ANNULUS );
+
+ compareApply< P2ElementwiseBlendingEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Apply",
+ thresholdOverMachineEpsApply );
+
+#ifdef TEST_DIAG
+ compareInvDiag< P2Function< real_t >, P2ElementwiseBlendingEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Inverse Diagonal",
+ thresholdOverMachineEpsInvDiag );
+#endif
+
+#ifdef TEST_ASSEMBLE
+ compareAssembledMatrix< P2ElementwiseBlendingEpsilonOperator, operatorgeneration::TestOpEpsilon >(
+ makeRefOp,
+ makeTestOp< operatorgeneration::TestOpEpsilon >,
+ level,
+ storageSetup,
+ storageSetup.description() + " Assembly",
+ thresholdOverMachineEpsAssembly );
+#endif
+
+ return EXIT_SUCCESS;
+}
--
GitLab
From 58f76ba13126bb33b793fa90ed6366098d2db4e6 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Fri, 31 May 2024 09:33:59 +0200
Subject: [PATCH 32/33] Re-added tests I accidentally removed.
---
hyteg_integration_tests/src/CMakeLists.txt | 30 ++++++++++++++++++++++
1 file changed, 30 insertions(+)
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index 03c3b8f..54837bb 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -37,6 +37,36 @@ endfunction()
# tests without blending
+add_operator_test(FILE CurlCurl.cpp DEF TEST_ASSEMBLE FORM CurlCurl ABBR noopts GEN_ARGS --quad-degree 0 0)
+add_operator_test(FILE CurlCurl.cpp FORM CurlCurl ABBR VIQ GEN_ARGS --quad-degree 0 0 -o MOVECONSTANTS VECTORIZE QUADLOOPS)
+add_operator_test(FILE CurlCurl.cpp FORM CurlCurl ABBR CVIP GEN_ARGS --quad-degree 0 0 --loop-strategy CUBES -o MOVECONSTANTS VECTORIZE POLYCSE)
+
+add_operator_test(FILE CurlCurlPlusMass.cpp DEF TEST_ASSEMBLE FORM CurlCurlPlusMass ABBR noopts GEN_ARGS --quad-degree 2 2)
+add_operator_test(FILE CurlCurlPlusMass.cpp FORM CurlCurlPlusMass ABBR VI GEN_ARGS --quad-degree 2 2 -o MOVECONSTANTS VECTORIZE)
+
+add_operator_test(FILE DiffusionP1.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 1 GEN_ARGS -s P1 --quad-degree 0 0)
+add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1VIP GEN_ARGS -s P1 --quad-degree 0 0 -o MOVECONSTANTS VECTORIZE POLYCSE)
+add_operator_test(FILE DiffusionP1.cpp FORM Diffusion ABBR 1CVI GEN_ARGS -s P1 --quad-degree 0 0 --loop-strategy CUBES -o MOVECONSTANTS VECTORIZE)
+add_operator_test(FILE DiffusionP2.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR 2 GEN_ARGS -s P2 --quad-degree 2 2)
+
+add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR noopts GEN_ARGS -s P2 --quad-degree 4 4)
+add_operator_test(FILE DivKGrad.cpp FORM DivKGrad ABBR VIQ GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS VECTORIZE QUADLOOPS )
+add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR IQT GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS QUADLOOPS TABULATE )
+
+add_operator_test(FILE Div.cpp DEF REF_OP=P2ToP1ElementwiseDivxOperator FORM Div_0 ABBR VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
+add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM DivT_2 ABBR VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
+add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
+add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
add_operator_test(FILE EpsilonVector.cpp DEF FORM Epsilon ABBR VecVIQT GEN_ARGS -s P2Vector -o MOVECONSTANTS QUADLOOPS TABULATE --quad-degree 2 2 --dimensions 2 LIBS constant_stencil_operator)
add_operator_test(FILE EpsilonVectorAnnulus.cpp DEF FORM Epsilon ABBR VecbVIQT GEN_ARGS -s P2Vector -b AnnulusMap -o MOVECONSTANTS QUADLOOPS TABULATE --quad-degree 2 2 --dimensions 2 LIBS constant_stencil_operator)
+
+# tests with blending
+
+add_operator_test(FILE DiffusionP1Annulus.cpp DEF TEST_ASSEMBLE FORM Diffusion ABBR b2 GEN_ARGS -s P1 -b AnnulusMap --quad-rule hillion_07 yu_3)
+add_operator_test(FILE DiffusionP1Annulus.cpp FORM Diffusion ABBR b2VIQT GEN_ARGS -s P1 -b AnnulusMap -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-rule hillion_07 yu_3)
+add_operator_test(FILE DiffusionP1IcosahedralShell.cpp FORM Diffusion ABBR b3VIQT GEN_ARGS -s P1 -b IcosahedralShellMap -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-rule hillion_07 yu_3)
+
+add_operator_test(FILE FullStokes.cpp DEF TEST_DIAG FORM=forms::p2_full_stokesvar_0_0_blending_q3 FORM FullStokes_0_0 ABBR 00b3VIQT GEN_ARGS -s P2 -b IcosahedralShellMap -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-rule yu_3 yu_3)
+add_operator_test(FILE FullStokes.cpp DEF FORM=forms::p2_full_stokesvar_2_1_blending_q3 FORM FullStokes_2_1 ABBR 21b3VIQT GEN_ARGS -s P2 -b IcosahedralShellMap -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-rule yu_3 yu_3)
+
--
GitLab
From a374af87b717f18fe9cd7cc0f60bae139fe1cd68 Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Fri, 31 May 2024 09:58:18 +0200
Subject: [PATCH 33/33] Fixing function spaces for div and divt.
---
hyteg_integration_tests/src/CMakeLists.txt | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/hyteg_integration_tests/src/CMakeLists.txt b/hyteg_integration_tests/src/CMakeLists.txt
index 54837bb..66b018d 100644
--- a/hyteg_integration_tests/src/CMakeLists.txt
+++ b/hyteg_integration_tests/src/CMakeLists.txt
@@ -53,8 +53,8 @@ add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR noopts
add_operator_test(FILE DivKGrad.cpp FORM DivKGrad ABBR VIQ GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS VECTORIZE QUADLOOPS )
add_operator_test(FILE DivKGrad.cpp DEF TEST_ASSEMBLE FORM DivKGrad ABBR IQT GEN_ARGS -s P2 --quad-degree 4 4 -o MOVECONSTANTS QUADLOOPS TABULATE )
-add_operator_test(FILE Div.cpp DEF REF_OP=P2ToP1ElementwiseDivxOperator FORM Div_0 ABBR VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
-add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM DivT_2 ABBR VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
+add_operator_test(FILE Div.cpp DEF REF_OP=P2ToP1ElementwiseDivxOperator FORM Div_0 ABBR VIQT GEN_ARGS -s P2ToP1 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
+add_operator_test(FILE DivT.cpp DEF REF_OP=P1ToP2ElementwiseDivTzOperator FORM DivT_2 ABBR VIQT GEN_ARGS -s P1ToP2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 2 2)
add_operator_test(FILE Epsilon.cpp DEF TEST_DIAG FORM=forms::p2_epsilonvar_0_0_affine_q4 FORM Epsilon_0_0 ABBR 00VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
add_operator_test(FILE Epsilon.cpp DEF FORM=forms::p2_epsilonvar_2_1_affine_q4 FORM Epsilon_2_1 ABBR 21VIQT GEN_ARGS -s P2 -o MOVECONSTANTS QUADLOOPS TABULATE VECTORIZE --quad-degree 4 4)
--
GitLab