diff --git a/hog/forms.py b/hog/forms.py
index 4a9439ad920f1796126b7f3085311ebe056adf36..85e93c2b02d67bbe5c189a48733b55c542675718 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -682,6 +682,144 @@ where
fe_coefficients={"mu": coefficient_function_space},
)
+def full_stokes_viscoplastic(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_trial: int = 0,
+ component_test: int = 0,
+ variable_viscosity: bool = True,
+ velocity_function_space: Optional[FunctionSpace] = None,
+ coefficient_function_space: Optional[FunctionSpace] = None,
+) -> Form:
+ docstring = f"""
+Implements the fully coupled viscous operator of the Stokes problem.
+The latter is the extension of the Epsilon operator to the case where
+the velocity field need not be divergence-free. This is e.g. the case
+in the (truncated) anelastic liquid approximation of mantle convection.
+
+The strong representation of the operator is given by:
+
+ - div[ μ (grad(u)+grad(u)ᵀ) ] + 2/3 grad[ μ div(u) ]
+
+Note that the factor 2/3 means that for 2D this is the pseudo-3D form
+of the operator.
+
+Component trial: {component_trial}
+Component test: {component_test}
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (vectorial space: {trial})
+ v: test function (vectorial space: {test})
+ μ: coefficient (scalar space: {coefficient_function_space})
+
+ ∫ μ {{ ( 2 ε(u) : ε(v) ) - (2/3) [ ( ∇ · u ) · ( ∇ · v ) ] }}
+
+where
+
+ ε(w) := (1/2) (∇w + (∇w)ᵀ)
+"""
+
+ from hog.recipes.integrands.volume.full_stokes_viscoplastic import integrand
+
+ return process_integrand(
+ integrand,
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ is_symmetric=trial == test,
+ docstring=docstring,
+ fe_coefficients={"mu_lin": coefficient_function_space, "ux": velocity_function_space, "uy": velocity_function_space},
+ )
+
+def evaluate_viscosity_viscoplastic(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_trial: int = 0,
+ component_test: int = 0,
+ variable_viscosity: bool = True,
+ velocity_function_space: Optional[FunctionSpace] = None,
+ coefficient_function_space: Optional[FunctionSpace] = None,
+) -> Form:
+ docstring = f"""
+Implements the fully coupled viscous operator of the Stokes problem.
+The latter is the extension of the Epsilon operator to the case where
+the velocity field need not be divergence-free. This is e.g. the case
+in the (truncated) anelastic liquid approximation of mantle convection.
+
+The strong representation of the operator is given by:
+
+ - div[ μ (grad(u)+grad(u)ᵀ) ] + 2/3 grad[ μ div(u) ]
+
+Note that the factor 2/3 means that for 2D this is the pseudo-3D form
+of the operator.
+
+Component trial: {component_trial}
+Component test: {component_test}
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (vectorial space: {trial})
+ v: test function (vectorial space: {test})
+ μ: coefficient (scalar space: {coefficient_function_space})
+
+ ∫ μ {{ ( 2 ε(u) : ε(v) ) - (2/3) [ ( ∇ · u ) · ( ∇ · v ) ] }}
+
+where
+
+ ε(w) := (1/2) (∇w + (∇w)ᵀ)
+"""
+
+ from hog.recipes.expressions.viscoplastic_rheology import expression
+
+ return process_integrand(
+ expression,
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ is_symmetric=trial == test,
+ docstring=docstring,
+ fe_coefficients={"mu_lin": coefficient_function_space, "ux": velocity_function_space, "uy": velocity_function_space},
+ dont_integrate=True
+ )
+
+def all_ones(
+ trial: TrialSpace,
+ test: TestSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ velocity_function_space: Optional[FunctionSpace] = None,
+ coefficient_function_space: Optional[FunctionSpace] = None,
+) -> Form:
+ docstring = f"""
+"""
+
+ from hog.recipes.expressions.ones import expression
+
+ return process_integrand(
+ expression,
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ is_symmetric=trial == test,
+ docstring=docstring,
+ dont_integrate=True
+ )
def shear_heating(
trial: TrialSpace,
diff --git a/hog/integrand.py b/hog/integrand.py
index 07313705041256cc2c0ffb8f2f35d62ee448dbbc..5e86bbe9621b0c721cabb83ea08d51816e59f98f 100644
--- a/hog/integrand.py
+++ b/hog/integrand.py
@@ -88,11 +88,15 @@ class Form:
symmetric: bool
free_symbols: List[sp.Symbol] = field(default_factory=lambda: list())
docstring: str = ""
+ dont_integrate: bool = False
def integrate(self, quad: Quadrature, symbolizer: Symbolizer) -> sp.Matrix:
"""Integrates the form using the passed quadrature directly, i.e. without tabulations or loops."""
mat = self.tabulation.inline_tables(self.integrand)
+ if self.dont_integrate:
+ return
+
for row in range(mat.rows):
for col in range(mat.cols):
if self.symmetric and row > col:
@@ -163,6 +167,13 @@ class IntegrandSymbols:
# The keys are specified by the strings that are passed to process_integrand.
grad_k: Dict[str, sp.Matrix] | None = None
+ k_dof_symbols: Dict[str, List[sp.Symbol]] | None = None
+
+ grad_k_dof_symbols: Dict[str, List[sp.Symbol]] | None = None
+
+ k_dof_symbol: Dict[str, sp.Symbol] | None = None
+ grad_k_dof_symbol: Dict[str, sp.Symbol] | None = None
+
# The geometry of the volume element.
volume_geometry: ElementGeometry | None = None
@@ -233,6 +244,7 @@ def process_integrand(
is_symmetric: bool = False,
rot_type: RotationType = RotationType.NO_ROTATION,
docstring: str = "",
+ dont_integrate: bool = False,
) -> Form:
"""
Constructs an element matrix (:class:`~Form` object) from an integrand.
@@ -293,7 +305,7 @@ def process_integrand(
finite-element function coefficients, they will be available to the callable as `k`
supply None as the FunctionSpace for a std::function-type coeff (only works for old forms)
:param is_symmetric: whether the bilinear form is symmetric - this is exploited by the generator
- :param rot_type: whether the operator has to be wrapped with rotation matrix and the type of rotation that needs
+ :param rot_type: whether the operator has to be wrapped with rotation matrix and the type of rotation that needs
to be applied, only applicable for Vectorial spaces
:param docstring: documentation of the integrand/bilinear form - will end up in the docstring of the generated code
"""
@@ -387,6 +399,13 @@ def process_integrand(
s.k = dict()
s.grad_k = dict()
+
+ s.k_dof_symbols = dict()
+ s.grad_k_dof_symbols = dict()
+
+ s.k_dof_symbol = dict()
+ s.grad_k_dof_symbol = dict()
+
for name, coefficient_function_space in fe_coefficients_modified.items():
if coefficient_function_space is None:
k = scalar_space_dependent_coefficient(
@@ -400,15 +419,15 @@ def process_integrand(
symbolizer,
domain="reference",
function_id=name,
- basis_eval=tabulation.register_phi_evals(
- coefficient_function_space.shape(volume_geometry)
- ),
+ # basis_eval=tabulation.register_phi_evals(
+ # coefficient_function_space.shape(volume_geometry)
+ # ),
)
if isinstance(k, sp.Matrix) and k.shape == (1, 1):
k = k[0, 0]
- grad_k, _ = fem_function_gradient_on_element(
+ grad_k, grad_dof_symbols = fem_function_gradient_on_element(
coefficient_function_space,
volume_geometry,
symbolizer,
@@ -419,6 +438,9 @@ def process_integrand(
s.k[name] = k
s.grad_k[name] = grad_k
+ s.k_dof_symbols[name] = dof_symbols
+ s.grad_k_dof_symbols[name] = grad_dof_symbols
+
##############################
# Jacobians and determinants #
##############################
@@ -482,16 +504,46 @@ def process_integrand(
mat = create_empty_element_matrix(trial, test, volume_geometry)
it = element_matrix_iterator(trial, test, volume_geometry)
- for data in it:
- s.u = data.trial_shape
- s.grad_u = data.trial_shape_grad
- s.hessian_u = data.trial_shape_hessian
+ # dont_integrate option only works well for P1
+ ref_coords = [(0, 0), (1, 0), (0, 1)]
+
+ if dont_integrate:
+ for data in it:
+ s.u = data.trial_shape
+ s.grad_u = data.trial_shape_grad
+ s.hessian_u = data.trial_shape_hessian
+
+ s.v = data.test_shape
+ s.grad_v = data.test_shape_grad
+ s.hessian_v = data.test_shape_hessian
+
+ if data.row == data.col:
+ for name in s.k_dof_symbols.keys():
+ s.k_dof_symbol[name] = s.k[name].subs(
+ [
+ (sp.Symbol("x_ref_0"), ref_coords[data.row][0]),
+ (sp.Symbol("x_ref_1"), ref_coords[data.row][1]),
+ ]
+ )
+ s.grad_k_dof_symbol[name] = s.grad_k[name].subs(
+ [
+ (sp.Symbol("x_ref_0"), ref_coords[data.row][0]),
+ (sp.Symbol("x_ref_1"), ref_coords[data.row][1]),
+ ]
+ )
+
+ mat[data.row, data.col] = integrand(**asdict(s))
+ else:
+ for data in it:
+ s.u = data.trial_shape
+ s.grad_u = data.trial_shape_grad
+ s.hessian_u = data.trial_shape_hessian
- s.v = data.test_shape
- s.grad_v = data.test_shape_grad
- s.hessian_v = data.test_shape_hessian
+ s.v = data.test_shape
+ 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))
free_symbols_sorted = sorted(list(free_symbols), key=lambda x: str(x))
@@ -603,4 +655,5 @@ where
symmetric=is_symmetric,
free_symbols=free_symbols_sorted,
docstring=docstring,
+ dont_integrate=dont_integrate,
)
diff --git a/hog/operator_generation/operators.py b/hog/operator_generation/operators.py
index 18dbfede7ce05b0f9b57714e3c3203fd8f7a8918..d666edddbbf1c1bc543e86d58b47b96ce61f3546 100644
--- a/hog/operator_generation/operators.py
+++ b/hog/operator_generation/operators.py
@@ -373,29 +373,30 @@ class HyTeGElementwiseOperator:
else:
mat = form.tabulation.inline_tables(mat)
- if Opts.QUADLOOPS in optimizations:
- quad_loop = QuadLoop(
- self.symbolizer,
- quad,
- mat,
- self._type_descriptor,
- form.symmetric,
- blending,
- )
- mat = quad_loop.mat
- else:
- with TimedLogger(
- f"integrating {mat.shape[0] * mat.shape[1]} expressions",
- level=logging.DEBUG,
- ):
- for row in range(mat.rows):
- for col in range(mat.cols):
- if form.symmetric and row > col:
- mat[row, col] = mat[col, row]
- else:
- mat[row, col] = quad.integrate(
- mat[row, col], self.symbolizer, blending
- )
+ if not form.dont_integrate:
+ if Opts.QUADLOOPS in optimizations:
+ quad_loop = QuadLoop(
+ self.symbolizer,
+ quad,
+ mat,
+ self._type_descriptor,
+ form.symmetric,
+ blending,
+ )
+ mat = quad_loop.mat
+ else:
+ with TimedLogger(
+ f"integrating {mat.shape[0] * mat.shape[1]} expressions",
+ level=logging.DEBUG,
+ ):
+ for row in range(mat.rows):
+ for col in range(mat.cols):
+ if form.symmetric and row > col:
+ mat[row, col] = mat[col, row]
+ else:
+ mat[row, col] = quad.integrate(
+ mat[row, col], self.symbolizer, blending
+ )
if volume_geometry.space_dimension not in self.integration_infos:
self.integration_infos[volume_geometry.space_dimension] = []
diff --git a/hog/recipes/expressions/ones.py b/hog/recipes/expressions/ones.py
new file mode 100644
index 0000000000000000000000000000000000000000..cdc9859da99015adf40688f2009530f5c31d85d5
--- /dev/null
+++ b/hog/recipes/expressions/ones.py
@@ -0,0 +1,34 @@
+# 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 expression(
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ jac_b_inv,
+ jac_b_abs_det,
+ grad_u,
+ grad_v,
+ k_dof_symbol,
+ grad_k_dof_symbol,
+ scalars,
+ tabulate,
+ **_,
+):
+
+ return 1.0
diff --git a/hog/recipes/expressions/viscoplastic_rheology.py b/hog/recipes/expressions/viscoplastic_rheology.py
new file mode 100644
index 0000000000000000000000000000000000000000..b19614ac1e80311647526d575cc9176fa83e2742
--- /dev/null
+++ b/hog/recipes/expressions/viscoplastic_rheology.py
@@ -0,0 +1,68 @@
+# 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 expression(
+ *,
+ x,
+ jac_a_inv,
+ jac_a_abs_det,
+ jac_b_inv,
+ jac_b_abs_det,
+ grad_u,
+ grad_v,
+ k,
+ grad_k,
+ k_dof_symbol,
+ grad_k_dof_symbol,
+ scalars,
+ tabulate,
+ **_,
+):
+
+ grad_u_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)
+ grad_v_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)
+
+ def symm_grad(w):
+ return 0.5 * (w + w.T)
+
+ mu_lin = k_dof_symbol["mu_lin"]
+
+ grad_ux = grad_k_dof_symbol["ux"]
+ grad_uy = grad_k_dof_symbol["uy"]
+
+ print("grad_ux = ", grad_ux)
+ # print("x = ", x)
+
+ uxx = (jac_b_inv.T * jac_a_inv.T * grad_ux)[0]
+ uxy = (jac_b_inv.T * jac_a_inv.T * grad_ux)[1]
+
+ uyx = (jac_b_inv.T * jac_a_inv.T * grad_uy)[0]
+ uyy = (jac_b_inv.T * jac_a_inv.T * grad_uy)[1]
+
+ print("uxx = ", uxx)
+
+ mu_star = scalars("mu_star")
+ sigma_y = scalars("sigma_y")
+
+ strain = uxx * uxx + 0.5 * (uxy + uyx) * (uxy + uyx) + uyy * uyy
+
+ mu_nonlin = mu_star + (sigma_y / (sp.sqrt(strain)))
+
+ mu = 2.0 / ((1.0 / mu_lin) + (1.0 / mu_nonlin))
+
+ return mu
diff --git a/hog/recipes/integrands/volume/full_stokes_viscoplastic.py b/hog/recipes/integrands/volume/full_stokes_viscoplastic.py
new file mode 100644
index 0000000000000000000000000000000000000000..e1bf51663575e10ae9ac2b206b0caf95c494a7e7
--- /dev/null
+++ b/hog/recipes/integrands/volume/full_stokes_viscoplastic.py
@@ -0,0 +1,89 @@
+# 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(
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ jac_b_inv,
+ jac_b_abs_det,
+ grad_u,
+ grad_v,
+ k,
+ grad_k,
+ scalars,
+ tabulate,
+ **_,
+):
+
+ grad_u_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)
+ grad_v_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)
+
+ def symm_grad(w):
+ return 0.5 * (w + w.T)
+
+ symm_grad_u = symm_grad(grad_u_chain)
+ symm_grad_v = symm_grad(grad_v_chain)
+
+ div_u = (jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)).trace()
+ div_v = (jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)).trace()
+
+ mu_lin = k["mu_lin"]
+
+ grad_ux = grad_k["ux"]
+ grad_uy = grad_k["uy"]
+
+ uxx = (jac_b_inv.T * jac_a_inv.T * grad_ux)[0]
+ uxy = (jac_b_inv.T * jac_a_inv.T * grad_ux)[1]
+
+ uyx = (jac_b_inv.T * jac_a_inv.T * grad_uy)[0]
+ uyy = (jac_b_inv.T * jac_a_inv.T * grad_uy)[1]
+
+ mu_star = scalars("mu_star")
+ sigma_y = scalars("sigma_y")
+
+ strain = uxx * uxx + 0.5 * (uxy + uyx) * (uxy + uyx) + uyy * uyy
+
+ mu_nonlin = mu_star + (sigma_y / (1e-20 + sp.sqrt(strain)))
+
+ mu = 2.0 / ((1.0 / mu_lin) + (1.0 / mu_nonlin))
+
+ return mu * (
+ (
+ double_contraction(2 * symm_grad_u, symm_grad_v)
+ * tabulate(jac_a_abs_det)
+ * jac_b_abs_det
+ )
+ - (2.0 / 3.0) * div_u * div_v * tabulate(jac_a_abs_det) * jac_b_abs_det
+ )
+
+# uxx = (jac_affine_inv.T * grad_ux)[0]
+# uxy = (jac_affine_inv.T * grad_ux)[1]
+
+# uyx = (jac_affine_inv.T * grad_uy)[0]
+# uyy = (jac_affine_inv.T * grad_uy)[1]
+
+# etaStar = 0.001
+# sigmaY = 1.0
+
+# epsContraction = uxx * uxx + 0.5 * (uxy + uyx) * (uxy + uyx) + uyy * uyy
+
+# etaNonlin = etaStar + (sigmaY / (sp.sqrt(epsContraction)))
+
+# eta = 2.0 / ((1.0 / etaLin[0]) + (1.0 / etaNonlin))
\ No newline at end of file