In [1]:
import os
f3d_directory = "/home/pablo/fargo3d-1.3/" #Change it to match your fargo3d directory
os.chdir(f3d_directory)                    #It changes the current working directory

We will now build the executable. This process is platform dependent. See the documentation for details.

  • In a Linux platform, the standalone version usually works without additional configurations.
  • In a Mac platform, the environment variable FARGO_ARCH must be exported.

e.g. If you have a Mac platform, write in your terminal export FARGO_ARCH=MacIntel (valid for a bash shell)

Another way to get the same result is declaring the variable just before calling make:

In [2]:
#!make FARGO_ARCH=MacIntel #Execute this line only if you are in a Mac platform
In [3]:
!make FARGO_ARCH=LINUX     #Execute this line only if you are in a Linux platform
Skipping --
make[1]: Entering directory '/home/pablo/fargo3d-1.3/bin'
make[1]: '../src/rescale.c' is up to date.
make[1]: Leaving directory '/home/pablo/fargo3d-1.3/bin'
make[1]: Entering directory '/home/pablo/fargo3d-1.3/bin'

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

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


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-1.3/bin'

We can check the content of the directory using !ls

In [4]:
!ls
arch	 fargo3d-1.2	     license.txt  planets  setups  test_suite
bin	 fargo3d-1.2.tar.gz  Makefile	  README   src	   utils
fargo3d  in		     outputs	  scripts  std

And can run the code using the following instruction:

In [5]:
!./fargo3d -o "ntot=100" 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.


===========================
FARGO3D Public version 1.3
SETUP: 'fargo'
===========================

The default output directory root is ./
Reading parameters file in/fargo.par
Warning : NTOT is redefined on the command line.
The output directory is ./outputs/fargo/
I do not output the ghost values
Warning: The Y spacing is linear (default).
Field2D Reduction2D has been created
Field gasvx has been created
Field Vx_temp has been created
Field Moment_Plus_X has been created
Field Moment_Minus_X has been created
Field2D VxMed has been created
Field2D Vxhy has been created
Field2D Vxhyr has been created
Field2D Vxhz has been created
Field2D Vxhzr has been created
Field2D Eta_xi has been created
Field2D Eta_xizi has been created
Field2D Eta_zi has been created
Field2D Nshift has been created
Field2D Nxhy has been created
Field2D Nxhz has been created
Field gasvy has been created
Field Vy_temp has been created
Field Moment_Plus_Y has been created
Field Moment_Minus_Y has been created
Field potential has been created
Field Slope has been created
Field DivRho has been created
Field DensStar has been created
Field Qs has been created
Field Pressure has been created
Field gasdens has been created
Field gasenergy has been created
Grids Pressure and QLeft share their storage
Field QRight has been created
Field2D rho0 has been created
Field2D e0 has been created
Field2D vx0 has been created
Field2D vy0 has been created
Field2D vz0 has been created
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
Standard version with no ghost zones in X
OUTPUTS 0 at date t = 0.000000 OK
............
..............
...............
.................
.................
...................
.................
.................
................
................
...............
...............
...............
..............
..............
..............
..............
..............
..............
..............
OUTPUTS 1 at date t = 6.283185 OK
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
OUTPUTS 2 at date t = 12.566371 OK
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
OUTPUTS 3 at date t = 18.849556 OK
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
OUTPUTS 4 at date t = 25.132741 OK
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
..............
OUTPUTS 5 at date t = 31.415927 OK
..............
End of simulation!

After the execution, the output files are stored in the directory outputs/fargo

In [6]:
os.chdir(f3d_directory+"/outputs/fargo") #We change the current working directory and move into outputs directory

Let's see the content of this directory

