Solve the nonlinear coupled-field problem with FEniCSX

1

The solution scheme of the nonlinear single field problem is well established in the tutorial Hyperelasticity — FEniCSx tutorial (jorgensd.github.io). But how to solve coupled-field problems monolithically (i.e., not the staggered one) with the Newton solver?

For example, suppose we have two coupled PDEs f(u,y)=0 and g(u,y)=0 controlling the evolution of two fields u and y. It is trivial to obtain the corresponding variational forms of the residuals of the two PDEs as F_u(u,y)=0 and F_y(u,y)=0. Then how to construct the monolithic residual F_uy(u,y) such that it can be substituted in the class variable NonlinearPDEProblem in Hyperelasticity — FEniCSx tutorial (jorgensd.github.io)?

As for the minimal example based on Hyperelasticity — FEniCSx tutorial (jorgensd.github.io), we can simply set the second field y same as the original one u, although these two fields u and y are independent in this case.

w = ufl.TestFunction(V) # test function of the second field
y = dolfinx.Function(V) # trial function of the second field

In addition, we set that the boundaries conditions, residual vectors, and other related quantities of the field y have the same form as the filed u.

# Deformation gradient of the second field
F_y = ufl.variable(I + ufl.grad(y))

# Right Cauchy-Green tensor of the second field
C_y = ufl.variable(F_y.T * F_y)

# Invariants of deformation tensors of the second field
Ic_y = ufl.variable(ufl.tr(C_y))
J_y = ufl.variable(ufl.det(F_y))

# Stored strain energy density of the second field (compressible neo-Hookean model)
psi_y = (mu / 2) * (Ic_y - 3) - mu * ufl.ln(J_y) + (lmbda / 2) * (ufl.ln(J_y))**2

# Stress of the second field
P_y = ufl.diff(psi_y, F_y)

# Define form F of the second field (we want to find u such that F_y(y) = 0)
F_y = ufl.inner(ufl.grad(w), P_y)*dx - ufl.inner(w, B)*dx - ufl.inner(w, T)*ds(2)

Then how to solve this “coupled” problem?

problem = NonlinearPDEProblem(F_uy, u, bcs) # how to set up the monolithic residual F_uy?

Theorectically, we should obtain the same field for u and y, and they coincide with the result of the single phase problem Hyperelasticity — FEniCSx tutorial (jorgensd.github.io).

Thanks!

2

Try making a mixed function space with u and y:

v_el = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
mixed_el = ufl.MixedElement([v_el, v_el])
VP = dolfinx.FunctionSpace(mesh, mixed_el)
(v, w) = ufl.TestFunctions(VP)
q = dolfinx.Function(VP)
u, y = ufl.split(q)
1 Like