In [1]:
%pylab inline
%cd /home/pablo/fargo3d-1.3/outputs/fargo/
Populating the interactive namespace from numpy and matplotlib
/home/pablo/fargo3d-1.3/outputs/fargo

This notebook contains some advanced material that can help you when working with the data produced by FARGO3D. You can find the code used in this notebook in "utils/advanced.py".

In this tutorial, we assume a 2D cylindrical mesh.

Let's define a class to store the cylindrical mesh

In [2]:
class Mesh():
    """
    Mesh class, for keeping all the mesh data.
    Input: directory [string] -> this is where the domain files are.
    """
    def __init__(self, directory=""):
        if len(directory) > 1:
            if directory[-1] != '/':
                directory += '/'
        try:
            domain_x = np.loadtxt(directory+"domain_x.dat")
        except IOError:
            print "IOError with domain_x.dat"
        try:
            #We avoid ghost cells
            domain_y = np.loadtxt(directory+"domain_y.dat")[3:-3]
        except IOError:
            print "IOError with domain_y.dat"
        self.xm = domain_x #X-Edge
        self.ym = domain_y #Y-Edge
        
        self.xmed = 0.5*(domain_x[:-1] + domain_x[1:]) #X-Center
        self.ymed = 0.5*(domain_y[:-1] + domain_y[1:]) #Y-Center
        
        #(Surfaces computed from edge to edge)
        #We first build 2D arrays for x & y (theta,r)
        T,R = meshgrid(self.xm, self.ym)
        R2  = R*R
        self.surf = 0.5*(T[:-1,1:]-T[:-1,:-1])*(R2[1:,:-1]-R2[:-1,:-1])

And we can now use this class to get the domain of the simulation

