Good question. 

I agree it's probably better to avoid evaluating a function that divides by zero. 

In the expression I have a conditional expression that uses less than or equal to something.  It occurs to me that if I use only less than (but not equal to) that should do the trick.

I would guess that I should change le to lt but not sure.  Do you think this is what I should try?

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

From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Andrew McRae [A.T.T.McRae@bath.ac.uk]
Sent: Sunday, March 05, 2017 11:15 AM
To: firedrake@imperial.ac.uk
Subject: Re: [firedrake] Nans and infs

Pardon the obvious question, but is there a reason you don't simply rewrite the expression to avoid generating NaN/Inf in the first place?

Replacing Inf with 0 also seems like a strange thing to do!

Andrew

On 3 March 2017 at 21:15, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
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