# Radiation module for the code PLUTO

## Introduction

This module implements the radiation transport equation as described in the publication from Commerçon at al. In detail we solve this equations in the flux limited diffusion approximation within the co-moving frame. To be able to solve even large problems including radiation we use the gray approximation instead of solving the radiation transport equation for a set of frequency bins. Scattering is not accounted directly in this treatment instead it is included in the effective isotropic absorption and emission coefficients.

The radiation module was implemented by Stefan Kolb as a part of his diploma thesis at the Eberhard Karls Universität Tübingen.

## Features

• Implemented fully parallel with MPI
• Tested up to 1024 processor cores but should work efficient with much more
• Usable with the PLUTO modules
• HD
• MHD
• Supported cordinate systems in 3 dimensions
• Cartesian
• Polar
• Spherical
• Opacity
• Solver
• Artificial dissipation

## Source Code

The source code for the radiation module can be downloaded as tar.gz and zip file via the following link radiation_module-1.1.tar.gz [zip]. The archive contains the files listed below. To use the source it must be integrated into PLUTO which is described in the following section.
File Comment
Src/Radiation/SemenovOpacity/generate_opacity.c opacity.f from opacity.zip converted with f2c form Fortran to C
Src/hdf5_io.c Modified PLUTO file
Tools/Python/make.py Modified PLUTO file

## Installation

This installation instructions are for the operating system Unix/Linux, but in general the instructions should be similar on Mac OS and Windows.

The first step is to download PLUTO version 4.0 (pluto-4.0.tar.gz) form http://plutocode.ph.unito.it/Download.html (free registration is needed) and save the file to a directory of your choice. We will use here the directory /home/user/programs/. To extract the source of PLUTO enter the following commands in a terminal and don't forget to change "/home/user/programs/" if you use a different directory.

$cd /home/user/programs/$ tar -xvzf pluto-4.0.tar.gz


Detailed instructions for the installation and initial setup of PLUTO can be found in the user guide shipped with PLUTO. The newest version of the user guide can also be found online: userguide.pdf

After extracting pluto-4.0.tar.gz the directory PLUTO is available. Now download the patch radiation_module-1.1.patch into the current directory here /home/user/programs and then execute the following commands to patch PLUTO to contain all files needed to use the radiation module.

$cd PLUTO$ patch -p1 -i ../radiation_module-1.1.patch


The output of these commands should look like:

patching file Src/hdf5_io.c
patching file Tools/Python/make.py


PLUTO 4.0 is now patched and contains the directory $PLUTO_DIR/Src/Radiation, the file$PLUTO_DIR/Tools/Python/radiation.py and a modified version of the file $PLUTO_DIR/Tools/Python/make.py as well as$PLUTO_DIR/Src/hdf5_io.c. $PLUTO_DIR is a shell variable used to specify the location of PLUTO's source code and is here /home/user/programs/PLUTO. ## Create a new project If you are unfamiliar with PLUTO please read the PLUTO user manual and setup your first simulation without radiation transport. There are some example simulations provided with PLUTO where you can verify that your installation of PLUTO is working correctly. ### The setup script A new project (problem setup) with enabled radiation can be creates by using the ncurses based setup script provided by PLUTO. If PLUTO is installed correctly and the shell variable$PLUTO_DIR is set as described in the user guide, the setup can be started with the command:

$python$PLUTO_DIR/setup.py

If this command is executed in an empty directory the setup creates all files necessary to use PLUTO with the radiation module. After the setup is started the radiation module is found in the HD and the MHD menu depending on which of the modules is used. The menu entry to enable the radiation transport is called RADIATION and if it is set to YES the next menu will show different options for the radiation module.

The PETSc library is not provided with PLUTO or the radiation module. If you want to use the PETSc solver you have to install and configure the PETSc library on your system.

All steps for enabling the radiation module are shown in the following pictures. In this example we have chosen the solver PETSc. The steps are slightly different for the SOR solver.

In the case that the SOR solver is selected (RADIATION_SOLVER set to SOLVER_SOR) an additional menu for the configuration of the SOR solver is shown after the radiation menu. This is shown in the following pictures.

### pluto.ini

If the setup script was run successfully a couple of file are created. One of these is the file pluto.ini with the following content:

[Grid]

X1-grid    1    0.0    100    u    1.0
X2-grid    1    0.0    100    u    1.0
X3-grid    1    0.0    1      u    1.0

