Error with complex numbers in dolfinx

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.

I’ve ported your code to v0.4.2, and it runs nicely in the real mode of DOLFINx

import dolfinx
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

# mesh and function space
mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
                     [10, 10], dolfinx.mesh.CellType.quadrilateral)
V = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 1))
S = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
u = dolfinx.fem.Function(V)  # test function (displacment field)
v = ufl.TestFunction(V)  # trial function

# kinematics
lambda_ = 0.5769
mu = dolfinx.fem.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 = dolfinx.fem.Function(V)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, left)
with u1.vector.localForm() as loc:
    loc.set(0)
bc1 = dolfinx.fem.dirichletbc(u1, dolfinx.fem.locate_dofs_topological(
    V, fdim, boundary_facets))
u2 = dolfinx.fem.Function(V)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, right)
with u2.vector.localForm() as loc:
    loc.set(1)
bc2 = dolfinx.fem.dirichletbc(u2, dolfinx.fem.locate_dofs_topological(
    V, fdim, boundary_facets))
bcs = [bc1, bc2]

# solve the problem
problem = dolfinx.fem.petsc.NonlinearProblem(R, u, bcs)
solver = dolfinx.nls.petsc.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(dolfinx.fem.form(J))
print(J_old)

# partial derivative of J w.r.t. mu
dJdmu = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.derivative(J, mu)))
dJdmu.assemble()

# partial derivative of R w.r.t. mu
dRdmu = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(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.petsc.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)

There is nothing in your code that suggests that you want to solve a complex problem, thus you should run DOLFINx in real mode.

1 Like

Thank you for your reply.

So, is there an easy way to swap to real mode?
I think I am running the complex version by default.

How did you install DOLFINx?
This is dependent on which PETSc you have linked DOLFINx to.

1 Like

It is installed from source in a cluster that I am using.
I think that they have linked the complex version of PETSc.

You would need to have a separate installation with real PETSc.

As a side note, as of Allow `ufl.real/imag/conj` on vectors and tensors by jorgensd · Pull Request #496 · FEniCS/ffcx · GitHub you can solve your problem in complex mode with only minor adaptations:

import dolfinx
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

# mesh and function space
mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
                     [10, 10], dolfinx.mesh.CellType.quadrilateral)
V = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 1))
S = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
u = dolfinx.fem.Function(V)  # test function (displacment field)
v = ufl.TestFunction(V)  # trial function

# kinematics
lambda_ = 0.5769
mu = dolfinx.fem.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 = dolfinx.fem.Function(V)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, left)
with u1.vector.localForm() as loc:
    loc.set(0)
bc1 = dolfinx.fem.dirichletbc(u1, dolfinx.fem.locate_dofs_topological(
    V, fdim, boundary_facets))
u2 = dolfinx.fem.Function(V)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, right)
with u2.vector.localForm() as loc:
    loc.set(1)
bc2 = dolfinx.fem.dirichletbc(u2, dolfinx.fem.locate_dofs_topological(
    V, fdim, boundary_facets))
bcs = [bc1, bc2]

# solve the problem
problem = dolfinx.fem.petsc.NonlinearProblem(R, u, bcs)
solver = dolfinx.nls.petsc.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(dolfinx.fem.form(J))
print(J_old)

# partial derivative of J w.r.t. mu
dJdmu = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(ufl.derivative(J, mu, ufl.conj(ufl.TestFunction(S)))))
dJdmu.assemble()

# partial derivative of R w.r.t. mu
dRdmu = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(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, ufl.conj(ufl.TestFunction(V)))
problem = dolfinx.fem.petsc.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)

returning

[-1.60101558e-02+0.j -8.07170926e-03+0.j -6.59236668e-03+0.j
 -6.09747688e-03+0.j -7.37735007e-03+0.j -1.04254831e-04+0.j
 -5.08384394e-03+0.j -6.07808867e-03+0.j -3.21288677e-03+0.j
  4.89571867e-03+0.j -4.28682548e-03+0.j -5.31427653e-03+0.j
 -4.31069887e-03+0.j  1.14728912e-03+0.j  9.27162015e-03+0.j
 -3.50512617e-03+0.j -4.67718435e-03+0.j -4.35275019e-03+0.j
 -9.92771691e-04+0.j  5.49700643e-03+0.j  1.30468353e-02+0.j
 -2.66141735e-03+0.j -4.07593919e-03+0.j -4.14835915e-03+0.j
 -1.92005777e-03+0.j  3.14440675e-03+0.j  9.80778967e-03+0.j
  1.61469775e-02+0.j -1.71410859e-03+0.j -3.41852988e-03+0.j
 -3.83267290e-03+0.j -2.23695099e-03+0.j  1.83958822e-03+0.j
  7.82439509e-03+0.j  1.39919754e-02+0.j  1.84796335e-02+0.j
 -6.40016329e-04+0.j -2.60156136e-03+0.j -3.39922091e-03+0.j
 -2.28828468e-03+0.j  1.12854284e-03+0.j  6.58059039e-03+0.j
  1.27863114e-02+0.j  1.78454254e-02+0.j  1.99052154e-02+0.j
  6.99988738e-04+0.j -1.43870388e-03+0.j -2.74070548e-03+0.j
 -2.16255193e-03+0.j  7.13844308e-04+0.j  5.71701053e-03+0.j
  1.18845757e-02+0.j  1.75886831e-02+0.j  2.08231466e-02+0.j
  1.95021802e-02+0.j -9.92959123e-05+0.j -1.76141968e-03+0.j
 -1.85505412e-03+0.j  4.41865196e-04+0.j  5.03808249e-03+0.j
  1.10815122e-02+0.j  1.71414558e-02+0.j  2.14353079e-02+0.j
  2.18793475e-02+0.j -7.32868115e-04+0.j -1.41493759e-03+0.j
  1.96393835e-04+0.j  4.38923445e-03+0.j  1.03301674e-02+0.j
  1.66054241e-02+0.j  2.16083266e-02+0.j  2.35577151e-02+0.j
 -1.15453581e-03+0.j -2.07519126e-04+0.j  3.54103059e-03+0.j
  9.51718314e-03+0.j  1.60646375e-02+0.j  2.16548710e-02+0.j
  2.48773160e-02+0.j -1.09881044e-03+0.j  2.14980785e-03+0.j
  8.30147535e-03+0.j  1.54482079e-02+0.j  2.16534755e-02+0.j
  2.59123500e-02+0.j -1.51274036e-04+0.j  6.06653986e-03+0.j
  1.45122435e-02+0.j  2.16190611e-02+0.j  2.67402940e-02+0.j
  2.24423290e-03+0.j  1.21732749e-02+0.j  2.14708842e-02+0.j
  2.74472083e-02+0.j  6.95759135e-03+0.j  2.12639968e-02+0.j
  2.81875974e-02+0.j  1.63457221e-02+0.j  2.96036734e-02+0.j
  3.65642028e-02+0.j]
