How to use bessel functions in Expression

General
1

I am using legacy FEniCS 2019.1.0 and I need to find the eigenvalues and eigenfunctions of a Matern covariance operator, which is defined using modified bessel functions. However, the following mwe, defining an expression using bessel functions, does not compile:

from fenics import *
import numpy as np

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20)
V = FunctionSpace(mesh, ‘P’, 1)

u = TrialFunction(V)
v = TestFunction(V)

k = Expression(‘pow(sqrt(pow(x[0]-x_i,2)+pow(x[1]-y_i,2)),nu)*cyl_bessel_k(nu, sqrt(pow(x[0]-x_i,2)+pow(x[1]-y_i,2)))’, degree=1, nu=1, x_i=0, y_i=0)

c = np.zeros((mesh.num_vertices(), mesh.num_vertices()))
coords = mesh.coordinates()
dof2v = dof_to_vertex_map(V)
for i in range(len(coords)):
k.x_i = Vertex(mesh, dof2v[i]).midpoint()[0]
k.y_i = Vertex(mesh, dof2v[i]).midpoint()[1]
c[i, :] = assemble(k * v * dx)

vals, vecs = np.linalg.eig(c)

The error message is the following:
/tmp/tmp6vvyh3xk/dolfin_expression_4a2c2db2129efa596eb787cd10ed666d.cpp: In member function ‘virtual void dolfin::dolfin_expression_4a2c2db2129efa596eb787cd10ed666d::eval(Eigen::Ref<Eigen::Matrix<double, -1, 1> >, Eigen::Ref<const Eigen::Matrix<double, -1, 1> >) const’:
/tmp/tmp6vvyh3xk/dolfin_expression_4a2c2db2129efa596eb787cd10ed666d.cpp:64:69: error: ‘cyl_bessel_k’ was not declared in this scope

Apparently cyl_bessel_k is not declared even though it is listed in the documentation of cmath: Standard library header <cmath> - cppreference.com

I have also looked into this (old) solution but the code provided there does not compile either, giving an extremely long error message, mostly repeating error: ‘namespace’ definition is not allowed here https://fenicsproject.org/qa/3498/use-bessel-functions-in-jit-compiled-expressions/

Is there a different way of approaching this problem? Or am I just using the wrong alias for bessel functions?

2

I would use the ufl formulation of bessel functions directly, instead of trying to generate code for it.
i.e.

from fenics import *

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20)

x = SpatialCoordinate(mesh)
nu = 1
x_i, y_i = 0,0
k = sqrt((x[0]-x_i)**2+(x[1]-y_i)**2)**nu*bessel_K(nu, sqrt((x[0]-x_i)**2+(x[1]-y_i)**2))
print(assemble(k*dx))

This gives you a more exact approximation of your function as it is evaluated at the quadrature points, and not interpolated into a Lagrange space of order degree (1 in your case).

3

I had already tried a similar approach using the ufl formulation:

from fenics import *
import numpy as np
import time

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20)
V = FunctionSpace(mesh,‘P’,1)
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)

def materncov(x, nu, x_i, y_i):
k = sqrt((x[0]-x_i)**2+(x[1]-y_i)**2)**nu*bessel_K(nu, sqrt((x[0]-x_i)**2+(x[1]-y_i)**2))
return k

c = np.zeros((mesh.num_vertices(), mesh.num_vertices()))
coords = mesh.coordinates()
dof2v = dof_to_vertex_map(V)
for i in range(len(coords)):
starttime = time.time()
x_i = Vertex(mesh, dof2v[i]).midpoint()[0]
y_i = Vertex(mesh, dof2v[i]).midpoint()[1]
c[i, :] = assemble(materncov(x, 1, x_i, y_i) * v * dx)
print(time.time() - starttime)

vals, vecs = np.linalg.eig(c)

While that does work for one iteration, in order to find the eigenfunctions I need to compare each pair of points in the mesh, so x_i,y_i have to loop over all vertices. If I have to use assemble in each iteration of this loop, this works quickly for the first few iterations (taking ~3ms/iteration), but after a handful of iterations this code starts calling the JIT compiler in every iteration, taking 2+ secs per iteration. Hence my original plan of using Expressions, since then I could just update the parameters x_i,y_i and not have to call the compiler every time. Is there a fix to this?

4

You should use Constant for x_i and y_i to avoid recompilation:

from fenics import *
import time
import numpy as np
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20)
V = FunctionSpace(mesh, "Lagrange", 1)
x = SpatialCoordinate(mesh)
nu = 1
x_i, y_i = Constant(0),Constant(0)
k = sqrt((x[0]-x_i)**2+(x[1]-y_i)**2)**nu*bessel_K(nu, sqrt((x[0]-x_i)**2+(x[1]-y_i)**2))
v = TestFunction(V)
res = k*v*dx 


c = np.zeros((mesh.num_vertices(), mesh.num_vertices()))
coords = mesh.coordinates()
dof2v = dof_to_vertex_map(V)
for i in range(len(coords)):
    starttime = time.time()
    mp =  Vertex(mesh, dof2v[i]).midpoint()
    x_i.assign(mp[0])
    y_i.assign(mp[1])
    c[i, :] = assemble(res)
    print(time.time() - starttime)

Reduces the runtime of assembly to ~0.003 seconds per vertex

EDIT
If you tested this with the previous approach, I would strongly suggest to clean your cache with dijitso clean to avoid having tons of pre-compiled forms to search through.