Hello, I am doing a stability calculation for different profiles. The test case that I have used is Gaussian, where the velocity and gradient of PV is defined using something like this, Ub = Function(V).interpolate(Uj*exp(-(x[0]-Ly/2)**2)) dQb = Function(V).interpolate(-2.*Uj*exp(-(x[0]-Ly/2)**2)*(2.*(x[0]-Ly/2)**2 - 1.)) I have compared this with the results of using a spectral method and the results agree very well. I am now considering a less trivial example using something called a Bump function. Yes, technical name is Bump. Thanks to previous help I have figured out how to define it piecewise (thanks again). Ub = Function(V).interpolate(conditional(le(abs(x[0]-Ly/2.), 1.0), exp(-1./(1. -(x[0]-Ly/2)**2)), 0.0)) dQb = Function(V).interpolate(conditional(le(abs(x[0]-Ly/2.), 1.0), -2.*exp(-1./(1. -(x[0]-Ly/2)**2))*(3.*(x[0]-Ly/2.)**4 - 1)/((x[0]-Ly/2.)**2-1)**4, 0.0)) However, when I try this I get zero growth rate, which does not agree with the spectral method. Actually, in the spectral method we need to zero out the Nans and Infs in expressions, and that works out nicely. Is there an easy way to convert Nans and Infs to zero in interpolate? I tried using both a Krylov method and Lapack. The first gives zero, which I don't believe, and the second fails, I presume because I am dividing by zero somewhere. Any help would be greatly appreciated! Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637