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.
File | Comment |
---|---|
Src/Radiation/assert.h | |
Src/Radiation/backtrace.c | |
Src/Radiation/backtrace.h | |
Src/Radiation/error.h | |
Src/Radiation/irradiation.c | |
Src/Radiation/irradiation.h | |
Src/Radiation/LICENSE | |
Src/Radiation/makefile | |
Src/Radiation/modify_pluto.c | |
Src/Radiation/modify_pluto.h | |
Src/Radiation/opacity.c | |
Src/Radiation/opacity.h | |
Src/Radiation/petsc.c | |
Src/Radiation/petsc.h | |
Src/Radiation/petsc_initialize_advanced.c | |
Src/Radiation/petsc_initialize_advanced.h | |
Src/Radiation/petsc_tools.c | |
Src/Radiation/petsc_tools.h | |
Src/Radiation/radiation.c | |
Src/Radiation/radiation.h | |
Src/Radiation/radiation_hd_rhs.c | Modified version of file Src/HD/rhs.c |
Src/Radiation/radiation_main.c | Modified version of file Src/main.c |
Src/Radiation/radiation_mhd_rhs.c | Modified version of file Src/MHD/rhs.c |
Src/Radiation/radiation_tools.c | |
Src/Radiation/radiation_tools.h | |
Src/Radiation/SemenovOpacity/generate_opacity.c | opacity.f from opacity.zip converted with f2c form Fortran to C |
Src/Radiation/SemenovOpacity/generate_opacity.h | |
Src/Radiation/SemenovOpacity/kP_h2001.dat | Originally from opacity.zip |
Src/Radiation/SemenovOpacity/kR_h2001.dat | Originally from opacity.zip |
Src/Radiation/SemenovOpacity/semenov_opacity.c | |
Src/Radiation/SemenovOpacity/semenov_opacity.h | |
Src/Radiation/Templates/init.c | |
Src/Radiation/Templates/pluto.ini | |
Src/hdf5_io.c | Modified PLUTO file |
Tools/Python/make.py | Modified PLUTO file |
Tools/Python/radiation.py |
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.
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.
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.
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
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
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
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
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.
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].
void RadiationOpacityInit() double RadiationRosselandOpacity(double density, double temperature) double RadiationPlanckOpacity(double density, double temperature) double IrradiationOpacity(double density, double temperature)
double RosselandOpacityLinPapaloizou1985(double density, double temperature)
double RosselandOpacityBellLin1994(double density, double temperature)
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.
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.
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: |