Dear all,
I am solving the adjoint Poisson equation to perform a sensitivity analysis. This is a minimal code:
import dolfinx
import ufl
from dolfinx.generation import RectangleMesh
from dolfinx.cpp.mesh import CellType
from dolfinx.fem import VectorFunctionSpace, FunctionSpace, Function, DirichletBC, form
from dolfinx.nls import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
# mesh and function space
mesh = RectangleMesh(MPI.COMM_WORLD, [[0, 0, 0], [1, 1, 0]],
[10, 10], CellType.quadrilateral)
V = VectorFunctionSpace(mesh, ("CG", 1))
S = FunctionSpace(mesh, ("DG", 0))
u = Function(V) # test function (displacment field)
v = ufl.TestFunction(V) # trial function
# kinematics
lambda_ = 0.5769
mu = Function(S) # in this case, we set the mu as a design variable
mu.vector.array = 0.3846
P = 2.0*mu*ufl.sym(ufl.grad(u)) + lambda_ * ufl.tr(ufl.sym(ufl.grad(u)))*ufl.Identity(len(u))
dx = ufl.Measure("dx", metadata={"quadrature_degree": 4})
R = ufl.inner(P, ufl.grad(v))*dx
# boundary conditions
def left(x): return np.isclose(x[0], 0)
def right(x): return np.isclose(x[0], 1)
fdim = mesh.topology.dim - 1
u1 = Function(V)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, left)
with u1.vector.localForm() as loc:
loc.set(0)
bc1 = DirichletBC(u1, dolfinx.fem.locate_dofs_topological(
V, fdim, boundary_facets))
u2 = Function(V)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, right)
with u2.vector.localForm() as loc:
loc.set(1)
bc2 = DirichletBC(u2, dolfinx.fem.locate_dofs_topological(
V, fdim, boundary_facets))
bcs = [bc1, bc2]
# solve the problem
problem = dolfinx.fem.NonlinearProblem(R, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.solve(u)
# print(u.vector.array)
# set an arbitrary functinal
J = mu*ufl.dot(u, u)*dx
J_old = dolfinx.fem.assemble_scalar(J)
print(J_old)
# partial derivative of J w.r.t. mu
dJdmu = dolfinx.fem.assemble_vector(ufl.derivative(J, mu))
dJdmu.assemble()
# partial derivative of R w.r.t. mu
dRdmu = dolfinx.fem.assemble_matrix(ufl.adjoint(ufl.derivative(R, mu))) # partial derivative
dRdmu.assemble()
# reset the boundary condition
with u2.vector.localForm() as loc:
loc.set(0.0)
# solve the adjoint vector
lhs = ufl.adjoint(ufl.derivative(R, u))
rhs = -ufl.derivative(J, u)
problem = dolfinx.fem.LinearProblem(lhs, rhs, bcs=bcs)
lmbda = problem.solve()
dJdmu += dRdmu*lmbda.vector
dJdmu.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
print(dJdmu.array)
The error that I get is:
Traceback (most recent call last):
File “Poisson.py”, line 59, in
dJdmu = dolfinx.fem.assemble_vector(ufl.derivative(J, mu))
File “/home/jl2269/anaconda3/lib/python3.8/functools.py”, line 875, in wrapper
return dispatch(args[0].class)(*args, **kw)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/dolfinx/fem/assemble.py”, line 91, in assemble_vector
_L = _create_cpp_form(L)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/dolfinx/fem/assemble.py”, line 27, in _create_cpp_form
return Form(form)._cpp_object
File “/usr/dolfinx/venv/lib/python3.8/site-packages/dolfinx/fem/form.py”, line 42, in init
self._ufc_form, module, self._code = jit.ffcx_jit(
File “/usr/dolfinx/venv/lib/python3.8/site-packages/dolfinx/jit.py”, line 61, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/dolfinx/jit.py”, line 215, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], parameters=p_ffcx, **p_jit)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ffcx/codegeneration/jit.py”, line 167, in compile_forms
impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ffcx/codegeneration/jit.py”, line 232, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, parameters=parameters)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ffcx/compiler.py”, line 98, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, parameters)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ffcx/analysis.py”, line 64, in analyze_ufl_objects
form_data = tuple(_analyze_form(form, parameters) for form in forms)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ffcx/analysis.py”, line 64, in
form_data = tuple(_analyze_form(form, parameters) for form in forms)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ffcx/analysis.py”, line 157, in _analyze_form
form_data = ufl.algorithms.compute_form_data(
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ufl/algorithms/compute_form_data.py”, line 407, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ufl/algorithms/check_arities.py”, line 177, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File “/usr/dolfinx/venv/lib/python3.8/site-packages/ufl/algorithms/check_arities.py”, line 170, 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
It seems that the variable J is complex, which can’t be possible because J is proportional to the dot product of u and u. Therefore, to see what is happening, I printed the value of J_old and it gives me a complex number with zero imaginary part. I have just realised that dolfinx always returns me a complex number in all the functions (if the output is supossed to be a real number, it gives me a complex number with zero imaginary part).
Does anybody know how to fix this issue?
The version that I am using is:
0.3.0 35d705001bdd42b859ac7f32f124cbf9a41d9b3b
Thank you very much in advance.
