Shape optimization in complex dolfinx

variational formulation
1

I need to adapt a shape optimization example for my use case. The unknown solution field is complex and hence I need to work with the complex build. While porting to the recent version, I run into errors that I do not encounter in the real mode. I would appreciate if someone could point out the mistake.

Below is an MWE evaluated on complex build of 0.10 installed from conda:

import dolfinx, ufl, gmsh, mpi4py

gmsh.initialize()
gmsh.clear()
gmsh.option.setNumber("General.Terminal", 0)
occ = gmsh.model.occ

gdim = 2
b1 = occ.addRectangle(0, 0, 0, 1, 1)
occ.synchronize()

# add_sub physical tags
all_doms = occ.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = occ.getEntities(gdim-1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)

mpi_rank = 0
mesh_data = dolfinx.io.gmsh.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim)

mesh, ct, ft = mesh_data.mesh, mesh_data.cell_tags, mesh_data.facet_tags

fem_order = 2
Omega = dolfinx.fem.functionspace(mesh, ("CG", fem_order))

u = ufl.TrialFunction(Omega)
ut = ufl.TestFunction(Omega)

dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)

R = ufl.inner(ufl.grad(u), ufl.grad(ut))*dx
zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))
R += ufl.inner(zero, ut)*dx

u_sol = dolfinx.fem.Function(Omega)

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
import dolfinx.fem.petsc
dbc = []
fwd_solver = dolfinx.fem.petsc.LinearProblem(ufl.lhs(R), ufl.rhs(R), u=u_sol, bcs=dbc, petsc_options_prefix="fwd_",petsc_options=petsc_options)
fwd_solver.solve()

lmbda = dolfinx.fem.Function(Omega)
J = ufl.inner(u_sol, u_sol) * dx
L = ufl.replace(ufl.action(R, u_sol), {ut: lmbda}) + J

dRdu = ufl.derivative(ufl.action(R, u_sol), u_sol, u)
dRdu_adj = ufl.adjoint(dRdu)
dJdu = ufl.derivative(J, u_sol, ufl.conj(ut))

dbc_adj = []
# This gives error
adj_solver = dolfinx.fem.petsc.LinearProblem(ufl.replace(dRdu_adj, {u_sol: ut}), -dJdu, u=lmbda, bcs=dbc_adj, 
                                              petsc_options_prefix="adj_", petsc_options=petsc_options)
adj_solver.solve()

X = ufl.SpatialCoordinate(mesh)
W = dolfinx.fem.functionspace(mesh, mesh.ufl_domain().ufl_coordinate_element())
dL = ufl.derivative(L, X, ufl.TestFunction(W))
dLf = dolfinx.fem.form(dL) # this gives error too in complex mode
2

You need to be really careful with usage of inner,dot, and conj when working in complex mode.
Here is a corrected code:

import dolfinx, ufl, gmsh, mpi4py

gmsh.initialize()
gmsh.clear()
gmsh.option.setNumber("General.Terminal", 0)
occ = gmsh.model.occ

gdim = 2
b1 = occ.addRectangle(0, 0, 0, 1, 1)
occ.synchronize()

# add_sub physical tags
all_doms = occ.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = occ.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(
        edge[0], [edge[1]], edge[1]
    )  # create the main group/node

gmsh.model.mesh.generate(gdim)

mpi_rank = 0
mesh_data = dolfinx.io.gmsh.model_to_mesh(
    gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim
)

mesh, ct, ft = mesh_data.mesh, mesh_data.cell_tags, mesh_data.facet_tags

fem_order = 2
Omega = dolfinx.fem.functionspace(mesh, ("CG", fem_order))

u = ufl.TrialFunction(Omega)
ut = ufl.TestFunction(Omega)

dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)

R = ufl.inner(ufl.grad(u), ufl.grad(ut)) * dx
zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))
R += ufl.inner(zero, ut) * dx

u_sol = dolfinx.fem.Function(Omega)

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
}
import dolfinx.fem.petsc

