Arity Errors when Running Code

Hello everyone,

I am working on a fluid-structure interface problem and I have managed to define my weak form as follows. I have the following problem:

# Solid weak form (on subdomain 1)
a_solid = ufl.inner(sigma_s(u_s), eps(v_s)) * dx(1)

# Fluid weak form (on subdomain 2)
a_fluid = (ufl.real(ufl.inner(sigma_f(u_f, p_f), eps(v_f))) * dx(2)
           - ufl.div(v_f) * p_f * dx(2)
           - q_f * ufl.div(u_f) * dx(2))

# Nitsche interface terms (on ds(3))
nitsche = (
    - ufl.inner(sigma_s(u_s)*n, jump_v) * ds(3)
    - ufl.inner(sigma_f(u_f, p_f)*n, jump_v) * ds(3)
    - ufl.inner(jump_u, sigma_s(v_s)*n) * ds(3)
    - ufl.inner(jump_u, sigma_f(v_f, 0)*n) * ds(3) 
    + gamma * ufl.inner(jump_u, jump_v) * ds(3)
)

a = a_solid + a_fluid + nitsche

# Right-hand side: zero source terms
zero_vector = fem.Constant(domain, (0.0, 0.0))
L = ufl.inner(zero_vector, v_s) * dx + ufl.inner(zero_vector, v_f) * dx

where each of the terms are defined as follows:

jump_u = ufl.as_vector([u_s[i] - u_f[i] for i in range(dim)])
jump_v = ufl.as_vector([v_s[i] - v_f[i] for i in range(dim)])

def eps(u): 
    return ufl.sym(ufl.grad(u))

def sigma_s(u): 
    return lmbda_s * ufl.tr(eps(u)) * ufl.Identity(dim) + 2 * mu_s * eps(u)

def sigma_f(u, p): 
    return 2 * mu_f * eps(u) - p * ufl.Identity(dim)

n = ufl.FacetNormal(domain)

I’m getting an Arity Error when I try to solve this; any suggestions on where to debug? (I’m using the zero vector for testing and running purposes).

Please make the code snippets into a reproducible example, ie, use a basic mesh (unit cube or square, define the relevant function spaces, functions, testfunctions, trialfunctions and constants required to run your code up to the point of calling dolfinx.fem.form.

Thank you for the response. Here is a simplified version of the code that I am attempting to run.

import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, io
from dolfinx.mesh import locate_entities, locate_entities_boundary
import ufl
import basix.ufl
from petsc4py import PETSc
import os

# Create mesh and tag subdomains
domain = mesh.create_rectangle(MPI.COMM_WORLD,
                               [np.array([0.0, 0.0]), np.array([1.0, 1.0])],
                               [32, 32],
                               cell_type=mesh.CellType.triangle)

# Subdomain tagging: 1 = solid, 2 = fluid
def solid_cells(x): 
    return x[0] <= 0.5

def fluid_cells(x): 
    return x[0] > 0.5

solid_cells_indices = locate_entities(domain, domain.topology.dim, solid_cells)
fluid_cells_indices = locate_entities(domain, domain.topology.dim, fluid_cells)

cell_tags = np.zeros(domain.topology.index_map(domain.topology.dim).size_local, dtype=np.int32)
cell_tags[solid_cells_indices] = 1
cell_tags[fluid_cells_indices] = 2
domain.topology.create_connectivity(domain.topology.dim, 0)
cell_tag = mesh.meshtags(domain, domain.topology.dim, np.arange(len(cell_tags)), cell_tags)

dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_tag)

# Mark interface facets for Nitsche terms
fdim = domain.topology.dim - 1

def interface(x): 
    return np.isclose(x[0], 0.5)

interface_facets = locate_entities_boundary(domain, fdim, interface)
interface_values = np.full(len(interface_facets), 3, dtype=np.int32)
facet_tag = mesh.meshtags(domain, fdim, interface_facets, interface_values)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)


# Define finite elements 
degree = 2
cell_type = domain.basix_cell()

lagrange_s = basix.ufl.element("Lagrange", cell_type, degree, shape=(domain.topology.dim,))
lagrange_f = basix.ufl.element("Lagrange", cell_type, degree, shape=(domain.topology.dim,))
# - For fluid pressure, a scalar element:
lagrange_p = basix.ufl.element("Lagrange", cell_type, degree - 1)

# Create a mixed element: ordering: [solid displacement, fluid velocity, fluid pressure]
mixed_element = basix.ufl.mixed_element([lagrange_s, lagrange_f, lagrange_p])


# Create function spaces

V_s = fem.functionspace(domain, lagrange_s)      # Solid displacement space
V_f = fem.functionspace(domain, lagrange_f)      # Fluid velocity space
Q_f = fem.functionspace(domain, lagrange_p)      # Fluid pressure space
W = fem.functionspace(domain, mixed_element)     # Mixed space


# Define trial and test functions, and split components

