Dolfin randomly crashes when using interpolated boundary data

[acham]

Hello. I keep encountering a weird bug that I cannot explain or resolve. I am trying to model Stokes flow in a rectangular domain with a circular cutout. The boundary condition is Dirichlet, and inferred from experimental data that needs to be interpolated. I do this by constructing a UserExpression that contains two SciPy RectBivariateSpline objects that return the interpolated velocity components at any point in the domain. It is necessary to perform cubic 2D interpolation of the data because in my use case, the position of the inner circular cut out is not known a priori, and when post processing the data I need to calculate derivatives so the boundary data needs to be smooth.

The problem I encounter is that dolfin 2019.1.0 randomly crashes with the “Missing eval() function (must be overloaded)” error when I do this. It most often (but not always) does so after restarting the kernel, and often also on repeated execution. Sometimes the problem resolves itself when I just try to re-run it, but often it doesn’t and I can’t reliably use it.

Based on THIS topic I suspect it has something to do with the way interpolation is handled, but I’m not really sure how to work around it given the constraints mentioned above.

I’d be really grateful for an explanation of the error and suggestions how to resolve it! MWE below, with an example data file to execute it HERE: EDIT see comment below for updated MWE without the need to download a file.

EDIT: condensed the MWE slightly. Also to be clear, the error occurs in the line ‘solve(a == L, w, bcs)’:

EDIT2: condensed the MWE even further, the error also occurs if the markers of the inner and outer boundary are not distinct.

*** Error:   Unable to evaluate expression.
*** Reason:  Missing eval() function (must be overloaded).
*** Where:   This error was encountered inside Expression.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f

#remove boundary distinctions
from dolfin import *
import mshr
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as interp
import pandas as pd

'''Classes'''

class Velocity_from_PIV(UserExpression):
    # Fenics object to evaluate/interpolate PIV velocity data anywhere
    def __init__(self, data, **kwargs):
        super().__init__(kwargs)
        self.x=data['x'].drop_duplicates().values #vector of x coords (ascending)
        self.y=data['y'].drop_duplicates().values #vector of y coords (ascending)
        self.dx=data['dx'].values.reshape(self.x.size,self.y.size) #matrix of values
        self.dy=data['dy'].values.reshape(self.x.size,self.y.size)
        #Splines (default: cubic, no smoothing)
        self.splx = interp.RectBivariateSpline(self.x,self.y,self.dx)
        self.sply = interp.RectBivariateSpline(self.x,self.y,self.dy)
    def eval_cell(self, value, x, ufc_cell):        
        value[0] = self.splx.ev(x[0],x[1])
        value[1] = self.sply.ev(x[0],x[1])
    def value_shape(self):
        return (2,)

'''Methods''' 

def main():
    #load data
    data = pd.read_csv('Data.txt',sep='\t',names=['i','j','x','y','dx','dy'])

    #Calculate mesh
    mesh, boundary_markers = meshing()
    
    plot(mesh)

    #Define function space and boundary conditions
    W, bcs = defineBCs(mesh,boundary_markers,data)

    #Perform FEM
    u, p = calc_fem(W,bcs) 
    
    return u, p

def meshing():
    """Create the mesh and determine inner and outer boundaries"""    
    # For sake of MWE define rectangular dimensions manually:
    xmin, xmax, ymin, ymax, PIVres = np.array([0.0, 0.7424643186232591, 0.0, 0.39774874211960304, 0.02651658280797354])
    # x,y coords and radius of circular cutout:
    centrex, centrey, radius = np.array([0.27,0.163,0.11])
    
    #define interior and exterior boundary
    outerb = mshr.Rectangle(Point(xmin,ymin),Point(xmax,ymax))
    innerb = mshr.Circle(Point(centrex,centrey),radius)
    
    mesh = mshr.generate_mesh(outerb-innerb,1/PIVres)
    
    #label the boundary elements
    boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    boundary_markers.set_all(0)
    
    return mesh, boundary_markers

def defineBCs(mesh,boundary_markers,data):
    '''define FEM space, BCs'''   
    
    #Create fenics UserExpression to access data
    data_fenics = Velocity_from_PIV(data)
    
    # Define function space, choose Taylor-Hood elements
    P2 = VectorElement("Lagrange", mesh.ufl_cell(), 3)
    P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
    TH = P2 * P1
    W = FunctionSpace(mesh, TH)

    #Define Dirichlet BC on each part of the boundary
    bcs = DirichletBC(W.sub(0), data_fenics, boundary_markers,0)
    
    return W, bcs

