FARGO3D TOUR

In [1]:
import os
In [2]:
home = "/home/" #Where your home is stored.
user = "pablo" #You should define your user here.
In [3]:
path = home+user+"/"

First for all, you have to download the code from the webpage. For simplicity, we assume you have a tar.gz file called "fargo3d-x.x.tar.gz" in the Downloads directory:

In [4]:
os.chdir(path+'Downloads') #Equivalent to cd "/home/pablo/Downloads" but more versatile with the path variable.
In [5]:
!tar zxf fargo3d-0.9.tar.gz

We move the fargo3d-0.9 folder into the main directory:

In [6]:
os.system("mv fargo3d-0.9 "+ path)
Out[6]:
256

We move the cursor inside the main directory of our fresh fargo3d distribution:

In [7]:
os.chdir(path+"fargo3d-0.9") #At this time we are in /home/pablo/fargo3d-0.9
In [8]:
!pwd
!ls
/home/pablo/fargo3d-0.9
Advanced.ipynb	FARGO3D_Tour.ipynb  matplotlibrc.b  scripts  test_suite
bin		in		    outputs	    setups   utils
doc		license.txt	    planets	    src
fargo3d		Makefile	    README	    std

Now, we have the code ready for the making process. The make process is platform dependent, if you have a standard Linux platform, normally the standalone version works without additional configurations, but instead if you have a Mac platform, you should define an environment variable called FARGO_ARCH. (It is recommended to read the documentation in both cases).

Suppose you have a Mac platform:

You should export the variable FARGO_ARCH in this way: type in your shell "export FARGO_ARCH=MacIntel" (bash shell)

In IPython we have not a true shell, so we can not export a variable simply calling !export FARGO_ARCH, instead, we can pass this variable as a make variable:

In [9]:
!make FARGO_ARCH=MacIntel
make[1]: Entering directory `/home/pablo/fargo3d-0.9/bin'
make[1]: `../src/rescale.c' is up to date.
make[1]: Leaving directory `/home/pablo/fargo3d-0.9/bin'
make[1]: Entering directory `/home/pablo/fargo3d-0.9/bin'

          FARGO3D SUMMARY:           
          ===============            

This built is SEQUENTIAL. Use "make para" to change that

This built can be launched on
a CPU with a GPU card (1 GPU only).


SETUP:      'fargo'         
(Use "make SETUP=[valid_setup_string]" to change set up)
(Use "make list" to see the list of setups implemented)
(Use "make info" to see the current sticky build options)

make[1]: Leaving directory `/home/pablo/fargo3d-0.9/bin'

This notebook is running under a Linux platform so it raises errors. The way to do this under linux is:

In [10]:
!make GPU=1 FARGO_ARCH=LINUX
make[1]: Entering directory `/home/pablo/fargo3d-0.9/bin'
make[1]: `../src/rescale.c' is up to date.
make[1]: Leaving directory `/home/pablo/fargo3d-0.9/bin'
make[1]: Entering directory `/home/pablo/fargo3d-0.9/bin'

          FARGO3D SUMMARY:           
          ===============            

This built is SEQUENTIAL. Use "make para" to change that

This built can be launched on
a CPU with a GPU card (1 GPU only).


SETUP:      'fargo'         
(Use "make SETUP=[valid_setup_string]" to change set up)
(Use "make list" to see the list of setups implemented)
(Use "make info" to see the current sticky build options)

make[1]: Leaving directory `/home/pablo/fargo3d-0.9/bin'

In [11]:
!ls #Exploring the fargo3d-9.0 directory
Advanced.ipynb	FARGO3D_Tour.ipynb  matplotlibrc.b  scripts  test_suite
bin		in		    outputs	    setups   utils
doc		license.txt	    planets	    src
fargo3d		Makefile	    README	    std

We have the FARGO3D executable file (called fargo3d) ready to run! (This file was created after the make process). As the build process was without special parallel flags, this build is gpu-sequential (if you do not have a gpu, remove "GPU=1" from the compilation line). Now, we run the code:

In [12]:
!./fargo3d -m in/fargo.par

       !!!! WARNING !!!!

