Missing bessel functions

[jorgensd]

We currently do not support bessel_I or bessel_K functions from ufl.

import ufl
import basix.ufl

c_el = basix.ufl.element("Lagrange", "triangle", 1, shape=(2,))
mesh = ufl.Mesh(c_el)
x = ufl.SpatialCoordinate(mesh)
J = ufl.bessel_K(1, x[0]) * ufl.dx

and

import ufl
import basix.ufl

c_el = basix.ufl.element("Lagrange", "triangle", 1, shape=(2,))
mesh = ufl.Mesh(c_el)
x = ufl.SpatialCoordinate(mesh)
J = ufl.bessel_I(1, x[0]) * ufl.dx

throws

root@dokken-XPS-9320:~/shared/adios4dolfinx/tests# python3 -m ffcx dd.py 
Traceback (most recent call last):
  File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/root/shared/ffcx/ffcx/__main__.py", line 14, in <module>
    sys.exit(main())
  File "/root/shared/ffcx/ffcx/main.py", line 75, in main
    code_h, code_c = compiler.compile_ufl_objects(
  File "/root/shared/ffcx/ffcx/compiler.py", line 111, in compile_ufl_objects
    code = generate_code(ir, options)
  File "/root/shared/ffcx/ffcx/codegeneration/codegeneration.py", line 53, in generate_code
    code_integrals = [integral_generator(integral_ir, options) for integral_ir in ir.integrals]
  File "/root/shared/ffcx/ffcx/codegeneration/codegeneration.py", line 53, in <listcomp>
    code_integrals = [integral_generator(integral_ir, options) for integral_ir in ir.integrals]
  File "/root/shared/ffcx/ffcx/codegeneration/C/integrals.py", line 40, in generator
    parts = ig.generate()
  File "/root/shared/ffcx/ffcx/codegeneration/integral_generator.py", line 164, in generate
    all_quadparts += self.generate_quadrature_loop(rule)
  File "/root/shared/ffcx/ffcx/codegeneration/integral_generator.py", line 266, in generate_quadrature_loop
    definitions, intermediates_0 = self.generate_varying_partition(quadrature_rule)
  File "/root/shared/ffcx/ffcx/codegeneration/integral_generator.py", line 309, in generate_varying_partition
    return self.generate_partition(arraysymbol, F, "varying", quadrature_rule)
  File "/root/shared/ffcx/ffcx/codegeneration/integral_generator.py", line 342, in generate_partition
    vexpr = L.ufl_to_lnodes(v, *vops)
  File "/root/shared/ffcx/ffcx/codegeneration/lnodes.py", line 1127, in ufl_to_lnodes
    raise RuntimeError(f"Missing lookup for expr type {optype}.")
RuntimeError: Missing lookup for expr type <class 'ufl.mathfunctions.BesselI'>.

Traceback (most recent call last):
  File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/root/shared/ffcx/ffcx/__main__.py", line 14, in <module>
    sys.exit(main())
  File "/root/shared/ffcx/ffcx/main.py", line 75, in main
    code_h, code_c = compiler.compile_ufl_objects(
  File "/root/shared/ffcx/ffcx/compiler.py", line 111, in compile_ufl_objects
    code = generate_code(ir, options)
  File "/root/shared/ffcx/ffcx/codegeneration/codegeneration.py", line 53, in generate_code
    code_integrals = [integral_generator(integral_ir, options) for integral_ir in ir.integrals]
  File "/root/shared/ffcx/ffcx/codegeneration/codegeneration.py", line 53, in <listcomp>
    code_integrals = [integral_generator(integral_ir, options) for integral_ir in ir.integrals]
  File "/root/shared/ffcx/ffcx/codegeneration/C/integrals.py", line 40, in generator
    parts = ig.generate()
  File "/root/shared/ffcx/ffcx/codegeneration/integral_generator.py", line 164, in generate
    all_quadparts += self.generate_quadrature_loop(rule)
  File "/root/shared/ffcx/ffcx/codegeneration/integral_generator.py", line 266, in generate_quadrature_loop
    definitions, intermediates_0 = self.generate_varying_partition(quadrature_rule)
  File "/root/shared/ffcx/ffcx/codegeneration/integral_generator.py", line 309, in generate_varying_partition
    return self.generate_partition(arraysymbol, F, "varying", quadrature_rule)
  File "/root/shared/ffcx/ffcx/codegeneration/integral_generator.py", line 342, in generate_partition
    vexpr = L.ufl_to_lnodes(v, *vops)
  File "/root/shared/ffcx/ffcx/codegeneration/lnodes.py", line 1127, in ufl_to_lnodes
    raise RuntimeError(f"Missing lookup for expr type {optype}.")
RuntimeError: Missing lookup for expr type <class 'ufl.mathfunctions.BesselK'>

[chrisrichardson]

The problem is that these functions are not generally supported in the C math library, also often not implemented in C++, either. We can add a call to something besselK(), but would rely on the user linking an appropriate implementation.

[jorgensd]

The problem is that these functions are not generally supported in the C math library, also often not implemented in C++, either. We can add a call to something besselK(), but would rely on the user linking an appropriate implementation.

They can be derived from the besselJ and besselY functions https://en.m.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions so we could implement them as expressions of those?

[chrisrichardson]

https://www.gnu.org/software/libc/manual/html_node/Special-Functions.html

Generally they only operate on real values, so I think that won't help. Probably the table for complex in ffcx is wrong, here.

[jorgensd]

Then the question is, should they exist in ufl? firedrake doesn’t seem to support them https://github.com/firedrakeproject/tsfc/blob/90c20c5b06c31beb388e0f4c874dbd2dc4f80fcf/tsfc/loopy.py#L436-L451 and we dont since they arent part of C. @dham @mscroggs your thoughts?

[jorgensd]

As the aim is to be backend (language) agnostic they would be supported if we generated C++: https://en.cppreference.com/w/cpp/numeric/special_functions/cyl_bessel_i

[dham]

I think they're very niche and could be done as an external operator if someone really needed them.

[jhale]

I agree with @dham, we're around 15 years in and no one has noticed they aren't implemented, implies they are very niche.

[jorgensd]

@jhale not that I disagree that they are niche, but the issue was reported here: https://fenicsproject.discourse.group/t/defining-a-fem-function-using-a-conditional-ufl-statement/13730/3?u=dokken and the bessel_i and bessel_k functions do exist in legacy dolfin.

I think for the reported problem, one could probably just use: https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.iv.html#scipy.special.iv in interpolation with a lambda function