diff --git a/generate_all_operators.py b/generate_all_operators.py
index 05c7ee2dabc57b02e0b3baf93805bfb13d833b73..e8b67fe6c4d0d32b4be143f3ec971280a13d8d4c 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -51,8 +51,9 @@ from hog.forms import (
supg_diffusion,
advection,
supg_advection,
- 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,
+ grad_rho_by_rho_dot_u,
+ nonlinear_diffusion_stiffness,
+ nonlinear_diffusion_newton_galerkin_approximated_mass_lumped,
)
from hog.forms_boundary import (
mass_boundary,
@@ -73,7 +74,7 @@ from hog.operator_generation.kernel_types import (
ApplyWrapper,
AssembleWrapper,
AssembleDiagonalWrapper,
- KernelWrapperType,
+ KernelWrapperType, LumpWrapper,
)
from hog.operator_generation.loop_strategies import (
LoopStrategy,
@@ -573,6 +574,14 @@ class OperatorInfo:
)
)
+ self.kernel_types.append(
+ LumpWrapper(
+ self.trial_space,
+ type_descriptor=self.type_descriptor,
+ dims=dims,
+ )
+ )
+
def all_operators(
symbolizer: Symbolizer,
@@ -620,20 +629,24 @@ def all_operators(
ops.append(OperatorInfo("P1", "NonlinearDiffusion", TrialSpace(P1), TestSpace(P1),
form=partial(nonlinear_diffusion, coefficient_function_space=P1),
type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
+
ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkin", TrialSpace(P1),
TestSpace(P1), form=partial(nonlinear_diffusion_newton_galerkin,
- coefficient_function_space=P1, only_newton_galerkin_part_of_form=False),
+ coefficient_function_space=P1),
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),
+
+ ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinStiffness", TrialSpace(P1),
+ TestSpace(P1), form=partial(nonlinear_diffusion_stiffness,
+ 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),
+
+ ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinStiffnessDifferentiated", TrialSpace(P1),
+ TestSpace(P1), form=partial(nonlinear_diffusion_stiffness,
+ coefficient_function_space=P1, use_a_derivative=True),
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,
+
+ ops.append(OperatorInfo("P1", "NonlinearDiffusionNewtonGalerkinApproximatedMassLumped", TrialSpace(P1),
+ TestSpace(P1), form=partial(nonlinear_diffusion_newton_galerkin_approximated_mass_lumped,
coefficient_function_space=P1),
type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
diff --git a/hog/forms.py b/hog/forms.py
index fdee07a0c06ea7cd4545f6c6cf998b4781543149..b3338d8c6b98736c66ba8d97cef869368b8ddd6e 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -247,7 +247,6 @@ def nonlinear_diffusion_newton_galerkin(
symbolizer: Symbolizer,
coefficient_function_space: FunctionSpace,
blending: GeometryMap = IdentityMap(),
- only_newton_galerkin_part_of_form: Optional[bool] = True,
) -> Form:
docstring = f"""
@@ -299,10 +298,7 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
* tabulate(jac_a_abs_det)
)
- if only_newton_galerkin_part_of_form:
- return newton_galerkin_term
- else:
- return diffusion_term + newton_galerkin_term
+ return diffusion_term + newton_galerkin_term
return process_integrand(
integrand,
@@ -316,19 +312,17 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
fe_coefficients={"k": coefficient_function_space},
)
-
-def nonlinear_diffusion_newton_galerkin_symmetric_part(
- trial: TrialSpace,
- test: TestSpace,
- geometry: ElementGeometry,
- symbolizer: Symbolizer,
- coefficient_function_space: FunctionSpace,
- blending: GeometryMap = IdentityMap(),
+def nonlinear_diffusion_stiffness(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ coefficient_function_space: FunctionSpace,
+ blending: GeometryMap = IdentityMap(),
+ use_a_derivative: bool = False,
) -> 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
+todo
Weak formulation
@@ -337,9 +331,12 @@ Weak formulation
k: FE coefficient function (space: {coefficient_function_space})
∫ a(k) ∇u · ∇v
+
+ or ∫ a'(k) ∇u · ∇v if use_a_derivative was set
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."
@@ -351,18 +348,24 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
)
def integrand(
- *,
- jac_a_inv,
- jac_a_abs_det,
- grad_u,
- grad_v,
- k,
- tabulate,
- **_,
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ grad_u,
+ grad_v,
+ k,
+ tabulate,
+ **_,
):
a = sp.Rational(1, 8) + k["k"] * k["k"]
+ a_prime = 2 * k["k"]
- diffusion_term = a * tabulate(
+ if use_a_derivative:
+ used_a = a_prime
+ else:
+ used_a = a
+
+ diffusion_term = used_a * tabulate(
dot(jac_a_inv.T * grad_u, jac_a_inv.T * grad_v) * jac_a_abs_det
)
@@ -380,20 +383,19 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
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(),
+def nonlinear_diffusion_newton_galerkin_approximated_mass_lumped(
+ 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
+Bi-linear form with approximated mass lumping on the asymmetric part for the solution of the non-linear diffusion equation by a Newton-Galerkin approach
+This just uses the approximation 9.25 from Larson & Bengzon. This is not the actual mass lumping approach.
-Asymmetric part of bi-linear form for the solution of the non-linear diffusion equation by a Newton-Galerkin approach
Weak formulation
@@ -401,7 +403,9 @@ Weak formulation
v: test function (space: {test})
k: FE coefficient function (space: {coefficient_function_space})
- ∫ a'(k) u ∇k · ∇v
+ ∫ a(k) ∇u · ∇v + δ(i,j) · ∫ a'(k) ∇k · ∇v
+ where δ(i, j) = if i == j then 1 else 0, and i, j are the indices of the corresponding matrix element
+ i.e. δ = 1 on the diagonal and 0 otherwise
Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended for :math:`k = u`.
"""
@@ -416,100 +420,36 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
)
def integrand(
- *,
- jac_a_inv,
- jac_a_abs_det,
- u,
- grad_v,
- k,
- grad_k,
- tabulate,
- **_,
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ grad_u,
+ grad_v,
+ k,
+ grad_k,
+ tabulate,
+ row,
+ col,
+ **_,
):
+ a = sp.Rational(1, 8) + k["k"] * k["k"]
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."
+ diffusion_term = a * tabulate(
+ dot(jac_a_inv.T * grad_u, jac_a_inv.T * grad_v) * jac_a_abs_det
)
- 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 = (
+ approximated_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)
- )
+ )
+
+ if row != col:
+ return diffusion_term
+ else:
+ return diffusion_term + approximated_newton_galerkin_term
- return newton_galerkin_term
return process_integrand(
integrand,
@@ -518,7 +458,7 @@ Note: :math:`a(k) = 1/8 + k^2` is currently hard-coded and the form is intended
geometry,
symbolizer,
blending=blending,
- is_symmetric=False,
+ is_symmetric=True,
docstring=docstring,
fe_coefficients={"k": coefficient_function_space},
)
diff --git a/hog/integrand.py b/hog/integrand.py
index 07313705041256cc2c0ffb8f2f35d62ee448dbbc..2637677ed97a940582f4bb9dc7836300f3da8118 100644
--- a/hog/integrand.py
+++ b/hog/integrand.py
@@ -491,7 +491,7 @@ def process_integrand(
s.grad_v = data.test_shape_grad
s.hessian_v = data.test_shape_hessian
- mat[data.row, data.col] = integrand(**asdict(s))
+ mat[data.row, data.col] = integrand(**asdict(s), row=data.row, col=data.col)
free_symbols_sorted = sorted(list(free_symbols), key=lambda x: str(x))