This is a sequential built of the ./fargo3d code
If you planned to run the MPI-parallel version,
then you MUST rebuild the executable. In this case,
 issue:
make PARALLEL=1 (or make para)

Any subsequent invocation of make will build an MPI version.

=========================
PROCESS NUMBER       0
RUNNING ON DEVICE NÂș 0
GEFORCE GT 640
COMPUTE CAPABILITY: 3.0
VIDEO RAM MEMORY: 2 GB
=========================

The default output directory root is ./
Reading parameters file in/fargo.par
The output directory is ./outputs/fargo/
I do not output the ghost values
CompPressFieldIso runs on the GPU
CompPressFieldAd runs on the GPU
Substep1 runs on the GPU
Substep2 runs on the GPU
Substep3 runs on the GPU
DivideByRho runs on the GPU
Vanleer runs on the GPU
momenta runs on the GPU
update runs on the GPU
NewVelocity runs on the GPU
Reduction runs on the GPU
AdvectShift runs on the GPU
ComputeResidual runs on the GPU
ChangeFrame runs on the GPU
Potential runs on the GPU
CorrectVtheta runs on the GPU
cfl runs on the GPU
ComputeForce runs on the GPU
VanLeerPPA runs on the GPU
Copy velocities runs on the GPU
Viscous tensor is computed on the GPU
addviscosity runs on the GPU
boundaries runs on the GPU
Warning: your version of MPI does not allow direct GPU-GPU communications.
Communications are done through the hosts (CPUs), and may not be very efficient.
Monitoring runs on the GPU
CheckMute runs on the GPU
ComputeSlopes runs on the GPU
ComputeStar runs on the GPU
ComputeEmf runs on the GPU
UpdateMagneticField runs on the GPU
LorentzForce runs on the GPU
Resistivity runs on the GPU
FargoMHD runs on the GPU
Stockholm Boundary runs on the GPU
Warning: The Y spacing is linear (default).
./outputs/fargo/dimensions.dat
./outputs/fargo/dims.dat
Field2D Reduction2D has been created
Field Reduction2D has created on the GPU
Pitch = 2560 bytes (320 elements)
Field gasvx has been created
Field gasvx has been created on the GPU
Field Vx_temp has been created
Field Vx_temp has been created on the GPU
Field Moment_Plus_X has been created
Field Moment_Plus_X has been created on the GPU
Field Moment_Minus_X has been created
Field Moment_Minus_X has been created on the GPU
Field2D VxMed has been created
Field VxMed has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D Vxhy has been created
Field Vxhy has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D Vxhyr has been created
Field Vxhyr has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D Vxhz has been created
Field Vxhz has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D Vxhzr has been created
Field Vxhzr has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D Eta_xi has been created
Field Eta_xi has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D Eta_xizi has been created
Field Eta_xizi has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D Eta_zi has been created
Field Eta_zi has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D Nshift has been created
Integer field Nshift has been created on the GPU
Field2D Nxhy has been created
Integer field Nxhy has been created on the GPU
Field2D Nxhz has been created
Integer field Nxhz has been created on the GPU
Field gasvy has been created
Field gasvy has been created on the GPU
Field Vy_temp has been created
Field Vy_temp has been created on the GPU
Field Moment_Plus_Y has been created
Field Moment_Plus_Y has been created on the GPU
Field Moment_Minus_Y has been created
Field Moment_Minus_Y has been created on the GPU
Field potential has been created
Field potential has been created on the GPU
Field Slope has been created
Field Slope has been created on the GPU
Field DivRho has been created
Field DivRho has been created on the GPU
Field DensStar has been created
Field DensStar has been created on the GPU
Field Qs has been created
Field Qs has been created on the GPU
Field Pressure has been created
Field Pressure has been created on the GPU
Field gasdens has been created
Field gasdens has been created on the GPU
Field gasenergy has been created
Field gasenergy has been created on the GPU
Grids Pressure and QLeft share their storage
Field QRight has been created
Field QRight has been created on the GPU
Field2D rho0 has been created
Field rho0 has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D e0 has been created
Field e0 has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D vx0 has been created
Field vx0 has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D vy0 has been created
Field vy0 has created on the GPU
Pitch = 2560 bytes (320 elements)
Field2D vz0 has been created
Field vz0 has created on the GPU
Pitch = 2560 bytes (320 elements)
1 planet found.
Planet number 0
---------------
x = 1.0000000000	y = 0.0000000000	z = 0.0000000000
vx = -0.0000000000	vy = 1.0004998751	vz = 0.0000000000
Non-accreting.
Doesn't feel the disk potential
Doesn't feel the other planets potential

