Torus

Preface

TORUS is a three-dimensional radiative transfer code which uses an adaptive mesh refinement scheme and a Monte-Carlo method to solve for the radiative equilibrium, hydrostatic equilibrium, and dust sublimation in circumstellar discs around both low and high-mass pre-main-sequence stars. TORUS is either an acronym for Transport of Radiation Under Sobolev, or Transport of Radiation Using Stokes.

Much of this is covered in detail on the TORUS wiki at Exeter, TorusWeb. To gain access to that resource, you must make an account; if you have SVN access to the TORUS distribution, that username/password should also work for the wiki (I think..). Here, I will discuss TORUS from an example-driven perspective. In this mini-manual, I'm not going to address any molecular capabilities of TORUS, and I am leaving out a lot of TORUS' functionality I don't personally use.

Building TORUS

I've found the building of TORUS to be somewhat system dependent. At Exeter, they test with each version of TORUS how it builds with various combinations of OS and different compilers. Here, I'll outline what works on my machine, a 12 core mac running OS X 10.6.8 (Snow Leopard). I've enabled openmp for parallel computing (which TORUS uses during a Lucy iteration, and possibly when ray-tracing to compute outputs).

Whenever updating TORUS modules, always do:

make clean

make

Environment Variables

Here are the environment variables, etc, set in my .cshrc file regarding TORUS:

setenv OMP_NUM_THREADS 8

setenv TORUS_DATA /Sisu1/aarnio/torus/data

setenv SYSTEM davesmac

limit stacksize 60000

alias visit /Sisu1/aarnio/visit/bin/visit

alias kvis /usr/local/karma/bin/kvis

Makefile

media: alicia_torus_makefile.txt My makefile (cropped down to just my SYSTEM selection). Within the TORUS folder, this is just called Makefile.

Version info

My system parameters:

/usr/local/intel/composerxe-2011.1.122/bin/intel64/ifort

Darwin sisu.astro.lsa.umich.edu 10.8.0 Darwin Kernel Version 10.8.0: Tue Jun 7 16:32:41 PDT 2011; root:xnu-1504.15.3~1/RELEASE_X86_64 x86_64

Auxiliary software

Tim and the Exeter group use VisIt and kvis to display TORUS outputs. VisIt will plot the .vtk output, and kvis plots fits outputs (dust continuum images as well as the velocity data cubes and line profiles therein). I posted a tutorial to the Exeter wiki regarding how to generate specific plots with these tools; pdf here: media:Plottingtorusdata.pdf.

Physics modules in TORUS

Generally, what happens in TORUS is that the disk begins cold. Photons are first absorbed, and then escape (~mm re-radiation, the disk is optically thin to it). The first Lucy iteration goes quickly, and each successive iteration takes longer and longer as the disk becomes optically thick to its own radiation. Usually, convergence (radiative equilibrium only, not hydrostatic) is reached in 5-10 Lucy iterations. The primary criterion for convergence is that the change in total emissivity is less than 1% (with respect to the previous iteration). If you include hydro calculations, a secondary convergence criterion is forced that TORUS run through at least 5 hydro iterations (within each hydro iteration, there are however many radiative equilibrium iterations necessary to reach convergence). These criteria can be relaxed by editing lucy_mod.F90 (not recommended unless you're very comfortable with what you're doing!!! See editing lucy_mod.F90 for an example). Between iterations, if convergence is not met and the fraction of bad cells (<50 photon packets through the cell) is nonzero, the number of source photons doubles. The initial number of photons is ~10x the number of cells, but this can be adjusted with the nlucy parameter.

atomicphysics

This module must be enabled to do line transfer calculations. This module can either perform a comoving frame calculation (computationally expensive) or use the Sobolev approximation.

Example from the literature of implementation of this function of TORUS- Kurosawa, Harries & Symington (2006). In this paper, the authors were distinguishing which parameters would generate Halpha line profiles like those in Bo Reipurth's classification.

Right now, only H is fully implemented (He is in "beta" mode, so to speak). We've lobbied Tim to add Ca, Na, O next.

photoionphysics

Pulling directly from the Exeter wiki:

Domain decomposed photoionization models are run using multiple processors (threads). The computational domain is split into equal sized subdomains, each of which is handled by an individual thread. An additional thread is also required to perform governing operations. At present models can be run using one of the following:

Dimensionality	Nthreads	Description
1D	3	grid split into 2 plus 1 control thread
2D	5	grid split into 4 plus 1 control thread
2D	17	grid split into 16 plus 1 control thread
3D	9	grid split into 8 plus 1 control thread
3D	65	grid split into 64 plus 1 control thread

Types of calculations done in TORUS

Which type of calculation TORUS does depends on the star+disk geometry and physics chosen by the user. Once those are defined, TORUS will do the necessary combinations of the following:

stateq

Statistical equilibrium calculation by solving the rate equations. If atomicphysics T then a comoving frame algorithm is employed.