In [7]:
!ls
bigplanet0.dat	 gasenergy11.dat  gasvx7.dat   summary0.dat
density0_2d.dat  gasenergy12.dat  gasvx8.dat   summary10.dat
dimensions.dat	 gasenergy13.dat  gasvx9.dat   summary11.dat
dims.dat	 gasenergy1.dat   gasvy0.dat   summary12.dat
domain_x.dat	 gasenergy2.dat   gasvy10.dat  summary13.dat
domain_y.dat	 gasenergy3.dat   gasvy11.dat  summary1.dat
domain_z.dat	 gasenergy4.dat   gasvy12.dat  summary2.dat
gasdens0.dat	 gasenergy5.dat   gasvy13.dat  summary3.dat
gasdens10.dat	 gasenergy6.dat   gasvy1.dat   summary4.dat
gasdens11.dat	 gasenergy7.dat   gasvy2.dat   summary5.dat
gasdens12.dat	 gasenergy8.dat   gasvy3.dat   summary6.dat
gasdens13.dat	 gasenergy9.dat   gasvy4.dat   summary7.dat
gasdens1.dat	 gasvx0.dat	  gasvy5.dat   summary8.dat
gasdens2.dat	 gasvx10.dat	  gasvy6.dat   summary9.dat
gasdens3.dat	 gasvx11.dat	  gasvy7.dat   torq_1d_Y_raw_planet_0.dat
gasdens4.dat	 gasvx12.dat	  gasvy8.dat   torq_planet_0.dat
gasdens5.dat	 gasvx13.dat	  gasvy9.dat   tqwk0.dat
gasdens6.dat	 gasvx1.dat	  grid000.inf  used_rad.dat
gasdens7.dat	 gasvx2.dat	  IDL.var      variables.par
gasdens8.dat	 gasvx3.dat	  mass.dat     vx0_2d.dat
gasdens9.dat	 gasvx4.dat	  momx.dat     vy0_2d.dat
gasenergy0.dat	 gasvx5.dat	  orbit0.dat
gasenergy10.dat  gasvx6.dat	  planet0.dat

Dealing with the parameters of a simulation (e.g: Nx, SigmaSlope, AspectRatio, OmegaFrame, OutputDir, etc.) can be simplified by defining a class to store them. These parameters are written in outputs/fargo/variables.par:

In [8]:
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

And this is how we use it:

In [9]:
P = Parameters("variables.par")
print P.nx, P.ny, P.omegaframe, P.outputdir, P.sigmaslope, P.sigma0
384 128 1.0005 ./outputs/fargo/ 0 0.00063661977237

You can use autocompletion, e.g. type P. and press "TAB". We can also print the value of all the parameters:

In [10]:
P._params
Out[10]:
{'ALPHA': '0',
 'ASPECT': '"auto"',
 'ASPECTRATIO': '0.05',
 'AUTOCOLOR': '1',
 'CFL': '0.44',
 'CMAP': '"magma"',
 'COLORBAR': '1',
 'COORDINATES': '"cylindrical"',
 'CS': '1',
 'DAMPINGZONE': '1.15',
 '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': '1',
 'KILLINGBCCOLATITUDE': '-0.2',
 'MASSTAPER': '0',
 'NINTERM': '20',
 'NOISE': '0',
 'NSNAP': '0',
 'NTOT': '100',
 'NU': '1e-05',
 'NX': '384',
 'NY': '128',
 'NZ': '1',
 'OMEGAFRAME': '1.0005',
 'OORTA': '-0.75',
 'ORBITALRADIUS': '0',
 'OUTPUTDIR': '"./outputs/fargo/"',
 'PERIODICY': '0',
 'PERIODICZ': '0',
 'PLANETCONFIG': '"planets/jupiter.cfg"',
 'PLANETHEATING': '0',
 'PLANETMASS': '0',
 'PLOTLINE': '"field[:,:,0]"',
 'PLOTLOG': '1',
 'REALTYPE': '"float64"',
 'RELEASEDATE': '0',
 'RELEASERADIUS': '0',
 'RESONANCE': '0.5',
 'ROCHESMOOTHING': '0',
 'SEMIMAJORAXIS': '0',
 'SETUP': '"fargo"',
 'SIGMA0': '0.00063661977237',
 'SIGMASLOPE': '0',
 'SPACING': '"default"',
 'TAUDAMP': '0.3',
 'THICKNESSSMOOTHING': '0.6',
 'VERTICALDAMPING': '0',
 'VMAX': '1',
 'VMIN': '0',
 'VTK': '0',
 'WRITEBX': '0',
 'WRITEBY': '0',
 'WRITEBZ': '0',
 'WRITEDENSITY': '0',
 'WRITEDIVERGENCE': '0',
 'WRITEENERGY': '0',
 'WRITEENERGYRAD': '0',
 'WRITETAU': '0',
 'WRITEVX': '0',
 'WRITEVY': '0',
 'WRITEVZ': '0',
 'XMAX': '3.14159265358979',
 'XMIN': '-3.14159265358979',
 'YMAX': '2.5',
 'YMIN': '0.4',
 'ZMAX': '1',
 'ZMIN': '0'}