[Chombo Refinement]

Levels           4
Ref_ratio        2 2 2 2 2
Regrid_interval  2 2 2 2
Refine_thresh    0.3
Tag_buffer_size  3
Block_factor     8
Max_grid_size    64
Fill_ratio       0.75

[Time]

CFL              0.4
CFL_max_var      1.1
tstop            1.0
first_dt         1.e-4

[Solver]

Solver         tvdlf

[Boundary]

X1-beg        outflow
X1-end        outflow
X2-beg        outflow
X2-end        outflow
X3-beg        outflow
X3-end        outflow

rX1-beg        reflective
rX1-end        reflective
rX2-beg        fixedvalueT 5.0
rX2-end        symmetric
rX3-beg        periodic
rX3-end        periodic

absolute-tolerance  1e-80
relative-tolerance  1e-8
max-iterations      10000
fpe-signal-haldler  disable
#petsc-options      -help -ksp_type=cgne -pc_type=icc

semenov-opacity-nT       2048
semenov-opacity-T0       1.000e+01
semenov-opacity-T1       9.999e+03
semenov-opacity-nrho     2048
semenov-opacity-rho0     1.000e-16
semenov-opacity-rho1     9.999e-09
semenov-opacity-top      h
semenov-opacity-model    nrm
semenov-opacity-shape    s

[Static Grid Output]

uservar    0
dbl        1.0  -1   single_file
flt       -1.0  -1   single_file
vtk       -1.0  -1   single_file
tab       -1.0  -1
ppm       -1.0  -1
png       -1.0  -1
log        1
analysis  -1.0  -1

[Chombo HDF5 output]

Checkpoint_interval  -1.0  0
Plot_interval         1.0  0

[Parameters]

SCRH   0


The options in this file can be places at an arbitrary position but for reasons of clarity the options are grouped. The options newly introduced by the radiation module are described in the following sections, all other options are described in detail by the user guide of PLUTO.

#### Boundary Conditions

This section specifies the boundary conditions in the different directions used by the radiation module. An example how this section can look like is show below:

[Radiation Boundary]

rX1-beg        reflective
rX1-end        reflective
rX2-beg        fixedvalueT 5.0
rX2-end        symmetric
rX3-beg        periodic
rX3-end        periodic

There are 6 different boundary conditions available for the radiation module:
• periodic
• symmetric
• reflecting
• outflow
• fixedvalue, fixedvalueT
• shearingbox (contributed by Moritz Stoll)
Note the symmetric, reflective and outflow boundary conditions are equal. They would be different for vectorial quantities like the magnetic filed but for a scalar quantity as the radiation energy density they are equal. The fixedvalue boundary conditions sets the radiation energy density on the whole specified boundary to the given value. The fixedvalueT boundary condition does the same but with the difference that a temperature must be specified. The given temperature is then converted to the radiation energy density by the equation $E=a_R T^4$ with $a_R=\frac{8 \pi^5 k_B^4}{15 c^3 h ^3}$.

#### Solver Options

General solver unspecific options for the absolute and relative tolerance as well as the maximum number of used iterations. Note it is possible that a set of values work with one solver but not with the other.

[Radiation Solver Options]

absolute-tolerance   1e-40
relative-tolerance   1e-5
max-iterations       10000


The following option enables or disables the handling of float point exceptions.

fpe-signal-haldler   enable

As the name suggest the option petsc-options is specific for the solver PETSc. All arguments after petsc-options are passed to the PETSc library. The following example sets the matrix solver to cgne and the preconditioner to icc. Additionally the option -help is passed to PETSc which causes PETSc to print a help message with all available options and the current used values for the options.

petsc-options        -help -ksp_type=cgne -pc_type=icc

#### Semenov Opacity Options

If the Semenov opacity is used the opacity table is generated at the start of the simulation and the following options are needed. For more information about the options visit the site http://www.mpia.de/homes/henning/Dust_opacities/Opacities/opacities.html or have look at the comments in the file opacity.f in the zip archive opacity.zip.

[Radiation Semenov Opacity]

semenov-opacity-nT       2048
semenov-opacity-T0       1.000e+01
semenov-opacity-T1       9.999e+03
semenov-opacity-nrho     2048
semenov-opacity-rho0     1.000e-16
semenov-opacity-rho1     9.999e-09
semenov-opacity-top      h
semenov-opacity-model    nrm
semenov-opacity-shape    s


