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).
@Andrash: 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!