Bessel function and complex number

General
1

Hello,

It took me quite a while to understand why one of my code was not working. I managed to find the source of the issue. It is quite unclear to me if the issue comes from ufl or the use of create_real_functionspace but it seems that bessel function can not be used with complex argument. Here is a MWE:

import dolfinx
import ufl
from scifem import create_real_functionspace
from mpi4py import MPI
from scipy.special import jv as jv_scipy

msh = dolfinx.mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0., 0.), (1., 1.)), n=(1, 1), cell_type=dolfinx.mesh.CellType.triangle)
R = create_real_functionspace(msh)
x = dolfinx.fem.Function(R)
f = dolfinx.fem.Function(R)

x_value = 12.0j


x = dolfinx.fem.Function(R)
x.x.array[:] = x_value
x.x.scatter_forward()

def my_function(x):
    jv = ufl.bessel_J if isinstance(x, ufl.core.expr.Expr) else jv_scipy
    return jv(0, x)


f.interpolate(dolfinx.fem.Expression(my_function(x), R.element.interpolation_points))

print(f.x.array[:],' / ',my_function(x_value)) # [1.+0.j]  /  (18948.925349296307+0j)

We see that the value is not the same as the one returned by the scipy implementation of the Bessel function.

I’m quite unsure whether I should post that to the ufl issue tracker, let me know if that is the case.

2

Unfortunately, it looks like the standard C library does not have an implementation of bessel functions of the first kind for complex arguments: IBM Documentation

You can try this with something like

#include <stdio.h>
#include <math.h>
#include <complex.h>

int main() {
    double complex x = 12.0*I;
    int n = 0;
    double complex result = jn(n, x);
    printf("J%d(%.1f + %.1fi) = %.1f + %.1fi\n", n, x, result);
    return 0;
}

Which gives me

J0(0.0 + 12.0i) = 1.0 + 0.0i
3

One workaround may be to use GitHub - a-latyshev/dolfinx-external-operator: Extension of DOLFINx implementing the concept of external operator · GitHub and have it point to the scipy implementation for complex numbers.

4

Thank you for this solution, it looks like it could solve my problem.

5

To conclude this post, here is an example using External operator, and comparing it to assigning data directly to a DG-0 function:

from mpi4py import MPI

import basix.ufl
import dolfinx
import ufl
from dolfinx_external_operator import (
    FEMExternalOperator,
    evaluate_external_operators,
    replace_external_operators,
    evaluate_operands,
)
import numpy as np
from scipy.special import jv as jv_scipy


def g_impl(uuh):
    """Identity operator"""
    output = jv_scipy(0, uuh)
    return output.reshape(-1)


def g_external(derivatives):
    if derivatives == (0,):  # no derivation, the function itself
        return g_impl
    else:
        return NotImplementedError


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

quadrature_degree = 1
Qe = basix.ufl.quadrature_element(mesh.basix_cell(), degree=quadrature_degree, value_shape=())
Q = dolfinx.fem.functionspace(mesh, Qe)
coefficient = dolfinx.fem.Function(Q, dtype=np.complex128)

x_value = np.complex128(12.0j)
x = dolfinx.fem.Constant(
    mesh,
    x_value,
)


f = FEMExternalOperator(x, function_space=Q, external_function=g_external, coefficient=coefficient)

dx = ufl.Measure("dx", domain=mesh)

Q = dolfinx.fem.functionspace(mesh, ("DG", 0))
g = dolfinx.fem.Function(Q, dtype=np.complex128)
g.x.array[:] = g_impl(x_value)

a = (f - g) * ufl.conj(f - g) * dx

a_replaced, a_external_operators = replace_external_operators(a)
# Pack coefficients for g
evaluated_operands = evaluate_operands(a_external_operators)
evaluate_external_operators(a_external_operators, evaluated_operands)

error = dolfinx.fem.assemble_scalar(dolfinx.fem.form(a_replaced, dtype=np.complex128))
print("Error:", np.sqrt(error))

output

Error: 0j