In [3]:
mesh = Mesh(directory="")
print mesh.xm
[-3.14159265 -3.12523019 -3.10886773 -3.09250527 -3.07614281 -3.05978034
 -3.04341788 -3.02705542 -3.01069296 -2.9943305  -2.97796804 -2.96160557
 -2.94524311 -2.92888065 -2.91251819 -2.89615573 -2.87979327 -2.8634308
 -2.84706834 -2.83070588 -2.81434342 -2.79798096 -2.7816185  -2.76525603
 -2.74889357 -2.73253111 -2.71616865 -2.69980619 -2.68344372 -2.66708126
 -2.6507188  -2.63435634 -2.61799388 -2.60163142 -2.58526895 -2.56890649
 -2.55254403 -2.53618157 -2.51981911 -2.50345665 -2.48709418 -2.47073172
 -2.45436926 -2.4380068  -2.42164434 -2.40528188 -2.38891941 -2.37255695
 -2.35619449 -2.33983203 -2.32346957 -2.3071071  -2.29074464 -2.27438218
 -2.25801972 -2.24165726 -2.2252948  -2.20893233 -2.19256987 -2.17620741
 -2.15984495 -2.14348249 -2.12712003 -2.11075756 -2.0943951  -2.07803264
 -2.06167018 -2.04530772 -2.02894526 -2.01258279 -1.99622033 -1.97985787
 -1.96349541 -1.94713295 -1.93077049 -1.91440802 -1.89804556 -1.8816831
 -1.86532064 -1.84895818 -1.83259571 -1.81623325 -1.79987079 -1.78350833
 -1.76714587 -1.75078341 -1.73442094 -1.71805848 -1.70169602 -1.68533356
 -1.6689711  -1.65260864 -1.63624617 -1.61988371 -1.60352125 -1.58715879
 -1.57079633 -1.55443387 -1.5380714  -1.52170894 -1.50534648 -1.48898402
 -1.47262156 -1.45625909 -1.43989663 -1.42353417 -1.40717171 -1.39080925
 -1.37444679 -1.35808432 -1.34172186 -1.3253594  -1.30899694 -1.29263448
 -1.27627202 -1.25990955 -1.24354709 -1.22718463 -1.21082217 -1.19445971
 -1.17809725 -1.16173478 -1.14537232 -1.12900986 -1.1126474  -1.09628494
 -1.07992247 -1.06356001 -1.04719755 -1.03083509 -1.01447263 -0.99811017
 -0.9817477  -0.96538524 -0.94902278 -0.93266032 -0.91629786 -0.8999354
 -0.88357293 -0.86721047 -0.85084801 -0.83448555 -0.81812309 -0.80176063
 -0.78539816 -0.7690357  -0.75267324 -0.73631078 -0.71994832 -0.70358585
 -0.68722339 -0.67086093 -0.65449847 -0.63813601 -0.62177355 -0.60541108
 -0.58904862 -0.57268616 -0.5563237  -0.53996124 -0.52359878 -0.50723631
 -0.49087385 -0.47451139 -0.45814893 -0.44178647 -0.42542401 -0.40906154
 -0.39269908 -0.37633662 -0.35997416 -0.3436117  -0.32724923 -0.31088677
 -0.29452431 -0.27816185 -0.26179939 -0.24543693 -0.22907446 -0.212712
 -0.19634954 -0.17998708 -0.16362462 -0.14726216 -0.13089969 -0.11453723
 -0.09817477 -0.08181231 -0.06544985 -0.04908739 -0.03272492 -0.01636246
  0.          0.01636246  0.03272492  0.04908739  0.06544985  0.08181231
  0.09817477  0.11453723  0.13089969  0.14726216  0.16362462  0.17998708
  0.19634954  0.212712    0.22907446  0.24543693  0.26179939  0.27816185
  0.29452431  0.31088677  0.32724923  0.3436117   0.35997416  0.37633662
  0.39269908  0.40906154  0.42542401  0.44178647  0.45814893  0.47451139
  0.49087385  0.50723631  0.52359878  0.53996124  0.5563237   0.57268616
  0.58904862  0.60541108  0.62177355  0.63813601  0.65449847  0.67086093
  0.68722339  0.70358585  0.71994832  0.73631078  0.75267324  0.7690357
  0.78539816  0.80176063  0.81812309  0.83448555  0.85084801  0.86721047
  0.88357293  0.8999354   0.91629786  0.93266032  0.94902278  0.96538524
  0.9817477   0.99811017  1.01447263  1.03083509  1.04719755  1.06356001
  1.07992247  1.09628494  1.1126474   1.12900986  1.14537232  1.16173478
  1.17809725  1.19445971  1.21082217  1.22718463  1.24354709  1.25990955
  1.27627202  1.29263448  1.30899694  1.3253594   1.34172186  1.35808432
  1.37444679  1.39080925  1.40717171  1.42353417  1.43989663  1.45625909
  1.47262156  1.48898402  1.50534648  1.52170894  1.5380714   1.55443387
  1.57079633  1.58715879  1.60352125  1.61988371  1.63624617  1.65260864
  1.6689711   1.68533356  1.70169602  1.71805848  1.73442094  1.75078341
  1.76714587  1.78350833  1.79987079  1.81623325  1.83259571  1.84895818
  1.86532064  1.8816831   1.89804556  1.91440802  1.93077049  1.94713295
  1.96349541  1.97985787  1.99622033  2.01258279  2.02894526  2.04530772
  2.06167018  2.07803264  2.0943951   2.11075756  2.12712003  2.14348249
  2.15984495  2.17620741  2.19256987  2.20893233  2.2252948   2.24165726
  2.25801972  2.27438218  2.29074464  2.3071071   2.32346957  2.33983203
  2.35619449  2.37255695  2.38891941  2.40528188  2.42164434  2.4380068
  2.45436926  2.47073172  2.48709418  2.50345665  2.51981911  2.53618157
  2.55254403  2.56890649  2.58526895  2.60163142  2.61799388  2.63435634
  2.6507188   2.66708126  2.68344372  2.69980619  2.71616865  2.73253111
  2.74889357  2.76525603  2.7816185   2.79798096  2.81434342  2.83070588
  2.84706834  2.8634308   2.87979327  2.89615573  2.91251819  2.92888065
  2.94524311  2.96160557  2.97796804  2.9943305   3.01069296  3.02705542
  3.04341788  3.05978034  3.07614281  3.09250527  3.10886773  3.12523019
  3.14159265]

We now define a class to store the parameters of the simulation

