External operator depending upon the amplitude of temperature flux #45
Description
Hi,
I want to modify the heat tutorial so as to make the external operator ufl.sqrt(ufl.dot(grad(T), grad(T))) rather than
Below is my attempt to modify the tutorial example
from mpi4py import MPI
import numpy as np
import basix
import ufl
import ufl.algorithms
from dolfinx import fem, mesh
from dolfinx_external_operator import (
FEMExternalOperator,
evaluate_external_operators,
evaluate_operands,
replace_external_operators,
)
from ufl import Identity, Measure, TestFunction, TrialFunction, derivative, grad, inner
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.functionspace(domain, ("CG", 1))
DG_vec = fem.functionspace(domain, ("DG", 0, (2,)))
DG = fem.functionspace(domain, ("DG", 0))
T = fem.Function(V)
T.interpolate(lambda x: x[0] ** 2 + x[1])
q_norm = ufl.sqrt(ufl.dot(grad(T), grad(T)))
quadrature_degree = 2
Qe = basix.ufl.quadrature_element(domain.topology.cell_name(), degree=quadrature_degree, value_shape=())
Q = fem.functionspace(domain, Qe)
dx = Measure(
"dx",
metadata={"quadrature_scheme": "default", "quadrature_degree": quadrature_degree},
)
k = FEMExternalOperator(q_norm, function_space=Q)
T_tilde = TestFunction(V)
F = inner(-k * grad(T), grad(T_tilde)) * dx
A = 1.0
B = 1.0
gdim = domain.geometry.dim
def k_impl(q_norm):
# The input `q_norm` is a `np.ndarray` and has the size equal to the number of
# DoFs of the space `DG`.
output = 1.0 / (A + B * q_norm)
# The output must be returned flattened to one dimension
return output.reshape(-1)
def dkdq_norm_impl(q_norm):
return -B * k_impl(q_norm) ** 2
# %%
def k_external(derivatives):
"""Defines behaviour of the external operator and its derivatives.
Args:
derivatives: a tuple of integers representing a multi-index. Each index
is associated with an operand and indicates whether we take derivative
at this operand.
Returns:
a callable Python function.
"""
if derivatives == (0,): # no derivation, the function itself
return k_impl
elif derivatives == (1,): # the derivative with respect to the operand `T`
return dkdq_norm_impl
else:
raise NotImplementedError(f"No external function is defined for the requested derivative {derivatives}.")
k.external_function = k_external
T_hat = TrialFunction(V)
J = derivative(F, T, T_hat)
J_expanded = ufl.algorithms.expand_derivatives(J)
F_replaced, F_external_operators = replace_external_operators(F)
J_replaced, J_external_operators = replace_external_operators(J_expanded)
evaluated_operands = evaluate_operands(F_external_operators)
_ = evaluate_external_operators(F_external_operators, evaluated_operands)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)