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:
except IOError:
print "IOError with domain_x.dat"
try:
#We avoid ghost cells
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.
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):
"""
"""
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))

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

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))

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

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)
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)))

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)