diff --git a/docs/physics/est_and_conv/convergence.rst b/docs/physics/est_and_conv/convergence.rst new file mode 100644 index 00000000000..b11a98dfb85 --- /dev/null +++ b/docs/physics/est_and_conv/convergence.rst @@ -0,0 +1,7 @@ +.. _convergence: + +*********** +Convergence +*********** + +Coming soon. \ No newline at end of file diff --git a/docs/physics/est_and_conv/index.rst b/docs/physics/est_and_conv/index.rst index ad323ee1169..a17e2cce827 100644 --- a/docs/physics/est_and_conv/index.rst +++ b/docs/physics/est_and_conv/index.rst @@ -4,6 +4,29 @@ Estimators and Convergence ************************** +As light travels through a real plasma, it has effects on the properties of the plasma due to light-matter +interactions as well as the presence of extra energy from the light. Additionally, as :ref:`previously discussed +`, properties of the plasma affect how light travels through it. Things like this (where two things +both affect each other on different ways) frequently occur in physics. The solution is finding a steady-state for +the plasma properties; that is, the actual plasma will be in a state such that the plasma state will not change as +light propagates through it, because the effects of the light on the plasma and the effects of the plasma on the +light have been "balenced out." + +One of TARDIS's main goals is to determine this plasma state (as we need the actual plasma properties in order to +get an accurate spectrum). This is done in an iterative process. After each :ref:`Monte Carlo iteration +` (which sends light through the supernova ejecta), TARDIS uses objects called estimators to determine +how the propagating light affects the plasma state, after which the plasma state is updated (as will be demonstrated +in the eatimators page linked below). We do this many times, and attempt to have the plasma state converge +to the steady-state we are looking for. In fact, all but the last Monte Carlo iteration is used for this purpose +(after which TARDIS will have the needed plasma state for its last iteration which calculates the spectrum). + +.. note:: + For all but the last iteration, TARDIS uses the number of Monte Carlo packets specified in the + :ref:`configuration file ` under the ``no_of_packets`` argument. This is because + a different number of packets may be necessary to calculate the spectrum as opposed to calculate the + plasma state. + .. toctree:: - estimators \ No newline at end of file + estimators + convergence \ No newline at end of file diff --git a/docs/physics/montecarlo/images/expansion_animation.gif b/docs/physics/images/expansion_animation.gif similarity index 100% rename from docs/physics/montecarlo/images/expansion_animation.gif rename to docs/physics/images/expansion_animation.gif diff --git a/docs/physics/montecarlo/images/optical_depth_summation.png b/docs/physics/images/optical_depth_summation.png similarity index 100% rename from docs/physics/montecarlo/images/optical_depth_summation.png rename to docs/physics/images/optical_depth_summation.png diff --git a/docs/physics/montecarlo/images/scatter_downbranch_ma.png b/docs/physics/images/scatter_downbranch_ma.png similarity index 100% rename from docs/physics/montecarlo/images/scatter_downbranch_ma.png rename to docs/physics/images/scatter_downbranch_ma.png diff --git a/docs/physics/montecarlo/images/spherical_symmetry.png b/docs/physics/images/spherical_symmetry.png similarity index 100% rename from docs/physics/montecarlo/images/spherical_symmetry.png rename to docs/physics/images/spherical_symmetry.png diff --git a/docs/physics/montecarlo/basicprinciples.rst b/docs/physics/montecarlo/basicprinciples.rst index 8c18c790376..adee47f05bd 100644 --- a/docs/physics/montecarlo/basicprinciples.rst +++ b/docs/physics/montecarlo/basicprinciples.rst @@ -4,25 +4,126 @@ Monte Carlo Radiative Transfer - Basic Principles ************************************************* -Radiative transfer describes how the properties of the electromagnetic -radiation field change as it propagates through an ambient material. Due to the -complex coupling induced by the many interactions the photons which constitute -the radiation field may perform with the surrounding material, solving the -radiative transfer problem is in most astrophysical cases only possible by -means of numerical calculations. While there are many different numerical -techniques available to tackle this, Monte Carlo techniques have become a -successful and elegant tool particularly for radiative transfer problems in -supernovae. - -Monte Carlo Radiative Transfer methods are probabilistic techniques and draw -inspiration from the microscopical interpretation of how photons propagate -through and interact with the ambient material. A sufficiently large number of -representative "machine photons" are considered and their propagation history -solved in a stochastic process. The initial properties of these photons are -randomly (in a probabilistic sense) assigned in accordance with the macroscopic -properties of the radiation field (see :ref:`Energy Packets `) -and in a similar manner the decisions about when, where and how the machine -photons interact with the surrounding material are made (see :ref:`Propagation -`). If this process is repeated for a large enough number of machine -photons, the ensemble behaviour and thus the macroscopic evolution of the -radiation field is recovered (see :ref:`Estimators `). +Radiative transfer describes how electromagnetic radiation (light) propagates through a medium. Since there are +a large number of light-matter interactions that have the possibility of affecting the light propagation, solving +the radiative transfer problem is in most astrophysical cases only possible by means of numerical calculations. +While there are many different numerical techniques available to tackle this, Monte Carlo techniques have become a +successful and elegant tool particularly for radiative transfer problems in supernovae. + +Monte Carlo Radiative Transfer methods track a sufficiently large number of photons (light particles) as they +propagate through the supernova ejecta. The initial properties of these photons are randomly (in a probabilistic +sense) assigned in accordance with the macroscopic properties of the radiation field (see :ref:`initialization`) +and in a similar manner the decisions about when, where and how the photons interact with the surrounding material +are made (see :ref:`Propagation `). Given a large enough sample, these photons behave as a microcosom +of all of the light traveling through the ejecta -- that is, based on the behavior of these photons, we can draw +conclusions about the propagation of light through the ejecta as a whole. This is eventually used to determine the +actual steady-state plasma properties (see :ref:`est_and_conv`) and the emitted spectrum (see :ref:`spectrum`). + + +.. _randomsampling: + +Random Sampling Basics +====================== + +During both the initialization of these photons and their propagation through the ejecta are modeled through +probabilistic processes. This involves assigning probabilities to the occurrence of certain events or properties. +For example, during isotropic scattering, finding a photon scattering into any given direction is equally likely. +During the Monte Carlo simulation, assignments +according to these probabilities have to be frequently performed. For this purpose, so-called Random +Number Generators are available. These produce (pseudo-) random numbers +:math:`z` uniformly distributed on the interval :math:`[0,1]`. The challenge +now lies in using these numbers to sample any physical process involved in the +Radiative transfer calculation. From a probability theory point of view, this +just implies finding a mapping between the probability distribution governing the +physical process and the one underlying the Random Number Generator. This +process is typically referred to as random sampling. + +Inverse transformation method +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. note:: + This is a very superficial and sloppy description of Random sampling. More + detailed and rigorous accounts are found in the standard literature, for + example in :cite:`Kalos2008`. + +The simplest and most-used technique in Monte Carlo radiative transfer +applications is referred to as the inverse transformation method and involves +the cumulative distribution function. In general, a random process captured by +the random variable :math:`X` is governed by a probability density +:math:`\rho_X(x)` (the continuous counterpart to discrete probabilities), with +:math:`\rho_X(x) \mathrm{d}x` describing the probability of the variable taking +values in the interval :math:`[x, x+\mathrm{d}x]`. The cumulative distribution +function in turn describes, as the name suggests, the probability of the +variable taking any value between :math:`-\infty` and :math:`x`: + +.. math:: + + f_X(x) = \int_{-\infty}^x \mathrm{d}x \rho_X(x) + +Since the probability density is by definition always positive, the cumulative +distribution function is monotonically increasing. This constitutes the basis +for the inverse transformation function. Consider two random variables, +:math:`X` and :math:`Y`. A mapping between those may be established by equating +their cumulative distribution functions. Numbers :math:`y` distributed +according to one of the probability densities (in this case :math:`\rho_Y(y)`) +may then be used to sample the other process by + +.. math:: + + x = f_X^{-1}\left(f_Y(y)\right). + +For the Random Number Generators described above, the cumulative distribution +function is trivial, namely :math:`f_Z(z) = z`. However, the inverse +distribution sampling method relies on finding the analytic inverse function of +the cumulative distribution function governing the physical processes to be +sampled. If this is not possible, other sampling methods, such as von-Neumann +rejection sampling techniques, have to be used. + +Examples +^^^^^^^^ + +A few examples are provided to illustrate the random sampling process using the +inverse transformation technique. + +Isotropic Scattering +-------------------- + +Consider the case of an isotropic scattering. +Here, all scattering angles are equally likely. Thus, the corresponding +(normalized) probability density and the cumulative distribution function are given by + +.. math:: + + \rho_{\mu}(\mu) &= \frac{1}{2}\\ + f_{\mu}(\mu) &= \frac{1}{2} (\mu - 1). + +This leads to the sampling rule + +.. math:: + + \mu = 2 z - 1. + +Next Interaction Location +------------------------- + +The probability of a photon interacting after covering an optical depth +:math:`\tau` is given by (see :ref:`propagation`) + +.. math:: + + \rho_{\tau}(\tau) &= \exp(-\tau)\\ + f_{\tau}(\tau) &= 1 - \exp(-\tau). + + +With the inverse transformation method, the optical depth to the next interaction location may then be sampled by + +.. math:: + + \tau = - \mathrm{ln}(1 - z) + + +which is equivalent to + +.. math:: + + \tau = - \mathrm{ln}z. \ No newline at end of file diff --git a/docs/physics/montecarlo/index.rst b/docs/physics/montecarlo/index.rst index 30ffbeb7fb2..d4c69971330 100644 --- a/docs/physics/montecarlo/index.rst +++ b/docs/physics/montecarlo/index.rst @@ -4,6 +4,12 @@ Monte Carlo Iteration ********************* +After setting up the simulation, TARDIS runs the simulation using the ``.run()`` method. This runs several Monte +Carlo iterations (which will be described in the links below), corresponding to the number of iterations specified +in the :ref:`Monte Carlo Configuration `. As will be decribed in :ref:`est_and_conv` and +:ref:`spectrum`, most of these iterations will eventually be used to calculate the steady-state plasma properties, +with the last iteration being used to determine the spectrum. + The following pages provide a very basic introduction to Monte Carlo radiative transfer techniques as they are used in TARDIS. All the information listed here can also be found in various papers by L. Lucy and in the main TARDIS publication @@ -12,11 +18,9 @@ can also be found in various papers by L. Lucy and in the main TARDIS publicatio :cite:`Kerzendorf2014`). .. toctree:: + :maxdepth: 2 + basicprinciples initialization propagation - lineinteraction - estimators - virtualpackets - sourceintegration - randomsampling \ No newline at end of file + lineinteraction \ No newline at end of file diff --git a/docs/physics/montecarlo/initialization.ipynb b/docs/physics/montecarlo/initialization.ipynb index 7957c96cd13..e749b6d77d6 100644 --- a/docs/physics/montecarlo/initialization.ipynb +++ b/docs/physics/montecarlo/initialization.ipynb @@ -20,7 +20,7 @@ "metadata": {}, "source": [ "While it is instructive to think about tracking the propagation history of\n", - "photons when illustrating the basic idea behind Monte Carlo radiative transfer\n", + "individual photons when illustrating the basic idea behind Monte Carlo radiative transfer\n", "techniques, there are important numerical reasons for using a different\n", "discretization scheme. Instead of thinking in the photon picture, it brings\n", "significant advantages to follow the idea of [] and\n", @@ -29,7 +29,7 @@ "are commonly referred to as \"energy packets\" or simply \"packets\", and are composed of many photons with the same frequency.\n", "\n", "During a Monte Carlo calculation, $N$ (a large number) packets, all with a certain\n", - "energy $\\varepsilon$, are created at the inner boundary of the computational domain (as will be discussed in [Packet Propagation](propagation.rst)) known as the photosphere. Currently, the photosphere is modeled as a spherical blackbody with a radius $R_\\mathrm{phot}$ and temperature $T_\\mathrm{phot}$. In TARDIS, all packets are assigned identical energies, and the total energy of the packets is 1 erg (and thus each packet has an energy of $\\frac{1}{N}$ ergs).\n", + "energy $\\varepsilon$, are created at the inner boundary of the computational domain (which is discussed in [Model of Supernova Domain](../setup/model.rst)) known as the photosphere. Currently, the photosphere is modeled as a spherical blackbody with a radius $R_\\mathrm{phot}$ and temperature $T_\\mathrm{phot}$. In TARDIS, all packets are assigned identical energies, and the total energy of the packets is 1 erg (and thus each packet has an energy of $\\frac{1}{N}$ ergs).\n", "\n", ".. note:: The indivisible energy packet scheme does not require that all packets have the same energy. This is just a convenient and simple choice adopted in TARDIS.\n", "\n", @@ -329,4 +329,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/docs/physics/montecarlo/lineinteraction.rst b/docs/physics/montecarlo/lineinteraction.rst index 3d32cb260a0..e24389b8142 100644 --- a/docs/physics/montecarlo/lineinteraction.rst +++ b/docs/physics/montecarlo/lineinteraction.rst @@ -96,5 +96,5 @@ shows the situation in the resonant scatter mode, the middle one for the downbranching scheme and the right one the macro atom results. .. image:: - images/scatter_downbranch_ma.png + ../images/scatter_downbranch_ma.png :width: 700 diff --git a/docs/physics/montecarlo/propagation.rst b/docs/physics/montecarlo/propagation.rst index 1f81317fb86..394a58c154a 100644 --- a/docs/physics/montecarlo/propagation.rst +++ b/docs/physics/montecarlo/propagation.rst @@ -14,21 +14,6 @@ spectrum, the main product of a TARDIS calculation. The different spectral features are simply a combined product of the changes in the packet properties induced in the radiation-matter interactions. -.. _expansion: - -Model of Supernova Domain -========================= - -TARDIS models supernovae as expanding homologously. This means that at the beginning of the explosion, the supernova starts at a single point and proceeds to expand radially outward such that the ratio of the velocity of the ejecta to the distance from the ejecta to the supernova's center is uniform throughout the supernova. As an example, if the outer edge of the ejecta moves outward at some velocity :math:`v_\mathrm{outer\_boundary}`, the velocity of the ejecta half way between the outer edge and the center would be :math:`\frac{v_\mathrm{outer\_boundary}}{2}`. The animation below demonstrates this type of expansion. - -TARDIS simulates radiative transfer between an inner boundary (the photosphere, as discussed above) and an outer boundary (the outer edge of the supernova ejecta). The velocity of the inner boundary :math:`v_\mathrm{inner\_boundary}` and the velocity of the outer boundary :math:`v_\mathrm{outer\_boundary}` are supplied in the configuration file (see :ref:`config-file`), as well as the time after the explosion for which TARDIS is calculating the spectrum (:math:`t_\mathrm{explosion}`). The radii of the inner and outer boundaries are therefore calcuated by :math:`r_\mathrm{inner\_boundary}=v_\mathrm{inner\_boundary}*t_\mathrm{explosion}` and :math:`r_\mathrm{outer\_boundary}=v_\mathrm{outer\_boundary}*t_\mathrm{explosion}`. Plasma at a distance :math:`r` from the center of the supernova would then be traveling outward at a speed :math:`v=\frac{r}{r_\mathrm{outer\_boundary}}v_\mathrm{outer\_boundary} = \frac{r}{t_\mathrm{explosion}}`. This is also shown in the animation. - -Additionally, TARDIS divides the space between the inner and outer computational boundaries into cells -- radial shells for which the plasma state is (spatially) constant. In the animation, 6 cells are shown, being divided by the light blue lines. As TARDIS is a time-independent code which calculates the spectra at an instant in time, the radii of the boundaries (either of the computational domain or of the cells) do not chage throughout the simulation. - -.. image:: - images/expansion_animation.gif - :width: 500 - Propagation in a Spherical Domain ================================= @@ -46,7 +31,7 @@ following sketch (taken from :cite:`Noebauer2014`): .. image:: - images/spherical_symmetry.png + ../images/spherical_symmetry.png :width: 500 @@ -62,39 +47,85 @@ Note that the propagation direction has also changed and now takes the value .. math:: \mu_f = \frac{l + r_i \mu_i}{r_f}. - + +Supernova Expansion +=================== + +.. note:: + This section is a summary of part of :ref:`model` which is included here for easy reference. For a complete + explanation, please refer back to that page. + +TARDIS models supernovae as expanding homologously, as shown in the animation below. This means that plasma at a +distance :math:`r` from the center of the supernova will be moving outwards at a velocity +:math:`v=\frac{r}{t_\mathrm{explosion}}`, where :math:`t_\mathrm{explosion}` is the time after the explosion for +which TARDIS is running (which is provided in the :ref:`model configuration `). This is also +shown in the animation. + +TARDIS simulates radiative transfer between an inner boundary (the photosphere) with a radius +:math:`r_\mathrm{inner\_boundary}`, and an outer boundary (the outer edge of the supernova ejecta) with a radius +:math:`r_\mathrm{outer\_boundary}`. Additionally, TARDIS divides the space between the inner and outer computational +boundaries into cells -- radial shells for which the plasma state is (spatially) constant. In the animation, 6 cells +are shown, being divided by the light blue lines. The boundaries of the computational domain and of these cells are +computed during the simulation setup (refer back to :ref:`model`). As TARDIS is a time-independent code, meeaning +that it calculates the spectra at an instant in time (namely at the time :math:`t_\mathrm{explosion}`), the radii of +the boundaries (both of the computational domain and of the cells) do not chage throughout the simulation. + +.. image:: + ../images/expansion_animation.gif + :width: 500 + + .. _referenceframes: Reference Frames ================ -In TARDIS, two reference frames are of particular importance: the lab frame and the co-moving frame. In the lab frame, the center of the supernova is at rest; for example, the animation above is shown in the lab frame. This is the frame for which the spectra are calculated. -The co-moving frame at some point in the supernova, however, has the plasma at that point be at rest. This is the frame of reference "according to the plasma." +Because ejecta in the supernva is moving, TARDIS must take reference frames into account. + +In TARDIS, two reference frames are of particular importance: the lab frame and the co-moving frame. In the lab +frame, the center of the supernova is at rest; for example, the animation above is shown in the lab frame. +This is the frame for which the spectra are calculated. + +The co-moving frame at some point in the supernova, however, has the plasma at that point be at rest. This is the +frame of reference "according to the plasma." -If a photon is propagating in the ejecta with a frequency :math:`\nu_\mathrm{lab}` in the lab frame and a propagation direction :math:`\mu`, the doppler effect says that in the co-moving frame at a distance :math:`r` from the center of the supernova, the photon's frequency is shifted to +If a photon is propagating in the ejecta with a frequency :math:`\nu_\mathrm{lab}` in the lab frame and a propagation +direction :math:`\mu`, the doppler effect says that in the co-moving frame at a distance :math:`r` from the center of +the supernova, the photon's frequency is shifted to .. math:: \nu_\mathrm{co-moving} = \nu_\mathrm{lab}\frac{1-\beta\mu}{\sqrt{1-\beta^2}} -where :math:`\beta = \frac{v}{c} = \frac{r}{ct_\mathrm{explosion}}` (note again that :math:`v` is the velocity of the plasma at a radius :math:`r` from the center of the supernova). The term :math:`\frac{1-\beta\mu}{\sqrt{1-\beta^2}}` is known as the doppler factor. In the nonrelativistic limit (as :math:`v << c`), we get +where :math:`\beta = \frac{v}{c} = \frac{r}{ct_\mathrm{explosion}}` (note again that :math:`v` is the velocity of the +plasma at a radius :math:`r` from the center of the supernova). The term :math:`\frac{1-\beta\mu}{\sqrt{1-\beta^2}}` +is known as the doppler factor. In the nonrelativistic limit (as :math:`v << c`), we get .. math:: \nu_\mathrm{co-moving} = \nu_\mathrm{lab}(1-\beta\mu). -Note that if the photon is propagating away from the center of the supernova (:math:`\mu>0`) it is red-shifted (:math:`\nu_\mathrm{co-moving}<\nu_\mathrm{lab}`), and if the photon is propagating towards the center of the supernova (:math:`\mu<0`) it is blue-shifted (:math:`\nu_\mathrm{co-moving}>\nu_\mathrm{lab}`). +Note that if the photon is propagating away from the center of the supernova (:math:`\mu>0`) it is red-shifted +(:math:`\nu_\mathrm{co-moving}<\nu_\mathrm{lab}`), and if the photon is propagating towards the center of the +supernova (:math:`\mu<0`) it is blue-shifted (:math:`\nu_\mathrm{co-moving}>\nu_\mathrm{lab}`). Numerical and Physical Events ============================= -While a packet is propagating through the computational domain, TARDIS calculates the distance the packet will propagate until it (i) crosses into a new cell and (ii) interacts with the plasma in the ejecta. If the former distance is shorter, the packet will be moved into the new cell (and the plasma properties will be recalculated), and if the latter distance is shorter, the packet will be moved to the location of the interaction, and the interaction will be performed. +While a packet is propagating through the computational domain, TARDIS calculates the distance the packet will +propagate until it (i) crosses into a new cell and (ii) interacts with the plasma in the ejecta. If the former +distance is shorter, the packet will be moved into the new cell (and the plasma properties will be recalculated), and +if the latter distance is shorter, the packet will be moved to the location of the interaction, and the interaction +will be performed. Distance to Next Cell --------------------- .. note:: - In this documentation, and in TARDIS as a whole, the subscripts "inner" and "outer" refer respectively to the inner an outer boundaries of a cell. The subscripts "inner_boundary" and "outer_boundary" refer respectively to the inner and outer boundaries of the computational domain. + In this documentation, and in TARDIS as a whole, the subscripts "inner" and "outer" refer respectively to the + inner an outer boundaries of a cell. The subscripts "inner_boundary" and "outer_boundary" refer respectively to + the inner and outer boundaries of the computational domain. -As previously mentioned, the physical properties of the plasma are stored in a discrete mesh of cells for which the plasma state is spatially constant. As a consequence, whenever a packet propagates into a +As previously mentioned, the physical properties of the plasma are stored in a discrete mesh of cells for which the +plasma state is spatially constant. As a consequence, whenever a packet propagates into a new cell, important quantities which are relevant for performing radiation-matter interactions have to be re-evaluated in accordance with the new state of the ambient material. Thus, during the packet propagation, the @@ -118,19 +149,31 @@ propagation is stopped and it is no longer considered. Instead, processing the next packet of the population is started. Similarly, the propagation is stopped if the packet escapes through the outer surface of the domain. However, in this case the packet contributes to the final emergent spectrum (see :ref:`Spectrum -Formation `). +Formation `). -When a packet is moved into a new cell, as mentioned before, it is moved to the location at which it crosses the boundary, the plasma properties are recalculated, and the propagation direction of the packet is updated (using :math:`\mu_f = \frac{l + r_i \mu_i}{r_f}`). +When a packet is moved into a new cell, as mentioned before, it is moved to the location at which it crosses the +boundary, the plasma properties are recalculated, and the propagation direction of the packet is updated (using +:math:`\mu_f = \frac{l + r_i \mu_i}{r_f}`). Physical Interactions --------------------- -As a packet propagates through the computational domain, physical radiation-matter interactions can trigger changes in the packet properties. The probability that a photon/packet will interact with matter is characterized by its optical depth :math:`\tau`; the probability that a packet will have interacted after going through an optical depth :math:`\Delta \tau` is :math:`1-e^{-\Delta \tau}`. To model this (see :ref:`Random Sampling `), the packet is assigned a random value of optical depth :math:`\tau_\mathrm{interaction} = -\log z` (for another random :math:`z` between 0 and 1), and upon reaching that optical depth, the packet will interact. +As a packet propagates through the computational domain, physical radiation-matter interactions can trigger changes +in the packet properties. The probability that a photon/packet will interact with matter is characterized by its +optical depth :math:`\tau`; the probability that a packet will have interacted after going through an optical depth +:math:`\Delta \tau` is :math:`1-e^{-\Delta \tau}`. To model this (see :ref:`Random Sampling `), the +packet is assigned a random value of optical depth :math:`\tau_\mathrm{interaction} = -\log z` (for another random +:math:`z` between 0 and 1), and upon reaching that optical depth, the packet will interact. -TARDIS considers two different radiation-matter interactions within the simulation: electron scattering and atomic line interactions. As packets propagate, they accumulate optical depth due to the possibility of going through either of these interations. Since the main focus of TARDIS is to calculate optical spectra, +TARDIS considers two different radiation-matter interactions within the simulation: electron scattering and atomic +line interactions. As packets propagate, they accumulate optical depth due to the possibility of going through either +of these interations. Since the main focus of TARDIS is to calculate optical spectra, electron-scatterings are treated in the elastic low-energy limit as classical -Thomson scatterings. In this case, the electron scattering process is frequency-independent. As a consequence to the frequency independence, the rate at which a packet accumulates electron scattering optical depth depends only on the free electron density :math:`n_e`. The optical depth that a Monte Carlo packet accumulates along a path of length :math:`l` due to +Thomson scatterings. In this case, the electron scattering process is frequency-independent. As a consequence to the +frequency independence, the rate at which a packet accumulates electron scattering optical depth depends only on the +free electron density :math:`n_e`. The optical depth that a Monte Carlo packet accumulates along a path of length +:math:`l` due to Thomson scattering is calculated by the formula .. math:: @@ -138,7 +181,9 @@ Thomson scattering is calculated by the formula \Delta \tau = \sigma_{\mathrm{T}} n_e l. The Thomson cross section :math:`\sigma_{\mathrm{T}}`, which is a constant, -appears here. This corresponds to the fact that a packet has a probability of :math:`1-e^{\sigma_{\mathrm{T}} n_e l}` of going through a Thomson scattering prior to traveling a distance :math:`l` (in other words, the probability of the packet making it across a distance :math:`l` without scattering is :math:`e^{\sigma_{\mathrm{T}} n_e l}`). +appears here. This corresponds to the fact that a packet has a probability of :math:`1-e^{\sigma_{\mathrm{T}} n_e l}` +of going through a Thomson scattering prior to traveling a distance :math:`l` (in other words, the probability of the +packet making it across a distance :math:`l` without scattering is :math:`e^{\sigma_{\mathrm{T}} n_e l}`). The situation is complicated by the inclusion of frequency-dependent bound-bound interactions, i.e. interactions with atomic line transitions. @@ -162,17 +207,25 @@ line. The location at which this occurs is referred to as the resonance or Sobolev point. This effectively reduces the line optical depth determination to a pure local problem. -If a packet with a frequency :math:`\nu_\mathrm{lab}` in the lab frame is at a radius :math:`r_i` with a propagation direction :math:`\mu_i`, the distance that the packet must travel to reach the next Sobolev point is calculated by setting the frequency of the packet in the co-moving frame at the Sobolev point equal to the resonant frequency that it will next hit, which we will label :math:`\nu_\mathrm{line}` (which is, of course, in the co-moving frame). Using the nonrelativistic doppler shift formula, we get +If a packet with a frequency :math:`\nu_\mathrm{lab}` in the lab frame is at a radius :math:`r_i` with a propagation +direction :math:`\mu_i`, the distance that the packet must travel to reach the next Sobolev point is calculated by +setting the frequency of the packet in the co-moving frame at the Sobolev point equal to the resonant frequency that +it will next hit, which we will label :math:`\nu_\mathrm{line}` (which is, of course, in the co-moving frame). Using +the nonrelativistic doppler shift formula, we get .. math:: \nu_\mathrm{line} = (1-\beta_f \mu_f)\nu_\mathrm{lab} -where the subscript :math:`f` refers to being at the sobolev point. Using :math:`\beta_f=\frac{r_f}{ct_\mathrm{explosion}}` and :math:`\mu_f = \frac{l + r_i \mu_i}{r_f}`, we get that the distance :math:`l` to the next Sobolev point is +where the subscript :math:`f` refers to being at the sobolev point. Using +:math:`\beta_f=\frac{r_f}{ct_\mathrm{explosion}}` and :math:`\mu_f = \frac{l + r_i \mu_i}{r_f}`, we get that the +distance :math:`l` to the next Sobolev point is .. math:: l = \left( 1-\beta_i\mu_i - \frac{\nu_\mathrm{line}}{\nu_\mathrm{lab}} \right)ct_\mathrm{explosion} = \frac{\nu_{\mathrm{CMF},i}-\nu_\mathrm{line}}{\nu_\mathrm{lab}}ct_\mathrm{explosion} where :math:`\nu_{\mathrm{CMF},i}` is the frequency of the packet in the co-moving frame at the initial position. -At a Sobolev point, the packet instantaneously accumulates optical depth, the value of which is called the Sobolev optical depth :math:`\tau_\mathrm{Sobolev}` (see :ref:`tau_sobolev`). This corresponds to a probability of :math:`1-e^{-\tau_\mathrm{Sobolev}}` that the packet interacts with the atomic line. +At a Sobolev point, the packet instantaneously accumulates optical depth, the value of which is called the Sobolev +optical depth :math:`\tau_\mathrm{Sobolev}` (see :ref:`tau_sobolev`). This corresponds to a probability of +:math:`1-e^{-\tau_\mathrm{Sobolev}}` that the packet interacts with the atomic line. Distance to Next Event ---------------------- @@ -186,26 +239,35 @@ electron-scattering. At the Sobolev point, the accumulated optical depth is instantaneously incremented by the Sobolev optical depth. Afterwards, the procedure is repeated, now with respect to the next transition in the frequency-ordered list of all possible atomic line transitions. The point at -which the accumulated optical depth reaches the randomly generated interaction optical depth :math:`\tau_\mathrm{interaction}` determines the type of interaction the packet performs and at which location in the spatial mesh, as shown with the example cases in the sketch below (taken from :cite:`Noebauer2014`, adapted from +which the accumulated optical depth reaches the randomly generated interaction optical depth +:math:`\tau_\mathrm{interaction}` determines the type of interaction the packet performs and at which location in +the spatial mesh, as shown with the example cases in the sketch below (taken from :cite:`Noebauer2014`, adapted from :cite:`Mazzali1993`): .. image:: - images/optical_depth_summation.png + ../images/optical_depth_summation.png :width: 500 -Three possible cases are highlighted in the above diagram, with the dotted lines showing the (randomly assigned) optical depth :math:`\tau_\mathrm{interaction}` at which the packet interacts. In case I, the interaction optical +Three possible cases are highlighted in the above diagram, with the dotted lines showing the (randomly assigned) +optical depth :math:`\tau_\mathrm{interaction}` at which the packet interacts. In case I, the interaction optical depth value is reached on one of the path segments between successive Sobolev points, where the packet is accumulating electron scattering optical depth. -Thus, the packet performs a Thomson scattering at the point at which its accumulated optical depth reaches :math:`\tau_\mathrm{interaction}`. In case II, the interaction +Thus, the packet performs a Thomson scattering at the point at which its accumulated optical depth reaches +:math:`\tau_\mathrm{interaction}`. In case II, the interaction optical depth is reached during the instantaneous increment by the line optical depth at one of the Sobolev points. As a consequence, the packet performs an -interaction with the corresponding atomic line transition. In both of these cases, the packet is moved to the interaction location, the interaction will be performed (as will be described in the next section), and the process of accumulating optical depth starts over. Finally, if the packet reaches the shell boundary before the optical depth value necessary for a physical interaction is achieved (as in case III), the packet will be moved to the next cell, the plasma properties will be updated, and the accumulation of optical depth will continue in the next cell. +interaction with the corresponding atomic line transition. In both of these cases, the packet is moved to the +interaction location, the interaction will be performed (as will be described in the next section), and the process +of accumulating optical depth starts over. Finally, if the packet reaches the shell boundary before the optical depth +value necessary for a physical interaction is achieved (as in case III), the packet will be moved to the next cell, +the plasma properties will be updated, and the accumulation of optical depth will continue in the next cell. Performing an Interaction ------------------------- -When a physical interaction occurs, whether it is a Thomson scattering or a line interaction, the packet is moved to the interaction location and a new propagation direction is assigned. Since this process is isotropic, the new direction is -sampled according to +When a physical interaction occurs, whether it is a Thomson scattering or a line interaction, the packet is moved to +the interaction location and a new propagation direction is assigned. Since this process is isotropic, the new +direction is sampled according to (see :ref:`Random Sampling `) .. math:: @@ -213,7 +275,8 @@ sampled according to using a new random :math:`z` (between 0 and 1). -For Thomson scattering, the energy of the packet in the co-moving frame is conserved, and thus the new energy and frequency of the packet in the lab frame (due to the doppler effect) is: +For Thomson scattering, the energy of the packet in the co-moving frame is conserved, and thus the new energy and +frequency of the packet in the lab frame (due to the doppler effect) is: .. math:: @@ -221,11 +284,15 @@ For Thomson scattering, the energy of the packet in the co-moving frame is conse \nu_f & = \nu_i \frac{1 - \beta \mu_i}{1 - \beta \mu_f} Here, the subscripts highlight the packet properties before (:math:`i` for -initial) and after (:math:`f` for final) the scattering. Note that :math:`\mu_i` is the propagation direction prior to the interaction **but at the interaction location.** +initial) and after (:math:`f` for final) the scattering. Note that :math:`\mu_i` is the propagation direction prior +to the interaction **but at the interaction location.** -For line interactions, the energy of the packet after the interaction is still given by the same formula (based on energy conservation in the co-moving frame). However, the post-interaction frequency depends on the selected line interaction treatment (see :ref:`Line Interaction Treatments `). +For line interactions, the energy of the packet after the interaction is still given by the same formula (based on +energy conservation in the co-moving frame). However, the post-interaction frequency depends on the selected line +interaction treatment (see :ref:`Line Interaction Treatments `). -The ratio :math:`\frac{1 - \beta \mu_i}{1 - \beta \mu_f}` can be visualized with the following graph for a plasma speed of :math:`1.1 \times 10^4` km/s: +The ratio :math:`\frac{1 - \beta \mu_i}{1 - \beta \mu_f}` can be visualized with the following graph for a plasma +speed of :math:`1.1 \times 10^4` km/s: .. plot:: physics/pyplot/plot_mu_in_out_packet.py diff --git a/docs/physics/montecarlo/randomsampling.rst b/docs/physics/montecarlo/randomsampling.rst deleted file mode 100644 index 21e81a63fd7..00000000000 --- a/docs/physics/montecarlo/randomsampling.rst +++ /dev/null @@ -1,111 +0,0 @@ -.. _randomsampling: - -********************** -Random Sampling Basics -********************** - -Monte Carlo radiative transfer calculations require a probabilistic -representation of the relevant physical processes. This involves assigning -probabilities to the occurrence of certain events or properties. For example, -finding a photon of an isotropic radiation field propagating into any given -direction is equally likely. During the Monte Carlo simulation, assignments -according to these probabilities have to be frequently performed, for example -when initializing Monte Carlo packets. For this purpose, so-called Random -Number Generators are available. These produce (pseudo-) random numbers -:math:`z` uniformly distributed on the interval :math:`[0,1]`. The challenge -now lies in using these numbers to sample any physical process involved in the -Radiative transfer calculation. From a probability theory point of view, this -just implies finding a mapping between the probability distribution governing the -physical process and the one underlying the Random Number Generator. This -process is typically referred to as random sampling. - -Inverse transformation method -============================= - -.. note:: - This is a very superficial and sloppy description of Random sampling. More - detailed and rigorous accounts are found in the standard literature, for - example in :cite:`Kalos2008`. - -The simplest and most-used technique in Monte Carlo radiative transfer -applications is referred to as the inverse transformation method and involves -the cumulative distribution function. In general, a random process captured by -the random variable :math:`X` is governed by a probability density -:math:`\rho_X(x)` (the continuous counterpart to discrete probabilities), with -:math:`\rho_X(x) \mathrm{d}x` describing the probability of the variable taking -values in the interval :math:`[x, x+\mathrm{d}x]`. The cumulative distribution -function in turn describes, as the name suggests, the probability of the -variable taking any value between :math:`-\infty` and :math:`x`: - -.. math:: - - f_X(x) = \int_{-\infty}^x \mathrm{d}x \rho_X(x) - -Since the probability density is by definition always positive, the cumulative -distribution function is monotonically increasing. This constitutes the basis -for the inverse transformation function. Consider two random variables, -:math:`X` and :math:`Y`. A mapping between those may be established by equating -their cumulative distribution functions. Numbers :math:`y` distributed -according to one of the probability densities (in this case :math:`\rho_Y(y)`) -may then be used to sample the other process by - -.. math:: - - x = f_X^{-1}\left(f_Y(y)\right). - -For the Random Number Generators described above, the cumulative distribution -function is trivial, namely :math:`f_Z(z) = z`. However, the inverse -distribution sampling method relies on finding the analytic inverse function of -the cumulative distribution function governing the physical processes to be -sampled. If this is not possible, other sampling methods, such as von-Neumann -rejection sampling techniques, have to be used. - -Examples -======== - -A few examples are provided to illustrate the random sampling process using the -inverse transformation technique. - -Isotropic Scattering --------------------- - -Consider the case of an isotropic scattering. -Here, all scattering angles are equally likely. Thus, the corresponding -(normalized) probability density and the cumulative distribution function are given by - -.. math:: - - \rho_{\mu}(\mu) &= \frac{1}{2}\\ - f_{\mu}(\mu) &= \frac{1}{2} (\mu - 1). - -This leads to the sampling rule - -.. math:: - - \mu = 2 z - 1. - -Next Interaction Location -------------------------- - -The probability of a packet interacting after covering an optical depth -:math:`\tau` is given by - -.. math:: - - \rho_{\tau}(\tau) &= \exp(-\tau)\\ - f_{\tau}(\tau) &= 1 - \exp(-\tau). - - -With the inverse transformation method, the optical depth to the next interaction location may then be sampled by - -.. math:: - - \tau = - \mathrm{ln}(1 - z) - - -which is equivalent to - -.. math:: - - \tau = - \mathrm{ln}z. - diff --git a/docs/physics/setup/plasma/plasma_plots/lte_ionization_balance.py b/docs/physics/pyplot/lte_ionization_balance.py similarity index 100% rename from docs/physics/setup/plasma/plasma_plots/lte_ionization_balance.py rename to docs/physics/pyplot/lte_ionization_balance.py diff --git a/docs/physics/setup/plasma/plasma_plots/nebular_ionization_balance.py b/docs/physics/pyplot/nebular_ionization_balance.py similarity index 100% rename from docs/physics/setup/plasma/plasma_plots/nebular_ionization_balance.py rename to docs/physics/pyplot/nebular_ionization_balance.py diff --git a/docs/physics/setup/index.rst b/docs/physics/setup/index.rst index 0e2763ea389..8a08f646a2c 100644 --- a/docs/physics/setup/index.rst +++ b/docs/physics/setup/index.rst @@ -4,5 +4,12 @@ Setting Up the Simulation ************************* +The first step executed when TARDIS runs is to call an instance of the ``Simulation`` class. This sets up a lot of +things that TARDIS will need during its run. The two main things that are set up are the supernova model and the +initial plasma state (which may be updated throughout the simulation, see :ref:`est_and_conv`). The pages linked +below explain how these are calculated, as well as showing this in action by calling an instance of the +``Simulation`` class. + .. toctree:: + model plasma/index \ No newline at end of file diff --git a/docs/physics/setup/model.rst b/docs/physics/setup/model.rst new file mode 100644 index 00000000000..7a126c12335 --- /dev/null +++ b/docs/physics/setup/model.rst @@ -0,0 +1,32 @@ +.. _model: + +************************* +Model of Supernova Domain +************************* + +TARDIS models supernovae as expanding homologously. This means that at the beginning of the explosion, the +supernova starts at a single point and proceeds to expand radially outward such that the ratio of the velocity of +the ejecta to the distance from the ejecta to the supernova's center is uniform throughout the supernova. As an +example, if the outer edge of the ejecta moves outward at some velocity :math:`v_\mathrm{outer\_boundary}`, the +velocity of the ejecta half way between the outer edge and the center would be +:math:`\frac{v_\mathrm{outer\_boundary}}{2}`. The animation below demonstrates this type of expansion. + +TARDIS simulates radiative transfer between an inner boundary (the photosphere) and an outer +boundary (the outer edge of the supernova ejecta). The velocity of the inner boundary +:math:`v_\mathrm{inner\_boundary}` and the velocity of the outer boundary :math:`v_\mathrm{outer\_boundary}` are +supplied in the configuration file (see :ref:`model-csvy-and-config`), as well as the time after the explosion for +which TARDIS is calculating the spectrum (:math:`t_\mathrm{explosion}`). The radii of the inner and outer boundaries +are therefore calcuated by :math:`r_\mathrm{inner\_boundary}=v_\mathrm{inner\_boundary}*t_\mathrm{explosion}` and +:math:`r_\mathrm{outer\_boundary}=v_\mathrm{outer\_boundary}*t_\mathrm{explosion}`. Plasma at a distance :math:`r` +from the center of the supernova would then be traveling outward at a speed +:math:`v=\frac{r}{r_\mathrm{outer\_boundary}}v_\mathrm{outer\_boundary} = \frac{r}{t_\mathrm{explosion}}`. This is +also shown in the animation. + +Additionally, TARDIS divides the space between the inner and outer computational boundaries into cells -- radial +shells for which the plasma state is (spatially) constant. In the animation, 6 cells are shown, being divided by the +light blue lines. As TARDIS is a time-independent code which calculates the spectra at an instant in time, the radii +of the boundaries (either of the computational domain or of the cells) do not chage throughout the simulation. + +.. image:: + ../images/expansion_animation.gif + :width: 500 \ No newline at end of file diff --git a/docs/physics/setup/plasma/lte_plasma.rst b/docs/physics/setup/plasma/lte_plasma.rst index ff18efaf52a..c7c831cdcbe 100644 --- a/docs/physics/setup/plasma/lte_plasma.rst +++ b/docs/physics/setup/plasma/lte_plasma.rst @@ -45,6 +45,6 @@ the quantities calculated here. Example Calculations ^^^^^^^^^^^^^^^^^^^^ -.. .. plot:: physics/plasma/plasma_plots/lte_ionization_balance.py +.. .. plot:: physics/pyplot/lte_ionization_balance.py :include-source: diff --git a/docs/physics/setup/plasma/nebular_plasma.rst b/docs/physics/setup/plasma/nebular_plasma.rst index bea2d111e72..94f014a26c8 100644 --- a/docs/physics/setup/plasma/nebular_plasma.rst +++ b/docs/physics/setup/plasma/nebular_plasma.rst @@ -86,6 +86,6 @@ the quantities calculated here. Example Calculations ^^^^^^^^^^^^^^^^^^^^ -.. .. plot:: physics/plasma/plasma_plots/nebular_ionization_balance.py +.. .. plot:: physics/pyplot/nebular_ionization_balance.py :include-source: diff --git a/docs/physics/spectrum/basic.rst b/docs/physics/spectrum/basic.rst new file mode 100644 index 00000000000..f193fb07855 --- /dev/null +++ b/docs/physics/spectrum/basic.rst @@ -0,0 +1,23 @@ +.. _basic-spectrum: + +************************* +Basic Spectrum Generation +************************* + +The first and most simple way in which TARDIS calculates spectra calculates it directly from the Monte Carlo +packets after the final :ref:`Monte Carlo iteration `. This simply requires knowledge of the each +packet's energy and frequency in the lab frame (see :ref:`referenceframes`) at the end of the iteration. The only +other quantity needed is the time duration of the simulation :math:`\Delta t`, which isc alculated based off of the +luminosity of the supernova's photosphere (see :ref:`initialization`). + +.. note:: + + The only packets which are used for this calculation are the packets which escape the outer boundary of the + computational domain -- those reabsorbed into the photosphere are not included (see :ref:`propagation`). + +The spectrum calculation is very straightforward. A packet of energy :math:`E_\mathrm{packet}` contributes a +luminosity + +.. math:: L_\mathrm{packet} = \frac{E_\mathrm{packet}}{\Delta t} + +to the spectrum at its frequency. \ No newline at end of file diff --git a/docs/physics/spectrum/index.rst b/docs/physics/spectrum/index.rst index 38e7a872654..aaa7a9b242d 100644 --- a/docs/physics/spectrum/index.rst +++ b/docs/physics/spectrum/index.rst @@ -4,6 +4,17 @@ Spectrum Generation ******************* +During the final :ref:`Monte Carlo iteration `, TARDIS calculates the emitted spectrum. We currently +employ three diiferent methods for doing this: a basic spectrum generation directly from the Monte Carlo packets, +a method using so-called "virtual packets," and by the formal intrgral method. These are all detailed in the links +below. + +.. note:: + For this last iteration, TARDIS uses the number of Monte Carlo packets specified in the :ref:`configuration file + ` under the ``last_no_of_packets`` argument. Because this last iteration is used for + calculating the actual spectrum, users may want to employ more packets than are used in the previous iterations. + .. toctree:: + basic virtualpackets sourceintegration \ No newline at end of file diff --git a/docs/quickstart/quickstart.ipynb b/docs/quickstart/quickstart.ipynb index 1205a7cbd53..bc0442a7df6 100644 --- a/docs/quickstart/quickstart.ipynb +++ b/docs/quickstart/quickstart.ipynb @@ -11,7 +11,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Every simulation run requires the atomic database (for more info refer to [atomic data](../using/components/atomic/atomic_data.rst) ) and a configuration file (more info at [configuration](../using/components/configuration/index.rst) ).\n", + "Every simulation run requires the atomic database (for more info refer to [atomic data](../io/configuration/components/atomic/atomic_data.rst)) and a configuration file (more info at [configuration](../io/configuration/index.rst)).\n", "You can obtain a copy of the atomic database from the\n", "(https://github.com/tardis-sn/tardis-refdata) repository\n", "(`atom_data` subfolder). We recommended to use the\n", @@ -138,4 +138,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file