2 Likes

Ive now also added a section regarding this at:
https://jorgensd.github.io/dolfinx-tutorial/chapter1/complex_mode.html

2 Likes

Great! It will be very useful.

Thanks @dokken. Your tutorial has a paragraph on activating the complex number build. That would be a good place to mention the Debian/Ubuntu method for those who installed dolfinx via apt-get. On these systems the complex build is activated using the PETSC_DIR/SLEPC_DIR environment variables, e.g

PETSC_DIR=/usr/lib/petscdir/petsc-complex
SLEPC_DIR=/usr/lib/slepcdir/slepc-complex

I met the same problem to activate the complex number build in ubuntu20.04. Briefly, I installed fenicsx and the dependent packages via apt-get, and followed the instruction above to add a line “PETSC_DIR=/usr/lib/petscdir/petsc-complex” to ./bashrc, which is supposed to modify the environment variables.

The code:

from petsc4py import PETSc
print(PETSc.ScalarType)

returning

<class 'numpy.float64'>

Is there anything I missed? Many thanks!

Could be one of a number of problems here. Does your code work without the PETSC_DIR (i.e. in real-number mode) ?

My guess is that you need to export the environment variable, to make sure it’s getting propagated. For good measure I suggest also including the SLEPc counterpart

export PETSC_DIR=/usr/lib/petscdir/petsc-complex
export SLEPC_DIR=/usr/lib/slepcdir/slepc-complex

If that doesn’t fix it then we’ll need to check package versions and make sure nothing got installed by pip.

By the way, I wouldn’t necessarily recommend setting PETSC_DIR in .bashrc. If you do then you’ll always be running in complex mode, which can interfere if you later try to run jobs in real mode. I guess it’s okay to do it if you know you always want to use complex mode by default (you can override it on the command line back to real mode when launching a job if needed)

1 Like

Thanks a lot! Dparsons. I tested with demo_poisson.py. and got the output as expected in real number mode. but in complex number mode(just adding the PETSC_DIR to .bashrc), I got an error as below,

Traceback (most recent call last):
  File "/usr/lib/python3.8/idlelib/run.py", line 559, in runcode
    exec(code, self.locals)
  File "/home/spm/fenics/test/001/demo_poisson.py", line 140, in <module>
    problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
  File "/usr/lib/petscdir/petsc-complex/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 518, in __init__
    self._A = create_matrix(self._a)
  File "/usr/lib/petscdir/petsc-complex/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 130, in create_matrix
    return _cpp.fem.petsc.create_matrix(a)
TypeError: create_matrix(): incompatible function arguments. The following argument types are supported:
    1. (a: dolfinx::fem::Form<std::complex<double> >, type: str = '') -> mat

Invoked with: <dolfinx.fem.forms.Form object at 0x7f52541b8720>

Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
when compiling your pybind11 module.

Attached are the version of the packages installed



and some python packages

enics-basix           0.4.2               
fenics-dolfinx         0.4.1               
fenics-ffcx            0.4.2               
fenics-ufl             2022.1.0   
petsc4py               3.12.0              
petsc4py-complex       3.12.0              
petsc4py-real          3.12.0              

Thanks again for helping debug the problem!

I can reproduce the problem. There seems to be an inconsistency in the complex build on Ubuntu focal. My recommended workaround is to upgrade your system to the latest Ubuntu release (jammy 22.04).

Hi,
i have the same problem. I installed fenicsx through apt with:

sudo apt update
sudo apt install fenicsx

and i could use the real-mode without any problems. If i change the PETSC_DIR via

export PETSC_DIR=/usr/lib/petscdir/petsc-complex

dolfinx cant be imported anymore in the python console.
I have inspected the petscdir directory and there is no petsc-complex folder/file present, only a petsc-real folder. It seems like, the complex-petsc version is not installed in the fenicsx-package. Do i need to install an other fenicsx package via apt to obtain the complex-mode or am i missing something here? Im using Ubuntu 22.

Thanks in Advance,
Tobias

Hi Tobias, the fenicsx package is set up to install the “default” build (i.e. real build). To install the complex build, add

sudo apt install python3-dolfinx-complex

After upgrade to Jammy 22.04 and reinstall FEniCSx, the complex build works eventually! Thanks a lot!