I’m building a script to minimize the free energy of a liquid crystal model. The mesh and energy functions are created via dolfinx and ufl for the flexibility to define any geometry. The minimization of the free energy is passed to scipy minimizer (L-BFGS-B) by converting the parameter back and forth between fem.Function and numpy array.
However, the minimization doesn’t work as expected. I suspect that the definition of eval_gradient(x) is not correct, but I’m not sure. Please let me know if you have any suggestion. Thanks a lot!
Here is the code:
from dolfinx import mesh, fem
from dolfinx.fem import functionspace
from dolfinx.mesh import locate_entities_boundary, meshtags
from basix.ufl import element, mixed_element
import ufl
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
— Domain setup —
Lx, Ly = 10.0, 2.0 # Domain dimensions
Nx, Ny = 50, 10 # Mesh resolution
Create 2D rectangular mesh
domain = mesh.create_rectangle(
MPI.COMM_WORLD,
points=[(0.0, 0.0), (Lx, Ly)],
n=(Nx, Ny),
cell_type=mesh.CellType.triangle
)
Define scalar and vector elements
scalar_el = element(“Lagrange”, domain.topology.cell_name(), 1)
vector_el = element(“Lagrange”, domain.topology.cell_name(), 1, shape=(2,))
Combine into mixed element
mixed_el = mixed_element([scalar_el, vector_el])
V = fem.functionspace(domain, mixed_el)
there is no trail or test functions, because the problem is non-linear
ArityMismatch error when applying nonlinear operator "Power"
u = fem.Function(V)
rho, Q_vec = ufl.split(u)
Free energy coefficients (adjust as needed)
a = fem.Constant(domain, PETSc.ScalarType(-1))
b = fem.Constant(domain, PETSc.ScalarType(0))
c = fem.Constant(domain, PETSc.ScalarType(10))
B = fem.Constant(domain, PETSc.ScalarType(10))
q = fem.Constant(domain, PETSc.ScalarType(1))
K = fem.Constant(domain, PETSc.ScalarType(1))
l = fem.Constant(domain, PETSc.ScalarType(1))
W = fem.Constant(domain, PETSc.ScalarType(10000))
a_val = float(a.value)
b_val = float(b.value)
q_val = float(q.value)
K_val = float(K.value)
B_val = float(B.value)
l_val = float(l.value)
dx = ufl.Measure(“dx”, domain)
#bulk ============================================================================
f_bulk = 0.5 * a * rho2 + (1/3) * b * rho3 + 0.25 * c * rho**4
F_bulk = f_bulk * dx
#layer ============================================================================
1. Compute the Hessian of ρ (symmetric 2nd derivative tensor)
rho_grad = ufl.grad(rho) # ∇ρ: (2,)
rho_hessian = ufl.sym(ufl.grad(rho_grad)) # Hessian: (2,2) symmetric
2. Reconstruct full Q-tensor from Q_vec = [Qxx, Qxy]
Qxx = Q_vec[0]
Qxy = Q_vec[1]
Q_tensor = ufl.as_tensor([[Qxx, Qxy],
[Qxy, -Qxx]]) # Enforces tr(Q) = 0
3. Identity tensor in 2D
I = ufl.Identity(2)
4. Coupling term: q² (Q + I/2) * ρ
Q_plus_I = Q_tensor + 0.5 * I
coupling_tensor = q**2 * Q_plus_I * rho # (2,2)
5. Total tensor to square
M = rho_hessian + coupling_tensor # (2,2)
6. Frobenius norm squared: sum over all M_ij²
layering_energy_density = B * ufl.inner(M, M)
7. Add to weak form (residual)
F_layer = layering_energy_density * dx
#elastic ============================================================================
Q_grad = ufl.grad(Q_tensor) # shape (2,2,2)
4. Compute energy density: sum over all (∂ₖ Q_{ij})²
Use inner product over 3 indices
elastic_energy_density = 0.5 * K * ufl.inner(Q_grad, Q_grad) # scalar
5. Add to weak form
F_elastic = elastic_energy_density * dx
#phase ============================================================================
tr_Q2 = ufl.inner(Q_tensor, Q_tensor)
Compute Q^2 first
Q_sq = ufl.dot(Q_tensor, Q_tensor) # Matrix multiplication
Compute tr(Q^3) = Q_ij Q_jk Q_ki = inner(Q^2, Q)
tr_Q3 = ufl.inner(Q_sq, Q_tensor)
Landau-de Gennes potential
f_LdG = l * (-0.5 * tr_Q2 - (1.0 / 3.0) * tr_Q3 + 0.5 * tr_Q2**2)
Add to weak form
F_LdG = f_LdG * dx
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
top_facets = locate_entities_boundary(domain, dim=domain.topology.dim - 1,
marker=lambda x: np.isclose(x[1], Ly)) # y = Ly
bottom_facets = locate_entities_boundary(domain, dim=domain.topology.dim - 1,
marker=lambda x: np.isclose(x[1], 0.0))
facet_indices = np.concatenate([top_facets, bottom_facets])
facet_markers = np.concatenate([
np.full_like(top_facets, 1), # 1 = top
np.full_like(bottom_facets, 2) # 2 = bottom
])
facet_tag = meshtags(domain, 1, facet_indices, facet_markers)
ds = ufl.Measure(“ds”, domain=domain, subdomain_data=facet_tag)
Preferred Q tensors for top and bottom surfaces
Q0_top_tensor = ufl.as_tensor([[-0.5, 0.0],
[ 0.0, 0.5]])
Q0_bot_tensor = ufl.as_tensor([[ 0.5, 0.0],
[ 0.0, -0.5]])
Anchoring energy densities (Frobenius norm squared)
anchoring_top = W * ufl.inner(Q_tensor - Q0_top_tensor, Q_tensor - Q0_top_tensor)
anchoring_bot = W * ufl.inner(Q_tensor - Q0_bot_tensor, Q_tensor - Q0_bot_tensor)
Integrate only over top and bottom boundaries
F_anchor = anchoring_top * ds(1) + anchoring_bot * ds(2)
Total energy functional from all contributions
F_total = F_bulk + F_layer + F_elastic + F_LdG + F_anchor
Create FEniCSx functional for energy
energy_form = fem.form(F_total)
— DOF I/O Utilities —
def extract_dofs(u):
return np.array(u.x.array)
def assign_dofs(u, dofs):
u.vector.array[:] = np.array(dofs)
Flattened initial guess
x0 = u.x.array.copy()
Assemble scalar energy
def eval_energy(x_flat):
assign_dofs(u, x_flat)
return fem.assemble_scalar(energy_form)
#=================================Here======================================
q = ufl.TrialFunction(V)
dFdu = ufl.derivative(F_total, u, q)
dFdu_compiled = fem.form(dFdu)
dFdu_value = fem.Function(V)
def eval_gradient(x):
assign_dofs(u, x)
#dFdu_value.x.array[:] = 0
fem.assemble_vector(dFdu_value.x.array, dfdu_compiled)
return dFdu_value.x.array
#=======================================================================
res = minimize(
fun=eval_energy,
x0=x0,
jac=eval_gradient,
method=‘L-BFGS-B’,
options={‘disp’: True, ‘maxiter’: 500}
)
Update solution
assign_dofs(u, res.x)
Extract rho component
rho_fem = u.sub(0).collapse()
Coordinates
coords = domain.geometry.x
Cell-to-vertex connectivity
domain.topology.create_connectivity(2, 0)
connectivity = domain.topology.connectivity(2, 0)
num_cells = domain.topology.index_map(2).size_local
cells = [connectivity.links(i) for i in range(num_cells)]
cells = np.array(cells)
rho values at nodes
rho_vals = rho_fem.x.array
Plot with tricontourf
plt.figure(figsize=(8, 2))
plt.tricontourf(coords[:, 0], coords[:, 1], cells, rho_vals, levels=100, cmap=“viridis”)
plt.colorbar(label=r"\rho")
plt.xlabel(“x”)
plt.ylabel(“y”)
plt.title(“Colormap of \\rho”)
plt.axis(“equal”)
plt.tight_layout()
plt.show()