In [4]:
class Parameters():
    """
    Class for reading the simulation parameters.
    input: string -> name of the parfile, normally variables.par
    """
    def __init__(self, directory=''):
        if len(directory) > 1:
            if directory[-1] != '/':
                directory += '/'
        try:
            params = open(directory+"variables.par",'r') #Opening the parfile
        except IOError:                  # Error checker.
            print  paramfile + " not found."
            return
        lines = params.readlines()     # Reading the parfile
        params.close()                 # Closing the parfile
        par = {}                       # Allocating a dictionary
        for line in lines:             #Iterating over the parfile
            name, value = line.split() #Spliting the name and the value (first blank)
            try:
                float(value)           # First trying with float
            except ValueError:         # If it is not float
                try:
                    int(value)         #                   we try with integer
                except ValueError:     # If it is not integer, we know it is string
                    value = '"' + value + '"'
            par[name] = value          # Filling the dictory
        self._params = par             # A control atribute, actually not used, good for debbuging
        for name in par:               # Iterating over the dictionary
            exec("self."+name.lower()+"="+par[name]) #Making the atributes at runtime

And a new class to store Fields

In [5]:
class Field(Mesh, Parameters):
    """
    Field class, it stores the mesh, parameters and scalar data 
    for a scalar field.
    Input: field [string] -> filename of the field
           staggered='c' [string] -> staggered direction of the field. 
                                      Possible values: 'x', 'y', 'xy', 'yx'
           directory='' [string] -> where filename is
           dtype='float64' (numpy dtype) -> 'float64', 'float32', 
                                             depends if FARGO_OPT+=-DFLOAT is activated
    """
    def __init__(self, field, staggered='c', directory='', dtype='float64'):
        if len(directory) > 1:
            if directory[-1] != '/':
                directory += '/'
        Mesh.__init__(self, directory)       #All the Mesh attributes inside Field
        Parameters.__init__(self, directory) #All the Parameters attributes inside Field

        #Now, the staggering:
        if staggered.count('x')>0:
            self.x = self.xm[:-1] #Do not dump the last element
        else:
            self.x = self.xmed
        if staggered.count('y')>0:
            self.y = self.ym[:-1]
        else:
            self.y = self.ymed
    
        self.data = self.__open_field(directory+field,dtype) #The scalar data is here.

    def __open_field(self, f, dtype):
        """
        Reading the data
        """
        field = fromfile(f, dtype=dtype)
        return field.reshape(self.ny, self.nx)

    def plot(self, log=False, cartesian=False, cmap='magma', **karg):
        """
        A layer for plt.imshow or pcolormesh function.
        if cartesian = True, pcolormesh is launched.
        """
        ax = gca()
        if log:
            data = np.log(self.data)
        else:
            data = self.data
        if cartesian:
            T,R = meshgrid(self.x,self.y)
            X = R*cos(T)
            Y = R*sin(T)
            pcolormesh(X,Y,data,cmap=cmap,**karg)
        else:
            ax.imshow(data, cmap = cmap, origin='lower',aspect='auto',
                      extent=[self.x[0],self.x[-1],self.y[0],self.y[-1]],
                      **karg)
            
    #The same but for contours
    def contour(self, log=False, cartesian=False, **karg):
        if log:
            data = np.log(self.data)
        else:
            data = self.data
        ax = gca()
        T,R = meshgrid(self.x,self.y)
        if cartesian:
            X = R*cos(T)
            Y = R*sin(T)
            contour(X,Y,data,**karg)
        else:
            contour(T,R,data,**karg)

And this is how we use these classes:

In [6]:
rho = Field("gasdens5.dat")

fig = figure(figsize=(14,7))
ax1 = fig.add_subplot(121)

rho.contour(cartesian=True,log=True,
            levels=linspace(log(rho.data).min(),log(rho.data).max(),20),
            colors='w', 
            alpha=0.3, 
            linestyles='solid')

rho.plot(cartesian=True, log=True, vmax=-5.0)
ax = gca()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.set_title("Cartesian coordinates")
t = linspace(-pi,pi,1000)
plot(cos(t),sin(t),'w')
xlim(0.7,1.3)
ylim(-.35,0.35)

#And in polar coordinates

ax2 = fig.add_subplot(122)

rho.plot(log=True,interpolation='bicubic')

rho.contour(log=True, 
            levels=linspace(log(rho.data).min(),log(rho.data).max(),30),
            colors='w', 
            alpha=0.3, 
            linestyles='solid')

ax = gca()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.set_title("Polar coordinates")
tmax = arctan2(0.35,1.3)
ylim(0.7,1.3)
xlim(-tmax,tmax)
Out[6]:
(-0.26299473168091947, 0.26299473168091947)