dbc = []
fwd_solver = dolfinx.fem.petsc.LinearProblem(
    ufl.lhs(R),
    ufl.rhs(R),
    u=u_sol,
    bcs=dbc,
    petsc_options_prefix="fwd_",
    petsc_options=petsc_options,
)
fwd_solver.solve()

lmbda = dolfinx.fem.Function(Omega)
J = ufl.dot(u_sol, u_sol) * dx
L = ufl.replace(ufl.action(R, u_sol), {ut: lmbda}) + J

dRdu = ufl.derivative(ufl.action(R, u_sol), u_sol, u)
dRdu_adj = ufl.adjoint(dRdu)
dJdu = ufl.derivative(J, u_sol, ufl.conj(ut))

dbc_adj = []
# This gives error
adj_solver = dolfinx.fem.petsc.LinearProblem(
    ufl.replace(dRdu_adj, {u_sol: ut}),
    -dJdu,
    u=lmbda,
    bcs=dbc_adj,
    petsc_options_prefix="adj_",
    petsc_options=petsc_options,
)
adj_solver.solve()

X = ufl.SpatialCoordinate(mesh)
W = dolfinx.fem.functionspace(mesh, mesh.ufl_domain().ufl_coordinate_element())
dL = ufl.derivative(L, X, ufl.conj(ufl.TestFunction(W)))
dLf = dolfinx.fem.form(dL)  # this gives error too in complex mode

Note that I use ufl.dot for the J, and ufl.conj in the last derivative for the test function.
Iā€™ve tried to highlight some of these things in: The Poisson problem with complex numbers ā€” FEniCSx tutorial

2 Likes
3

Thanks a lot for the quick fix!

I might have guessed or seen the requirement on ufl.derivative, but to be forced to use J = u_sol**2*dx instead of J = ufl.inner(u_sol,u_sol)*dx is very unexpected. Luckily, I donā€™t need it for my current assignment. But quite often optimization based on the energy functional |u_sol|^2 will be necessary. How can we work with that?

4

Complex-valued optimization is not trivial.
UFL enforces some strict (yet reasonable) rules on what it can allow to compile.

If you have your functional as

J = \frac{1}{2}\int_\Omega u \cdot \bar{u}~\mathrm{d}x

You will get something a bit strange when computing \frac{\partial J}{\partial u}[v]:

\frac{\partial J}{\partial u}[v]= \frac{1}{2}\left(\int_{\Omega} v \cdot \bar{u} + u \cdot \bar{v}~\mathrm{d}x\right)

which UFL wouldnā€™t like, as illustrated below (only using DOLFINx to get the finite element spaces, everything is really on the symbolic analysis level)

from mpi4py import MPI
import dolfinx
import ufl

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a0 =u * ufl.conj(v) * ufl.dx
a1 = v * ufl.conj(u) * ufl.dx 
a = a0 + a1

for ai in [a0, a1, a]:
    pulled_back_ai = ufl.algorithms.compute_form_data(
        ai,
        do_apply_function_pullbacks=True,
        do_apply_integral_scaling=True,
        do_apply_geometry_lowering=True,
        preserve_geometry_types=(ufl.classes.Jacobian,),
        complex_mode=True
    )
    print(pulled_back_ai.integral_data[0])

yielding:

