Wiki
Clone wikiFLASH PP Pipeline / Home
Welcome
Welcome to the wiki for the FLASH PP Pipeline.
Table of contents:
1. Requirements
This pipeline requires the following packages to work:
-
Python 2.7 / 3.6+ with following packages:
-
RADMC-3D (http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/)
2. Data extraction
In this chapter I present the different routines which are used to extract the main data from the FLASH simulations. There are two main different ways to extract the data from the plot or checkpoint files of FLASH. The first one directly uses the output files from FLASH to extract the data while the other one uses QuickFLASH write out an oct-tree file which contains all data which we want to write out.
2.1 Simulation data
For writing out the input files directly from simulation we have two routines. The first one uses YT and writes out the data in a layered adaptive mesh refinement (AMR) grid while the second one uses an oct-tree AMR grid where we can also specify a region.
2.1.1 HDF routine
As stated above this routine is based on the YT package which maps the FLASH data on a layered AMR grid and writes out all available datasets. If you want to extract a specific dataset which is not included in the pipeline you can define it in the file scripts.py in the scripts folder. In order to do this you need to define a function:
#!python def _dens(field, data): return data['dens']
and add the following lines to the yt2radmc function:
#!python unit_dens = pf.field_info.get(('flash', 'dens')).units pf.add_field("Density", function=_dens, units=unit_dens, sampling_type='cell')
as well as:
#!python writer.write_line_file("Density", "gas_density.inp")
2.1.2 FLASH routine
This routine uses the oct-tree structure which can save some space and is also able cut out a region. For this the routine finds all blocks which lie within the given region and expands these to a cubic box. This can cause shifts to the center which should be kept in mind if you want to compare this with other data. To specify another dataset to be written out, you can just add the line:
#!python writeout_gas( data=pf_raw['dens'][()], fname='gas_density.inp', block_list=blist_maxref, leaf_nr=leafs, clist=counter_list, block_size=block_size, coords=coords, xmin=xmin, xmax=xmax, reduc=reduc, ftype=ftype )
where you only need to modify the field name (or a combination of different fields) in data and the filename of the file which is written out.
2.1.3 Oct-tree routine
The oct-tree routine is based on the QuickFlash routine map_FLASH_data. This writes out all given fields for a region around a given center. To add more fields you'll have to modify the QuickFlash routine as well as the silcc2radmc function in the scripts.py file.
3. Functions
In the following we present all functions included in this pipeline.
add_high_temp(fpath, dpath)
The function add_high_temp
modifies the input file molecules_%s.inp
so that the collision rate at high temperatures (in this case it is set
to 1E30 K) is given by the last given entry for the different
transitions. This gives a constant collision rate from T_max up
to 1E30 K.
The input parameters for this function are:
-
fpath
which is the path to the originalmolecules_%s.inp
file and -
dpath
which is the destination path of the modified file.
create_para_ortho(fpath)
The function create_para_ortho
calculates the number density for
para-/ortho-H2 [numberdens_p-h2.inp
, numberdens_o-h2.inp
] from
the gas temperature [gas_temperature.inp
], number density of H2
using a routine written by Dr. Annika Franeck.
At first we calculate the number density of para- and ortho-H2
depending on the temperature.
This is done accordingly to Rachford et al. (2009), Gerlich (1990).
If the ratio between ortho- and para-H2 exceeds the 3:1 limit we set the number density for para- and ortho-H2 to the 3:1 ratio.
The input parameters are:
fpath
which is the path to the folder containinggas_temperature.inp
,numberdens_h2.inp
.
create_cp_e(fpath, script_folder, threshold)
This function we calculate the modified number density for C+ ,
12C+ [numberdens_12cp.inp
] and number density for electrons
[numberdens_e.inp
] using the conversion factor presented in
Sutherland & Dopita (1993)
for cells with temperatures higher than the
threshold
temperature T_th.
The input parameters are:
-
fpath
which is the path to the folder containinggas_temperature.inp
,numberdens_hp.inp
andnumberdens_cp.inp
, -
script_folder
which is the script folder contains the file with correction factors according to @1993ApJS...88..253S and -
threshold
which describes the lower temperature threshold for which the @1993ApJS...88..253S correction is applied. The threshold temperature has to be higher than 1E4 K, as the conversion factor is only defined for temperature higher than this.
create_microturbulence_mol(fpath, mol_mass)
This function creates the microturbulence.inp
file from the
gas_temperature.inp
file by using the following equation:
Hereby we set M_mol to 2.3 as our standard value but
this can be changed with mol_mass
.
The input parameters are:
-
fpath
as the path to the folder which containsgas_temperature.inp
and -
mol_mass
which is the molar mass of the particle for which we want to calculate the microturbulence.
create_co(fpath, freezeout, crir=1.3E17)
The create_co
function calculates the number densities of CO,
13CO and C18O with or without the freezeout effect from the
number densities of H2, H, and CO as well as the dust temperature.\
Hereby the number densities of 13CO and C18O are calculated in
the following way, where the ratios are taken from Wilson & Rood (1994)
and Wilson (1999):
The freezeout effect is given by a factor which is multiplied to the number densities if the effect is to be considered.
The input parameters are:
-
fpath
is the path to folder contain all the different number density and gas temperature files, -
freezeout
is an option which can be set to either0
or1
to exclude or include the freezeout effect respectively and -
crir
is the cosmic ray ionization rate which is set to by default.
create_spintemp_wf(fpath, wf)
For HI observation we can also calculate the spin temperature with or
without the Wouthuysen-Field effect. This needs the gas temperature
(gas_temperature.inp
), number density for atomic hydrogen
(numberdens_ha.inp
) and electrons (numberdens_e.inp
) as input files.
Based on the paper from C.-G. Kim, Ostriker, W.-T. Kim (2014) we calculate the
spin-temperature with the following equations:
-
fpath
which gives the path to the folder containing thegas_temperature.inp
, (numberdens_ha.inp
) and (numberdens_e.inp
) file and -
wf
which enables/disables the Wouthuysen-Field effect.
create_dust_dens(fpath, dust_to_gas)
This function creates the from the gas_density.inp
file according to
the given dust_to_gas
-ratio. The input parameters are:
-
fpath
which gives the path to the folder containing thegas_density.inp
file and -
dust_to_gas
which gives the conversion factor and is set to0.01
by default.
make_wavelength_file(lam, delv, nlam)
This function creates the necessary camera_wavelength_micron.inp
file
which contains all wavelengths for which RADMC3D will create images.
Here fore we assume that we want to observe a spectrum centered around
the resting frequency of the molecule transition with a range given in
km s^-1. The input parameters are:
-
lam
which is the resting wavelength of the molecule transition which we want to observe, -
delv
is the width of the spectrum given in km s^-1 and -
nlam
is the number of channels we want to observe. This will be always set to an odd number so that the resting wavelength is in his channel.
4. Input parameters and variables
To configure your simulation for your needs you have a range of input parameters which you have to set or change. These are grouped in different block which we show in the following. Each section in the settings file has to start with the section name in caps and end in section name plus 'END'. For example:
FLAGS
...
FLAGSEND
All these settings need to be written in a file which path you have to give to the run_scripts.py as an argument. The path can also be relative as the scripts converts it to an absolute path internally.\ In the following sections these blocks are filled with example settings which you can modify to your needs.
Each flag setting has the following form:
flag_name [parameter1,parameter2,...]
i.e. if several parameters for one flag type are used, they need to be separated by a comma.
Flags
In the flags section we set all main settings:
-
mol_type
has to contain the different molecule type for which we want to create synthetic observations (e.g. mol_type ['co','13co']). The spelling has to be consistent with themolecule_%s.inp
file. It can also be set to a string directly if you only observe one molecule. -
wf
represents the Wouthuysen-Field effect which can be turned off or on (e.g. wf [0,1]). This will always be set to 0 if the molecule type is not atomic hydrogen ('ha') and the spin temperature is not calculated. -
spin_temp
represent the spin temperature which we can calculate for atomic hydrogen which is an effect that can be included or not (e.g. spin_temp [0,1]). For all other molecules this is effect will always be set to 0. -
dust_to_gas
is the dust to gas ratio which by default is set to 0.01 but can be changed with this flag (e.g. dust_to_gas [0.01, 0.1]). -
chemnet
tells us which chemical network was used in the simulations. For now only Nelson & Langer 1997 ('NL97') and 1999 ('NL99') are included. This has to be set otherwise the routine will stop (e.g. chemnet 'NL97'). -
dtype
gives the data type of the given data files. For now we only can read in octtree (i.e. dtype 'octtree') and FLASH (hdf5) data (i.e. dtype 'flash'). -
coords
gives the coordinates of the region we are interested in. For now this is only implemented for the data type octtree and flash. The coordiantes have to be aranged in the following order: for octtree [dx,dy,dz,x1,x2,y1,y2,z1,z2] and for flash [xmin,ymin,zmin,xmax,ymax,zmax]. If one is missing the routine will stop. You can also set coords to [] if no region should be cut out. -
mode
can be set to 'parallel' or 'serial'. In parallel mode the job is divided multiple steps which makes it possible to run multiple instance at once, but these will need to be started by hand. in serial mode you the start the routine once and it goes through all simulations one by one which makes this process very time consuming. -
mol_mass
gives the molecular mass for calculating the microturbulence. This is set to 2.3 by default but can be changed with this flag (e.g. mol_mass [2.1]). -
freezeout
tells the script if we want to include freeze-out or not (e.g. freezeout [0,1]). -
crir
is the cosmic ray ionization rate which can be changed to any value and is set to 1.3E-17 by default. This will only be considered if freeze-out is included (e.g. crir [1.3e-17,2.5e-17]) -
ftype
is format in which the output will be written. Options are 'ascii'. -
cuboid
allows to cut out a non-cubic region as specified by the coordinates. Set with 'True' -
reduc
When set to 'True' it specifies for dtype 'flash' that everything outside of the given coordinates is set to 1e-60. -
threshold
is a CII specifies flag which specifies that above the temperature given by its value all carbon is put into the form of CIII, i.e. double-ionized.
For all flags which show a list in its example (excluding the coordinates flag) you an insert also a single value instead of the list. If an array with length bigger than 1 is given, we create a list of settings over all possible and scientific useful permutations.
Wavelength settings
Another import settings block is the one for the wavelength settings.
WAVELENGTHSETTINGS
nlam [...]
widthkms [...]
lam_0 [...]
WAVELENGTHSETTINGSEND
This one only includes three flags which are:
-
nlam
which gives the number of wavelength bins and will always be increased to an odd number if one is not given, so that the resting wavelength is in the central bin and not split over 2 bins. -
widthkms
gives the width of the spectrum in km s^-1. -
lam_0
is the rest wavelength which is always at the center of the spectrum.
Structure
The structure block is tells the routine where the data files are located within the data folder. For example: This would need the following section in the settings files, to include all folders;
STRUCTURE
MC1
MC2
STRUCTUREEND
You can also omit folders for which you don't want to run the simulations.
Orientations
To determine the direction from which you want to observe your simulation you have to define a block like this.
ORIENTATIONS
incl 10 phi 0
ORIENTATIONSEND
If you want to just keep the standard orientation just leave the block empty. It should be noted that these names will be modified for folder creation and spaces are replaced with underscores.
RADMC3D call
This section contains all the RADMC3D commands you have to use when
calling RADMC3D. For example, if you want to make an image over all the
wavelengths defined in the camera_wavelength.inp
file with a
resolution of 512 pixel in each direction:
RADMC3DCALL
image
loadlambda
npix 512
RADMC3DCALLEND
You can also write all of that in one line, just like calling RADMC3D in the console, if you prefer it that way.
RADMC3D settings
Instead of parsing the options all over the command line RADMC3D also
has a separate input file radmc3d.inp
which can contain settings. To
parse these information you have write the section as follows:
RADMC3DSET
lines_mode = 3
lines_nonlte_maxiter = 100
lines_nonlte_convcrit = 0.001
tgas_eq_tdust = 0
scattering_mode_max = 0
incl_dust = 0
incl_lines = 1
catch_doppler_resolution = 0.05
RADMC3DSETEND
This is of course just an example and can be configured to your needs.
Lines settings
At last there are the lines settings which contains the information which molecules need to modeled for the simulation and radiative transfer. An example for this would be:
LINESSETTINGS
2
1
c18o leiden 0 0 1
h2
LINESSETTINGSEND
5. Using the pipeline
To use this pipeline the python packages listed in chapter 1 need to be installed and you need to create/ modify the following files:
-
create settings file (name can be chosen by yourself but needs to conform to layout given in chapter 4
- If you are using the example settings file please check if the flag dtype corresponds to your data files (see chapter 2). If your simulation data are FLASH output files you can leave the parameter at 'flash'.
-
modify header of run_data_parallel.sh and run_data_parallel.sh to fit your job submission standards
The pipeline is supposed to be used in parallel mode to be able to run multiple instances of RADMC-3D at once. In this case there are two steps to run the pipeline.
-
create all data in parallel
-
run all simulations in parallel
For the first step you use the command (in one line):
python <path to run_script.py> <path to data folder> <path to scripts folder>
<path to settings file>
All of these paths can be relative as the scripts will convert them to absolute path internally. This will create all shell scripts for the data creation of the different simulations which will be written to the folder where the scripts is executed. They have names in the following form run_data_parallel_i.sh where i is the corresponding folder number.
After these are executed and have finished all data creation, another shell script is created for each of the different simulations. These are called run_radmc_parallel_i_j.sh where i is the folder number and j a reference number to the direction of that simulation in the same order given in the settings file.
In case that for some reason the files run_radmc_parallel_i_j.sh have to be re-written but the input data themselves not, the pipeline can be called with the add-on radmc_parallel at the end of the call. This will create an run_radmc_parallel_i_j.sh file for each file located in the data folder (note that this can lead to many more *sh files than required if the data folder contains additional files like output during previous calculations.
python <path to run_script.py> <path to data folder> <path to scripts folder>
<path to settings file> radmc_parallel
Updated