Found 0 communicators
OUTPUTS 0 at date t = 0.000000 OK
TotalMass = 0.0121800000 
!................
.........................
.............................
..............................
..............................
.............................
.............................
.............................
............................
...........................
..........................
........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
OUTPUTS 1 at date t = 6.283185 OK
TotalMass = 0.0122521810 
!........................
.........................
.........................
.........................
.........................
........................
........................
........................
........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
..........................
OUTPUTS 2 at date t = 12.566371 OK
TotalMass = 0.0123551955 
!.........................
..........................
..........................
.........................
.........................
.........................
........................
........................
........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
OUTPUTS 3 at date t = 18.849556 OK
TotalMass = 0.0123872556 
!........................
.........................
.........................
.........................
.........................
........................
........................
........................
........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
OUTPUTS 4 at date t = 25.132741 OK
TotalMass = 0.0123874059 
!........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
.........................
........................
........................
........................
........................
.........................
.........................
.........................
.........................
.........................
........................
OUTPUTS 5 at date t = 31.415927 OK
TotalMass = 0.0123994503 
!.......................
........................
........................
........................
........................
........................
........................
........................
.......................
.......................
.......................
.......................
........................
........................
........................
........................
........................
........................
........................
........................
OUTPUTS 6 at date t = 37.699112 OK
TotalMass = 0.0123772379 
!.......................
........................
........................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
........................
........................
.......................
.......................
OUTPUTS 7 at date t = 43.982297 OK
TotalMass = 0.0123925723 
!......................
.......................
.......................
.......................
.......................
......................
......................
......................
......................
......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
.......................
OUTPUTS 8 at date t = 50.265482 OK
TotalMass = 0.0123915308 
!......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
.......................
.......................
.......................
.......................
......................
......................
......................
......................
OUTPUTS 9 at date t = 56.548668 OK
TotalMass = 0.0123948336 
!.....................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
......................
OUTPUTS 10 at date t = 62.831853 OK
TotalMass = 0.0123974276 
!.....................
End of simulation!

After the execution, we have the data files stored in outputs/fargo directory.

In [13]:
os.chdir(path+"fargo3d-0.9/outputs/fargo") #At this time we are in /home/pablo/fargo3d-0.9
In [14]:
!ls
density0_2d.dat  gasdens9.dat	  gasvx13.dat  gasvy4.dat
dimensions.dat	 gasenergy0.dat   gasvx14.dat  gasvy5.dat
dims.dat	 gasenergy10.dat  gasvx1.dat   gasvy6.dat
domain_x.dat	 gasenergy11.dat  gasvx2.dat   gasvy7.dat
domain_y.dat	 gasenergy12.dat  gasvx3.dat   gasvy8.dat
domain_z.dat	 gasenergy13.dat  gasvx4.dat   gasvy9.dat
gasdens0.dat	 gasenergy14.dat  gasvx5.dat   grid000.inf
gasdens10.dat	 gasenergy1.dat   gasvx6.dat   IDL.var
gasdens11.dat	 gasenergy2.dat   gasvx7.dat   mass.dat
gasdens12.dat	 gasenergy3.dat   gasvx8.dat   momx.dat
gasdens13.dat	 gasenergy4.dat   gasvx9.dat   orbit0.dat
gasdens14.dat	 gasenergy5.dat   gasvy0.dat   planet0.dat
gasdens1.dat	 gasenergy6.dat   gasvy10.dat  torq_1d_Y_raw_planet_0.dat
gasdens2.dat	 gasenergy7.dat   gasvy11.dat  torq_planet_0.dat
gasdens3.dat	 gasenergy8.dat   gasvy12.dat  used_rad.dat
gasdens4.dat	 gasenergy9.dat   gasvy13.dat  variables.par
gasdens5.dat	 gasvx0.dat	  gasvy14.dat  vx0_2d.dat
gasdens6.dat	 gasvx10.dat	  gasvy1.dat   vy0_2d.dat
gasdens7.dat	 gasvx11.dat	  gasvy2.dat
gasdens8.dat	 gasvx12.dat	  gasvy3.dat

