How to define conjugates for variational formulation?

variational formulation
1

Hi, I am trying to solve a problem which requires me to form a variational formulation with a term involving the conjugate of a function, as in the following minimum working code example:

  from petsc4py import PETSc
  import dolfinx
  import dolfinx.fem.petsc
  import dolfinx.nls.petsc
  import ufl
  from mpi4py import MPI  
  mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
  V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
  u = dolfinx.fem.Function(V)
  v = ufl.TestFunction(V)
  F = ufl.inner(ufl.grad(ufl.conj(u)), ufl.grad(v)) * ufl.dx
  mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
  boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
  boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
  bc = dolfinx.fem.dirichletbc(dolfinx.fem.Function(V), boundary_dofs)
  ns_problem = dolfinx.fem.petsc.NonlinearProblem(F, u, bcs=bc)
  newton_solver = dolfinx.nls.petsc.NewtonSolver(mesh.comm, ns_problem) 
  newton_solver.solve(u)

The error message I am having with this code is:

Traceback (most recent call last):
  File "/app/fx4f/run_simulation.py", line 40, in <module>
    run_simulation(simulation.main)
  File "/app/fx4f/run_simulation.py", line 25, in run_simulation
    main( logger, **kwargs )
  File "/app/simulation.py", line 35, in main
    ns_problem = dolfinx.fem.petsc.NonlinearProblem(F, u, bcs=bc)
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 907, in __init__
    self._a = _create_form(
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 249, in form
    return _create_form(form)
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 244, in _create_form
    return _form(form)
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 186, in _form
    ufcx_form, module, code = jit.ffcx_jit(
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/jit.py", line 51, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/jit.py", line 201, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 256, in compile_forms
    impl = _compile_objects(
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 383, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
  File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 108, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 94, in analyze_ufl_objects
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 94, in <genexpr>
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 180, in _analyze_form
    form_data = ufl.algorithms.compute_form_data(
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/compute_form_data.py", line 427, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py", line 211, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py", line 205, in check_integrand_arity
    raise ArityMismatch(f"Argument {arg} is spuriously conjugated in complex Form")
ufl.algorithms.check_arities.ArityMismatch: Argument v_1 is spuriously conjugated in complex Form

My fenicsx version is 0.8.0, and I installed it using Docker. I am using my PETSc in complex mode, which was again set from the Dockerfile.
The code above is created only to represent the problem I am facing with my original code.

I don’t know how I am supposed to include the conjugate of a function in the variational formulation, and why I am getting this error. I would appreciate it if you could inform me about the usage of conjugates. Thank you in advance!

2

ufl.inner already does the conj for you, see ufl/ufl/operators.py at 004a678320e5c65015d56071b968ee08314aecc1 · FEniCS/ufl · GitHub
Furthermore, ufl expects that conjugation happens on the test function, not on the trial, see ufl/ufl/algorithms/check_arities.py at 004a678320e5c65015d56071b968ee08314aecc1 · FEniCS/ufl · GitHub

The following

import dolfinx
import dolfinx.fem.petsc
import ufl
from mpi4py import MPI

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)

print("Case 1")
try:
    a1 = ufl.inner(ufl.grad(ufl.conj(u)), ufl.grad(v)) * ufl.dx
    dolfinx.fem.form(a1)
except BaseException as e:
    print(e)
else:
    print("ok")

print("Case 2")
a2 = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
dolfinx.fem.form(a2)
print("ok")

print("Case 3: same as case 2")
a3 = ufl.dot(ufl.grad(u), ufl.grad(ufl.conj(v))) * ufl.dx
dolfinx.fem.form(a3)
print("ok")


print("ok")
print("Case 4: same as case 1, but different error")
try:
    a4 = ufl.dot(ufl.grad(ufl.conj(u)), ufl.grad(v)) * ufl.dx
    dolfinx.fem.form(a4)
except BaseException as e:
    print(e)
else:
    print("ok")

prints

Case 1
Argument v_1 is spuriously conjugated in complex Form
Case 2
ok
Case 3: same as case 2
ok
ok
Case 4: same as case 1, but different error
Failure to conjugate test function in complex Form
3

Hi Francesco, thank you for your answer.

In my case, u is not a TrialFunction but a fem.Function. Also, my problem’s variational form includes the conjugate of a function explicitly, as I provided in the code above. Is there any chance to implement that?

4

Is the form sesquilinear?
This is a requirement: Form language — Unified Form Language (UFL) 2024.1.0 documentation

And please note that inner also incorporates conj on the second argument, i.e.
ufl.inner(ufl.grad(ufl.conj(u)), ufl.grad(v))*ufl.dx corresponds to
ufl.dot(ufl.grad(ufl.conj(u)), ufl.conj(ufl.grad(v)))*ufl.dx.
You are getting an error when the jacobian is formed, i.e. when
J=ufl.derivative(F, u, du) is formed.
The error stems from: ufl/ufl/algorithms/check_arities.py at 004a678320e5c65015d56071b968ee08314aecc1 · FEniCS/ufl · GitHub
where the assumption is that a TrialFunction cannot be conjugated.
The TrialFunction itself is real valued: Form language — Unified Form Language (UFL) 2024.1.0 documentation meaning that having a ufl.conj on it doesn’t really make sense.

So the big question is, what is the Jacobian of your complex variational form?

5

Hi dokken, thank you for your answer.

Consider the example in here. It would seem that taking a conjugate of such a solution should be a valid operation.

Suppose in that tutorial, if the PDE were not ‘Delta u_c = f’, but rather ‘Delta conj(u_c) = f’, then how would one solve that?

6

That equation doesn’t really make sense to me, as you are searching for complex valued coefficients a+ib, so you could solve the non-conjugate version of this and just get a negative b coeff.

7

I see, for that example your suggestion is indeed the easiest solution. What would you suggest for a variational formulation like the one in the below code?

  import dolfinx
  import dolfinx.fem.petsc
  import dolfinx.nls.petsc
  import ufl
  from mpi4py import MPI
  import basix
  mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
  Ve = basix.ufl.element("Lagrange", mesh.topology.cell_name(), Pu, shape=(mesh.geometry.dim,)) 
  Qe = basix.ufl.element("Lagrange", mesh.topology.cell_name(), Pp)
  W_element = basix.ufl.mixed_element([Ve, Qe])
  W = dolfinx.fem.functionspace(mesh, W_element)
  V, _ = W.sub(0).collapse()
  vq = ufl.TestFunction(W)
  (v, q) = ufl.split(vq)
  dup = ufl.TrialFunction(W)    
  up = dolfinx.fem.Function(W)
  (u, p) = ufl.split(up)
  F = ufl.inner( ufl.dot(ufl.conj(u),ufl.nabla_grad(u))  , v) * ufl.dx
  J = ufl.derivative(F, up, dup)
  mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
  boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
  boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
  bc = dolfinx.fem.dirichletbc(dolfinx.fem.Function(V), boundary_dofs)
  ns_problem = dolfinx.fem.petsc.NonlinearProblem(F, up, bcs=[bc], J =J)
  newton_solver = dolfinx.nls.petsc.NewtonSolver(mesh.comm, ns_problem) 
  newton_solver.solve(up)

The code above yields the error below:

Traceback (most recent call last):
  File "/app/fx4f/run_simulation.py", line 40, in <module>
    run_simulation(simulation.main)
  File "/app/fx4f/run_simulation.py", line 25, in run_simulation
    main( logger, **kwargs )
  File "/app/simulation.py", line 66, in main
    ns_problem = dolfinx.fem.petsc.NonlinearProblem(F, up, bcs=[bc], J =J)
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 907, in __init__
    self._a = _create_form(
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 249, in form
    return _create_form(form)
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 244, in _create_form
    return _form(form)
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 186, in _form
    ufcx_form, module, code = jit.ffcx_jit(
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/jit.py", line 51, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/jit.py", line 201, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 256, in compile_forms
    impl = _compile_objects(
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 383, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
  File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 108, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 94, in analyze_ufl_objects
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 94, in <genexpr>
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 180, in _analyze_form
    form_data = ufl.algorithms.compute_form_data(
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/compute_form_data.py", line 427, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py", line 211, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py", line 192, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
  File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
    (result,) = map_expr_dags(
  File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 103, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py", line 57, in sum
    raise ArityMismatch(
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments ('conj(v_1)',) vs ('v_1',).

It again raises an ArityMismatch, but of a different kind. Is it again because of the way I am forming the Jacobian? I didn’t understand how the code takes the conjugate of the TrialFunction, although I am not taking it anywhere in the code. I also tried using inner to avoid using conjugates explicitly, and changed the F to:

F = ufl.inner( ufl.inner(ufl.nabla_grad(u), u)  , v) * ufl.dx

Which resulted in the following error:

Traceback (most recent call last):
  File "/app/fx4f/run_simulation.py", line 25, in run_simulation
    main( logger, **kwargs )
  File "/app/simulation.py", line 62, in main
    F = ufl.inner( ufl.inner(ufl.nabla_grad(u), u)  , v) * ufl.dx
  File "/usr/local/lib/python3.10/dist-packages/ufl/operators.py", line 213, in inner
    return Inner(a, b)
  File "/usr/local/lib/python3.10/dist-packages/ufl/tensoralgebra.py", line 166, in __new__
    raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <NablaGrad id=139896952732160> and <ListTensor id=139896712508928>

I don’t understand why I cannot use conjugate explicitly in the variational formulation above. Could you tell me how I should form the variational formulation / Jacobian for the code here? Thank you in advance.

8

Again, could you please write out the strong formulation of your problem, to me I do not understand what PDE would give rise to \int_{\Omega} (\bar{u}\cdot(\nabla u)^T) \cdot \bar{v}~\mathrm{d}x.

9

I am trying to mimic the paper here.

You can see that in the example on the 3rd page of the article, due to the nonlinear convection operator, I need to use the conjguate of my velocity component in the variational formulation.

10

I would start by writing out the analytical expression of the Jacobian that will be used in Newton’s method.