Bug? Interpolated expression gives wrong result for nth root of a function

1

EDIT: should I perhaps flag this as a bug? If I compute the nodal values by acting on the function vector and using numpy ** operation for the power, I get the correct result. It would be good to know another’s opinion before I file a bug report, in case I’m doing something wrong.

Hi all,

I’m trying to compute the nrth root of a function, and get the wrong result at the mesh edge.

Even though the actual function I use is more complicated, I can reproduce the issue by using this function:

w( r ) = { exp(-r/tau) if 0 < r < 10^9; 0 elsewhere

with tau=10^8, of which I compute the nth root:

y( r ) = w^(1/n)

To compute the nth root, I define a dolfin Expression and then interpolate it onto a continuous function space. Below is a minimum working example.

If I run this example, I find that even though w(10^9) = 0 (correctly), I get y(10^9) = w(10^9)^(1/5) = 6.60502302e-05, which is wrong. Why is that?

My understanding of Expressions is that they are evaluated at all degrees of freedom using the string code provided, with the degrees of freedom subsequently used as coefficients for the basis functions upon interpolation. If this understanding is correct, than I don’t understand why I get a wrong result.

Can you please help? Thanks in advance!

P.s.: I’m using FEniCS 2017.2.0.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import dolfin as d

def get_y( w, n, degree):
    if n==2:
         code = "sqrt(fabs(w))"
    elif n==3:
        code = "cbrt(w)"
    else:
        if n%2==0: 
            code = "pow(fabs(w),1./n)"
        else: 
            code = "pow(w,1./n)"
    y_expression = d.Expression( code, w=w, n=n, degree=degree )
    y = d.interpolate( y_expression, V )

    return y


# build non-linear mesh (power law of uniform mesh, cutting the points close to 0)
A = 3.163e-50
gamma = 6.5
c = 3.
N = 4083
num_cells = 4000
r_min = 1e-2
r_max = 1e9

mesh = d.IntervalMesh( num_cells, r_min, r_max )
unif_mesh = np.linspace( 0., r_max, N+1, endpoint=True )
refined_mesh = A * (unif_mesh + c)**gamma
mesh.coordinates()[:,0] = refined_mesh[N-num_cells : ]

# define function space
degree=4
Pn = d.FiniteElement( 'CG', d.interval, degree )
V = d.FunctionSpace( mesh, Pn )

# get w and y=w^(1/n)
w = d.Expression('x[0] < r_max ? exp(-x[0]/1e8) : 0.', degree=degree, r_max=r_max )
w = d.interpolate( w, V )
y = get_y( w, 5, degree )
y.compute_vertex_values()