IntegralData over domain(cell, ('otherwise',)) with integrals:
weight * |(J[0, 0] * J[1, 1] + -1 * J[0, 1] * J[1, 0])| * (conj((reference_value(v_0)))) * (reference_value(v_1)) * dx(<Mesh #0>[('otherwise',)], {'estimated_polynomial_degree': 2}, {})
and metadata:
{metadata}
Traceback (most recent call last):
  File "/root/shared/mwe.py", line 14, in <module>
    pulled_back_ai = ufl.algorithms.compute_form_data(
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/compute_form_data.py", line 520, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py", line 213, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py", line 205, in check_integrand_arity
    raise ArityMismatch("Failure to conjugate test function in complex Form")
ufl.algorithms.check_arities.ArityMismatch: Failure to conjugate test function in complex Form

The most appropriate place to raise issues regarding complex optimization with UFL would be in the UFL issue tracker.

1 Like
5

I am facing another issue while adapting the original shape optimization problem. I notice that volume (or area) derivative is always zero in the complex mode. What is the right way of getting it to work? The MWE below calculates the derivative correctly when evaluated with the real build (0.10, conda) but gives zero with complex.

import dolfinx, ufl, mpi4py, gmsh
import numpy as np
import matplotlib.pyplot as plt
import dolfinx.fem.petsc

gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

occ = gmsh.model.occ
gdim = 2

box = occ.addRectangle(0, 0, 0, 1, 1)
occ.synchronize()

all_doms = occ.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

all_edges = occ.getEntities(gdim-1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)

mpi_rank = 0

mesh_data = dolfinx.io.gmsh.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim)
mesh = mesh_data.mesh

X = ufl.SpatialCoordinate(mesh)
W = dolfinx.fem.functionspace(mesh, mesh.ufl_domain().ufl_coordinate_element())

volume = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1.0))*ufl.dx
volume_form = dolfinx.fem.form(volume)
target_volume = mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(volume_form), op=mpi4py.MPI.SUM)
dvol = dolfinx.fem.form(ufl.derivative(volume, X, ufl.conj(ufl.TestFunction(W))))

dvol_func = dolfinx.fem.Function(W)
dolfinx.fem.petsc.assemble_vector(dvol_func.x.petsc_vec, dvol)

print("Current volume: ", target_volume)
print("min, max of dvol:", dvol_func.x.array.min(), dvol_func.x.array.max())
6

This seems to be a bug in UFL, as the following, ā€œpureā€ UFL example does not give the right result:

import ufl
from utils import LagrangeElement
import numpy as np

scalar_type=np.float64

mesh = ufl.Mesh(LagrangeElement(ufl.triangle, 1, shape=(2, )))

X = ufl.SpatialCoordinate(mesh)
W = ufl.FunctionSpace(mesh, mesh.ufl_coordinate_element())

volume = ufl.Constant(mesh)*ufl.dx
#volume_form = dolfinx.fem.form(volume)
#target_volume = mesh.comm.allreduce(
#    dolfinx.fem.assemble_scalar(volume_form), op=mpi4py.MPI.SUM)
ddvol = ufl.derivative(volume, X, ufl.conj(ufl.TestFunction(W)))

complex_mode = np.issubdtype(scalar_type, np.complexfloating)

fd = ufl.algorithms.compute_form_data(ddvol, 
                                       do_apply_function_pullbacks=True,
                                        do_apply_integral_scaling=True,
                                        do_apply_geometry_lowering=True,
                                        preserve_geometry_types=(ufl.classes.Jacobian,),
                                        do_apply_restrictions=True,
                                        do_append_everywhere_integrals=False, 
                                                                                                                complex_mode=complex_mode)


print(complex_mode, str(fd.preprocessed_form))
scalar_type=np.complex64
complex_mode = np.issubdtype(scalar_type, np.complexfloating)

fd = ufl.algorithms.compute_form_data(ddvol, 
                                       do_apply_function_pullbacks=True,
                                        do_apply_integral_scaling=True,
                                        do_apply_geometry_lowering=True,
                                        preserve_geometry_types=(ufl.classes.Jacobian,),
                                        do_apply_restrictions=True,
                                        do_append_everywhere_integrals=False, 
                                                                                                                complex_mode=complex_mode)


print(complex_mode, str(fd.preprocessed_form))

Issue created at:

Edit
Attempt at fixing this has been pushed to: Extend shape derivative to complex mode by jorgensd Ā· Pull Request #462 Ā· FEniCS/ufl Ā· GitHub

1 Like
7

Thanks for the clarification.

Could you share your opinion on how far this bug might be affecting complex mode shape derivatives in UFL/dolfinx? I am planning to do some experiments in frequency domain electromagnetics. If the issue likely goes beyond what I chanced upon then I will happily pause the activity until the problem is sorted out.

8

I think the fix above should resolve the issue. However, Iā€™d have to cook up some tests to verify it

1 Like