Non-convergence when using mixed element

General
1

Hello,
I wanted to add another filed to a working minimization problem. I added a simple quadratic term, so it’s minimum will be a zero tensor. However, the solver does not converge, even though my initial guess was close to the actual solution. What did I do wrong?

I provide a MWE:

import ufl
from dolfinx import mesh, fem, nls
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import numpy as np
from ufl import as_tensor, indices, inner, derivative

# Create the computational domain
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

# Define finite element spaces for vector and tensor fields
VFS = ufl.VectorElement('CG', domain.ufl_cell(), 1, dim=2)  # Vector field
TFS = ufl.TensorElement('CG', domain.ufl_cell(), degree=2)  # Tensor field
mel = ufl.MixedElement([VFS, TFS])  # Mixed element combining vector and tensor

# Create function spaces
V = fem.FunctionSpace(domain, VFS)
Q = fem.FunctionSpace(domain, TFS)
Y = fem.FunctionSpace(domain, mel)

# Define functions and split the mixed function
y = fem.Function(Y)
u, q = ufl.split(y)  # Split into vector and tensor components
yt = ufl.TestFunction(Y)
v, p = ufl.split(yt)

# Define initial guesses for the mixed function components
y.sub(0).interpolate(lambda x: [x[0], x[1]])  # Initial guess for the vector component
y.sub(1).interpolate(lambda x: np.zeros((4, x.shape[1])))  # Initial guess for the tensor component

# Define the Eucledean metric tensor `eta` as a UFL tensor
eta = as_tensor(((1, 0), (0, 1)))

# Define the elastic tensor `A`
i, j, k, l = indices(4)
A = as_tensor((eta[i, j] * eta[k, l] + (eta[i, k] * eta[j, l] + eta[i, l] * eta[j, k])), (i, j, k, l))

# Define a function to compute the actual metric tensor `g` from a vector field `u`
def g(u):
    gradu = ufl.grad(u)
    t1 = gradu[:, 0]
    t2 = gradu[:, 1]
    return as_tensor([[ufl.dot(t1, t1), ufl.dot(t1, t2)], [ufl.dot(t1, t2), ufl.dot(t2, t2)]])

# Define elastic strain and stress
def epsilon(u):
    return 0.5 * (g(u) - eta)
def sigma(u):
    return as_tensor(A[i, j, k, l] * epsilon(u)[k, l], (i, j))
def sigmaq(q):
    return as_tensor(A[i, j, k, l] * q[k, l], (i, j))

# Define boundary conditions
def lower_boundary(x):
    return np.isclose(x[1], 0)
def upper_boundary(x):
    return np.isclose(x[1], 1)

fdim = domain.topology.dim - 1
boundary_facets_up = mesh.locate_entities_boundary(domain, fdim, upper_boundary)
dofs_up = fem.locate_dofs_topological((Y.sub(0), V), fdim, boundary_facets_up)
u_up = fem.Function(V)
u_up.interpolate(lambda x: np.broadcast_to(ScalarType([[0], [0.1]]), (2, x.shape[1])))
bc_up = fem.dirichletbc(u_up, dofs_up, Y)

boundary_facets_down = mesh.locate_entities_boundary(domain, fdim, lower_boundary)
dofs_down = fem.locate_dofs_topological((Y.sub(0), V), fdim, boundary_facets_down)
u_down = fem.Function(V)
u_down.interpolate(lambda x: np.broadcast_to(ScalarType([[0], [0.0]]), (2, x.shape[1])))
bc_down = fem.dirichletbc(u_down, dofs_down, Y)

bc = [bc_up, bc_down]

# Define the total elastic energy
Energy = inner(sigma(u), epsilon(u)) * ufl.dx + inner(sigmaq(q), q) * ufl.dx

# Define the residual and solve the problem
Res = derivative(Energy, y, yt)
problem = fem.petsc.NonlinearProblem(Res, y, bcs=bc)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.solve(y)

# Output the total energy
print("Energy = ", fem.assemble_scalar(fem.form(Energy)))