Radiation module for the code PLUTO

• Introduction
• Features
• Source Code
• Installation
• Create a new project
  • The setup script
  • pluto.ini
    • Boundary Conditions
    • Solver Options
    • Semenov Opacity Options
  • init.c
    • Mean molecular weight
    • Initial radiation energy density
    • Opacity
• Example Setup
  • Plot

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
  • Rosseland mean opacity after Lin & Papaloizou (1985)
  • Rosseland mean opacity after Bell & Lin (1994)
  • Rosseland and Planck mean opacity after Semenov et al. (2003) (implementation contributed by Markus Flaig)
  • Easy interface to implement custom opacities
• Solver
  • SOR (successive overrelaxation)
  • Solver based on the PETSc Library (list of available matrix solvers supported by PETSc)
• Additional physics
  • Irradiation (in spherical coordinates)
  • Artificial dissipation

Source Code

Installation

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 Src/Radiation/assert.h
patching file Src/Radiation/backtrace.c
patching file Src/Radiation/backtrace.h
patching file Src/Radiation/error.h
patching file Src/Radiation/irradiation.c
patching file Src/Radiation/irradiation.h
patching file Src/Radiation/LICENSE
patching file Src/Radiation/makefile
patching file Src/Radiation/modify_pluto.c
patching file Src/Radiation/modify_pluto.h
patching file Src/Radiation/opacity.c
patching file Src/Radiation/opacity.h
patching file Src/Radiation/petsc.c
patching file Src/Radiation/petsc.h
patching file Src/Radiation/petsc_initialize_advanced.c
patching file Src/Radiation/petsc_initialize_advanced.h
patching file Src/Radiation/petsc_tools.c
patching file Src/Radiation/petsc_tools.h
patching file Src/Radiation/radiation.c
patching file Src/Radiation/radiation.h
patching file Src/Radiation/radiation_hd_rhs.c
patching file Src/Radiation/radiation_main.c
patching file Src/Radiation/radiation_mhd_rhs.c
patching file Src/Radiation/radiation_tools.c
patching file Src/Radiation/radiation_tools.h
patching file Src/Radiation/SemenovOpacity/generate_opacity.c
patching file Src/Radiation/SemenovOpacity/generate_opacity.h
patching file Src/Radiation/SemenovOpacity/kP_h2001.dat
patching file Src/Radiation/SemenovOpacity/kR_h2001.dat
patching file Src/Radiation/SemenovOpacity/semenov_opacity.c
patching file Src/Radiation/SemenovOpacity/semenov_opacity.h
patching file Src/Radiation/Templates/init.c
patching file Src/Radiation/Templates/pluto.ini
patching file Tools/Python/make.py
patching file Tools/Python/radiation.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

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.

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

[Radiation-Boundary]

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

[Radiation Option]

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

[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

[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

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

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
    radiation module.

    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"

/**
    The Init() function can be used to assign initial conditions as
    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;

    RadiationSetMu( 0.6 );

    #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
}

#if RADIATION == YES
/**
    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
    attribute Erad of the global struct radiation_data.

    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)
    {
        radiation_data.Erad[k][j][i] = 1.0;
    }
    \endcode
*/
void RadiationSetInitialErad( Grid *grid, Data *data )
{
    RadiationSetInitialEradBlackBody( grid, data );
}

/**
    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
*/
void RadiationOpacityInit()
{

}

/**
    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 );
}

#if IRRADIATION == YES
/**
    Same as \ref RadiationRosselandOpacity and \ref RadiationPlanckOpacity but for the opacity used for irradiation.

    \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.

Initial radiation energy density

void RadiationSetInitialErad(Grid* grid, Data* data)
{
    int k,j,i;
    
    TOT_LOOP(k,j,i)
    {
        radiation_data.Erad[k][j][i] = 1.0;
    }
}

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

void RadiationOpacityInit()
double RadiationRosselandOpacity(double density, double temperature)
double RadiationPlanckOpacity(double density, double temperature)
double IrradiationOpacity(double density, double temperature)

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 $math$ \rho=8.08 \cdot 10^{-4} \; \left[\frac{\text{g}}{\text{cm}^3}\right] $math$, the pressure to $math$ p = \frac{T \rho k_B}{\mu m_u}$math$ with the temperature $math$ T = 450\; [\text{K}] $math$. The mean molecular wight is set to $math$ \mu = 0.6$math$ and the ratio of specific heat is set to $math$\gamma = 1.43$math$.

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 $math$1000\;[\text{K}]$math$ and at the end of the domain to $math$100\;[\text{K}]$math$.

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:

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.

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.