### init.c

The radiation module must be compiled in a parallel environment. This means an implementation of MPI like OpenMPI must be installed and configured. To compile PLUTO it is necessary to use a makefile like Linux.mpicc.defs, where PARALLEL is set to YES.

Another file created by the setup script is the file init.c. The default content, in the case that the radiation module is enabled, is the following:

/**
\file
\brief Contains basic functions for problem initialization.

The file was modified to contain all needed function to use the

The init.c file collects most of the user-supplied functions useful
for problem configuration.
It is automatically searched for by the makefile.

\author A. Mignone (mignone@ph.unito.it)
\author S. Kolb

\date   March 14, 2013
*/

#include "pluto.h"

/**
as a function of spatial position.

\param [out] v   a pointer to a vector of primitive variables
\param [in] x1   coordinate point in the 1st dimension
\param [in] x2   coordinate point in the 2nd dimension
\param [in] x3   coordinate point in the 3rdt dimension

The meaning of x1, x2 and x3 depends on the geometry:
\f[ \begin{array}{cccl}
x_1  & x_2    & x_3  & \mathrm{Geometry}    \\ \noalign{\medskip}
\hline
x    &   y    &  z   & \mathrm{Cartesian}   \\ \noalign{\medskip}
R    &   z    &  -   & \mathrm{cylindrical} \\ \noalign{\medskip}
R    & \phi   &  z   & \mathrm{polar}       \\ \noalign{\medskip}
r    & \theta & \phi & \mathrm{spherical}
\end{array}
\f]

Variable names are accessed by means of an index v[nv], where
nv = RHO is density, nv = PRS is pressure, nv = (VX1, VX2, VX3) are
the three components of velocity, and so forth.

*/
void Init( double *v, double x1, double x2, double x3 )
{
v[RHO] = 1.0;
v[VX1] = 0.0;
v[VX2] = 0.0;
v[VX3] = 0.0;
v[PRS] = 1.0;
v[TRC] = 0.0;

#if PHYSICS == MHD || PHYSICS == RMHD
v[BX1] = 0.0;
v[BX2] = 0.0;
v[BX3] = 0.0;

v[AX1] = 0.0;
v[AX2] = 0.0;
v[AX3] = 0.0;
#endif
}

/**
This function is called at the initialization of the radiation module to set the initial values of the
radiation energy density \f$E \f$. The default is to call the function \ref RadiationSetInitialEradBlackBody
which sets \f$E \f$ to \f$E=a_R T^4\f$.

To customize this function iterate over all array element of the local domain (ghost cells included). And
set \f$E \f$ to the value you want. To access the radiation energy density data you have to use the

To set all cells in \f$E \f$ to 1 in code units you can write:

\code{.c}
int k,j,i;

TOT_LOOP(k,j,i)
{
}
\endcode
*/
{
}

/**
This function can be used to initialize the opacity as needed for example to use the Semenov opacities.

To use the Semenov opacities call the function \ref SemenovOpacityInit() here
*/
{

}

/**
This function is called every time the Rosseland opacity is needed in the radiation module.
\Note It is essential implement this function.

There are different available routines which can be used here.

Rosseland mean opacity after paper Lin & Papaloizou (1985): \ref RosselandOpacityLinPapaloizou1985(density,temperature)

Rosseland mean opacity after paper Bell & Lin (1994): \ref RosselandOpacityBellLin1994(density,temperature)

Rosseland mean opacity after paper Semenov et al. (2003): \ref SemenovRosselandOpacity(density,temperature)
\Note It is necessary to call \ref SemenovOpacityInit() in \ref RadiationOpacityInit() to use the Semenov opacities.

\param[in] density      density in cgs units
\param[in] temperature  temperature in cgs units
\return the Rosseland mean opacity in cgs units
*/
double RadiationRosselandOpacity( double density, double temperature )
{
return RosselandOpacityLinPapaloizou1985( density, temperature );
}

/**
Same as \ref RadiationRosselandOpacity but for the Planck opacity.

\Note Except for the Semenov opacities there is currently no Planck mean opacity routine
implemented. A solution to this is to use here for simplicity also the Rosseland mean opacity.

If you want to use the Semenov opacities use the function \ref SemenovPlanckOpacity(density,temperature)

\param[in] density      density in cgs units
\param[in] temperature  temperature in cgs units
\return the Planck mean opacity in cgs units
*/
double RadiationPlanckOpacity( double density, double temperature )
{
return RosselandOpacityLinPapaloizou1985( density, temperature );
}