Let's check some staggered fields

In [7]:
vx = Field("gasvx5.dat",staggered='x')
vy = Field("gasvy5.dat",staggered='y')

Might be useful to have some method to center staggered fields along some direction:

In [8]:
def shift_field(Field,direction):
    import copy
    """
    Half-cell shifting along the direction ('x','y', 'xy', 'yx').
    
    Note that after a call to this function, Field.xm/xmed does not represent 
    the coodinates of the field anymore (it is not hard to improve).
    """
    F = copy.deepcopy(Field)
    if direction.count('x')>0:
        F.data = 0.5*(Field.data[:,1:]+Field.data[:,:-1])
        F.x = 0.5*(Field.x[1:]+Field.x[:-1])
    if direction.count('y')>0:
        F.data = 0.5*(F.data[1:,:]+F.data[:-1,:])
        F.y = 0.5*(F.y[1:]+F.y[:-1])

    F.nx = len(F.x)
    F.ny = len(F.y)
    
    return F
In [9]:
vxc = shift_field(vx,'x')
vyc = shift_field(vy,'y')

We can check if the field was centered correcly along X by comparing its x coordinates, vx.x, with rho.x:

In [10]:
print any(vxc.x-rho.x[:-1])
print any(vyc.y-rho.y[:-1])
False
False

Unfortunately, after the centering, the dimensions of vx, vy and rho are not the same. We need a function to trim the arrays:

In [11]:
def cut_field(Field, direction, side):
    """                                                                                                  
    Cutting a field:
    Input: field --> a Field class                                                                       
           axis  --> 'x', 'y' or 'xy'                                                                    
           side  --> 'p' (plus), 'm' (minnus), 'pm' (plus/minnus)                                        
    """                                                                                                  
    import copy
    
    _cut_field = copy.deepcopy(Field)                                                                     
    ny,nx = Field.ny, Field.nx
    mx = my = px = py = 0
    
    if direction.count('x')>0:
        if side.count('m')>0:
            mx = 1                                                                                       
        if side.count('p')>0:
            px = 1
    if direction.count('y')>0:
        if side.count('m')>0:
            my = 1
        if side.count('p')>0:
            py = 1
                                                                                                         
    _cut_field.data = Field.data[my:ny-py,mx:nx-px]                                                       
    _cut_field.x = _cut_field.x[mx:nx-px]                                                                  
    _cut_field.y = _cut_field.y[my:ny-py]                                                                  
                                                                                                         
    return _cut_field

Let's check how it works...

In [12]:
vxcc = cut_field(vxc,direction='y',side='p')
vycc = cut_field(vyc,direction='x',side='p')
rhoc = cut_field(rho,direction='xy', side='p')
print vxcc.data.shape
print vycc.data.shape
print rhoc.data.shape
(127, 383)
(127, 383)
(127, 383)

Let's now compute the total velocity

In [13]:
import copy
v = copy.deepcopy(vxcc)
v.data = sqrt(vxcc.data**2+vycc.data**2)
In [14]:
figure(figsize=(9,7))
v.plot(log=True, cartesian=True, vmin=-4.5); colorbar()
ax = gca()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
title("Velocity", size=25)
xlim(0.8,1.2)
ylim(-.2,0.2)
Out[14]:
(-0.2, 0.2)

If we want a quick plot of the velocity field, we can a function like this:

In [15]:
def vector_field(vx,vy, **karg):
    nsx = nsy = 3
    T,R = meshgrid(vx.x[::nsx],vx.y[::nsy])
    X = R*cos(T)
    Y = R*sin(T)
    vx = vx.data[::nsy,::nsx]
    vy = vy.data[::nsy,::nsx]
    U = vy*cos(T) - vx*sin(T)
    V = vy*sin(T) + vx*cos(T)
    ax = gca()
    ax.quiver(X,Y,U,V,scale=3,pivot='middle', **karg)
In [16]:
figure(figsize=(8,7))
rho.plot(log=True, cartesian=True); colorbar()
vector_field(vxcc, vycc,color='w')
ax = gca()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
title("Velocity", size=15)
xlim(0.7,1.3)
ylim(-.3,0.3)
Out[16]:
(-0.3, 0.3)

As a real-life example, we show how to compute the vortencity $$\omega/\rho = \frac{\nabla\times\vec{v}|_z}{\rho}$$

