BUG: Fix streamplot when velocity component is exactly zero.

[tonysyu]

When one component of the velocity field is zero, streamplot fails to draw streamlines. Here's a short example of the issue:

import matplotlib.pyplot as plt
import numpy as np

N = 10
Y, X = np.mgrid[:N, :N]
U = np.ones((N, N))
V = np.zeros((N, N))

plt.streamplot(X, Y, U, V)
plt.show()

Explanation of fix:
Negative and positive steps are treated differently by the code. Before, a zero gradient (cx, cy) was getting treated as a negative step. This resulted in a step size of -inf; treating zero as a positive step results in a step size of (positive) inf. This inf value gets filtered out since the integration routine chooses the smaller of the two (dsx, dsy) step sizes.

[travisbot]

This pull request fails (merged 5639974 into 9448ff2).

[efiring]

@tonysyu, I'm sorry, but I really don't like either the original form or the fix. One should not divide by zero. Relying on it to return inf is not good. At the very least, it raises a RuntimeWarning. Depending on the types of the arguments, it could fail completely. The zero-denominator case should be trapped and handled appropriatedly.

[tonysyu]

@efiring Yeah, it felt a little flakey to me too, but I was just trying to make a minimal change.

Do you think checking for zero and setting the step size to inf would be appropriate? This fix would still rely on min removing the inf later on; if both steps are inf then it would fail. This shouldn't happen, though, since _euler_step is only called when the local velocity causes a trajectory to go beyond the domain boundaries (which wouldn't happen if both velocity components are zero).

[efiring]

@tonysyu, I think your suggestion above is OK. I don't see anything better, offhand.

[tonysyu]

I've added the checks to prevent divide-by-zero warnings. I could explicitly check if both values are zero, but like I said, this shouldn't happen.

I also made a couple of changes to suppress other warnings. I usually don't like changes that venture beyond the PR's topic, but it seemed minor enough. I can revert these, if necessary.

[travisbot]

This pull request fails (merged 1d4682c into 9448ff2).

[dmcdougall]

@tonysyu Perhaps I'm missing something but I have a potentially silly question. If the velocity field is close to zero, passive tracers (basically streamlines) don't move in the fluid. Is there any reason you're setting the step-size to inf as opposed to just not doing an Euler step?

[tonysyu]

Not a silly question at all: the idea here is that only one of the two (x, y) values should be set to inf. E.g. in a horizontal flow, it'd take an infinite step to reach the upper or lower boundary, but the distance to the left or right boundary would be finite.

This routine shouldn't get called if both velocity components are zero, which is why I'm not inclined to add an explicit check for that case. As far as being "close to zero", I'm not sure; were you thinking of some sort of threshold?

[dmcdougall]

Well, I was more thinking that, in your example of a horizontal flow, you would simply ignore the vertical component.

[tonysyu]

That's essentially what happens. The inf stuff in this PR is used to calculate the step size along the stream line. Then, this step size gets scaled to x and y components (for horizontal flow, the y-scale is exactly zero). Is this equivalent to what you had in mind?

[dmcdougall]

By 'stepsize' I meant the Euler step in the Runge Kutta scheme for the odes:

\dot{x} = u(x)
\dot{y} = v(x)

I think we're talking about different things.

[tonysyu]

There's definitely a lot of room for misunderstanding since there are 4 step sizes: dx, dy, ds, and dt. By your last comment, I assume you mean dt. Basically, the scheme calculates the ds (or similarly dt, since this takes a single step based on the velocities at x, y) required to reach the closest horizontal and vertical boundaries and chooses the shorter of the two.

When you say "we're talking about different things", do you mean that your original suggestion no longer applies? Or am I missing something?

[dmcdougall]

My suggestion was to enforce that ds (or dy) is zero when the velocity at x (or y) is zero. Then you don't need to deal with inf.