diff --git a/generate_all_operators.py b/generate_all_operators.py
index 6ce4a12dcab8ddddd60fb2e05107591d8415a84a..05c7ee2dabc57b02e0b3baf93805bfb13d833b73 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -51,7 +51,8 @@ from hog.forms import (
supg_diffusion,
advection,
supg_advection,
- grad_rho_by_rho_dot_u,
+ grad_rho_by_rho_dot_u, nonlinear_diffusion_newton_galerkin_symmetric_part,
+ nonlinear_diffusion_newton_galerkin_asymmetric_part, nonlinear_diffusion_newton_galerkin_asymmetric_lumped_part,
)
from hog.forms_boundary import (
mass_boundary,
@@ -623,6 +624,18 @@ def all_operators(
TestSpace(P1), form=partial(nonlinear_diffusion_newton_galerkin,
coefficient_function_space=P1, only_newton_galerkin_part_of_form=False),
type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+ ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinSymmetricPart", TrialSpace(P1),
+ TestSpace(P1), form=partial(nonlinear_diffusion_newton_galerkin_symmetric_part,
+ coefficient_function_space=P1),
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+ ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinAsymmetricPart", TrialSpace(P1),
+ TestSpace(P1), form=partial(nonlinear_diffusion_newton_galerkin_asymmetric_part,
+ coefficient_function_space=P1),
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+ ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinAsymmetricLumpedPart", TrialSpace(P1),
+ TestSpace(P1), form=partial(nonlinear_diffusion_newton_galerkin_asymmetric_lumped_part,
+ coefficient_function_space=P1),
+ type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
ops.append(OperatorInfo("P1Vector", "Diffusion", TrialSpace(P1Vector), TestSpace(P1Vector),
form=diffusion, type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
diff --git a/hog/forms.py b/hog/forms.py
index f67f0781d698489310f163fcd996231aecee66e0..fdee07a0c06ea7cd4545f6c6cf998b4781543149 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -317,6 +317,213 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
)
+def nonlinear_diffusion_newton_galerkin_symmetric_part(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ coefficient_function_space: FunctionSpace,
+ blending: GeometryMap = IdentityMap(),
+) -> Form:
+ docstring = f"""
+todo: heavy code duplication, this can be implemented with small changes to nonlinear_diffusion_newton_galerkin
+
+Symmetric part of the bi-linear form for the solution of the non-linear diffusion equation by a Newton-Galerkin approach
+
+Weak formulation
+
+ u: trial function (space: {trial})
+ v: test function (space: {test})
+ k: FE coefficient function (space: {coefficient_function_space})
+
+ ∫ a(k) ∇u · ∇v
+
+Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended for :math:`k = u`.
+"""
+ if trial != test:
+ raise HOGException(
+ "Trial space must be equal to test space to assemble diffusion matrix."
+ )
+
+ if blending != IdentityMap():
+ raise HOGException(
+ "The nonlinear_diffusion_newton_galerkin form does currently not support blending."
+ )
+
+ def integrand(
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ grad_u,
+ grad_v,
+ k,
+ tabulate,
+ **_,
+ ):
+ a = sp.Rational(1, 8) + k["k"] * k["k"]
+
+ diffusion_term = a * tabulate(
+ dot(jac_a_inv.T * grad_u, jac_a_inv.T * grad_v) * jac_a_abs_det
+ )
+
+ return diffusion_term
+
+ return process_integrand(
+ integrand,
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ is_symmetric=True,
+ docstring=docstring,
+ fe_coefficients={"k": coefficient_function_space},
+ )
+
+
+def nonlinear_diffusion_newton_galerkin_asymmetric_part(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ coefficient_function_space: FunctionSpace,
+ blending: GeometryMap = IdentityMap(),
+) -> Form:
+ docstring = f"""
+
+todo: heavy code duplication, this can be implemented with small changes to nonlinear_diffusion_newton_galerkin
+
+Asymmetric part of bi-linear form for the solution of the non-linear diffusion equation by a Newton-Galerkin approach
+
+Weak formulation
+
+ u: trial function (space: {trial})
+ v: test function (space: {test})
+ k: FE coefficient function (space: {coefficient_function_space})
+
+ ∫ a'(k) u ∇k · ∇v
+
+Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended for :math:`k = u`.
+"""
+ if trial != test:
+ raise HOGException(
+ "Trial space must be equal to test space to assemble diffusion matrix."
+ )
+
+ if blending != IdentityMap():
+ raise HOGException(
+ "The nonlinear_diffusion_newton_galerkin form does currently not support blending."
+ )
+
+ def integrand(
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ u,
+ grad_v,
+ k,
+ grad_k,
+ tabulate,
+ **_,
+ ):
+ a_prime = 2 * k["k"]
+
+ newton_galerkin_term = (
+ a_prime
+ * u
+ * dot(jac_a_inv.T * grad_k["k"], tabulate(jac_a_inv.T * grad_v))
+ * tabulate(jac_a_abs_det)
+ )
+
+ return newton_galerkin_term
+
+ return process_integrand(
+ integrand,
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ is_symmetric=False,
+ docstring=docstring,
+ fe_coefficients={"k": coefficient_function_space},
+ )
+
+
+def nonlinear_diffusion_newton_galerkin_asymmetric_lumped_part(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ coefficient_function_space: FunctionSpace,
+ blending: GeometryMap = IdentityMap(),
+) -> Form:
+ docstring = f"""
+
+todo: heavy code duplication, this can be implemented with small changes to nonlinear_diffusion_newton_galerkin
+
+Asymmetric part of bi-linear form for the solution of the non-linear diffusion equation by a Newton-Galerkin approach. Optimized with mass lumping.
+
+Weak formulation
+
+ u: trial function (space: {trial})
+ v: test function (space: {test})
+ k: FE coefficient function (space: {coefficient_function_space})
+
+ ∫ a'(k) u ∇k · ∇v
+
+Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended for :math:`k = u`.
+"""
+ if trial != test:
+ raise HOGException(
+ "Trial space must be equal to test space to assemble diffusion matrix."
+ )
+
+ if blending != IdentityMap():
+ raise HOGException(
+ "The nonlinear_diffusion_newton_galerkin form does currently not support blending."
+ )
+
+ def integrand(
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ u,
+ v,
+ grad_v,
+ k,
+ grad_k,
+ tabulate,
+ **_,
+ ):
+ # todo this turns into a diagonal matrix, can this further get optimized?
+ # todo can you get the indices of u and v and compare those instead of the whole expressions?
+ if not u.equals(v):
+ return 0
+ else:
+ a_prime = 2 * k["k"]
+
+ newton_galerkin_term = (
+ a_prime
+ * dot(jac_a_inv.T * grad_k["k"], tabulate(jac_a_inv.T * grad_v))
+ * tabulate(jac_a_abs_det)
+ )
+
+ return newton_galerkin_term
+
+ return process_integrand(
+ integrand,
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ is_symmetric=False,
+ docstring=docstring,
+ fe_coefficients={"k": coefficient_function_space},
+ )
+
+
def epsilon(
trial: TrialSpace,
test: TestSpace,
@@ -372,9 +579,11 @@ where
is_symmetric=trial == test,
docstring=docstring,
fe_coefficients={"mu": coefficient_function_space},
- rot_type=RotationType.PRE_AND_POST_MULTIPLY
- if rotation_wrapper
- else RotationType.NO_ROTATION,
+ rot_type=(
+ RotationType.PRE_AND_POST_MULTIPLY
+ if rotation_wrapper
+ else RotationType.NO_ROTATION
+ ),
)
@@ -580,9 +789,9 @@ Weak formulation
blending=blending,
is_symmetric=False,
docstring=docstring,
- rot_type=RotationType.POST_MULTIPLY
- if rotation_wrapper
- else RotationType.NO_ROTATION,
+ rot_type=(
+ RotationType.POST_MULTIPLY if rotation_wrapper else RotationType.NO_ROTATION
+ ),
)
@@ -621,9 +830,9 @@ def gradient(
blending=blending,
is_symmetric=False,
docstring=docstring,
- rot_type=RotationType.PRE_MULTIPLY
- if rotation_wrapper
- else RotationType.NO_ROTATION,
+ rot_type=(
+ RotationType.PRE_MULTIPLY if rotation_wrapper else RotationType.NO_ROTATION
+ ),
)
@@ -906,18 +1115,20 @@ Weak formulation
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,
- },
+ 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,
+ }
+ ),
)
@@ -1097,7 +1308,7 @@ Weak formulation
∫ ((∇ρ / ρ) · u) v
"""
-
+
from hog.recipes.integrands.volume.frozen_velocity import integrand
return process_integrand(
@@ -1114,6 +1325,7 @@ Weak formulation
},
)
+
def zero_form(
trial: TrialSpace,
test: TestSpace,