diff --git a/generate_all_operators.py b/generate_all_operators.py
index 50a66bcd0e613cad016ff67423ef3ebfba7e3e31..6ce4a12dcab8ddddd60fb2e05107591d8415a84a 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -49,6 +49,8 @@ from hog.forms import (
nonlinear_diffusion,
nonlinear_diffusion_newton_galerkin,
supg_diffusion,
+ advection,
+ supg_advection,
grad_rho_by_rho_dot_u,
)
from hog.forms_boundary import (
@@ -628,6 +630,18 @@ def all_operators(
ops.append(OperatorInfo("P2", "SUPGDiffusion", TrialSpace(P2), TestSpace(P2),
form=partial(supg_diffusion, velocity_function_space=P2, diffusivityXdelta_function_space=P2),
type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+
+ ops.append(OperatorInfo("P2", "AdvectionConstCp", TrialSpace(P2), TestSpace(P2),
+ form=partial(advection, velocity_function_space=P2, coefficient_function_space=P2, constant_cp = True),
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+
+ ops.append(OperatorInfo("P2", "AdvectionVarCp", TrialSpace(P2), TestSpace(P2),
+ form=partial(advection, velocity_function_space=P2, coefficient_function_space=P2),
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+
+ ops.append(OperatorInfo("P2", "SUPGAdvection", TrialSpace(P2), TestSpace(P2),
+ form=partial(supg_advection, velocity_function_space=P2, coefficient_function_space=P2),
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo("P2", "BoundaryMass", TrialSpace(P2), TestSpace(P2), form=None,
form_boundary=mass_boundary, type_descriptor=type_descriptor, geometries=list(geometries),
diff --git a/hog/forms.py b/hog/forms.py
index b3f87b4d9d92f2f43a7797b67904bd7339af486c..774752bb4c99066c8a46b8157c2bb01c7bc11209 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -812,6 +812,97 @@ Weak formulation
)
+def advection(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ velocity_function_space: FunctionSpace,
+ coefficient_function_space: FunctionSpace,
+ constant_cp: bool = False,
+ blending: GeometryMap = IdentityMap(),
+) -> Form:
+ docstring = f"""
+advection operator which needs to be used in combination with SUPG
+
+Geometry map: {blending}
+
+Weak formulation
+
+ T: trial function (scalar space: {trial})
+ s: test function (scalar space: {test})
+ u: velocity function (vectorial space: {velocity_function_space})
+
+ ∫ cp ( u · ∇T ) s
+"""
+
+ from hog.recipes.integrands.volume.advection import integrand
+
+ return process_integrand(
+ integrand,
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ is_symmetric=False,
+ docstring=docstring,
+ fe_coefficients={
+ "ux": velocity_function_space,
+ "uy": velocity_function_space,
+ "uz": velocity_function_space,
+ } if constant_cp else {
+ "ux": velocity_function_space,
+ "uy": velocity_function_space,
+ "uz": velocity_function_space,
+ "cp": coefficient_function_space,
+ },
+ )
+
+
+def supg_advection(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ velocity_function_space: FunctionSpace,
+ coefficient_function_space: FunctionSpace,
+ blending: GeometryMap = IdentityMap(),
+) -> Form:
+ docstring = f"""
+advection operator which needs to be used in combination with SUPG
+
+Geometry map: {blending}
+
+Weak formulation
+
+ T: trial function (scalar space: {trial})
+ s: test function (scalar space: {test})
+ u: velocity function (vectorial space: {velocity_function_space})
+
+ ∫ cp ( u · ∇T ) 𝛿(u · ∇s)
+"""
+
+ from hog.recipes.integrands.volume.supg_advection import integrand
+
+ return process_integrand(
+ integrand,
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ is_symmetric=False,
+ docstring=docstring,
+ fe_coefficients={
+ "ux": velocity_function_space,
+ "uy": velocity_function_space,
+ "uz": velocity_function_space,
+ "cp_times_delta": coefficient_function_space,
+ },
+ )
+
+
def supg_diffusion(
trial: TrialSpace,
test: TestSpace,
diff --git a/hog/recipes/integrands/volume/advection.py b/hog/recipes/integrands/volume/advection.py
new file mode 100644
index 0000000000000000000000000000000000000000..5bdaeacf294a68dd6795471a625b78235811375e
--- /dev/null
+++ b/hog/recipes/integrands/volume/advection.py
@@ -0,0 +1,51 @@
+# HyTeG Operator Generator
+# Copyright (C) 2024 HyTeG Team
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+from hog.recipes.common import *
+
+
+def integrand(
+ *,
+ v,
+ jac_a_inv,
+ jac_a_abs_det,
+ jac_b_inv,
+ jac_b_abs_det,
+ grad_u,
+ grad_v,
+ k,
+ scalars,
+ volume_geometry,
+ tabulate,
+ **_,
+):
+ if volume_geometry.dimensions > 2:
+ u = sp.Matrix([[k["ux"]], [k["uy"]], [k["uz"]]])
+ else:
+ u = sp.Matrix([[k["ux"]], [k["uy"]]])
+
+ if "cp" in k.keys():
+ coeff = k["cp"]
+ else:
+ coeff = scalars("cp")
+
+ return (
+ coeff
+ * dot(jac_b_inv.T * tabulate(jac_a_inv.T * grad_u), u)
+ * v
+ * tabulate(jac_a_abs_det)
+ * jac_b_abs_det
+ )
diff --git a/hog/recipes/integrands/volume/supg_advection.py b/hog/recipes/integrands/volume/supg_advection.py
new file mode 100644
index 0000000000000000000000000000000000000000..560bf45b97216dafd68e49679f6e9fdffd3955ee
--- /dev/null
+++ b/hog/recipes/integrands/volume/supg_advection.py
@@ -0,0 +1,45 @@
+# HyTeG Operator Generator
+# Copyright (C) 2024 HyTeG Team
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+from hog.recipes.common import *
+
+
+def integrand(
+ *,
+ v,
+ jac_a_inv,
+ jac_a_abs_det,
+ jac_b_inv,
+ jac_b_abs_det,
+ grad_u,
+ grad_v,
+ k,
+ volume_geometry,
+ tabulate,
+ **_,
+):
+ if volume_geometry.dimensions > 2:
+ u = sp.Matrix([[k["ux"]], [k["uy"]], [k["uz"]]])
+ else:
+ u = sp.Matrix([[k["ux"]], [k["uy"]]])
+
+ return (
+ k["cp_times_delta"]
+ * dot(jac_b_inv.T * tabulate(jac_a_inv.T * grad_u), u)
+ * dot(jac_b_inv.T * tabulate(jac_a_inv.T * grad_v), u)
+ * tabulate(jac_a_abs_det)
+ * jac_b_abs_det
+ )