Before to read and plot some sacalar field, we will see how to deal with the simulation parameters (eg: nx, SigmaSlope, AspectRatio, OmegaFrame, OutputDir, etc). For comfort we dump all these parameters (standard & user-defined) in two files, called IDL.var or variables.par. The first one has a fancy format for IDL. The second one is a normal .par file that could be used to start (restart) a normal simulation, but includes all the parameters, stored in the std/standard.par & SETUPNAME.par file, stored in the SETUP directory.

We will load all the variables inside Python in a fancy way. We wrote the next Class (it is not hard to understand it, the magic is the function "exec", please read the python documentation for a good explanation):

In [15]:
class Parameters(object):
    """
    Class for reading the simulation parameters.
    input: string -> name of the parfile, normally variables.par
    """
    def __init__(self, paramfile):
        try:
            params = open(paramfile,'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

example on how to use this Class:

In [16]:
P = Parameters("variables.par")

At this time, if you type "P." and then press "Enter", you will see all the variables stored in variables.par as atributes of the Object P:

In [17]:
print P.nx, P.ny, P.omegaframe, P.outputdir, P.sigmaslope, P.sigma0
512 256 1.0005 ./outputs/fargo/ 0 0.00063661977237

When you develop python scripts for dealing with your simulations, this class should be very useful for developing general scripts.

In [18]:
P._params
Out[18]:
{'ASPECT': '"auto"',
 'ASPECTRATIO': '0.05',
 'AUTOCOLOR': '1',
 'CMAP': '"cubehelix"',
 'COLORBAR': '1',
 'COORDINATES': '"cylindrical"',
 'CS': '1',
 'DISK': '0',
 'DT': '0.314159265359',
 'ECCENTRICITY': '0',
 'ETA': '0',
 'EXCLUDEHILL': '0',
 'FIELD': '"gasdens"',
 'FLARINGINDEX': '0',
 'FRAME': '"G"',
 'FUNCARCHFILE': '"std/func_arch.cfg"',
 'GAMMA': '1.66666667',
 'INCLINATION': '0',
 'INDIRECTTERM': '0',
 'LOG': '0',
 'MASSTAPER': '0',
 'NINTERM': '20',
 'NOISE': '0',
 'NTOT': '200',
 'NU': '1e-05',
 'NX': '512',
 'NY': '256',
 'NZ': '1',
 'OMEGAFRAME': '1.0005',
 'OORTA': '-0.75',
 'OUTPUTDIR': '"./outputs/fargo/"',
 'PERIODICY': '0',
 'PERIODICZ': '0',
 'PLANETCONFIG': '"planets/jupiter.cfg"',
 'PLOTLINE': '"field[:,:,0]"',
 'REALTYPE': '"float64"',
 'ROCHESMOOTHING': '0',
 'SETUP': '"fargo"',
 'SIGMA0': '0.00063661977237',
 'SIGMASLOPE': '0',
 'SPACING': '"default"',
 'THICKNESSSMOOTHING': '0.6',
 'VERTICALDAMPING': '0',
 'VMAX': '1',
 'VMIN': '0',
 'VTK': '0',
 'XMAX': '3.14159265358979',
 'XMIN': '-3.14159265358979',
 'YMAX': '2.5',
 'YMIN': '0.4',
 'ZMAX': '1',
 'ZMIN': '0'}

It is time to make some plots of the scalar fields. First, we need to load the data:

In [19]:
n = 10
rho = fromfile("gasdens{0:d}.dat".format(n),dtype='float64') #If the simulation was performed in simple precision, dtype should be 'float32'

If we try to plot the density, the result is:

In [20]:
figure(figsize=(15,5))
plot(rho)
Out[20]:
[<matplotlib.lines.Line2D at 0x270c050>]

Really it is not a good plot! The fields (density, velocity, magnetic field, etc) are stored in a 1D array of doubles/floats, and we should perform a reshape before to plot the colormap. In this case, the simulation is 2D, so, the reshape is from a 1D array --> 2D array (in 3D is 1D-->3D).

In [21]:
rho = rho.reshape(P.ny,P.nx)
In [22]:
figure(figsize=(10,10))
imshow(log10(rho), origin='lower', cmap='cubehelix', aspect='auto', extent=[P.xmin,P.xmax,P.ymin,P.ymax])
Out[22]:
<matplotlib.image.AxesImage at 0x2edb4d0>

Much better! Supose you want to study the evolution of the planet gap in the disk:

In [23]:
figure(figsize=(15,5))
for i in range(10):
    rho = fromfile("gasdens{0:d}.dat".format(i)).reshape(P.ny,P.nx).mean(axis=1)
    plot(linspace(P.ymin,P.ymax,P.ny),rho)
xlim(P.ymin,P.ymax)
Out[23]:
(0.4, 2.5)

Please, be careful with the "extent" argument in imshow(); it should be used for making quick plots, and not for precise studies. For example, P.xmin & P.xmax are refered to the edges of the density field, but for the center of the x-velocity, and the same for P.ymin & P.ymax for the y-velocity. In the documentation you can read about the staggering of the fields. Maybe the best way to deal with the staggered fields is computing all the coordinates for each field, using the domain_i.dat files:

In [24]:
xmin = loadtxt("domain_x.dat") #x-edges
ymin = loadtxt("domain_y.dat")[3:-3] #y-edges, dropping the ghost cells included. For z direction is the same.
xc   = 0.5*(xmin[1:]+xmin[:-1]) #x-center
yc   = 0.5*(ymin[1:]+ymin[:-1]) #y-center

The number of elements of xc, is the number of cells in the outputs. For y & z directions, is the same except for some additional cells (ghost cells) that must be dropped:

In [25]:
print xc.shape, yc.shape
(512,) (256,)

In [26]:
print xmin[0], xmin[-1], ymin[0], ymin[-1]
print P.xmin, P.xmax, P.ymin, P.ymax
-3.14159265359 3.14159265359 0.4 2.5
-3.14159265359 3.14159265359 0.4 2.5

You can see that P.xmin, P.xmax, P.ymin and P.ymax are actually the edges of the cells.

As a final example, we will compute the velocity on the mesh: \(\displaystyle{v = \sqrt{v_x^2+v_y^2}}\)

Fargo3d dump both the velocity-x and the velocity-y with the same number of cells that the density, but actually the x-velocity should be an array of (nx+1,ny) and the y-velocity should be (nx,ny+1). For simplicity, we cut the final element (x in x & y in y) of these field. Maybe it is not the best choise, but it makes the read process more compact. The price is the lost of elements when we want to recenter the fields:

In [27]:
vx = fromfile("gasvx{0:d}.dat".format(n)).reshape(P.ny,P.nx) #X-satggered Y-centered field
vy = fromfile("gasvy{0:d}.dat".format(n)).reshape(P.ny,P.nx) #X-centered Y-staggered field

We can calculate v as (not precise):

In [28]:
figure(figsize=(10,10))
v = sqrt(vx**2+vy**2)
imshow(log10(v), origin='lower', cmap='cubehelix', aspect='auto', extent=[P.xmin,P.xmax,P.ymin,P.ymax])
xlim(-1,1)
ylim(0.6,1.4)
Out[28]:
(0.6, 1.4)

Note in this case extent is an approximation! The right way to do that is:

In [29]:
figure(figsize=(10,10))
vxc = 0.5*(vx[:-1,1:]+vx[:-1,:-1])
vyc = 0.5*(vy[1:,:-1]+vx[:-1,:-1])
v = sqrt(vxc**2+vyc**2)
imshow(log10(v), origin='lower', cmap='cubehelix', aspect='auto', extent=[xc.min(),xc.max(),yc.min(),yc.max()])
xlim(-1,1)
ylim(0.6,1.4)
Out[29]:
(0.6, 1.4)

If you want to see more examples, go to the notebook about advanced scripting.