radeq

This requires that dustphysics T and uses the Monte-Carlo radiative equilibrium method of Lucy (1999)

photoion eq

Solves photoionization equilibrium including thermal equilibrium using a method similar to that of Ercolano xxxx. Dust may be included.

hydrodynamics

hydrodynamics T	! vertical hydrostatic eq. puffs up inner disk
nhydrothreads 17	! number of hydro threads.. ehhh? 16+1 control.. see how this goes.
nhydro 4	! max number of hydro iterations. i think default is 5
splitovermpi T	! for the flux limiting algorithm, this must be set to T

My usage of this module is somewhat limited to just making sure the inner edge of the disk puffs up; for more detailed discussion of the defaults (artificial viscosity, the flux limiting algorithm used, etc), see the Exeter hydrodynamics wiki page.

radiationhydro

This module combines photoionization and hydrodynamics; the evolution of the radiation field and the material dynamics are coupled. TORUS begins with a radiative transfer iteration, then the following iterations alternate between hydrodynamics and radiative transfer. To use this physics, the following modules must be enabled:

radiationhydro T

hydrodynamics T

splitovermpi T

photoionphysics T

Troubleshooting

In the most common application of TORUS, performing a radiative equilibrium calculation, output is generated after each Lucy iteration. The file convergence_lucy.dat documents each Lucy iteration- the iteration number, mean temperature change in the disk, minimum/maximum temperature changes within the grid cells, percentage of bad cells (<50 photon packets through a given cell), dust emissivity/stellar emissivity, total emissivity, maximum fractional change in temperature across all cells, and nMonte, the number of photon packets.

AMR grid

Grid setup

amrgridsize

amrgridcentre[x/y]

amr2d

maxdepthamr

The volume of the smallest grid cell (your finest resolution in the grid) is:

Vol = (grid size / 2^(maxdepthamr))^3 So, for a grid 2000 AU across (1000 AU, radially speaking) and a max cell depth of 20, the smallest cell is 0.0019 AU to a side.

Here is a sample AMR mesh setup:

