| mocassin(1) | mocassin man page | mocassin(1) |
mocassin - Monte Carlo Simulations of Ionised Nebulae
mocassin, mocassinWarm, mocassinOutput, mocassinPlot
mocassin is a fully 3D or 2D photoionisation and dust radiative transfer code which employs a Monte Carlo approach to the transfer of radiation through media of arbitrary geometry and density distribution. It was originally developed for the modelling of photoionised regions like HII regions and planetary nebulae and has since expanded and been applied to a variety of astrophysical problems, including modelling clumpy dusty supernova envelopes, star forming galaxies, protoplanetary disks and inner shell fluorence emission in the photospheres of stars and disk atmospheres. The code can deal with arbitrary Cartesian grids of variable resolution, it has successfully been used to model complex density fields from SPH calculations and can deal with ionising radiation extending from Lyman edge to the X-ray. The dust and gas microphysics is fully coupled both in the radiation transfer and in the thermal balance.
mocassin is the main MOCASSIN driver that is used to start a new simulation. mocassinWarm resumes an interrupted simulation. mocassinOutput runs the output routines using the current grid files in the output subdirectory. mocassinPlot uses the current grid files in the output/ subdirectory to create 3d-emission maps.
The benchmark problems
Once you have downloaded, unpacked and compiled mocassin successfully you can first of all try to run one or more of the available benchmark problems. These were designed at a series of workshops on photoionised plasma codes held at Meudon (France) and Lexington (USA). The latest version of these benchmarks are included in Pequignot et al. (2001) in Spectroscopic Challenges of Photoionised Plasmas, ASP Conference Series Vol 247; G. Ferland and D.W. Savina eds (http://adsabs.harvard.edu/abs/2001ASPC..247.....F). Should these conference proceedings not be available to you, the same benchmarks are also reported in Ercolano et al., 2003, MNRAS, 340, 1136 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E). Here you will also find some more general info about mocassin. (There is also a link to this paper on the publications page (http://mocassin.nebulousresearch.org/publications)).
Please note that the mocassin results included in the Pequignot et al. (2001) (http://cdsads.u-strasbg.fr/abs/2001ASPC..247..533P) were only from a very early version of the code. Best results are those listed in the Ercolano et al., 2003, MNRAS, 340, 1136 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E) paper.
Input files for benchmarks
Copy the input file of the benchmark you wish to run (you should have found these in a subdirectory called benchmarks, also included in the tar ball) into the mocassin.X.XX/input subdirectory and rename it input.in or link the input.in file to the benchmark file.
e.g. for the Meudon standard HII region (Tstar=40kK) (assume mocassin is in ~/mocassin.X.XX/ ) type
cp ~/mocassin.X.XX/benchmarks/gas/HII40/meudonHII40.in
~/mocassin.X.XX/input/input.in
or
ln -s ~/mocassin.X.XX/benchmarks/gas/HII40/meudonHII40.in
~/mocassin.X.XX/input.in
also copy the respective nebular abundance file into the same directory (no need to rename!). In the case of the Meudon standard HII region benchmark this will be abunHII40.in (you can check the name of the abundance file, which is given in the input.in file with keyword nebComposition)
cp ~/mocassin.X.XX/benchmarks/gas/HII40/abunHII40.in
~/mocassin.X.XX/input/abunHII40.in
Run a mocassin pure photoionisation benchmark
The command you need to use to actually run mocassin will depend on how your system is set up, and whether you are required to run your models through a queueing system, in which case you will probably need some sort of shell script to do that (ask your network administrator if unsure).
If you want to run mocassin outside any queueing system then you can try one of the following commands and hope that one works, if not ask your network administrator about how to interactively run MPI programs on your system.
For the SUN:
mprun -np N ./mocassin (where N=number of processors you wish to use)
for the IBM
./mocassin -procs N (where N=number of processors you wish to use)
For the SGI, Beowulf and single Linux PCs using the standard MPI as above:
mpirun -np N ./mocassin (where N=number of processors you wish to
use)
Checking that everything is going OK
If the 'output' keyword is included in the input file, at the end of each iteration mocassin will write out a number of output files in the output/ directory. Some of them are written out regardless of whether the 'output' keyword is included, these are grid files that are needed by the mocassinWarm driver (should you decide to stop the simulation and re-start it from the same iteration); see the 'writeGrid' keyword to see how to control the output grid files. Also mocassin will keep you updated on what he is doing by outputting stuff to the screen (if that annoys you, just pipe it to a file e.g. ./mocassin -procs N > log; and then you can check the log at any time you wish).
After each iteration it will tell you the percentage of grid-cells converged so far, as well as the number of 'dark cells' (i.e. those cells not reached by any of the photon packets) and also will provide a summary of eventual convergence problems generally this is the last bit of info written before the beginning of the new iteration, so, if you have decided to pipe the screen output into a file called log, then you can easily check the progress so far by typing tail -f log, or grep Summary log, or more simply, the convergence history is summarised in a file called 'output/summary.out'.
For the benchmark cases you can start having an idea of whether mocassin is working properly by waiting until it's about 60-70% converged and then checking the 'output/lineFlux.out' file. A description of this file is given in output files. Also check the 'output/temperature.out' and the 'output/ionratio.out' files (also described in output files).
Completing the benchmark runs
The benchmark input files included are set up to run a simulation of 13*13*13 grid cells (this is sufficient to reproduce the benchmark results; see Ercolano et al., 2003 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E)) and a given starting number of energy packets (e.g. this can be set by the keyword nPhotons 100000); one iteration of 100000 packet should run in a few minutes on most systems.
The benchmarks input files make use of the autoPackets option, which automatically increases the number of energy packets to be used as soon as a convergence plateau is reached.
This could also be done manually, by excluding the autoPackets option from the input file and starting off with about 1000000 packets. However 1000000 photons are not enough to reach full convergence in this case and you will have to stop the simulation when the percentage convergence reaches 30-40%, then edit the grid3.out file changing 1000000 to 3000000, and restart the simulation using the previously compiled mocassinWarm driver, e.g.
mprun -np N ./mocassinWarm
This will basically tell mocassin to increase the number of energy packets to be used in the simulation from 1000000 to 3000000. It will, of course, now take longer to complete each iteration, but the convergence percentage should increase quite quickly. You can stop the simulation when you reach convergence > 90%. For simplicity, until you familiarise yourself with the program, it is advised that you keep the autoPackets option as set in the input files, in which case it will not be needed to manually stop and restart the simulation.
Comparing with the benchmark tables
The description of the output files is given in output files. The output files contain all the info (and much more) you need to compare your simulation with the benchmark tables. Remember that mocassin employs a statistical method, so do NOT expect to see exactly the same figures as those quote in the benchmark paper as some of the differences will be due to the statistical error and also to the convergence level/limit employed. Furthermore some of the atomic data has been updated since the publication of the paper.
Version 2.0 was designed for the modelling of regions where dust grains and photoionised gas coexist in the same volume; however extreme care has been taken so that the code could also be run efficiently when only one or the other process is included. In this section we will see briefly how to run dust-only models by attempting to reproduce pure dust 1D and 2D models as shown by Ercolano, Barlow and Storey (2005, MNRAS, 362, 1038) (http://cdsads.u-strasbg.fr/abs/2005MNRAS.362.1038E).
1D Benchmarks
A set of spherically symmetric benchmark models and solutions are described by Ivezic et al. (1997, MNRAS, 291, 121) (http://cdsads.u-strasbg.fr/abs/1997MNRAS.291..121I). The input files used by mocassin for some of these benchmarks are included in the directory ~/mocassin.X.XX/benchmarks/dust/1D.
Copy the input file of the benchmark you wish to run into the mocassin.X.XX/input subdirectory and rename it input.in or link the input.in file to the benchmark file. Then run the code (see Section 3.1.3) and finally compare the output files (SED.out dustGrid.out, see output files) with the results published by Ercolano, Barlow and Storey (2005, MNRAS, 362, 1038) (http://cdsads.u-strasbg.fr/abs/2005MNRAS.362.1038E)
2D Benchmarks
A set of 2D disk benchmark models and solutions are described by Pascucci et al. (2004, A&A, 417, 793) (http://cdsads.u-strasbg.fr/abs/2004A%26A...417..793P). Sample input files used by mocassin for some of these benchmarks are included in the directory ~/mocassin.X.XX/benchmarks/dust/2D.
As for the 1D benchmarks you will have to copy the input files of the benchmark you wish to run into the mocassin.X.XX/input subdirectory and rename it input.in or link the input.in file to the benchmark file. Then run the code (see Section 2.1.3) and finally compare the output files (SED.out dustGrid.out, see Section 4) with the results published by Ercolano, Barlow and Storey (2005, MNRAS, 362, 1038) (http://cdsads.u-strasbg.fr/abs/2005MNRAS.362.1038E).
Once you have convinced yourself that mocassin is behaving as it should, you can run your own simulations. If you are not confident that this is the case, please contact me (http://mocassin.nebulousresearch.org/contact) before going any further.
The first step is to define your model. This is done via a set of keywords included in the input.in file. YOU SHOULD NEVER HAVE TO CHANGE ANYTHING IN THE SOURCE CODE (I hope). If you find you have to, please let me know. The rest of this section will list and describe all the keywords that can be included in the input file. The number of keywords you decide to use and the order you decide to list them in is entirely up to you. If you have missed out some required and fundamental keyword the simulation will stop and tell you what you have missed out. However, you should be aware that, for most of these keywords, default values are defined, so make sure that the default value is actually what you want before you decide to leave a given keyword out. Default values to each keyword are given below, enclosed by square brackets. The fundamental keywords, which must be specified in every input file are marked by an asterisk in the square brackets and have no default value.
List of keywords
The number of energy packets to be used in the next iteration will automatically be increased by a factor of real2 whenever a convergence plateau (defined by real1) is reached, i.e. if the convergence level increase is less than the value specified by the user in real1. For example autoPackets 0.10 5. 1e8 will cause the total number of energy packets to increase by a factor of 5. whenever the convergence level increase from one iteration to the other is less than 10%. This only occurs when the convergence level is less than 85% and is the maximum number of packets defined by integer3 (1e8 in this example) has not been reached.
Default value: .false.
This keyword creates a 3D cube of the escaped continuum radiation in direction given by the inclination keyword and over all directions. If no inclination is specified in the input the continuum cube will include packets escaping in all directions. The continuum band wavelength limits are defined by real1 and real2 (given in μm). A negative value or the omission of this keyword will result in no continuum cube being written out.
Default value: -1.
The shape of the ionising continuum. The default value is 'blackbody', in which case the ionising stellar continuum is approximated by the Planck function for the stellar temperature defined by the keyword TStellar. If 'string' == 'powerlaw' then this must be followed by a real number indicating the index of the power law distribution such that Fnu ∝ nu^(-index). E.g. contShape powerlaw 1. will generate an input spectrum following a nu^(-1) distribution. Note that the final result will be normalised to the luminosity entered. The power law keyword is not compatible with the LPhot keyword, but only with the LStar keyword. If a stellar atmosphere data file is to be used, the 'string' must specify the path of the external file containing the data. For example contShape NLTE140lg65 tells the program to look for the nLTe140lg65 data file in the current directory. The stellar atmosphere files must be in a format consisting of two columns: the first listing the wavelength points in units of [A] and the second containing the corresponding wavelength-dependent stellar Eddington fluxes in units of [erg/cm^2/s/A/sr]. A set of stellar atmosphere flux tables have been compiled by Dr T. Rauch in a mocassin-friendly format and are available from http://astro.uni-tuebingen.de/~rauch/ (but please manually remove the header from the flux tables.
Default value: blackbody
This is the convergence limit for the variation of a given parameter, in each grid cell from one Monte Carlo iteration to the next (e.g. 0.05 means changes of 5% maximum). In the case of a pure-gas (or gas+dust) simulation the criterion is based on the rate of change of neutral hydrogen. In the case of dust-only the criterion is based on the rate of change of the dust temperatures.
Other convergence criteria can also be used, at the moment, this would require a simple editing of some source modules. If you would like to use a different convergence criterion please email me (http://mocassin.nebulousresearch.org/contact) and I can do the editing for you.
Default value: 0.05
Logical switch to enable the debugging mode. When this keyword is included mocassin will calculate a number of extra quantities (see Section 5.), which will, of course, slow the process down and also require more memory.
Default value: .false.
The density structure of the nebula can be defined cell by cell by using an external density file. mocassin knows that a density file is to be used when the densityFile 'string' is included in the input file, where 'string' contains the name and path of the file where the data is stored. This file must consist of four, five or six columns, with the first three columns containing the x-, y-, and z- coordinates of the grid cell in [cm] and the fourth columns containing the value of the hydrogen density by number in [cm^{-3}] at the particular grid cell. The x, y and z axis do not to be equally spaced - irregular grids are perfectly acceptable by mocassin and also the extent of each axis can vary (as long as this is consistent with the values given in the nx, ny and nz fields). The fifth column is optional. If the multiChemistry keyword is specified the fifth column must contain an integer number in the range [1, Ncomponents] which indicates what component this cell belongs to (so that mocassin can assign the chemical abundances for this component).
It is possible to specify nx, ny and nz from within 'string2' - the first row of the file should then be
# nx ny nz
and the keywords can be omitted from input.in
Default value: none
This keyword is usually followed by a set of parameters which are to be fed into the density law routine, included in the grid_mod module. Any density law can be specified by editing the code in the setGrid subroutine. If the nebula is homogeneous, this keyword must be omitted and the Hdensity keyword included instead. Note that if neither of the two keywords is included, and an external density file is not specified with the densityFile keyword, the nebular density distribution is left undefined and the simulation halted with an error message being produced.
Default value: 0. 0. etc..
This keyword can be used in the case of a non-central source such as the heating by radioactive decay of 56Co in supernova remnants. real1 is the total luminosity of the diffuse source in [1e36 erg/s], and real2 is the temperature of the diffuse source. string1 is the shape of the source spectrum, with the same syntax as contShape.
integer1 is the number of diffuse photons to use at the start of the
model (same use as nPhotons).
integer2 is the index of the grid in which the diffuse source emits,
which allows control of the extent of the diffuse source. For example, a
medium with clumps embedded could be represented by a mother grid (index
0) and subgrids (indices 1..n). To restrict emission to the inter-clump
medium only, integer2 would be set to 0. If there are no subgrids, set
integer2 to 0.
Currently, for "legacy" reasons, a value of TStellar must still be specified (i.e. TStellar=0.) but nPhotons, LStar and contShape can be removed from the input.in file.
diffuseSource should not be used at the same time as symmetricXYZ.
Default value: *
names dust data files -
string1 = grain species file
string2 = grain size info file
Default value: "none", "none"
By default, mocassin calculates the total dust mass from the user-specified number densities, dust species and grain size distributions. If this keyword is specified, the specified number densities are scaled when read in, such that the total dust mass is equal to real1. real1 has units of [Msun]
Default value: .false.
If one needs to account for light travel time effects, for example if the source luminosity is changing, then this keyword will cause the results in SED.out to be integrated only over grid cells lying between ellipsoids corresponding to light travel times of real1 and real2. real1 and real2 are in light days, and real1 must be less than real2. The ellipsoids open out along the z-axis, and so one should also specify 'inclination 1 0.0 -1' in input.in when using this keyword. A review of the geometries of light echoes can be found in Sugerman, 2003, AJ, 126, 1939 (http://cdsads.u-strasbg.fr/abs/2003AJ....126.1939S), and a case of mocassin modelling including light travel time considerations is described by Wesson et al, 2010, MNRAS, 403, 474 (http://cdsads.u-strasbg.fr/abs/2010MNRAS.403..474W).
Default value: .false.
Defines the grid edges (only to be used with automatic grids). real1 real2 and real3 are respectively the x, y and z edges in cm
Default value: -1., -1., -1. *must be given if automatic grids are used
real1 can assume all values from 0. to 1. to defines the gas volume (and/or dust) filling factor
Default value: 1.
Logical switch to enable equivalent optical depth calculations (see Ercolano et al. 2007 (http://cdsads.u-strasbg.fr/abs/2007MNRAS.375..753E)) - useful for diffuseSource and multiPhotoSources cases. The last column in the output file equivalentTau.out gives the source SED in the same units as SED.out, which may be useful to the user as well.
Default value: .false.
Logical switch to enable optical depth calculations and output - may be time consuming for large grids
Default value: .false.
This keyword specifies a constant hydrogen density, by number, throughout the nebular region The command Hdensity 300 will e.g. set the hydrogen density, by number, to the constant value of 300 cm-3.
Default value: 0.
This keyword controls the viewing angles at which the SED will be calculated as it will appear in the SED.out file. The number of viewing angles is given first (integer; in this version integer <= 2) and then the θ and φ inclination in radians. To turn off the φ dependence, real_2 must be set to to -1. (θ=0. when the line of sight coincides with the z-axis)
Default value: 0
This indicates that the density distribution of the grid is in terms of electron densities instead of H density. This will cause the code to calculate at each iteration the values of H density from the local Ne values by taking into account the local ionisation structure.
Default value: .false.
This keyword activates the logical switch that turns off unisotropic scattering (implemented via the Heyney-Greenstein phase function - with <cos θ> calculated from the dielectric constants via MIE theory).
Default value: .false.
This is the number of hydrogen-ionizing photons emitted by the source per unit time, which is generally referred to as Q(H0), in the literature, with units of [E36 sec-1]. If this is given then the stellar luminosity, LStar is automatically derived from it.
Default value: * if LStar not given
This is the stellar luminosity in units of [E36 erg sec-1]. If this is given as an input, then the number of hydrogen- ionizing photons, Q(H0), is automatically derived from it and from the input spectrum.
Default value: * if LPhot not given
integer1 is the maximum number of Monte Carlo iterations to be performed
in the simulation. real2 is the minimum convergence level to be achieved
before the end of a simulation. The program will however stop after
integer1 iteration even if real1% of convergence has yet to be
reached.
Default value: 30 95.
Dust to Gas ratio by mass. If string1 = 'constant' then it must be followed by real1, containing the value of MdMg to be applied homogeneously to all cells in the grid. If string1=file then it must be followed by string2, the name of the file defining the MdMg at each location. Note that MdMg, MdMh and Ndust are mutually exclusive.
This file must consists of four columns, with the first three columns containing the x-, y-, and z- coordinates of the grid cell in [cm] and the fourth columns containing the dust to gas mass ratio at the particular grid cell. The x, y and z axis do not to be equally spaces, irregular grids are perfectly acceptable by mocassin and also the extent of each axis can vary (as long as this is consistent with the values given in the nx, ny and nz fields).
It is possible to specify nx, ny and nz from within 'string2' - the first row of the file should then be
# nx ny nz
and the keywords can be omitted from input.in
Default value: .false. 0./'none'
Dust to Hydrogen ratio by mass. If string1 = 'constant' then it must be followed by real1, containing the value of MdMh to be applied homogeneously to all cells in the grid. If string1=file then it must be followed by string2, the name of the file defining the MdMh at each location. Note that MdMg, MdMh and Ndust are mutually exclusive.
This file must consists of four columns, with the first three columns containing the x-, y-, and z- coordinates of the grid cell in [cm] and the fourth columns containing the dust to hydrogen mass ratio at the particular grid cell. The x, y and z axis do not to be equally spaces, irregular grids are perfectly acceptable by mocassin and also the extent of each axis can vary (as long as this is consistent with the values given in the nx, ny and nz fields).
It is possible to specify nx, ny and nz from within 'string2' - the first row of the file should then be
# nx ny nz
and the keywords can be omitted from input.in
Default value: .false. 0./'none'
This keyword must be used when a chemically inhomogeneous model is being performed. The integer1 value defined the number of components to be used. string(1), string(2) ... string( integer1) are the names of the files describing the abundances in each component. When the multiChemistry keyword is included the density distribution MUST be specified via the densityFile keyword. The fifth column of the density file must contain the index for the abundance file describing the chemical composition at each location.
Default value: .false.
This defines a multiple grid environment. The integer1 is the total number of grids to be used (mother + subgrids) and string1 is the name of the file containing the subgrid information. Please see the Running multiple spatial grids section in this document.
Default value: .false. 'none'
This keyword is used to define multiple (or single) ionising sources, that can be placed at any locations. 'string1' is the name of the file containing the central star parameters. This file must contain in the first line the number of sources to be included and then each source should be specified on successive lines providing (in this order) Teff [K], L* [E36erg/s], ContShape, x-, y-, z-position. Note that the location of the source must be given normalised to the length of the relative axis. As example is included in the example/ directory.
The total number of points to be used in the frequency mesh.
Default value: 600
Number density [cm-3] of dust grains. If string1='constant' then it must be followed by real1, containing the value of Ndust to be applied homogeneously to all cells in the grid. If string1=file then it must be followed by string2, the name of the file defining Ndust at each location. Note that MdMg and Ndust are mutually exclusive. This file must consists of four columns, with the first three columns containing the x-, y-, and z- coordinates of the grid cell in [cm] and the fourth columns containing the number density of dust in [cm-3] at the particular grid cell. The x, y and z axis do not to be equally spaced; irregular grids are perfectly acceptable by mocassin and also the extent of each axis can vary (as long as this is consistent with the values given in the nx, ny and nz fields).
It is possible to specify nx, ny and nz from within 'string2'
- the first row of the file should then be
# nx ny nz
and the keywords can be omitted from input.in
Default value: .false. 0./'none'
This keyword specifies the path of the nebular composition data file. If the default 'solar' composition (defined in the file solar.dat) is to be used, this keyword can be omitted. However the solar.dat file includes all elements with Z<30: this will result in a more memory expensive simulation. It is therefore advised to set to zero those elements which are not needed in the simulation. Otherwise the nebular composition can be specified in the user-defined 'string' file to be found in the current directory. All composition input files must be in a format consisting of one column containing the abundances by number and relative to hydrogen for the first thirty elements in order of ascending atomic number. The abundances of elements which are not to be included in the simulations must be set to zero (this will automatically exclude them by flagging them out throughout the program). If the multiChemistry keyword is specified the nebComposition keyword must be omitted.
Default value: * if gas is present and no multiChemistry
Initial guess for the nebular electron density.
Default value: 0.
When this keyword is included in the input file all procedures associated with the calculation of the grain charges and photoelectric yield are switched off as well as gas-grain collision processes.
Default value: .false.
This is the number of energy packets to be used in the Monte Carlo simulation and it has to be specified for each model.
Default value: *
This is the number of ionisation stages to be used in the model. Max allowed is currently 10. If you have the atomic data necessary and would like to use more than 10 ionisation stages please contact me (http://mocassin.nebulousresearch.org/contact), or if you are confident you can edit the data/fileNames.dat and include the new data files to the pool.
Default value: 7
High limit of the frequency mesh in units of [Ryd]
Default value: 15.
Low limit of the frequency mesh in units of [Ryd]
Default value: 0.
Number of axial points in the x-direction. Can also be specified on the first row of densityFile, Ndust, MdMg and MdMh files, in the form of a row containing
# nx ny nz
The keywords nx, ny and nz can then be omitted.
Default value: 30
Number of axial points in the y-direction. Can also be specified on the first row of densityFile, Ndust, MdMg and MdMh files, in the form of a row containing
# nx ny nz
The keywords nx, ny and nz can then be omitted.
Default value: 30
Number of axial points in the z-direction. Can also be specified on the first row of densityFile, Ndust, MdMg and MdMh files, in the form of a row containing
# nx ny nz
The keywords nx, ny and nz can then be omitted.
Default value: 30
when this keyword is included the output files will be written at the end of each iteration. If it is omitted no output will be created, however if the grid files are being created the output files can easily be recovered using the mocassinOutput driver.
Default value: .false.
This keyword is used when the ionisation source is from a plane and not from a point source. real1 must contain the value of the incident ionizing flux above υ = real2 [Ryd] (constant at each point on the plane) in units of [phot/s/cm2]. If real2 = 0. then the real1 is assumed to be the bolometric flux.
When this keyword is specified the density distribution must be defined via the densityFile option. The ionizing plane is the x-z plane. the energy packets are reflected when they hit the y-z and the y-x planes and can escape from the x-z planes. (This can be easily changed to suit. please contact me (http://mocassin.nebulousresearch.org/contact) if your work requires it to be different)
Default value: .false.
This keyword activated the temperature spiking routines (quantum grain heating). It is only valid for simulations including dust. real1 is the limiting size of grain radii [μm] that will be allowed to spike (i.e. grains with a < real1 will spike). real2 is the minimum convergence level for the spiking to occur. Please read the notes on Grain temperature spiking procedures.
Default value: 1.e-3, 99.
This keyword provides controls over some internal parameters in the quantum grain heating procedures. They should not be modified unless you really know what you are doing. real1 is the max temperature to be considered in the grain temperature mesh of the spiking. integer1 is the number of temperature (and enthalpy) bins considered. logical1 switches on and off the writing out of a file containing all the probability functions for all the grains in every cell. The resulting file can be really gigantic and so this value should be set to .true. only for debugging purposes and used with care. Please read the notes on Grain temperature spiking procedures contained in this manual.
Default value: 700., 300, .false.
Inner radius of the ionised region, in units of [cm].
Default value: *
Outer radius of the ionised region in units of [cm].
Default value: 0.
If this keyword is introduced, recombination line intensity of astrophysically abundant ions will be computed and appended to the lineFlux.out file
Default value: .false.
real tells at what level of grid convergence the resonance line photons escape fractions should be calculated. This should be included when both dust and gas are present.
Default value: 101.
This keyword will cause the results in lineFlux.out, temperature.out and ionratio.out to be integrated over only those cell that fall under the projection of a slit aligned along the z-axis of the grid. The slit x and y dimensions (in [cm]) are defined by real1 and real2 respectively. If real1 and real2 have value 0. or if they are omitted, no slit is used and the results are integrated over the whole active volume.
Default value: 0., 0.
When the nebula to be modelled shows axial symmetry in the x- y- and z-directions, this keyword can be used to enable the symmetric grid procedures. This will result in the ionizing source being put in a corner of the grid, instead of being put in the centre, meaning that only one eighth of the nebula will have to be computed.
symmetricXYZ should not be used in models illuminated by a diffuse source.
Default value: .false.
This switch enables the verbose version of the program.
Default value: .false.
Initial guess for the nebular temperature.
Default value: 10000.
Logical variable to switch on the thermal balance channel tracing. When this is included in the input file a file called thermalBalance.out will be written to the output/ directory. Be aware that depending on the size of your grid this may be quite a large file and time-consuming in the I/O phase.
Default value: .false.
Temperature in [K] of the ionizing stellar source.
Default value: *
real indicates the minimum grid convergence percentage after which the grid files will be written out.
Default value: 0.
The source files are contained in a subdirectory called source/. mocassin looks for the input.in file from a subdirectory called input/. The atomic data files should all be contained in a subdirectory named data/. Most of the atomic data files should not need to be changed at all. Unless you decide to update some of them, in which case (under the GPL agreement) you should also email me (http://mocassin.nebulousresearch.org/contact) with the changes so that I can include them in the public version of the code. The dust optical data library and other dust related data files are contained in a subdirectory named dustData/.
The user's input files may be a combination of the following files, depending on the processes included in a given simulation.
input.in
This is the main input file where you can specify all the keywords to define your simulation. Some example input files are given for the Meudon/Lexington benchmark cases (see pure photoionisation benchmarks).
gas abundances file
This is the nebular abundances file which should have the name specified by the user in the nebComposition field of the input.in file. Some sample files are given for the Meudon/Lexington benchmarks.
gas density file
This the nebular density structure file which should have the name specified by the user in the densityFile field of the input.in file. The format of this file is given in Section 4 (see densityFile).
stellar atmosphere file
This is the stellar atmosphere file which should have the name specified by the user in the contShape field of the input.in file. The format of this file is given in Section 4 (see contShape).
dust number density file
This should contain the dust number density distribution across the grid. It's name and path should be specified in the input.in file by the keyword Ndust.
dust to hydrogen or dust to gas ratio
This contains the dust to hydrogen or dust to gas ratio distribution across the grid. It's name and path should be specified in the input.in file by the keyword MdMh or MdMg, respectively.
dust species and grain size distribution files
The names of these two files must be specified in the input.in file following the keyword dustFile. The dust species file should contain a first line specifying how many species are to be included, and then successive lines containing the names of the optical data (n,k or Qs) file and the relative abundance of the species.
e.g. for a pure silicate model (using the Draine and Lee 1984 data (http://cdsads.u-strasbg.fr/abs/1984ApJ...285...89D)) :
1
'dustData/sil-dl.nk' 1.0
The grain size distribution file should contain a first line specifying how many grain sizes are to be included, the rest of the file should consist of three columns : index, radius (in um), weight. Grain size distribution files can be created using the makeGrainSizeDistribution.f90 program included in the accessories/ subdirectory.
e.g. for a single grain size
1 size
1 0.16 1.0
Input files for multigrid simulations are described in [[running multiple spatial grids.
plot.in
The plot.in file is used by the mocassinPlot driver in order to create 3D grids of line emission. This file must be place (or linked to) the input/ subdirectory. The plot.out and the grid4.out files are written out to output/ and can then be used to create emission line maps by integration along any given line of sight.
Monochromatic grids are created using the mono keyword, and individual lines using the line keyword.
For example:
mono
line 2 4861. 4861.
line 93 4686. 4686.
line 1529 5007. 5007.
line 1540 4363. 4363.
line 2407 6733. 6733.
line 2408 6718. 6718.
line 929 6583. 6583.
The integer following the keyword line is the line index number as given in the lineFlux.out file. NOTE that the line indices will be different for different simulations as they depend on which elements are included and on the number of ionisation stages accounted for. The 2nd and 3rd indices contain the central wavelength of the line (these are redundant for monochromatic plots, however they must be included).
mocassin produces several output files at various times during the simulation. This will be contained in a subdirectory named mocassin.X.XX/output/ The files ionratio.out, lineFlux.out, temperature.out, (tau.out), ionDen.out and SED.out are all produced by the output_mod module.
The plot.out and grid4.out files are produced by the mocassinPlot file
ionratio.out
Ionratio.out contains the volume averaged ionic fractions. Different authors is the past have used slightly different definitions of this quantity in their models. Please refer to Ercolano et al. (2003) (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E) for further information on the description used by mocassin.
The first two columns of the ionratio.out file give the atomic number of the element and its ionisation stage (1 for neutral, 2 for singly ionised etc.), and the third column gives the required quantity. If a multiChemistry model is being run, then the results will be given for each individual component.
lineFlux.out
The file lineFlux.out contains the volume integrated intensities of all the emission lines calculated by mocassin. These are all given relative to Hβ, which is in absolute units.
The first two columns give the element and ion number, these are followed by mocassin codes for the levels of the transition; these are followed by the wavelength in [A]. The wavelength column is followed by the analytical and Monte Carlo line intensities relative to Hβ, which is given in absolute units at the top of each region. Finally the last column gives the ratio of the two previous columns. NOTE that the Monte Carlo line intensities are only calculated if the debug keyword is included in the input file. In normal mode only the intensities calculated using the formal solution (which are in general more accurate) will be available.
If a multiChemistry model is being run, then the results will be given for each individual component.
temperature.out
The file temperature.out contains the mean electronic temperatures weighted by the ionic species (see Ercolano et al., 2003 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E) for definition). This file has the same structure as the ionratio.out file.
If a multiChemistry model is being run, then the results will be given for each individual component.
equivalentTau.out
Please see Section 2.3 of Ercolano, Barlow and Sugerman (2007) (http://cdsads.u-strasbg.fr/abs/2007MNRAS.375..753E) for a description of this quantity. The first column contains Energy [Ryd] the second column contains wavelength [μm] the third column contains equivalent tau [see paper] and the fourth column contains Fλ0 [see paper]. Note that in the case of diffuse illumination this is the only meaningful quantity - see discussion in paper.
tauNu.out
The file tauNu.out contains the frequency dependent optical depth tau(nu) [which includes ALL opacity sources] at the edge of the grid in the three positive axial directions starting from the origin of the axes.
tau.out
THIS FILE IS NOT PRODUCED ANY MORE IN VERSIONS >= 2.02.49 PLEASE CONTACT ME (http://mocassin.nebulousresearch.org/contact) IF YOU NEED IT. The file tau.out contains the run of the optical depth from the centre of the nebula to the outer edge along the three axial directions. The optical depths are calculated at the neutral hydrogen ionisation threshold, nu = 1.0 Ryd (13.6 eV), at the neutral Helium ionisation threshold, nu = 1.8 Ryd, and at the singly ionised Helium ionisation threshold, nu = 4.0 Ryd. The first column of the file gives the distance in [cm] from the centre of the nebula and the second column gives the optical depth from the centre to that point. This file is not always calculated correctly if the internal switches are not set up properly. Please use equivalentTau. See description above.
ionDen.out
The file ionDen.out contains the ionic fractions at each grid cell. The first three columns give the x-, y- and z-axis indices of the cell, the fourth and fifth columns give the atomic number and the ionisation stage of the element (as above, 1 for atom, 2 for singly ionised etc.) and, finally, the sixth column gives the corresponding ionic fraction.
SED.out
This file contains the emerging spectral energy distribution from the grid. As indicated in the files' header column 1 and column 2 contain the frequency [Ryd] and wavelength [μm] grid, column 3 contains the direction averaged SED per unit direction (must multiply by π to obtain the total overall directions); the following columns contain the SED emerging in any given line of sight as requested by the user in the input.in file with the inclination keyword.
The grid files and photoSource.out
As we have already mentioned elsewhere, grid files are also written out after each iteration by routines contained in the grid_mod module. These are needed by the warm start driver (mocassinWarm) to re-initialise an interrupted simulation. These files are formatted such that they can be written out and read back in quickly and therefore they may not be very clear to the human eye. However, most of the information they contain is also given in a more intelligible form in the other output files listed above, dust temperatures are an exception as discussed in plotting dust temperatures. Gas-only simulations will result in 4 grid files being written out: grid0.out, grid1.out, grid2.out and grid3.out.
The first line of the grid0.out file gives the number of grids included (mother+ subgrids), on the next line are the x-, y- and z-axes points in the mother grid, followed by the outer radius in [cm]. The next few lines list the x-axis points, then the y-axis points and, finally, the z-axis points. The rest of the file contains the convergence info for each grid cell. The active index of the cell in the first column, whether it has converged (1=yes,0=no), and whether it is a black cell (i.e. if could not be reached by any photons, 1=black, 0=normal). Cells that have a 0 index are inactive cells; cells with a negative index are cells that have been replaced by a subgrid, whose index is equal to the absolute value of the negative mother-cell index. In the case of multiGrid simulations, the file will loop around all subgrids included.
The grid1.out file contains the electron temperatures, electron densities and hydrogen densities for each grid cell in each grid. As for the grid0.out file, this information is given for each grid cell, with the last index varying the fastest (i.e. (1,1,1), (1,1,2), etc.).
The file grid2.out contains the ionic fractions at each grid cell for the ions included in the simulation. These are given in order of increasing atomic number and ionisation stage, with each element occupying one line. The grid cell indices vary in the same fashion as in the grid1.out file.
The file grid3.out contains a list of the specified simulation parameters in a fixed order (the keywords are indicated on the right of each value). NOTE that it is not possible any more to change nuStepSize from the input.in file. If you wish to change this parameter (you shouldn't need to), this is defined in the set_input_mod module.
Dust-only simulations only produce grid0.out, grid3.out and dustGrid.out
The file dustGrid.out contains the dust number density at each cell followed by the grain temperatures for each grain size of each grain species. For each cell Ndust is written on one line and this is followed by n_size+1 lines of n_species+1 columns containing the individual grain temperatures for each size and species, where n_size=number of size bins and n_species=number of grain species included. The average temperature for the grain mixture weighted by the size distribution and the species abundances is given at (0,0).
Dust+gas simulations will produce all the files above.
grid4.out is written out by mocassinPlot and it just contains the volume of each gridcell in [e45 cm3], which is needed for visualisation purposes.
The ionising source(s) input parameters are written out to a file named photoSource.out.
plot.out
The mocassinPlot driver produces also an output file, containing the luminosities of each individual grid volume element in the required emission lines. This file, named plot.out, is written in a format which should be easily readable into a data visualisation software, such as IDL or PDL. A grid4.out is also written out, containing the volume of each gridcell in [e45 cm3].
The total luminosity of the nebula in various emission lines longward of the Lyman limit can be obtained by using two methods. The first method, which is only available to Monte Carlo codes, consists of summing up the number of energy packets in the given line over all the grid cells. From this, the power emitted in that line can be readily obtained (see e.g. Ercolano et al., 2003 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E)). The second method consists of using the values of the local electron temperatures and ionic abundances given by the final model solution to obtain the line emissivities for each grid cell. The luminosity of the nebula in any given line can then be calculated easily by integrating the emissivity of the required line over the volume of the nebula.
A comparison of the results obtained using the two methods described above, provides an indication of the level of statistical accuracy achieved during the simulation, as the two methods will give consistent results only if enough energy packets are used in order to yield good statistics for every line. However, in general, the second method (formal solution) yields the most accurate results, particularly for weak lines, which only emit a few photons.
Calculation of the line emissivities using the first method, although straightforward, requires extra book-keeping which can be expensive for larger simulations. For this reason, this calculation will only be carried out when the keyword debug is included in the input file, otherwise the more speedy (and accurate) formal solution will only be employed.
WARNING: This section is out of date -- new help is in the process of being written as we speak -- almost ;-) . Please email me (http://mocassin.nebulousresearch.org/contact) if you want to use multiple grids and need some help
From Version 2.00 onwards you can choose to run simulations including multiple spatial sampling. A typical example of when this would be need is for the modelling of knots embedded in a gaseous nebula. The density enhancement in the knot requires a finer spatial grid than one that may be sufficient to model the main nebula. In such cases the subgrid describing the knot can be modelled simultaneously and self consistently by mocassin's multigrid routines. The radiation field will be safely transferred from mother to sub grids with (hopefully) no error being introduced by the process. The overheads involved only reflect the eventual introduction of extra grid cells.
A typical example of how such a model would be set up is briefly described below. All files mentioned should be included in the examples/multigrid subdirectory.
The example below shows the input.in file for a plane parallel model consisting of a cubic blob embedded into a cubic grid:
autoPackets 0.20 2. 1000000000
output
contShape blackbody
nebComposition 'examples/abunHII40.in'
maxIterateMC 30 95.
nPhotons 10000000
nx 10
ny 10
nz 10
planeIonization 3.0
Rin 0.
Rout 2.1e19
TStellar 45000.
writeGrid 5.
convLimit 0.03
densityFile 'examples/cubenew.dat'
multiGrids 2 'examples/subGrids.in'
The density distribution of the mother grid (containing the main nebula) is provided provided by an external file using the densityFile keyword as described in Section 3.1; e.g.
densityFile 'examples/cube.dat'
mocassin knows that it will have to deal with more than one grid (in this case 2) thanks to the line
multiGrids 2 'examples/subGrids.in'
the integer after the multiGrids keyword is the total number of grids including the mother grid; the file 'example/subGrids.in' contains all the information regarding the position of the subgrids and their specification. This file is created automatically by using the makeSubGrids.f90 program included in the 'accessories/' directory.
The makeSubGrids.f90 program reads the makeSubGrids.in input file, that has the following form
1
cube.dat 10 10 10 1
cubenew.dat
1 9 9 9 knot_den.in
7.e18 14.e18 7.e18 14.e18 7.e18 14.e18
100. -1
line 1: number of subgrids to be included line 2: name of mother grid file; its x,y,z dimensions; multiChemistry? (1=no; >1=yes) line 3: name of the modified mother grid file (to be created by makeSubGrids.f90) line 4: knot index; its x,y,z extent; name of knot file (to be created by makeSubGrids.f90) line 5: xmin,xmax,ymin,ymax,zmin,zmax in cm, defining the region occupied by knot 1 line 6: H number density in [cm-3] of knot 1, abundance index (negative if no multiChemistry) The maximum number of subgrids that may be included in a simulation is controlled by the parameter maxGrids in the source/constants_mod.f90 module; in the example above only one subgrid is included, additional subgrids could be included by simply repeating lines 4,5 and 6 for each subgrid.
The knot_den.in file above has the following format
0.0000000E+00 0.0000000E+00 0.0000000E+00 100.0000
0.0000000E+00 0.0000000E+00 0.1250000 100.0000
0.0000000E+00 0.0000000E+00 0.2500000 100.0000
0.0000000E+00 0.0000000E+00 0.3750000 100.0000
0.0000000E+00 0.0000000E+00 0.5000000 100.0000
0.0000000E+00 0.0000000E+00 0.6250000 100.0000
0.0000000E+00 0.0000000E+00 0.7500000 100.0000
0.0000000E+00 0.0000000E+00 0.8750000 100.0000
0.0000000E+00 0.0000000E+00 1.000000 100.0000
0.0000000E+00 0.1250000 0.0000000E+00 100.0000
0.0000000E+00 0.1250000 0.1250000 100.0000
etc ........
Columns 1, 2 and 3 contain the x, y, and z positions in the subgrid, normalised to 1. The fourth column contains the h number density. (when running a multiChemistry model a fifth column would appear; containing the abundance file index).
It is worth noticing that since the x y and z positions in the knot_den.in files are given normalised to unity, we can include knots with the same geometry and density at several grid locations without the need to create a density file for each knot. this can be done by specifying the same filename for the knots in the makeSubgrids.in file.
The makeSubgrids.f90 file will also process the cube.dat file (which is the original mother grid density file provided by the user) turning off all the cells that are to be replaced by subgrids. The new file will be called cubenew.dat in the above example and this is the file that must be specified in the input.in file of mocassin.
The makeSubgrids.f90 file can only create subgrids of homogeneous density and chemical indices. It should be very easy for the user to customise the grids obtained to included their 'favourite' density/chemistry distribution.
The output from multigrid models will be summarised in the usual files and an overall, as well as grid by grid result will be provided.
Christophe Morisset of UNAM is currently in the process of creating an IDL visualisation tool that is able to deal with multiGrids and reconstruction for a quick and efficient analysis of the results. It is intended to distribute the new tool as soon as it becomes available.
The gas and dust radiative transfer routines in mocassin are now fully integrated. It is therefore now possible to run mocassin in its original gas-only mode, as well as dust-only and of course dust+gas. The basics on how to run the code for dust-only or gas-only cases have already been given in pure photoionisation benchmarks and pure dust benchmarks; here we will concentrate on an example where both dust and gas are present.
Below is the input.in file used for the dust and gas model of NGC 3918 as described by Ercolano, Barlow and Storey (2005, MNRAS, 362, 1038) (http://cdsads.u-strasbg.fr/abs/2005MNRAS.362.1038E).
autoPackets 20. 2. 500000000
densityFile 'ngc3918/ngc3918den.dat'
symmetricXYZ
multiChemistry 2 'ngc3918/ngc3918.dat' 'ngc3918/ngc3918.dat
contShape 'ngc3918/nLTe140lg65'
maxIterateMC 30 95.
nPhotons 1000000
nx 23
ny 23
nz 23
nbins 700
LStar 27.64
nuMax 23.7
nuMin 3.1e-4
Rin 0.
Rout 3.27142e17
TStellar 140000.
MdMh constant 0.0011
dustFile 'ngc3918/primary_grainspecies.dat'
'ngc3918/primary_grainsizes.dat'
writeGrid 50.
convLimit 0.03
resLinesTransfer 90.
In summary, it should be clear from the example above that the only keywords that differentiate the file above from the input of a gas-only model are : MdMh (dust to hydrogen mass ratio) dustFile (files containing the grain size distribution and species information) resLinesTransfer (which tells at what level of grid convergence the resonance line photons escape fractions should be calculated). Note that the first two keywords that define how much, what type and what size grains are to be used is also needed to run dust-only models (although MdMh could be substituted by Ndust or by MdMg). resLinesTransfer is only needed if there is gas also in the simulation (which would then be emitting resonance lines capable of heating the grains).
A number of useful (or not) FORTRAN and IDL programs are included in the accessories/ subdirectory. Some words of warning: the programs are very basic and poorly commented, as they were developed for personal use. Anyone is welcome to use them at their own risk!! The IDL programs, in particular, are specific to the simulation they were developed for and some editing will be necessary to customise them to the user's needs.
Depending on the complexity of the dust model employed in a given model and the geometrical complexity of the grid, it may be quite challenging to explore the dust temperature distribution calculated by mocassin and written out to dustGrid.out.
Sometimes the only way is to plot out the results or create 3D maps using a visualisation program such as IDL, PDL etc..
The accessories/ subroutine contains an example (dustTemperatures2.pro) of how such a grid may visualised using IDL. This was written for a grid containing 2 grain species and 20 grain sizes. The simulation was a multiChemistry dust and gas one, so the results are also split by sector.
The grain temperature spiking routines included in mocassin are based on the Guhathakurta & Draine (1989 ApJ 345, 230) (http://cdsads.u-strasbg.fr/abs/1989ApJ...345..230G) method. This is a very powerful method and allows the time-efficient computation of the time-dependent grain temperatures due to quantum heating. For the limitations of the method also see Siebenmorgen et al. (1992 A&A 266, 501) (http://cdsads.u-strasbg.fr/abs/1992A%26A...266..501S). The temperature spikes only affect the output SED from dust grains, it is therefore advisable not to include these time consuming procedures until the model has almost converged in the case of gas+grain simulations (keep a high value of real2 - see quantumHeatGrain). In the gas of dust only models, it is worth to have the procedures working right from the start since the convergence criterion is then based on dust temperatures and therefore one must take this effect into account right from the start in order to avoid convergence fluctuations. It is in general not worth running the quantum heating routines on large grains that are unlikely to spike. For a discussion of the general cases when quantum heating routines must be considered please see Siebenmorgen et al. (1992 A&A 266, 501) (http://cdsads.u-strasbg.fr/abs/1992A%26A...266..501S).
At present the grain temperature spiking routines are only implemented for carbonaceous or silicate grains. mocassin will expect to be told what type of grain he is dealing with when calculating the spiking. This is done by adding a -capital- 'S' or 'C' at the beginnig of the species label in the optical constants file. e.g. for the Draine & Laor (1993) (http://cdsads.u-strasbg.fr/abs/1993ApJ...402..441L) silicates data (dustData/sil-dlaor.nk):
nk
Ssil_dl 1400. 3.3 0.588 20.077
0.10000E-02 0.99956E+00 0.97380E-04
0.10120E-02 0.99955E+00 0.10160E-03
0.10230E-02 0.99954E+00 0.10610E-03
0.10350E-02 0.99953E+00 0.11060E-03
0.10470E-02 0.99952E+00 0.11520E-03
0.10590E-02 0.99951E+00 0.12000E-03
.........
Please email me (http://mocassin.nebulousresearch.org/contact) if you would like to include T-spiking effects for other species.
mocassin's principal limitation is imposed by the computer power available. The great volume of data which has to be handled in a three-dimensional simulation, implies the need for a system with multi-processing capabilities in order to accelerate the computational time. However, the fast development of Beowulf Linux clusters is making parallel computing more affordable, and this is also a reason why the MPI formalism was chosen (as opposed to openMP), as this allows information to be passed from one processor to another and, therefore, it does not necessarily require a system with shared memory facilities. Such systems, which include the Silicon Graphics Origin 2000 machine used for this work, are generally much more expensive than Beowulf clusters.
Pure-dust models are less computationally expensive than gas or dust and gas ones and reasonably sizes grids can be feasibly run on single processor machines.
mocassin was designed for the modelling of the photoionised region of planetary nebulae and H II regions and it does not, at present, include the high energy physical processes which are needed, for example, for the modelling of AGNs. However the inclusion of processes such as inner shell photoionisation and Compton heating is straightforward and this is intended to be one of the developments of the near future. Moreover, current efforts are geared toward the inclusion of a self-consistent treatment for PDR processes, this requires the incorporation of a chemical network.
No known bugs.
Barbara Ercolano
| 31 Dec 2015 | 2.02.70 |