/**

\param[in] density      density in cgs units
\param[in] temperature  temperature in cgs units
\return opacity in cgs units
*/
double IrradiationOpacity( double density, double temperature )
{
return RosselandOpacityLinPapaloizou1985( density, temperature );
}
#endif

#endif

/**
Perform runtime data analysis.

\param [in] d the PLUTO Data structure
\param [in] grid   pointer to array of Grid structures
*/
void Analysis( const Data *d, Grid *grid )
{

}

#if PHYSICS == MHD
/**
Define the component of a static, curl-free background
magnetic field.

\param [in] x1  position in the 1st coordinate direction \f$x_1\f$
\param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
\param [in] x3  position in the 3rd coordinate direction \f$x_3\f$
\param [out] B0 array containing the vector components of the background magnetic field
*/
void BackgroundField( double x1, double x2, double x3, double *B0 )
{
B0[0] = 0.0;
B0[1] = 0.0;
B0[2] = 0.0;
}
#endif

/**
Assign user-defined boundary conditions.

\param [in,out] d  pointer to the PLUTO data structure containing
cell-centered primitive quantities (d->Vc) and staggered
magnetic fields (d->Vs, when used) to be filled.
\param [in] box    pointer to a RBox structure containing the lower
and upper indices of the ghost zone-centers/nodes or edges at
which data values should be assigned.
\param [in] side   specifies the boundary side where ghost zones need
to be filled. It can assume the following pre-definite values:
X1_BEG, X1_END, X2_BEG, X2_END, X3_BEG, X3_END. The special value
side == 0 is used to control a region inside the computational domain.
\param [in] grid  pointer to an array of Grid structures.
*/
void UserDefBoundary( const Data *d, RBox *box, int side, Grid *grid )
{
int   i, j, k, nv;
double  *x1, *x2, *x3;

x1 = grid[IDIR].x;
x2 = grid[JDIR].x;
x3 = grid[KDIR].x;

if( side == 0 )     /* -- check solution inside domain -- */
{
DOM_LOOP( k, j, i ) {};
}

if( side == X1_BEG )  /* -- X1_BEG boundary -- */
{
if( box->vpos == CENTER )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X1FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X2FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X3FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
}

if( side == X1_END )  /* -- X1_END boundary -- */
{
if( box->vpos == CENTER )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X1FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X2FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X3FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
}

if( side == X2_BEG )  /* -- X2_BEG boundary -- */
{
if( box->vpos == CENTER )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X1FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X2FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X3FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
}

if( side == X2_END )  /* -- X2_END boundary -- */
{
if( box->vpos == CENTER )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X1FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X2FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X3FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
}

if( side == X3_BEG )  /* -- X3_BEG boundary -- */
{
if( box->vpos == CENTER )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X1FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X2FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X3FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
}

if( side == X3_END )  /* -- X3_END boundary -- */
{
if( box->vpos == CENTER )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X1FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X2FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
else if( box->vpos == X3FACE )
{
BOX_LOOP( box, k, j, i ) {  }
}
}
}

#if BODY_FORCE != NO
/**
* Prescribe the acceleration vector as a function of the coordinates
* and the vector of primitive variables *v.

\param [in] v  pointer to a cell-centered vector of primitive variables
\param [out] g acceleration vector
\param [in] x1  position in the 1st coordinate direction \f$x_1\f$
\param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
\param [in] x3  position in the 3rd coordinate direction \f$x_3\f$

*/
void BodyForceVector( double *v, double *g, double x1, double x2, double x3 )
{
g[IDIR] = 0.0;
g[JDIR] = 0.0;
g[KDIR] = 0.0;
}

/**
Return the gravitational potential as function of the coordinates.

\param [in] x1  position in the 1st coordinate direction \f$x_1\f$
\param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
\param [in] x3  position in the 3rd coordinate direction \f$x_3\f$

\return The body force potential \f$\Phi(x_1,x_2,x_3) \f$.
*/
double BodyForcePotential( double x1, double x2, double x3 )
{
return 0.0;
}
#endif


#### Mean molecular weight

To set the value used for the mean molecular wight μ to an user defined value you have to use the function RadiationSetMu. This function must be called from the Init function. And it doesn't matter that RadiationSetMu is called more than once. But note, the value set by the last call is used for the computation.