readgrid F	! we aren't reading a grid, we will set one up from scratch
! inputfile grid_out.dat	! if we did read in a grid, this is how to call it
writegrid T	! write the grid to file (this includes the grid cells, and the EOS in each cell)
outputfile grid_out.dat	! name of the output grid
amrgridsize 2.0e6	! units of 10^10cm (here, I've made grid a little bigger than disk itself- the outermost cells will be huge and empty)
amrgridcentrex 1.0e6	! the linear size of the top-level AMR mesh in units of 10^10 cm. This is useful if you use multiple sources
amr2d T	! this is a 2d (cylindical) model
maxdepthamr 22	! capping the AMR mesh depth helps TORUS to converge faster, saves some CPU. Set this based on how fine a resolution you need in final model.

Special grid options

TORUS can smooth the grid, reducing large cell-to-cell variations in cell refinement and optical depth.

smoothgridtau T	! smooths the grid for optical depth, in order to resolve disc photosphere
dosmoothgrid T	! smooth the grid for jumps in cell refinement
smoothfactor 3.0	! make sure that neighboring cells are not only one AMR depth apart
lambdasmooth 5500.0
taumax 1.
taumin 0.01

TORUS output

Once your TORUS model has reached the criteria for convergence (more on this later), the final grid is written to a temporary output file, and then TORUS will calculate SEDs, images, and/or line profiles by sending however many photons you specify through the converged grid. If you don't calculate the SEDs/images/line profiles while running TORUS, fear not; change the lucy_grid_tmp.dat file to a different name and use it as an input to TORUS, turning off all the other physics modules. In this way, you can 'hot start' TORUS and calculate output data without re-running your model.

Sample output calls

SEDs

nphotons 100000

! the number of photon packets in SED

! Output SEDs

spectrum T

! produce a spectrum

! SED parameters

ninc 2	! number of inclinations
firstinc 1.0	! the first inclination (degrees)
lastinc 48.0	! the last inclination (degrees)
filename MWC275	! the root of the output filename
sised T	! Write spectrum as lambda vs F lambda in SI units
sedlammin 0.12	! minimum wavelength in SED file
sedlammax 2000	! maximum wavelength in SED file
sedwavlin F	! Linear spacing in SED file?
sednumlam 1000	! number of wavelength points in SED

The comments make this fairly self-explanatory, but note: you must set nphotons for an SED or an image. In general, the SEDs need fewer photons than the images to get decent signal to noise (in the SED, divide the total number of photons by the number of wavelength points). I've found 50,000 photons or more for the SEDs works well. Also, if you don't specify sednumlam, the default is 200 wavelength points, and that can produce a jagged, noisy SED.

Images

nphotons 10000000	! the number of photon packets in image
image T	! produce images
nimage 4	! how many images?
imageaxisunits AU
imagesize 20	! Size of your image across each side in AU. divide this by your npixels to get desired AU/pixel resolution.
imagefile1 kband_20AU.fits	! name of first image file
lambdaimage1 21590.	! monochromatic wavelength in Angstroms
npixels1 256	! number of pixels. I generally don't adjust this (the more pixels, the more photons you need)
inclination1 0	! inclination of the system; if unknown, try a few different values in multiple images
imagetype1 dustonly	! choose freefree, forbidden, recombination, or dustonly
imagefile2 nband_1_20AU.fits
lambdaimage2 77000	! 7.7um
npixels2 256
inclination2 0
imagetype2 dustonly
imagefile3 nband_2_20AU.fits
lambdaimage3 99000	! 9.9um
npixels3 256
inclination3 0
imagetype3 dustonly
imagefile4 nband_3_20AU.fits
lambdaimage4 126000	! 12.6um
npixels4 256
inclination4 0
imagetype4 dustonly

Line profiles

To calculate line profiles, you must run the TORUS model with the comoving frame option enabled.

cmf T

! comoving frame (vs sobolev approx)</nowiki>

also, specify the atomic physics option and its parameters:

atomicphysics T	! Include atomic physics
natom 1	! One model atom
atom1 H.atm	! Hydrogen
xabundance 1.0	! Pure hydrogen
yabundance 0.	! no helium
vturb 20.	! microturbulence in km/s

TORUS will output a velocity space data cube per your input parameters:

! output a datacube

datacube T	! produce a fits datacube
inclination 1.	! viewing angle
positionangle 0.	! position angle
datacubefile cmf_i1_ha.fits	! title of fits output
imageside 1500.	! size of image in 10^10cm
npixels 200	! number of pixels
nv 200	! number of velocity bins
maxVel 800.d0	! -800 to +800 km/s
distance 140.	! distance to object in pc
lamline 6563.	! wavelength in Angstroms</nowiki>

The nitty-gritty: input parameters

In general, TORUS is pretty forgiving in terms of parameters- if any two inputs are redundant, or if you forget something, TORUS will halt and give an intelligible error message. If you include extra parameters TORUS doesn't need, there's no catastrophic fail- just a message saying those parameters will be ignored.

Source parameters

Awesomely, TORUS can handle multiple input sources; to add additional sources, simply change the numerical value after each parameter, and set nsource to whatever value greater than 1 applies.

nsource 1	! there is just one source
radius1 2.0	! it has a radius of 1 solar radius
teff1 10000.	! the source effective temperature
contflux1 kurucz	! the continuum flux (other option is blackbody)
mass1 2.5	! the source has a mass of one solar mass
sourcepos1 0. 0. 0	! it is located at the grid cntre
distance 150.	! Distance to observer, pc

Disk geometries

TORUS can calculate very simple disk configurations, or as complicated a disk (plus magnetosphere, in the ttauri case) as you want. If you talk to Ilse, she can tell you about her own disk geometry modules. Since TORUS is written in a modular way, creating whatever geometry you'd like is more straightforward than one might think.

benchmark

A Pascucci benchmark disk

shakara

This is a sample pulled from an MWC 275 parameter file.

geometry shakara	! flared protostellar disk
rinner 0.22	! inner disc radius (AU) --this is from ajay's paper
router 200.	! outer disc radius (AU)
height 10.	! disc scaleheight at 100 AU (in AU)
mdisc 0.01	! Msun
alphadisc 1.0	! this is from ajay's paper
betadisc 1.125	! disc scaleheight goes as r^beta</nowiki>

Options that apply to any disk geometry

smoothinneredge T	! exponential density decay at inner disk edge
gasopacity T
vardustsub T

The gasopacity switch is, I believe, the only way that gas is applied in TORUS: the gas is a source of opacity, and it is a radiation source (probably mostly emission), but it isn't coupled to the gas (i.e., doesn't play a role in the hydrostatic equilibrium of the disk).

'Important!' If you set variable dust sublimation, vertical hydrostatic equilibrium must be on (hydro T) to allow the inner disk to puff up!! This will add the secondary convergence criterion that convergence can only be reached if the number of Lucy iterations is greater than 9- this is an excerpt from lucy_mod.F90:

if (variableDustSublimation) then

if (iIter_grand < 9) then

converged = .false.

endif

endif

As far as I know, you can only change this within the lucy_mod.F90 file (don't forget to re-build TORUS after editing that file!). You can limit the number of hydro iterations on top of these Lucy iterations by setting the nhydro keyword. The default is 5 hydro iterations.

Dust options

These dust parameters apply when dustphysics T.

iso_scatter T	! Assume isotropic scattering (assumed by benchmark)
ndusttype 1
graintype1 sil_dl	! Drain and Lee silicates
amin1 0.01	! minimum grain size (microns)
amax1 0.25	! maximum grain size (microns)
qdist1 1.5	! power law index (a^-qdist)
dusttogas 0.01	! torus assumes this value, even if you don't explicitly define it.