In [17]:
def vortencity(rho, vx, vy):
    """
    Computing the vortencity from a staggered velocity field.
    """
    Tx,Rx = np.meshgrid(vx.x,vx.y)
    Ty,Ry = np.meshgrid(vy.x,vy.y)

    rvx = Rx*(vx.data)

    curldata = ((rvx[1:,1:]-rvx[:-1,1:])/(Rx[1:,1:]-Rx[:-1,1:]) -
                (vy.data[1:,1:]-vy.data[1:,:-1])/(Ty[1:,1:]-Ty[1:,:-1]))

    curl = copy.deepcopy(vx)
    curl.nx = curl.nx-1
    curl.ny = curl.ny-1
    curl.x    = vx.x[:-1]
    curl.y    = vy.y[:-1]

    rho_corner = 0.25*(rho.data[1:,1:]+rho.data[:-1,:-1]+
                       rho.data[1:,:-1]+rho.data[:-1,1:])

    T,R = np.meshgrid(curl.x,curl.y)
    curl.data = curldata/R
    
    curl.data = (curl.data + 2.0*rho.omegaframe)/rho_corner
    
    return curl
In [18]:
vort = vortencity(rho, vx, vy)
In [19]:
fig = figure(figsize=(14,7))

ax1 = fig.add_subplot(121)
vort.plot(cartesian=True)
title("Cartesian coordinates", size=15)
xlim(0.6,1.4)
ylim(-.4,0.4)

ax2 = fig.add_subplot(122)
vort.plot()
title("Polar coordinates", size=15)
ylim(0.7,1.3)
Out[19]:
(0.7, 1.3)

Let's now compute some streamlines. We will use a first order integrator, an interpolator and a recursive stream mannager:

In [20]:
def euler(vx, vy, x, y, reverse):
    """
    Euler integrator for computing the streamlines.
    Parameters:
    ----------

    x,y: Floats.
         Initial condition
    reverse: Boolean.
             If reverse is true, the integration step is negative.

    Output
    ------

    (dx,dy): (float,float).
             Are the azimutal and radial increments.
             Only works for cylindrical coordinates.
    """
    sign = 1.0
    if reverse:
        sign = -1
    vphi = get_v(vx,x,y)
    vrad = get_v(vy,x,y)
    if vphi == None or vrad == None: #Avoiding problems...                                                                                                                                                                                    
        return None,None
    l = np.min((((vx.xmax-vx.xmin)/vx.nx),((vx.ymax-vx.ymin)/vx.ny)))
    h = 0.5*l/np.sqrt((vphi**2+vrad**2))

    return sign*h*np.array([vphi/y,vrad])


def bilinear(x,y,f,p):
    """
    Bilinear interpolation.
    Parameters
    ----------
    x = (x1,x2); y = (y1,y2)
    f = (f11,f12,f21,f22)
    p = (x,y)
    where x,y are the interpolated points and
    fij are the values of the function at the
    points (xi,yj).
    Output
    ------
    f(p): Float.
          The interpolated value of the function f(p) = f(x,y)
    """

    xp  = p[0]; yp   = p[1]; x1  = x[0]; x2  = x[1]
    y1  = y[0]; y2  = y[1];  f11 = f[0]; f12 = f[1]
    f21 = f[2]; f22 = f[3]
    t = (xp-x1)/(x2-x1);    u = (yp-y1)/(y2-y1)

    return (1.0-t)*(1.0-u)*f11 + t*(1.0-u)*f12 + t*u*f22 + u*(1-t)*f21

def get_v(v, x, y):
    """
    For a real set of coordinates (x,y), returns the bilinear            
    interpolated value of a Field class.
    """

    i = int((x-v.xmin)/(v.xmax-v.xmin)*v.nx)
    j = int((y-v.ymin)/(v.ymax-v.ymin)*v.ny)

    if i<0 or j<0 or i>v.data.shape[1]-2 or j>v.data.shape[0]-2:
        return None

    f11 = v.data[j,i]
    f12 = v.data[j,i+1]
    f21 = v.data[j+1,i]
    f22 = v.data[j+1,i+1]
    try:
        x1  = v.x[i]
        x2  = v.x[i+1]
        y1  = v.y[j]
        y2  = v.y[j+1]
        return bilinear((x1,x2),(y1,y2),(f11,f12,f21,f22),(x,y))
    except IndexError:
        return None
    
    
    
def get_stream(vx, vy, x0, y0, nmax=1000000, maxlength=4*np.pi, bidirectional=True, reverse=False):
    """
    Function for computing a streamline.
    Parameters:
    -----------

    x0,y0: Floats.
          Initial position for the stream
    nmax: Integer.
          Maxium number of iterations for the stream.
    maxlength: Float
               Maxium allowed length for a stream
    bidirectional=True
                  If it's True, the stream will be forward and backward computed.
    reverse=False
            The sign of the stream. You can change it mannualy for a single stream,
            but in practice, it's recommeneded to use this function without set reverse
            and setting bidirectional = True.

    Output:
    -------

    If bidirectional is False, the function returns a single array, containing the streamline:
    The format is:

                                      np.array([[x],[y]])

    If bidirectional is True, the function returns a tuple of two arrays, each one with the same
    format as bidirectional=False.
    The format in this case is:

                            (np.array([[x],[y]]),np.array([[x],[y]]))

    This format is a little bit more complicated, and the best way to manipulate it is with iterators.
    For example, if you want to plot the streams computed with bidirectional=True, you can do:

    stream = get_stream(x0,y0)
    ax.plot(stream[0][0],stream[0][1]) #Forward
    ax.plot(stream[1][0],stream[1][1]) #Backward

    """
    
    if bidirectional:
        s0 = get_stream(vx, vy, x0, y0, reverse=False, bidirectional=False, nmax=nmax,maxlength=maxlength)
        s1 = get_stream(vx, vy, x0, y0, reverse=True,  bidirectional=False, nmax=nmax,maxlength=maxlength)
        return (s0,s1)

    l = 0
    x = [x0]
    y = [y0]

    for i in xrange(nmax):
        ds = euler(vx, vy, x0, y0, reverse=reverse)
        if ds[0] == None:
            if(len(x)==1):
                print "There was an error getting the stream, ds is NULL (see get_stream)."
            break
        l += np.sqrt(ds[0]**2+ds[1]**2)
        dx = ds[0]
        dy = ds[1]
        if(np.sqrt(dx**2+dy**2)<1e-13):
            print "Warning (get_stream): ds is very small, check if you're in a stagnation point."
            print "Try selecting another initial point."
            break
        if l > maxlength:
            print "maxlength reached: ", l
            break
        x0 += dx
        y0 += dy
        x.append(x0)
        y.append(y0)
        
    return np.array([x,y])

def get_random_streams(vx, vy, xmin=None, xmax=None, ymin=None, ymax=None, n=30, nmax=100000):
    if xmin == None:
        xmin = vx.xmin
    if ymin == None:
        ymin = vx.ymin
    if xmax == None:
        xmax = vx.xmax
    if ymax == None:
        ymax = vx.ymax
    X = xmin + np.random.rand(n)*(xmax-xmin)
    Y = ymin + np.random.rand(n)*(ymax-ymin)
    streams = []
    counter = 0
    for x,y in zip(X,Y):
        stream = get_stream(vx, vy, x, y, nmax=nmax, bidirectional=True)
        streams.append(stream)
        counter += 1
    return streams

def plot_random_streams(streams, cartesian=False, **kargs):
    ax = plt.gca()
    for stream in streams:
        for sub_stream in stream:
            if cartesian:
                ax.plot(sub_stream[1]*cos(sub_stream[0]),sub_stream[1]*sin(sub_stream[0]),**kargs)
            else:
                ax.plot(sub_stream[0],sub_stream[1],**kargs)
In [21]:
random_streams = get_random_streams(vx,vy, n=100)
maxlength reached:  12.57056309127447
maxlength reached:  12.567712208934692
There was an error getting the stream, ds is NULL (see get_stream).
There was an error getting the stream, ds is NULL (see get_stream).
In [22]:
figure(figsize=(9,9))
rho.plot(log=True, cartesian=True, cmap='gray')
vector_field(vx, vy, color='r', alpha=0.5)
plot_random_streams(random_streams,cartesian=True,color='w', linewidth=2, alpha=0.5)
ax = gca()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
title("Velocity", size=15)
xlim(0.5,1.4)
ylim(-0.45,0.45)
Out[22]:
(-0.45, 0.45)