Welcome

Welcome to the wiki for the FLASH PP Pipeline.

1. Requirements

This pipeline requires the following packages to work:

Python 2.7 / 3.6+ with following packages:
- Astropy (http://www.astropy.org/)
- Cython (https://cython.org/)
- H5py (https://www.h5py.org/)
- mpmath (http://mpmath.org/)
- NumPy (http://www.numpy.org/) (>=v1.10.4)
- SciPy (https://www.scipy.org/)
- Setuptools (https://pypi.org/project/setuptools/) (>=v30.0.0)
- yt (https://yt-project.org/)
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 original molecules_%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 containing gas_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 containing gas_temperature.inp, numberdens_hp.inp and numberdens_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 contains gas_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 either 0 or 1 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 the gas_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 the gas_density.inp file and
dust_to_gas which gives the conversion factor and is set to 0.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 the molecule_%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

Wiki

FLASH PP Pipeline / Home

Welcome

Table of contents:

1. Requirements

2. Data extraction

2.1 Simulation data

2.1.1 HDF routine

2.1.2 FLASH routine

2.1.3 Oct-tree routine

3. Functions

`add_high_temp(fpath, dpath)`

`create_para_ortho(fpath)`

`create_cp_e(fpath, script_folder, threshold)`

`create_microturbulence_mol(fpath, mol_mass)`

`create_co(fpath, freezeout, crir=1.3E17)`

`create_spintemp_wf(fpath, wf)`

`create_dust_dens(fpath, dust_to_gas)`

`make_wavelength_file(lam, delv, nlam)`

4. Input parameters and variables

Flags

Wavelength settings

Structure

Orientations

RADMC3D call

RADMC3D settings

Lines settings

5. Using the pipeline