def calc_fem(W,bcs):
    """obtain u and p from Stokes FEM"""
    
    # Define variational problem
    (u, p) = TrialFunctions(W)
    (v, q) = TestFunctions(W)
    
    #Stokes equations in weak form
    a = (inner(grad(u), grad(v)) - div(v)*p-div(u)*q -(1E-10)*p*q)*dx
    L = -Constant(0)*q*dx
    
    # Compute solution
    w = Function(W)
    solve(a == L, w, bcs)   

    #extract velocity and pressure fields
    (u, p) = w.split()
    return u, p

u, p = main()

[Emanuel_Gebauer]

Hey,

I think the problem lies here:

@acham

 def eval_cell(self, value, x, ufc_cell):        
        value[0] = self.splx.ev(x[0],x[1])
        value[1] = self.sply.ev(x[0],x[1])

I could not test your problem as it loads a file and I currently dont have pandas installed. However I had a similar problem before: I think the DCBoundray calls the function eval(value, x) and not eval_cell(value, x, ufc_cell). Therefore chanigng your UserExpression to

class Velocity_from_PIV(UserExpression):
    # Fenics object to evaluate/interpolate PIV velocity data anywhere
    def __init__(self, data, **kwargs):
        super().__init__(kwargs)
        self.x=data['x'].drop_duplicates().values #vector of x coords (ascending)
        self.y=data['y'].drop_duplicates().values #vector of y coords (ascending)
        self.dx=data['dx'].values.reshape(self.x.size,self.y.size) #matrix of values
        self.dy=data['dy'].values.reshape(self.x.size,self.y.size)
        #Splines (default: cubic, no smoothing)
        self.splx = interp.RectBivariateSpline(self.x,self.y,self.dx)
        self.sply = interp.RectBivariateSpline(self.x,self.y,self.dy)
    def eval(self, value, x):        
        value[0] = self.splx.ev(x[0],x[1])
        value[1] = self.sply.ev(x[0],x[1])
    def value_shape(self):
        return (2,)

could solve your problem

Vage explanation: Internally (if not overloaded) the function eval_cell(value, x, ufc_cell) calls eval(value, x). However not the other way around, therefore eval(value, x) does not call eval_cell(value, x, ufc_cell) and I also have not found a way to do that (You also linked the topic). In your case in your userExpression you don’t need the ufc_cell variable in your eval_cell function, therefore I’d just leave it out. If you want to evaluate a MeshFunction in your eval_cell function (then you need the ufc_cell argument), then (I think) you first have to interpolate your MeshFunction, therefore creating a function, which then can be called with a coordintate. See: Error: Problem with Expression inside an Expression: Missing eval() function (must be overloaded) - #2 by Emanuel_Gebauer

Best Regards

Emanuel

[acham]

Hi Emanuel thanks for your suggestion Unfortunately the error occurs (again, randomly) regardless of whether I define eval, eval_cell or both… it’s really annoying, sometimes the identical piece of code works fine 10 times in a row and sometimes it will crash forever.

I think the file may not be necessary in the end. Here is a modified MWE that just uses a randomly generated matrix as input.

Thanks for giving my issue thought!

from dolfin import *
import mshr
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as interp
#import pandas as pd

'''Classes'''

class Velocity_from_PIV(UserExpression):
    # Fenics object to evaluate/interpolate PIV velocity data anywhere
    def __init__(self, **kwargs):
        super().__init__(kwargs)
        self.x=np.linspace(0,1,20) #vector of x coords (ascending)
        self.y=np.linspace(0,1,20) #vector of y coords (ascending)
        self.dx=np.random.rand(20,20) #matrix of values
        self.dy=np.random.rand(20,20)
        #Splines (default: cubic, no smoothing)
        self.splx = interp.RectBivariateSpline(self.x,self.y,self.dx)
        self.sply = interp.RectBivariateSpline(self.x,self.y,self.dy)
    def eval_cell(self, value, x, ufc_cell):        
        value[0] = self.splx.ev(x[0],x[1])
        value[1] = self.sply.ev(x[0],x[1])
    def value_shape(self):
        return (2,)

'''Methods''' 

def main():
    #load data
    #data = pd.read_csv('Data.txt',sep='\t',names=['i','j','x','y','dx','dy'])

    #Calculate mesh
    mesh, boundary_markers = meshing()
    
    plot(mesh)

    #Define function space and boundary conditions
    W, bcs = defineBCs(mesh,boundary_markers)

    #Perform FEM
    u, p = calc_fem(W,bcs) 
    
    return u, p

def meshing():
    """Create the mesh and determine inner and outer boundaries"""    
    # For sake of MWE define rectangular dimensions manually:
    xmin, xmax, ymin, ymax, PIVres = np.array([0.0, 1.0, 0.0, 1.0, 0.1])
    # x,y coords and radius of circular cutout:
    centrex, centrey, radius = np.array([0.27,0.163,0.11])
    
    #define interior and exterior boundary
    outerb = mshr.Rectangle(Point(xmin,ymin),Point(xmax,ymax))
    innerb = mshr.Circle(Point(centrex,centrey),radius)
    
    #Define classes to let fenics know whether you are on a certain section of the boundary
    class Boundary_outer(SubDomain):
            def inside(self, x, on_boundary):
                tol = 1E-10
                return on_boundary and (near(x[1],ymin,tol) or near(x[1],ymax,tol) or near(x[0],xmin,tol) or near(x[0],xmax,tol))
    
    mesh = mshr.generate_mesh(outerb-innerb,2/PIVres)
    
    #label the boundary elements
    boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    boundary_markers.set_all(1) #default to 1 (if not on outer, then on inner)
    bx0 = Boundary_outer()
    bx0.mark(boundary_markers, 0) #outer line element ds(0)   
    
    return mesh, boundary_markers

def defineBCs(mesh,boundary_markers):
    '''define FEM space, BCs'''   
    
    #Create fenics UserExpression to access data
    data_fenics = Velocity_from_PIV()
    
    # Define function space, choose Taylor-Hood elements
    P2 = VectorElement("Lagrange", mesh.ufl_cell(), 3)
    P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
    TH = P2 * P1
    W = FunctionSpace(mesh, TH)

    #Define Dirichlet BC on each part of the boundary
    bc0 = DirichletBC(W.sub(0), data_fenics, boundary_markers,0)
    bc1 = DirichletBC(W.sub(0), data_fenics, boundary_markers,1)

    #collect BCs
    bcs = [bc0, bc1]
    
    return W, bcs

def calc_fem(W,bcs):
    """obtain u and p from Stokes FEM"""
    
    # Define variational problem
    (u, p) = TrialFunctions(W)
    (v, q) = TestFunctions(W)
    
    #Stokes equations in weak form
    a = (inner(grad(u), grad(v)) - div(v)*p-div(u)*q -(1E-10)*p*q)*dx
    L = -Constant(0)*q*dx
    
    # Compute solution
    w = Function(W)
    solve(a == L, w, bcs)   

    #extract velocity and pressure fields
    (u, p) = w.split()
    return u, p

u, p =main()

plot(u)

[Emanuel_Gebauer]

Interesting…
no Idea why copying the definition of the bcs inside the main function helped, but after that i did not come accross this issue…

Consider the folllowing variation of your code:

from dolfin import *

'''Classes'''

class Velocity_from_PIV(UserExpression):
    # Fenics object to evaluate/interpolate PIV velocity data anywhere
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
    def eval_cell(self, value, x, ufc_cell):        
        value[0] = 1.
        value[1] = 2.
    def value_shape(self):
        return (2,)

'''Methods''' 

def main():
    #Calculate mesh
    mesh = UnitSquareMesh(10,10)
    
    plot(mesh)

    #Define function space and boundary conditions
    
    data_fenics = Velocity_from_PIV()
    
    W, bcs=defineBCs(mesh, data_fenics)

    #Perform FEM
    u, p = calc_fem(W,bcs) 
    
    return u, p


def defineBCs(mesh, data_fenics):
    '''define FEM space, BCs'''   
    
    #Create fenics UserExpression to access data
    if data_fenics is None:
        data_fenics = Velocity_from_PIV()


    # Define function space, choose Taylor-Hood elements
    P2 = VectorElement("Lagrange", mesh.ufl_cell(), 3)
    P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
    TH = P2 * P1
    W = FunctionSpace(mesh, TH)

    #Define Dirichlet BC on each part of the boundary
    bc0 = DirichletBC(W.sub(0), data_fenics, 'on_boundary')

    bcs = [bc0]
    
    return W, bcs

