Merge branch 'ponsuganth/advection-form' into 'main'

Add advection and SUPG advection forms

See merge request !23

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),

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,

hog/recipes/integrands/volume/advection.py

+# 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
+    )

hog/recipes/integrands/volume/supg_advection.py

+# 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
+    )