Shape optimization in complex dolfinx

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())