I think the best way to check it would be to reproduce a manufactured solution, or solve a problem where you know the exact solution and then calculate the error.
Hello All, as I understand it, the fix achieved in this thread to insert an external function into the variational formulation is to update the diffu values after each time-step. If I understand it correctly, this amounts to treating the diffu term explicitly.
My question is: Does the delfinx-external-operator (GitHub - a-latyshev/dolfinx-external-operator: Extension of DOLFINx implementing the concept of external operator) allow a python function (taking an ndarray as input and output) to be treated implicitly in time?
As a first try, I attempted to solve the cahn-hilliard equation shown in the dolfinx demo using python written code for the chemical potential. However, the output is not as expected (the solver runs without error, but the output is erroneous). I suspect the values pulled from the python code are not updating each time step (?) Thank you for any insight you can lend on this question and whether dolfinx-external-operator is the right approach to solving the presented issue. Otherwise, is the only way to treat source terms like this explicitly by interpolating the values from a python function at the beginning of each time step. Thanks!
import os
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element, quadrature_element
from dolfinx import default_real_type, log, plot
from dolfinx.fem import Function, functionspace, form, assemble_vector, assemble_matrix
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_unit_interval
from dolfinx.nls.petsc import NewtonSolver
from dolfinx_external_operator import (
FEMExternalOperator,
evaluate_external_operators,
evaluate_operands,
replace_external_operators,
)
from ufl import dx, grad, inner, derivative, Measure
lmbda = 1.0e-02 # surface parameter
dt = 1.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
N= 1024
msh = create_unit_interval(MPI.COMM_WORLD, N)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
V = ME.sub(0)
Q = functionspace(msh, P1)
# Trial and test functions of the space `ME` are now defined:
q, v = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)
# Interpolate initial condition
rng = np.random.default_rng(42)
u.sub(0).interpolate(lambda x: 0.63 + 0.02 * (0.5 - rng.random(x.shape[1])))
u.x.scatter_forward()
mu_mid = (1.0 - theta) * mu0 + theta * mu
quadrature_degree = 1
Qe = quadrature_element(msh.topology.cell_name(), degree=quadrature_degree, value_shape=())
H = functionspace(msh, Qe)
dfdc_ext = FEMExternalOperator(c, function_space=H)
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc_ext, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F = F0 + F1
def dfdc_ext_impl(c):
output = 200 * c * (1 - c)**2 - 200 * c**2 * (1 - c)
return(output.reshape(-1))
def dfdc_ext_impl_jac(c):
output = 200 * (1 - c)**2 - 800 * c * (1 - c) + 200 * c**2
return(output.reshape(-1))
def k_external(derivatives):
if derivatives == (0,): # no derivation, the function itself
return(dfdc_ext_impl)
elif derivatives == (1,): # the derivative with respect to the operand `T`
return(dfdc_ext_impl_jac)
else:
return(NotImplementedError)
dfdc_ext.external_function = k_external
J = derivative(F, u)
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)
F_compiled = form(F_replaced)
J_compiled = form(J_replaced)
b_vector = assemble_vector(F_compiled)
A_matrix = assemble_matrix(J_compiled)
problem = NonlinearProblem(F_compiled, u, J=J_compiled)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2
ksp = solver.krylov_solver
opts = PETSc.Options() # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
sys = PETSc.Sys() # type: ignore
# For factorisation prefer MUMPS, then superlu_dist, then default.
if sys.hasExternalPackage("mumps"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
ksp.setFromOptions()
file = XDMFFile(MPI.COMM_WORLD, "demo_ch/output.xdmf", "w")
file.write_mesh(msh)
# Step in time
t = 0.0
T = 200 * dt
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
V0, dofs = ME.sub(0).collapse()
c = u.sub(0)
u0.x.array[:] = u.x.array
while t < T:
t += dt
r = solver.solve(u)
print(f"Step {int(t / dt)}: num iterations: {r[0]}")
u0.x.array[:] = u.x.array
u.x.scatter_forward()
file.write_function(c, t)
file.close()
Hi Thomas, I just ran your code and the output seems to make sense. Would you mind sharing why do you think the output is erroneous?
Dear Pavan,
Thank you for for your shared interest in the question. The code runs to completion but the output is nonsensical. Something seems to be getting “stuck”. I cleaned up the code a bit and reposted it below with 1) the original formulation of the chemical potential using internal ufl operators and 2) using a new formulation of the chemical potential that uses an external operator; you just need to uncomment/comment the appropriate sections. If you run and visualize the output using the two approaches, you will see that the concentration field differs dramatically.
import os
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from basix.ufl import element, mixed_element, quadrature_element
from dolfinx import default_real_type, log, plot
from dolfinx.fem import Function, functionspace, form, assemble_vector, assemble_matrix
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_unit_interval
from dolfinx.nls.petsc import NewtonSolver
from dolfinx_external_operator import (
FEMExternalOperator,
evaluate_external_operators,
evaluate_operands,
replace_external_operators)
from ufl import grad, inner, derivative, Measure
lmbda = 1.0e-02 # surface parameter
dt = 1.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
N= 1000
msh = create_unit_interval(MPI.COMM_WORLD, N)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
SE = functionspace(msh, P1)
quadrature_degree = 1
QE = quadrature_element(msh.topology.cell_name(), degree=quadrature_degree, value_shape=())
Q = functionspace(msh, QE)
dx = Measure("dx", metadata={"quadrature_scheme": "default", "quadrature_degree": quadrature_degree})
me = ufl.TrialFunction(ME)
(c_hat, mu_hat) = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
rng = np.random.default_rng(42)
u.sub(0).interpolate(lambda x: 0.63 + 0.02 * (0.5 - rng.random(x.shape[1])))
u.x.scatter_forward()
(c, mu) = ufl.split(u)
(c0, mu0) = ufl.split(u0)
mu_mid = (1.0 - theta) * mu0 + theta * mu
################################################################################################
######## ORIGINAL IMPLEMENTATION OF THE CHEMICAL POTENTIAL USING UFL OPERATORS #################
################################################################################################
"""
mu_original = 200 * c * (1 - c)**2 - 200 * c**2 * (1 - c)
F0 = inner(c, c_hat) * dx - inner(c0, c_hat) * dx + dt * inner(grad(mu_mid), grad(c_hat)) * dx
F1 = inner(mu, mu_hat) * dx - inner(mu_original, mu_hat) * dx - lmbda * inner(grad(c), grad(mu_hat)) * dx
F = F0 + F1
problem = NonlinearProblem(F, u)
"""
################################################################################################
######## ATTEMPTED IMPLEMENTATION OF THE CHEMICAL POTENTIAL USING EXTERNAL OPERATOR ############
################################################################################################
mu_external = FEMExternalOperator(c, function_space=Q)
F0 = inner(c, c_hat) * dx - inner(c0, c_hat) * dx + dt * inner(grad(mu_mid), grad(c_hat)) * dx
F1 = inner(mu, mu_hat) * dx - inner(mu_external, mu_hat) * dx - lmbda * inner(grad(c), grad(mu_hat)) * dx
F = F0 + F1
def mu_external_impl(c):
output = 200 * c * (1 - c)**2 - 200 * c**2 * (1 - c)
return(output.reshape(-1))
def dmudc_external_impl(c):
output = 200 * (1 - c)**2 - 800 * c * (1 - c) + 200 * c**2
return(output.reshape(-1))
def mu_external_function(derivatives):
if derivatives == (0,):
return(mu_external_impl)
elif derivatives == (1,):
return(dmudc_external_impl)
else:
return(NotImplementedError)
mu_external.external_function = mu_external_function
J = derivative(F, u, me)
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)
F_compiled = form(F_replaced)
J_compiled = form(J_replaced)
problem = NonlinearProblem(F_compiled, u, J=J_compiled)
################################################################################################
################################################################################################
################################################################################################
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2
ksp = solver.krylov_solver
opts = PETSc.Options() # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
sys = PETSc.Sys() # type: ignore
# For factorisation prefer MUMPS, then superlu_dist, then default.
if sys.hasExternalPackage("mumps"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
ksp.setFromOptions()
file = XDMFFile(MPI.COMM_WORLD, "demo_ch/output.xdmf", "w")
file.write_mesh(msh)
# Step in time
t = 0.0
T = 100 * dt
u0.x.array[:] = u.x.array
while t < T:
t += dt
r = solver.solve(u)
print(f"Step {int(t / dt)}: num iterations: {r[0]}")
file.write_function(u.sub(0), t)
u0.x.array[:] = u.x.array[:]
file.close()
Screenshots of the output running either option of the code.
I took a quick look at your code examples.
The ‘issue’ with your code is that you are using the built-in Newton Solver which does not evaluate the external operators prior to assembling the Newton step system.
You can see an example around here: dolfinx-external-operator/doc/demo/demo_plasticity_von_mises.py at main · a-latyshev/dolfinx-external-operator · GitHub
I think this solution is not hugely satisfactory though, and we should improve how we show the use and implementation of non-linear solvers in our examples.
2 Likes
I personally believe that the external operator should be part of FFCx/DOLFINx.
Very little code injection is neeeded, only a update operator coefficients at the same time as pack coefficients.
3 Likes
Dear @Thomas_Petersen,
Thanks for your interest in the dolfinx-external-operator project.
I would like to complement the @jackhale reply here above in case if you need a solution right now.
In the tutorials: Plasticity of von Mises — External operators within FEniCSx/DOLFINx, we implemented the Newton solver manually via Python to demonstrate how the framework works. Inside the Newton solver, we call
evaluated_operands = evaluate_operands(F_external_operators)
to update the input for external operators (which may be an UFL expression) and then call
((_, sigma_new, dp_new),) = evaluate_external_operators(J_external_operators, evaluated_operands)
sigma.ref_coefficient.x.array[:] = sigma_new
that actually updates the operators. So, instead of using NonlinearProblem and NewtonSolver, you could manually implement the Newton method and control every step of your modelling.
If the use of PETSc solvers is preferable I may propose the following “hack”. I suggest changing the method F of the class NonlinearProblem dolfinx/python/dolfinx/fem/petsc.py at main · FEniCS/dolfinx · GitHub by putting the update of the external operators inside. For example here, I put a callable field self.inside_Newton() convex-plasticity/src/plasticity_framework.py at main · a-latyshev/convex-plasticity · GitHub that can wrap additional computations which should be made on every iteration of the Newton solver. In this notebook, you can find an application of such a hack to a convex plasticity problem within PETSc quasi-Newton SNES: convex-plasticity/tutorials/convex_plasticity/convex_plasticity_SNESQN.ipynb at main · a-latyshev/convex-plasticity · GitHub.
def inside_Newton():
fs.interpolate_quadrature(eps_vec(Du), deps) # eps_xy * sqrt(2.)!
for q in range(N_patches):
return_mapping.deps.value[:] = deps_values[q,:].T
return_mapping.sig_old.value[:] = sig_old_values[q,:].T
return_mapping.p_old.value = p_old_values[q,:]
return_mapping.solve(**conic_solver_params)
sig_values[q,:] = return_mapping.sig.value[:].T
p_values[q,:] = return_mapping.p.value
if residue_size != 0: #how to improve ?
return_mapping_residue.deps.value[:] = deps_values_residue[0,:].T
return_mapping_residue.sig_old.value[:] = sig_old_values_residue[0,:].T
return_mapping_residue.p_old.value = p_old_values_residue[0,:]
return_mapping_residue.solve(**conic_solver_params)
sig_values_residue[0,:] = return_mapping_residue.sig.value[:].T
p_values_residue[0,:] = return_mapping_residue.p.value
residual = ufl.inner(as_3D_tensor(sig), eps(u_))*dx - F_ext(u_)
J = ufl.derivative(ufl.inner(sigma(eps(Du)), eps(u_))*dx, Du, v_)
snes_problem = pf.SNESProblem(residual, Du, J, bcs, petsc_options=petsc_options_SNESQN, inside_Newton=inside_Newton)
I hope this will help for now 
Let us know if you have more questions concerning external operators within DOLFINx.
All the best
1 Like
Thank you Andrash and Jack,
Both of your comments are extremely helpful. I implemented the customized Newton solver by iterating the LinearProblem and the simulation runs as expected now. The updated code for the interested reader is below (note that I had to copy the solver.py file from the plasticity von Mises demo into my current directory).
@a-latyshev: I will be very curious to try your hack of placing the updates to the external operators as a method in the source code for the NonlinearProblem.
Thanks again!
import numpy as np
from mpi4py import MPI
import ufl
from basix.ufl import element, mixed_element, quadrature_element
# from dolfinx import default_real_type, default_scalar_type, log, plot
from dolfinx.fem import Function, functionspace
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_interval
from dolfinx_external_operator import (
FEMExternalOperator,
evaluate_external_operators,
evaluate_operands,
replace_external_operators)
from solvers import LinearProblem
from ufl import grad, inner, derivative, Measure
lmbda = 1.0e-02 # surface parameter
dt = 1.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
N= 1000
msh = create_unit_interval(MPI.COMM_WORLD, N)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
SE = functionspace(msh, P1)
quadrature_degree = 1
QE = quadrature_element(msh.topology.cell_name(), degree=quadrature_degree, value_shape=())
Q = functionspace(msh, QE)
dx = Measure("dx", metadata={"quadrature_scheme": "default", "quadrature_degree": quadrature_degree})
me = ufl.TrialFunction(ME)
(c_hat, mu_hat) = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
du = Function(ME) # change in solution
rng = np.random.default_rng(42)
u.sub(0).interpolate(lambda x: 0.4 + 0.02 * (0.5 - rng.random(x.shape[1])))
u.x.scatter_forward()
(c, mu) = ufl.split(u)
(c0, mu0) = ufl.split(u0)
mu_mid = (1.0 - theta) * mu0 + theta * mu
mu_external = FEMExternalOperator(c, function_space=Q)
F0 = inner(c, c_hat) * dx - inner(c0, c_hat) * dx + dt * inner(grad(mu_mid), grad(c_hat)) * dx
F1 = inner(mu, mu_hat) * dx - inner(mu_external, mu_hat) * dx - lmbda * inner(grad(c), grad(mu_hat)) * dx
F = F0 + F1
def mu_external_impl(c):
output = 200 * c * (1 - c)**2 - 200 * c**2 * (1 - c)
return(output.reshape(-1))
def dmudc_external_impl(c):
output = 200 * (1 - c)**2 - 800 * c * (1 - c) + 200 * c**2
return(output.reshape(-1))
def mu_external_function(derivatives):
if derivatives == (0,):
return(mu_external_impl)
elif derivatives == (1,):
return(dmudc_external_impl)
else:
return(NotImplementedError)
mu_external.external_function = mu_external_function
J = derivative(F, u, me)
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)
problem = LinearProblem(J_replaced, -F_replaced, du)
# Using customized, self-built solver
# Step in time
t = 0.0
T = 100 * dt
max_iter = 100
rel_tol = 1e-8
u0.x.array[:] = u.x.array
file = XDMFFile(MPI.COMM_WORLD, "demo_ch/output.xdmf", "w")
file.write_mesh(msh)
while t < T:
t += dt
problem.assemble_vector()
residual_0 = problem.b.norm()
if MPI.COMM_WORLD.rank == 0:
print(f"Step {int(t / dt)}:")
for iteration in range(max_iter):
if iteration > 0:
if rel_error < rel_tol:
break
problem.assemble_matrix()
problem.solve(du)
du.x.scatter_forward()
u.x.petsc_vec.axpy(1.0, du.x.petsc_vec)
evaluated_operands = evaluate_operands(F_external_operators)
_ = evaluate_external_operators(F_external_operators, evaluated_operands)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)
problem.assemble_vector()
residual = problem.b.norm()
rel_error = np.abs(residual / residual_0)
# residual_0 = residual
if MPI.COMM_WORLD.rank == 0:
print(f" it# {iteration} relative error: {rel_error}")
file.write_function(u.sub(0), t)
u0.x.array[:] = u.x.array[:]
file.close()
2 Likes
@Thomas_Petersen You are a star! The code you wrote is super helpful! It seems that the LinearSolver class from the demo (convex-plasticity/src/plasticity_framework.py at main · a-latyshev/convex-plasticity · GitHub) needs a little bit of tweaking to make it compatible with dolfinx 0.9.0.
It does seem the current implementation for interfacing external operators with coupled nonlinear problems is a bit clunky, but its a good start!
We’ve updated the von Mises plasticity demo to show use of the built-in NewtonSolver for problems with internal state.
2 Likes
I just saw this. Thank you for demonstrating how to implement the external-operator into the built-in NewtonSolver!
Dear @Thomas_Petersen and @pavaninguva,
I wanted to inform you that we have implemented an interface for Newton solvers that facilitates the use of external operators with two of the most common solvers: the one from DOLFINx and PETSc.SNES. This is achieved via the NonlinearProblemWithCallback and PETScNonlinearProblem classes, respectively (source). Both classes require defining an external_callback, where external operator evaluations will be executed.
You can find application examples in our tutorials: NonlinearProblemWithCallback is used in von Mises plasticity and PETScNonlinearProblem is applied in Mohr-Coulomb plasticity.
I hope this simplifies the interaction with external operators in your nonlinear problems.
Please feel free to reach out if you have any questions!
The demo refers to solvers,
from solvers import NonlinearProblemWithCallback
but it’s not clear what the solvers are.
I guess you don’t mean solvers · PyPI
It’s the file solvers.py in the ˋdemo` folder.
I see it. Thanks.
Hello a-latyshev,
This is amaizing, thank you! I still have the work-around solution using the external_operators implemented in my code, but am certain to start using the new NonlinearProblemWithCallback class in the next iteration of the code.
TP
Hi all, I came across the dolfinx-external-operator module recently and was trying to use the NonlinearProblemWithCallback class that @a-latyshev said that they implemented a little while back. I am trying to use it for a case where multiple material properties in my system’s residual would need to be implemented as FEMExternalOperator objects. Specifically, I am coupling together 1D instances of heat and momentum equations, where my independent variables are v_theta and T, and then both my viscosity and thermal conductivity are functions of T.
Going off the workflow shown in the demo pages of the documentation, I first define
def make_external_from_scalar_fun(python_function, h=1e-9):
"""Take a python function and create the external function for FEMExternalOperator."""
print(f"DEBUG: Creating dispatcher for '{python_function}'", flush=True) # Check creation
def f(T):
print(f"DEBUG: Calling function {python_function}"),
return python_function(T).ravel()
def f_prime(T):
print(f"DEBUG: Calling derivative of {python_function}")
return ((python_function(T + h) - python_function(T - h)) / (2.0 * h)).ravel()
def dispatcher(derivatives):
print(f"DEBUG: Dispatcher for '{python_function}' called with derivatives={derivatives}")
if derivatives == (0,):
return f
elif derivatives == (1,):
return f_prime
else:
raise NotImplementedError
return dispatcher
where I use a finite difference approximation to compute the derivative, and the print() statements are there purely for debugging purposes. As mentioned above, I then define my material property functions—i.e. viscosity and thermal conductivity—as functions of the temperature:
def mu(T):
# Viscosity fit for hydrogen
return -8.22278*10**(-6) + 5.03195*10**(-7) * T**0.58505
def kappa(T):
# Thermal conductivity fit for hydrogen
return -0.475048 + 0.0191871 * T**0.569845
Lastly, I create the actual dolfinx framework and attempt to solve the problem using an instance of NonlinearProblemWithCallback:
mesh_size = 1000
r_in, r_out = 0.0075, 0.06
Itot, L, B = 200.0, 0.25, 0.5
# Construction of domain and function spaces
domain = mesh.create_interval(MPI.COMM_WORLD, mesh_size, points=(r_in, r_out))
V = fem.functionspace(domain, element("Lagrange", domain.topology.cell_name(), 1, shape=(2,)))
u = fem.Function(V)
v_theta, T = ufl.split(u)
v_theta_test, T_test = ufl.TestFunctions(V)
# Quadrature space where the coefficients will live
quad_deg = 2
Qe = basix.ufl.quadrature_element(domain.topology.cell_name(), degree=quad_deg, value_shape=())
Q = fem.functionspace(domain, Qe)
dx = ufl.Measure("dx", domain=domain, metadata={"quadrature_degree": quad_deg})
# Build the material property objects
kappa_sym = FEMExternalOperator(T, function_space=Q)
mu_sym = FEMExternalOperator(T, function_space=Q)
sigma_perp_sym = fem.Constant(domain, 10.77) # assumed constant
# Weak form
x = ufl.SpatialCoordinate(domain)
r = x[0]
j = - Itot / (2 * np.pi * r * L)
F = ufl.inner(r * kappa_sym * T.dx(0), T_test.dx(0)) * dx
F += - ufl.inner(mu_sym * (r * (v_theta / r).dx(0))**2, T_test) * r * dx
F += - ufl.inner(j**2 / sigma_perp_sym, T_test) * r * dx
F += - ufl.inner((mu_sym / r) * (r * v_theta).dx(0), (v_theta_test * r).dx(0)) * dx
F += - ufl.inner(j * B, v_theta_test) * r * dx
# Assign the external functions to the symbolic objects
kappa_sym.external_function = make_external_from_scalar_fun(kappa)
mu_sym.external_function = make_external_from_scalar_fun(mu)
u_trial = ufl.TrialFunction(V)
J = ufl.derivative(F, u, u_trial)
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)
def constitutive_update():
all_ops = J_external_operators + F_external_operators
evaluated_operands = evaluate_operands(all_ops)
new_values = evaluate_external_operators(all_ops, evaluated_operands)
for op, arr in zip(all_ops, new_values):
op.ref_coefficient.x.array[:] = arr
problem = NonlinearProblemWithCallback(F_replaced, u, J=J_replaced, external_callback=constitutive_update)
solver = NewtonSolver(domain.comm, problem)
solver.solve(u)
However, it seems like only the first material property instance defined in the residual, in this case kappa, is actually updated as the solver proceeds, as this is what the print statements from make_external_from_scalar_fun() look like when I try to run the solver:
DEBUG: Creating dispatcher for '<function kappa at 0x7d8419afaca0>'
DEBUG: Creating dispatcher for '<function mu at 0x7d8419afade0>'
DEBUG: Dispatcher for '<function kappa at 0x7d8419afaca0>' called with derivatives=(0,)
DEBUG: Calling function <function kappa at 0x7d8419afaca0>
DEBUG: Dispatcher for '<function kappa at 0x7d8419afaca0>' called with derivatives=(1,)
DEBUG: Calling derivative of <function kappa at 0x7d8419afaca0>
DEBUG: Dispatcher for '<function kappa at 0x7d8419afaca0>' called with derivatives=(0,)
DEBUG: Calling function <function kappa at 0x7d8419afaca0>
DEBUG: Dispatcher for '<function kappa at 0x7d8419afaca0>' called with derivatives=(0,)
...
As you can see, the dispatcher for mu is created, but is never updated during the solver iterations. Am I doing something incorrectly here with my setup?
I am rather novice to the dolfinx framework in general so it’s entirely possible I’m just setting up the problem poorly 
Thanks in advance!
For those who will face the same issue, it was fixed (see PR: Fix external operators detection by a-latyshev · Pull Request #25 · a-latyshev/dolfinx-external-operator · GitHub) and tested for dolfinx v0.9.0 and the current nightly.
3 Likes
Hi @a-latyshev Thanks so much for all the awesome work!
I do have a follow up question re possibly extending things to cases where my external operator has the form f(a,b), i.e., it depends on two variables.
Using your example for nonlinear diffusion from the external operators repo, but with a modified k(T1, T2), I am not sure how this would look like regarding the code because I would need to specify the partial derivatives dk/dT1, dk/dT2 in the k_external function. Would appreciate your help to see how that can be implemented (my guess is that the indexing in the if/elif statements might change).
Copied the k_external function for ease of reference!
def k_external(derivatives):
if derivatives == (0,): # no derivation, the function itself
return k_impl
elif derivatives == (1,): # the derivative with respect to the operand `T`
return dkdT_impl
That’s easy, just follow “More complex case” within the nonlinear heat tutorial: More complex case — External operators within FEniCSx/DOLFINx. It has all you need.
In short, the case of multiple operands with high-order differentiation is governed by the multi-index, which is the tuple of integers in Python.
For example, in the tutorial mentioned above you will find the following external function
def q_external(derivatives):
if derivatives == (0, 0):
return q_impl
elif derivatives == (1, 0):
return dqdT_impl
elif derivatives == (0, 1):
return dqdsigma_impl
where the multi-index (0, 1) is associated with the partial derivative with respect to the second operand.
1 Like