To use the radiation module it is necessary to implement the function
as done in the file init.c shown above. The default is to initialize the radiation energy density to the energy density emitted by a black body at the temperature T of the gas $E_{i,j,k} = a_R T^4$. This is done for each grid cell in the domain individually by calling the function
To set the initial radiation energy density for example to a constant value of 1.0 in code units. The following implementation can be used.
void RadiationSetInitialErad(Grid* grid, Data* data)
{
int k,j,i;

TOT_LOOP(k,j,i)
{
}
}


In this example we access the radiation energy density through the global struct radiation_data defined in radiation.h with its attribute Erad. Erad is a 3D array like for example data->Vc[RHO].

#### Opacity

Additional to the function RadiationSetInitialErad it is necessary to implement the following four functions for setting up and defining the opacities which are used in the computation.
void RadiationOpacityInit()

The function RadiationOpacityInit can be used to initialize the opacity. It is called at the initialization of the radiation module. And it is needed if the Semenov opacities are used. In this case it is necessary to call the function SemenovOpacityInit() inside RadiationOpacityInit().

The other three functions are used to define the opacities. The function RadiationRosselandOpacity is called every time the Rosseland opacity is needed in the computation. It expects as the other functions the density and the temperature in cgs units and the functions must return the opacity in cgs units. The function RadiationPlanckOpacity is similar to RadiationRosselandOpacity but it is called every time the Planck opacity is needed in the computation.

The function IrradiationOpacity is only needed if irradiation is used. It is called every time the opacity is used for the computation of the optical depth.
There are different available routines which can be used in the described functions.
• Rosseland mean opacity after paper Lin & Papaloizou (1985):
double RosselandOpacityLinPapaloizou1985(double density, double temperature)
• Rosseland mean opacity after paper Bell & Lin (1994):
double RosselandOpacityBellLin1994(double density, double temperature)
• Rosseland and Planck mean opacity after paper Semenov et al. (2003):
double SemenovRosselandOpacity(double density, double temperature)
double SemenovPlanckOpacity(double density, double temperature)

Note: It is necessary to call the function SemenovOpacityInit() in RadiationOpacityInit() to use the Semenov opacities.
In the template for init.c shown above all opacity function call the function RosselandOpacityLinPapaloizou1985 to compute the opacity.

## Example Setup

This example is without hydrodynamics and very simple. The domain is quasi one dimensional with a length of 300 cm and a width and height of 3 cm. The grid was chose to have 300x3x3 grid cells and the simulation was performed in Cartesian coordinates. The initial density was set to $\rho=8.08 \cdot 10^{-4} \; \left[\frac{\text{g}}{\text{cm}^3}\right]$, the pressure to $p = \frac{T \rho k_B}{\mu m_u}$ with the temperature $T = 450\; [\text{K}]$. The mean molecular wight is set to $\mu = 0.6$ and the ratio of specific heat is set to $\gamma = 1.43$.

The boundary conditions are also very simple. All except of the boundaries in z direction are set to periodic. Note even without solving the equations of hydrodynamics it is important to set the boundary conditions in the periodic case for both the hydrodynamics and the radiation module to periodic.
The beginning boundary condition in z direction is set to fixedvalueT with a value of $1000\;[\text{K}]$ and at the end of the domain to $100\;[\text{K}]$.

To disable the hydrodynamics solver add the following line

CFLAGS += -DDISABLE_HYDRODYNAMICS=1

to the file local_make. This will add -DDISABLE_HYDRODYNAMICS=1 to each call of the compiler.

All files for this setup can be found here:

init.c
local_make
definitions.h
pluto.ini

To generate the makefile it is necessary to call the setup script as described above. Then enter the entry "Change makefile" and select the makefile template which fits your needs. Now you can leave the setup script and run the command make.

If you want to run this setup with more as one processor you have to use the following command to disable the parallelization in the x and y direction. This is necessary because of the very view grid cells in the dirextions which we used in this example. To run the simulation for example on 8 cores use a command like:
mpirun -n 8 ./pluto -no-x1par -no-x2par

### Plot

Here we show the results of a simulation which was run up to 2000 seconds. With the slider below the plot it is possible to select the time at which the graph is shown and it is also possible to chose the plotted variable with the shown buttons.

 Select Variable: