Interpolated expression gives wrong result for nth root of a function
Note: this issue was reported in the community forums here:
It received no answers, so unfortunately I’m not 100% sure it’s a bug. I’m reporting it in case it is: I apologise in advance if turn out to be wrong.
I get an incorrect result when computing the nrth root of a function. The error is at the mesh edge.
I can reproduce the issue by using this function, defined on a 1D (radial) mesh:
w( r ) = { exp(-r/tau) if 0 < r < 10^9; 0 elsewhere
where tau=10^8. To compute the nth root of this function:
y( r ) = w^(1/n)
I define a dolfin Expression and then interpolate it onto a continuous function space ('CG'). 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.
If I compute the nodal values by casting the function vector as a numpy array, and then using numpy’s ** operation for 1/n, I get the correct result.
Minimum working example:
import numpy as np
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()