Let's plot some scalar fields:

In [11]:
%pylab inline
n   = 5
rho = fromfile("gasdens{0:d}.dat".format(n),dtype='float64').reshape(P.ny, P.nx) #change dtype to 'float32' if your simulation is single precision
Populating the interactive namespace from numpy and matplotlib

The reshape instruction is needed because the scalar fields gasdens, gasvx, gasvy, etc are stored in a 1D array of doubles but the data is 2D in this particular example.

In [12]:
figure(figsize=(10,5))
imshow(log10(rho), origin='lower', aspect='auto',extent=[P.xmin,P.xmax,P.ymin,P.ymax])
xlabel(r"$\varphi$")
ylabel(r"$r$")
Out[12]:
Text(0,0.5,'$r$')

To study the temporal evolution of the fields, you need to load different outputs:

In [13]:
figure(figsize=(10,5))
for i in range(5):
    mean_rho = fromfile("gasdens{0:d}.dat".format(i)).reshape(P.ny,P.nx).mean(axis=1)
    #We plot the azimuthally averaged surface density
    plot(linspace(P.ymin,P.ymax,P.ny),mean_rho)
xlim(P.ymin,P.ymax)
Out[13]:
(0.4, 2.5)

We can avoid using linspace for the domain and use the real domain of the fields by loading the "domain_x/y/z.dat" files:

In [14]:
xmin = loadtxt("domain_x.dat")       #face-centered coodinates
ymin = loadtxt("domain_y.dat")[3:-3] #face-centered coordinates. Ghosts are included, so we remove them from the array.
xc   = 0.5*(xmin[1:]+xmin[:-1]) #cell-centered
yc   = 0.5*(ymin[1:]+ymin[:-1]) #cell-centered

We can compare the corrdinates of these arrays with those stored in the parameters file

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

As a final example, let's compute the velocity: $$\displaystyle{v = \sqrt{v_x^2+v_y^2}}$$

FARGO3D writes both the x-velocity (which is the azimuthal velocity) and y-velocity (radial velocity). These two fields have the same number of cells. However, the x- and y-velocities are staggered. A centering is required before summing them:

In [16]:
n    = 5
vphi = fromfile("gasvx{0:d}.dat".format(n)).reshape(P.ny,P.nx)
vrad = fromfile("gasvy{0:d}.dat".format(n)).reshape(P.ny,P.nx)
#And we center them
vphic = 0.5*(vphi[:-1,1:]+vphi[:-1,:-1])
vradc = 0.5*(vrad[1:,:-1]+vrad[:-1,:-1])
In [17]:
figure(figsize=(10,5))
v = sqrt(vphic**2+vradc**2)
imshow(log10(v), origin='lower', aspect='auto', extent=[xc.min(),xc.max(),yc.min(),yc.max()])
xlim(-1,1)
ylim(0.6,1.4)
Out[17]:
(0.6, 1.4)