Hi Anton,
One possible solution would be to write your own Newton solver, in which you interpolate the material properties at each Newton iteration based on the current value of uh. See below, where I have adapted this solver to DOLFINx. On my machine, the same max temperature is achieved by both the “polynomial” and the “interpolation” method, although the number of Newton steps is different.
Note that this only works in serial, however.
import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp
mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [50, 50])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
# Define initial conditions
u_n_p = dolfinx.Function(V)
with u_n_p.vector.localForm() as loc:
loc.set(1000)
u_n_i = dolfinx.Function(V)
with u_n_i.vector.localForm() as loc:
loc.set(1000)
#Define solution variable
uh_p = dolfinx.Function(V)
uh_i = dolfinx.Function(V)
with uh_p.vector.localForm() as loc:
loc.set(1000)
with uh_i.vector.localForm() as loc:
loc.set(1000)
v = ufl.TestFunction(V)
f = dolfinx.Constant(mesh, 0)
## Material properties
# Definition using polynom
def heat_cond_polynom(u):
return 0.00004*u + 0.02
# Definition using interpolation
def heat_cond_interpolation(u):
heat_conduction_data = sp.interp1d([0,1000],[0.02,0.06])
return heat_conduction_data(u)
### NewtonSolver implementation adapted from javjiv4
# https://fenicsproject.discourse.group/t/custom-newton-solver-problem-with-dirichlet-conditions/2329/5
def newtonSolve(G, dG, u, bcs, q, q_interp):
"""
Solves the nonlinear system of equations G(u; v) == 0 using Newton's method:
G'(u; du, v) = -G(u; v)
where the nonlinearity occurs due to the coefficient `q` that depends on
`u`.
Parameters
==========
`G` :: semilinear form
`dG` :: Jacobian of G
`u` :: `Function` to be solved for
`bcs` :: `DirichletBC`s for the problem
`q` :: `Function` or list of functions that occurs as a coefficient(s) in
`G` and `dG`, which depends on `u` but has to be evaluated
numerically by interpolation
`q_interp` :: Python function or list of functions of a vector that
interpolates for q
"""
if isinstance(q, list):
if isinstance(q_interp, list):
if not len(q) == len(q_interp):
raise ValueError("q and q_interp must have same len()")
else:
raise ValueError("Mismatched types for q and q_interp")
elif isinstance(q, dolfinx.Function):
if not callable(q_interp):
raise ValueError("If q is a Function, q_interp must be callable")
else:
raise ValueError("q must be a Function or list of Functions")
tol = 1e-3
E = 1.
it = 0
max_it = 10
for bc in bcs:
bc.apply(u.vector())
bc.homogenize()
while E > tol and it < max_it:
# Interpolate material properties using current value of `u`
if isinstance(q, list):
for qi, q_int_i in zip(q, q_interp):
qi.vector.array[:] = q_int_i(np.array(u.vector.array))
else:
q.vector.array[:] = q_interp(np.array(u.vector.array))
x = dolfinx.Function(u.function_space)
problem = dolfinx.fem.LinearProblem(dG, -G, bcs, x)
problem.solve()
error = (x)**2 * ufl.dx
E = np.sqrt(dolfinx.fem.assemble_scalar(error))
# print('Iteration:', str(it), '; Error: %3.3e' % E)
if E < tol:
break
# Update
u.vector.array[:] = u.vector.array[:] + x.vector.array[:]
it += 1
return it, it < max_it
t = 0
dt=0.1
#Time loop
for i in range(1,6):
t += dt
print("time: " + str(np.round(t, 2)))
q_p = heat_cond_polynom(uh_p)
q_i = dolfinx.Function(V)
q_i.vector.array[:] = heat_cond_interpolation(np.array(uh_i.vector.array))
Fp = 0.54 * 0.0078 * uh_p * v * ufl.dx + q_p*dt*ufl.dot(ufl.grad(uh_p), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n_p + dt*f) * v * ufl.dx
Fi = 0.54 * 0.0078 * uh_i * v * ufl.dx + q_i*dt*ufl.dot(ufl.grad(uh_i), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n_i + dt*f) * v * ufl.dx
# Boundary conditions (Robin on all sides)
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
marker = 1
facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 2) | np.isclose(x[1], 0) | np.isclose(x[1], 2))
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
Fp += dt * 0.03 * ufl.inner(uh_p-20, v) * ds(marker)
Fi += dt * 0.03 * ufl.inner(uh_i-20, v) * ds(marker)
#Problem definition - polynomial
problem = dolfinx.fem.NonlinearProblem(Fp, uh_p)
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(uh_p)
assert(converged)
print(f"Polynomial number of iterations: {n:d}")
#Problem definition - interpolation
u = ufl.TrialFunction(V)
J = ufl.derivative(Fi, uh_i, u)
n, converged = newtonSolve(Fi, J, uh_i, [], q_i, heat_cond_interpolation)
assert(converged)
print(f"Interpolation number of iterations: {n:d}")
#Print max temperature
print("Max temp, polynomial:",np.max(np.round(np.array(uh_p.vector.array),2)))
print("Max temp, interpolation:",np.max(np.round(np.array(uh_i.vector.array),2)))
# Update solution at previous time step (u_n)
with uh_p.vector.localForm() as loc, u_n_p.vector.localForm() as loc_n:
loc.copy(result=loc_n)
with uh_i.vector.localForm() as loc, u_n_i.vector.localForm() as loc_n:
loc.copy(result=loc_n)
EDIT: updated my code to deal with the case of multiple nonlinear coefficients q.