ArityMismatch: Adding expressions with non-matching form arguments ('v_1',) vs ('conj(v_1)',)

General
1

` Hello,

I have a similar issue with the one in this post Problem with use of a ufl function in a NonlinearProblem

I’m trying to solve an eigenvalue problem in FEniCSx with complex-valued functions, and I’m encountering an ArityMismatch error when trying to assemble the discrete matrices from the variational form.

ArityMismatch: Adding expressions with non-matching form arguments (‘v_1’,) vs (‘conj(v_1)’,).

Here, I post the respective code snippet:

v=ufl.TrialFunction(V)
w=ufl.TestFunction(V)   
  
#Terms of left-hand side
grad_term= ufl.dot(grad_r(v, omega, domain), ufl.conj(grad_r(w, omega, domain))) 

pot_term= ufl.dot((potential_r + beta * ufl.inner(phi, phi)) * v,ufl.conj(w))

beta_expression=(ufl.dot(ufl.conj(phi),v))
beta_expression_real= 0.5 * (beta_expression + ufl.conj(beta_expression)) # line that produces the error
beta_term=2*beta*ufl.dot(beta_expression_real*phi,ufl.conj(w))

#Left-hand side
lhs = (grad_term + pot_term + beta_term)*ufl.dx(degree=4)

# Right-hand side
rhs = ufl.dot(v, ufl.conj(w)) * ufl.dx(degree=4)

# Assemble the matrices
form_lhs = dolfinx.fem.form(lhs)
form_rhs = dolfinx.fem.form(rhs)

A=assemble_matrix(form_lhs,[bc])
A.assemble()

B=assemble_matrix(form_rhs,[bc])
B.assemble()

Here, phi is a precomputed wavefunction, beta and omega some constants, and potential_r another interpolated function on the mesh.

The problem is with the computation of beta_expression_real, which causes incompatibility with the form arguments. Alternatively, when I tried intead to extract the real part of beta_expression by using ufl.real(), I encouter the following error:

ArityMismatch: Applying nonlinear operator Real to expression depending on form argument v_1.

Is there any recommended way to ompute the real part of beta expression?

I didn’t post the entire code, as I thought it might be too overwhelming, but in case it is needed let me know.

2

Could you provide a minimal working example up till that line with the error.

For instance, I do not know what is grad_r in your code.

3

Of course. Thank you for your answer. I tried to simplified the code and make it more accesible.

import dolfinx
import math
import numpy as np
import ufl
from dolfinx import mesh, fem
from dolfinx.fem import Function
from mpi4py import MPI
from ufl import SpatialCoordinate, grad
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import  assemble_matrix


#Create the mesh
xmin,ymin=-6,-6
xmax,ymax=6,6
h=0.025
divisions_x=math.ceil((xmax - xmin)/h)
divisions_y=math.ceil((ymax - ymin)/h)
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                               points=[[xmin, ymin], [xmax, ymax]],
                               n=(divisions_x, divisions_y),
                               cell_type=mesh.CellType.triangle)
topology=domain.topology
geometry=domain.geometry


# Get the number of nodes (vertices)
num_nodes = topology.index_map(0).size_local
coords = domain.geometry.x

#Create a function space
V = fem.functionspace(domain, ("Lagrange", 1))

#Insert boundary conditions (homogenous)
phi_D=Function(V,dtype=np.complex128)
phi_D.x.array[:]=ScalarType(0)
tdim = topology.dim
fdim = tdim - 1
topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(phi_D, boundary_dofs)

#Define L2 norm
def l2_norm(phi,quadrature_degree):
    l2_norm = dolfinx.fem.assemble_scalar(
            dolfinx.fem.form(ufl.inner(phi,phi) * ufl.dx(degree=quadrature_degree))
        )
    l2_norm = np.sqrt(l2_norm)
    return l2_norm

#Covariant Gradient
def grad_r(phi,omega,domain):
    x=SpatialCoordinate(domain)
    R_T=ufl.as_vector((-x[1], x[0]))
    return ufl.grad(phi)+1j*(omega/2)*R_T*phi  

#Rotational momentum operator
def Lz(phi,domain): 
    x = SpatialCoordinate(domain)
    Lz_phi = -1j* (x[0] * grad(phi)[1] - x[1] * grad(phi)[0])
    return Lz_phi


#Harmonic trapping Potential
def Vtr(gamma_x, gamma_y,V,omega ):
    trap_potential = Function(V)
    trap_potential.interpolate(
        lambda x: 
        ((gamma_x**2 * x[0]**2 + gamma_y**2 * x[1]**2))-1/4*(omega**2)*(x[0]**2 +  x[1]**2)
    )    
    return trap_potential

#Solution function
def phi_init(omega,V,quadrature_degree):
    phi_0=Function(V,dtype=np.complex128)

    phi_0.interpolate(
        lambda x:
            (1 - omega) * (1/np.sqrt(np.pi)*np.exp(-(x[0]**2+x[1]**2)/2))\
            + omega * (1/np.sqrt(np.pi) * (x[0] * np.exp(-(x[0]**2 + x[1]**2)/2) + 1j * x[1] * np.exp(-(x[0]**2 + x[1]**2)/2)))
    )
    l2_norm_phi_0=l2_norm(phi_0,quadrature_degree)
    phi_0.x.array[:] /= l2_norm_phi_0    

    return phi_0

#Set parameters
omega=1.4
beta=200
gamma_x=1.1
gamma_y=1.0
potential_r= Vtr(gamma_x, gamma_y,V,omega)
potential= Vtr(gamma_x, gamma_y,V,0)

#Set solution function
phi=phi_init(omega,V,quadrature_degree=4)

#Define variational problem
v=ufl.TrialFunction(V)
w=ufl.TestFunction(V)   
  
#Terms of left-hand side
grad_term= ufl.dot(grad_r(v, omega, domain), ufl.conj(grad_r(w, omega, domain))) 

pot_term= ufl.dot((potential_r + beta * ufl.inner(phi, phi)) * v,ufl.conj(w))

beta_expression=(ufl.dot(ufl.conj(phi),v))
beta_expression_real= 0.5 * (beta_expression + ufl.conj(beta_expression)) # line that produces the error
beta_term=2*beta*ufl.dot(beta_expression_real*phi,ufl.conj(w))

#Left-hand side
lhs = (grad_term + pot_term + beta_term)*ufl.dx(degree=4)

# Right-hand side
rhs = ufl.dot(v, ufl.conj(w)) * ufl.dx(degree=4)

# Assemble the matrices
form_lhs = dolfinx.fem.form(lhs)
form_rhs = dolfinx.fem.form(rhs)

A=assemble_matrix(form_lhs,[bc])
A.assemble()

B=assemble_matrix(form_rhs,[bc])
B.assemble()
4

I tried running your code and got the following error:

tmp: for the --non-global-value-max-name-size option: may only occur zero or one times!
Traceback (most recent call last):
  File "/home/marazzato/PostDoc/Coding/arity.py", line 97, in <module>
    phi=phi_init(omega,V,quadrature_degree=4)
  File "/home/marazzato/PostDoc/Coding/arity.py", line 83, in phi_init
    l2_norm_phi_0=l2_norm(phi_0,quadrature_degree)
  File "/home/marazzato/PostDoc/Coding/arity.py", line 47, in l2_norm
    dolfinx.fem.form(ufl.inner(phi,phi) * ufl.dx(degree=quadrature_degree))
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 249, in form
    return _create_form(form)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 244, in _create_form
    return _form(form)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 229, in _form
    f = ftype(
TypeError: __init__(): incompatible function arguments. The following argument types are supported:
    1. __init__(self, spaces: list[dolfinx.cpp.fem.FunctionSpace_float64], integrals: dict[dolfinx.cpp.fem.IntegralType, list[tuple[int, int, ndarray[dtype=int32, writable=False, shape=(*), order='C']]]], coefficients: list[dolfinx.cpp.fem.Function_float64], constants: list[dolfinx.cpp.fem.Constant_float64], need_permutation_data: bool, entity_maps: dict[dolfinx.cpp.mesh.Mesh_float64, ndarray[dtype=int32, writable=False, shape=(*), order='C']], mesh: Optional[dolfinx.cpp.mesh.Mesh_float64]) -> None
    2. __init__(self, form: int, spaces: list[dolfinx.cpp.fem.FunctionSpace_float64], coefficients: list[dolfinx.cpp.fem.Function_float64], constants: list[dolfinx.cpp.fem.Constant_float64], subdomains: dict[dolfinx.cpp.fem.IntegralType, list[tuple[int, ndarray[dtype=int32, writable=False, shape=(*), order='C']]]], entity_maps: dict[dolfinx.cpp.mesh.Mesh_float64, ndarray[dtype=int32, writable=False, shape=(*), order='C']], mesh: Optional[dolfinx.cpp.mesh.Mesh_float64]) -> None

Invoked with types: dolfinx.cpp.fem.Form_float64, _CDataBase, list, list, list, dict, dict, dolfinx.cpp.mesh.Mesh_float64

Could you simplify your code further? Only putting in what creates the error you are asking about?

5

Yes, I can do this. Here’s the minimal code to reproduce the error:

from dolfinx import mesh, fem
from mpi4py import MPI
import ufl
import numpy as np
import dolfinx
from dolfinx.fem.petsc import  assemble_matrix

# Create mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[-6, -6], [6, 6]], [480, 480], cell_type=mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))

# Define complex function
phi = fem.Function(V, dtype=np.complex128)
phi.interpolate(lambda x: 1 + 1j * x[1])

# Variational problem
v, w = ufl.TrialFunction(V), ufl.TestFunction(V)
expression = ufl.dot(v, ufl.conj(phi))
expression_real = 0.5 * (expression + ufl.conj(expression))  # Error source
lhs = ufl.dot(expression_real * phi, ufl.conj(w)) * ufl.dx

# Assemble the matrix
form_lhs = dolfinx.fem.form(lhs)
A=assemble_matrix(form_lhs)
A.assemble()
6

Currently you end up mixing conj(v) and v in the same form (due to expression_real), this is not supported in UFL

What kind of strong form of a PDE are you considering as you need this?

7

This is the strong form of the eigenvalue problem, I am trying to solve:

-\Delta v + \left(V_R + \beta |\phi|^2\right)v + 2\beta \operatorname{Re}(\phi v)\phi = \lambda v ,

where:
\beta a real constant, \phi a pre-computed complex valued function on the mesh and V_R another real valued interpolated function.

The problem is with the computation of

\operatorname{Re}(\phi v).

In the weak form with v being the trial function, ufl.real() results in an error, as nonlinear operator on a form argument and if I extract the real part as

\frac{\phi v + \overline{\phi v}}{2}

I end up mixing conj(v) and v.

8

The problem here doesn’t fall under the natural complex number abstraction in UFL.

Consider the following:

\int_\Omega u \cdot \bar{v}~\mathrm{d}x=\int_\Omega f \cdot \bar{v}~\mathrm{d}c

If we want to solve this with «complex numbers» in dolfinx, we have the following anzats

u=\sum_{i=1}^N c_i \phi_i(x) where c_i=a_i+1j*b_i, a_i, b_i\in\mathbb{R}
and \phi_k:\Omega\mapsto\mathbb{R}.

Inserting this into the above equation, you would get

Ac=B

where

A_{ji}=\int_\Omega \phi_i(x)\phi_j(x)~\mathrm{d}x
b_j = \int_\Omega f \cdot \phi_j~\mathrm{d}x

If you would introduce a term \int_\Omega\mathrm{Re}(u)\cdot\bar{v}~\mathrm{d}x it wouldn’t fit into a row of the aforementioned A as you would only have coefficients a_i, not b_i.

To solve your problem, you could:

  1. split your problem into two real valued function spaces, one for a_i and another one for b_i.
  2. Write out your variational for with this split, and gather all real valued terms in one block of the matrix, all complex valued in another. You can then get a real valued matrix system coupling a and b