Use of symmetric tensor in dolfinx/ufl results in error

Errors
1

Hi,
Currently I am porting an old application code from Fenics (2019.1.0) to new Fenicsx.

The versions of the libraries I used are as follows:
dolfinx : 0.7.2
ufl: 2023.2.0
ffcx: 0.7.0
I compiled and built these programs from source which I obtained from Fenics git repository.

In my code I would like to use a symmetric trace-free tensor to represent stress tensor. UFL library does allow symmetric tensor flag. When I execute the following MWE code

import dolfinx as dfx
from mpi4py import MPI
import ufl as ufl

mesh_dfx = dfx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 8, 8,
        dfx.mesh.CellType.tetrahedron
        ) 

# create 2 finite elements
# Velocity u (vector),
u_elem = ufl.VectorElement("Lagrange", mesh_dfx.ufl_cell(), 1)
# Stress sigma (Tensor)
# Variable sigma should be symmetric and trace-free. 
# UFL library allows me to create symmetric tensor with a flag `symmetric`. 
#NOTE: With symmetry - CODE RETURNS AN ERROR
sigma_elem = ufl.TensorElement("Lagrange", mesh_dfx.ufl_cell(),1, symmetry=True)
# Without symmetry - SEEMS TO RUN WITHOUT ERROR
#sigma_elem = ufl.TensorElement("Lagrange", mesh_dfx.ufl_cell(),1)

# create mixed finite element vector space
mxd_fspace = dfx.fem.FunctionSpace(mesh_dfx,
        ufl.MixedElement([u_elem, sigma_elem ])
        )

# setup trial and test functions
(u, sigma) = ufl.TrialFunctions(mxd_fspace)
(v, psi ) = ufl.TestFunctions(mxd_fspace)

# Convert symmetric tensor to symmetric-and-trace-free tensor 
# using the following custom function
def gen3DTFdim3(sigma):
    return ufl.as_tensor([
        [sigma[0, 0], sigma[0, 1], sigma[0, 2]],
        [sigma[1, 0], sigma[1, 1], sigma[1, 2]],
        [sigma[2, 0], sigma[2, 1], -sigma[0, 0] - sigma[1, 1]]]
        )
sigma = gen3DTFdim3(sigma)
psi = gen3DTFdim3(psi)

# For the purpose of demonstration I include only one term 
# from the actual weak formulation of my problem.
a = ufl.dot(ufl.div(sigma), v) * ufl.dx 
print(dfx.fem.assemble_scalar(dfx.fem.form(a)))

results in the following error

Traceback (most recent call last):
  File "/home/sysadmin/py/fenicsxR13/demo_symTensor_unitCube.py", line 43, in <module>
    print(dfx.fem.assemble_scalar(dfx.fem.form(a)))
  File "/home/sysadmin/.local/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 188, in form
    return _create_form(form)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 183, in _create_form
    return _form(form)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 141, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
  File "/home/sysadmin/.local/lib/python3.10/site-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 199, in compile_forms
    raise e
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 190, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/compiler.py", line 102, in compile_ufl_objects
    ir = compute_ir(analysis, object_names, prefix, options, visualise)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/ir/representation.py", line 197, in compute_ir
    irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/ir/representation.py", line 197, in <listcomp>
    irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/ir/representation.py", line 492, in _compute_integral_ir
    integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/ir/integral.py", line 72, in compute_integral_ir
    S = build_scalar_graph(expression)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/ir/analysis/graph.py", line 80, in build_scalar_graph
    scalar_expressions = rebuild_with_scalar_subexpressions(G)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/ir/analysis/graph.py", line 125, in rebuild_with_scalar_subexpressions
    V_symbols = value_numberer.compute_symbols()
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/ir/analysis/valuenumbering.py", line 68, in compute_symbols
    symbol = f(expr)
  File "/home/sysadmin/.local/lib/python3.10/site-packages/ffcx/ir/analysis/valuenumbering.py", line 184, in _modified_terminal
    raise RuntimeError("Internal error in value numbering.")
RuntimeError: Internal error in value numbering.

I have read few posts in this forum regarding the erroneous behaviour of ufl when symmetric flag used. However, the error message mentioned in each post seems to be different from one another. Therefore, I’ve posted mine here.

I would like to know how can I fix this problem? Can I use any other version of dolfinx/ufl/ffcx combination to circumvent this error?

Thanks for reading my post!

2

Since you are building from source anyways, can you try building basix, ffcx and dolfinx from the current main branch?

Fix __hash__ and __eq__ for basix.ufl elements by conpierce8 · Pull Request #718 · FEniCS/basix · GitHub might have fixed this (see also [BUG]: Symmetric spaces broken · Issue #2750 · FEniCS/dolfinx · GitHub). Even if it that does not actually fix this, you can use the updated code to open an issue at Issues · FEniCS/basix · GitHub (please do not open an issue with this code at your currently installed version, since you would immediately be asked to try with main, because we are aware that this part was changed very recently)

While updating you will need to replace ufl.VectorElement/ufl.TensorElement creation with basix elements, since former ufl ones have been removed from main. Something like

basix.ufl.element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
basix.ufl.element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,mesh.geometry.dim), symmetry=True_or_False)

should be the way to do it.

3

Thank you for your prompt reply!
Following your recommendation, I built the dolfinx/ffcx/basix/ufl sources from each of their main repository.
The version numbers are:
dolfinx: 0.8.0.0
ffcx: 0.8.0.dev0
basix: 0.8.0.0
ufl: 2023.3.0.dev0
Thanks to your suggestion I modified the code. When I run the following code

import dolfinx as dfx
from mpi4py import MPI
import ufl as ufl
import basix as basix
# create mesh
mesh_dfx = dfx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 8, 8,
        dfx.mesh.CellType.tetrahedron) 
# define two finite elments: a vector and a symmetric tensor
u_elem = basix.ufl.element("Lagrange", mesh_dfx.basix_cell(), 1,
        shape=(mesh_dfx.geometry.dim,))
sigma_elem = basix.ufl.element("Lagrange", mesh_dfx.basix_cell(),1, 
        shape=(mesh_dfx.geometry.dim, mesh_dfx.geometry.dim), 
        symmetry=True)
# define mixed finite element function space
mxd_fspace = dfx.fem.functionspace(mesh_dfx,
        basix.ufl.mixed_element([u_elem, sigma_elem ]))
# define trial and test functions
(u, sigma) = ufl.TrialFunctions(mxd_fspace)
(v, psi ) = ufl.TestFunctions(mxd_fspace)
# define a form expression
a = ufl.dot(ufl.div(sigma), v) * ufl.dx 
# cast it explicitly to form type
compiled_form = dfx.fem.form(a)
# assemble  --> NOTE: the following line produces error
print(dfx.fem.assemble_scalar(compiled_form))

, it produces the following error:

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: ---------------------  Stack Frames ------------------------------------
[0]PETSC ERROR: No error traceback is available, the problem could be in the main program. 
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 59.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

It seems to me that because of the new modifications there is something going on under
the hood which I am not able to understand. Could you please verify if the above code snippet compiles correctly for you?

4

This does not make sense as a command.
The line that produces an error tries to assemble a matrix into a single scalar value.

Could you try with
dolfinx.fem.petsc.assemble_matrix(compiled_form)?

5

Hi Dokken,
Thanks for pointing out the error. Indeed, with dolfinx.fem.petsc.assemble_matrix() the code does run without error.