def calc_fem(W,bcs):
    """obtain u and p from Stokes FEM"""
    
    # Define variational problem
    (u, p) = TrialFunctions(W)
    (v, q) = TestFunctions(W)
    
    #Stokes equations in weak form
    a = (inner(grad(u), grad(v)) - div(v)*p-div(u)*q -(1E-10)*p*q)*dx
    L = -Constant(0)*q*dx
    
    # Compute solution
    w = Function(W)
    solve(a == L, w, bcs)   

    #extract velocity and pressure fields
    (u, p) = w.split()
    return u, p

for i in range(1,100):
    print(str(i))
    u, p =main()

if you change the line W, bcs=defineBCs(mesh, data_fenics) to W, bcs=defineBCs(mesh, None) it crashes similar to your problem.

Interesting… maybe someone else can spot the underlying issue

[conpierce8]

I’d recommend removing the subfunctions and implementing this in the global namespace of the script, i.e.:

from dolfin import *
import mshr
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as interp
#import pandas as pd

'''Classes'''

class Velocity_from_PIV(UserExpression):
    # Fenics object to evaluate/interpolate PIV velocity data anywhere
    def __init__(self, **kwargs):
        super().__init__(kwargs)
        self.x=np.linspace(0,1,20) #vector of x coords (ascending)
        self.y=np.linspace(0,1,20) #vector of y coords (ascending)
        self.dx=np.random.rand(20,20) #matrix of values
        self.dy=np.random.rand(20,20)
        #Splines (default: cubic, no smoothing)
        self.splx = interp.RectBivariateSpline(self.x,self.y,self.dx)
        self.sply = interp.RectBivariateSpline(self.x,self.y,self.dy)
    def eval_cell(self, value, x, ufc_cell):        
        value[0] = self.splx.ev(x[0],x[1])
        value[1] = self.sply.ev(x[0],x[1])
    def value_shape(self):
        return (2,)

'''Methods''' 

"""Create the mesh and determine inner and outer boundaries"""    
# For sake of MWE define rectangular dimensions manually:
xmin, xmax, ymin, ymax, PIVres = np.array([0.0, 1.0, 0.0, 1.0, 0.1])
# x,y coords and radius of circular cutout:
centrex, centrey, radius = np.array([0.27,0.163,0.11])

#define interior and exterior boundary
outerb = mshr.Rectangle(Point(xmin,ymin),Point(xmax,ymax))
innerb = mshr.Circle(Point(centrex,centrey),radius)

#Define classes to let fenics know whether you are on a certain section of the boundary
class Boundary_outer(SubDomain):
        def inside(self, x, on_boundary):
            tol = 1E-10
            return on_boundary and (near(x[1],ymin,tol) or near(x[1],ymax,tol) or near(x[0],xmin,tol) or near(x[0],xmax,tol))

mesh = mshr.generate_mesh(outerb-innerb,2/PIVres)

#label the boundary elements
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary_markers.set_all(1) #default to 1 (if not on outer, then on inner)
bx0 = Boundary_outer()
bx0.mark(boundary_markers, 0) #outer line element ds(0)   

#Create fenics UserExpression to access data
data_fenics = Velocity_from_PIV()

# Define function space, choose Taylor-Hood elements
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 3)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

#Define Dirichlet BC on each part of the boundary
bc0 = DirichletBC(W.sub(0), data_fenics, boundary_markers,0)
bc1 = DirichletBC(W.sub(0), data_fenics, boundary_markers,1)

#collect BCs
bcs = [bc0, bc1]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

#Stokes equations in weak form
a = (inner(grad(u), grad(v)) - div(v)*p-div(u)*q -(1E-10)*p*q)*dx
L = -Constant(0)*q*dx

# Compute solution
w = Function(W)
solve(a == L, w, bcs)   

#extract velocity and pressure fields
(u, p) = w.split()

plot(mesh)
plot(u)

This ran 50 times in a row for me without encountering an error, although I’m not sure why this fixes the problem. My only guess is that it has something to do with namespaces, garbage collection, pybind, old variables remaining from previous executions of the script, or a combination of some or all of these.

[acham]

Hi and sorry for the very late reply. Thanks for giving your thought to this. It is very strange that this solves the problem but I will not complain Maybe something the devs should keep at the back of their mind though.