u = ufl.TrialFunction(W)
v = ufl.TestFunction(W)
u_s, u_f, p_f = ufl.split(u)
v_s, v_f, q_f = ufl.split(v)

# Build explicit jump terms so that components are formed as standard UFL vectors.
dim = domain.topology.dim
jump_u = ufl.as_vector([u_s[i] - u_f[i] for i in range(dim)])
jump_v = ufl.as_vector([v_s[i] - v_f[i] for i in range(dim)])

# Material parameters and constitutive functions
mu_s = 1.0
lmbda_s = 1.0
mu_f = 1.0
gamma = 10 * degree**2  # Nitsche penalty

def eps(u): 
    return ufl.sym(ufl.grad(u))

def sigma_s(u): 
    return lmbda_s * ufl.tr(eps(u)) * ufl.Identity(dim) + 2 * mu_s * eps(u)

def sigma_f(u, p): 
    return 2 * mu_f * eps(u) - p * ufl.Identity(dim)

n = ufl.FacetNormal(domain)


# Define variational forms

# Solid weak form (on subdomain 1)
a_solid = ufl.inner(sigma_s(u_s), eps(v_s)) * dx(1)

# Fluid weak form (on subdomain 2)
a_fluid = (ufl.real(ufl.inner(sigma_f(u_f, p_f), eps(v_f))) * dx(2)
           - ufl.div(v_f) * p_f * dx(2)
           - q_f * ufl.div(u_f) * dx(2))

# Nitsche interface terms (on ds(3))
nitsche = (
    - ufl.inner(sigma_s(u_s)*n, jump_v) * ds(3)
    - ufl.inner(sigma_f(u_f, p_f)*n, jump_v) * ds(3)
    - ufl.inner(jump_u, sigma_s(v_s)*n) * ds(3)
    - ufl.inner(jump_u, sigma_f(v_f, 0)*n) * ds(3) 
    + gamma * ufl.inner(jump_u, jump_v) * ds(3)
)

a = a_solid + a_fluid + nitsche

# Right-hand side: zero source terms
zero_vector = fem.Constant(domain, (0.0, 0.0))
L = ufl.inner(zero_vector, v_s) * dx + ufl.inner(zero_vector, v_f) * dx

# Apply boundary conditions
def left(x): 
    return np.isclose(x[0], 0.0)
def right(x): 
    return np.isclose(x[0], 1.0)
def walls(x): 
    return np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)

left_facets = locate_entities_boundary(domain, fdim, left)
right_facets = locate_entities_boundary(domain, fdim, right)
wall_facets = locate_entities_boundary(domain, fdim, walls)

# Map BCs onto subspaces of the mixed space
left_dofs = fem.locate_dofs_topological((W.sub(0), V_s), fdim, left_facets)
wall_dofs = fem.locate_dofs_topological((W.sub(1), V_f), fdim, wall_facets)
right_p_dofs = fem.locate_dofs_topological((W.sub(2), Q_f), fdim, right_facets)

u_solid_bc = fem.Function(V_s)
u_solid_bc.x.array[:] = 0.0
u_fluid_bc = fem.Function(V_f)
u_fluid_bc.x.array[:] = 0.0
p_out_bc = fem.Function(Q_f)
p_out_bc.x.array[:] = 0.0

bc_solid_fixed = fem.dirichletbc(u_solid_bc, left_dofs, W.sub(0))
bc_noslip = fem.dirichletbc(u_fluid_bc, wall_dofs, W.sub(1))
bc_pressure_outlet = fem.dirichletbc(p_out_bc, right_p_dofs, W.sub(2))

bcs = [bc_solid_fixed, bc_noslip, bc_pressure_outlet]


# Assemble and solve
problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "minres", "pc_type": "lu"})
U = problem.solve()

I am not getting any arity-errors when running your example, but:

raceback (most recent call last):
  File "/root/shared/mwe.py", line 185, in <module>
    U = problem.solve()
        ^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/petsc.py", line 856, in solve
    self._solver.solve(self._b, self._x)
  File "petsc4py/PETSc/KSP.pyx", line 1782, in petsc4py.PETSc.KSP.solve
petsc4py.PETSc.Error: error code 73
[0] KSPSolve() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:1075
[0] KSPSolve_Private() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:826
[0] KSPSetUp() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:415
[0] PCSetUp() at /usr/local/petsc/src/ksp/pc/interface/precon.c:1071
[0] PCSetUp_LU() at /usr/local/petsc/src/ksp/pc/impls/factor/lu/lu.c:87
[0] MatLUFactorSymbolic() at /usr/local/petsc/src/mat/interface/matrix.c:3253
[0] MatLUFactorSymbolic_SeqAIJ() at /usr/local/petsc/src/mat/impls/aij/seq/aijfact.c:305
[0] Object is in wrong state
[0] Matrix is missing diagonal entry 2317

To me, it seems like you should use sub-meshes, as you do not define your variational form over the whole domain for u_s and u_f. See for instance: Multiphysics: Solving PDEs on subdomains — FEniCS Workshop

1 Like