diff --git a/docs/index.rst b/docs/index.rst index 8bea5336233..c6db0e50d26 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -77,7 +77,7 @@ Mission Statement physics/intro/index physics/setup/index physics/montecarlo/index - physics/est_and_conv/index + physics/update_and_conv/update_and_conv physics/spectrum/index physics/energy_input/index diff --git a/docs/io/configuration/components/models/index.rst b/docs/io/configuration/components/models/index.rst index 7fb87af0e18..6189ce2b72a 100644 --- a/docs/io/configuration/components/models/index.rst +++ b/docs/io/configuration/components/models/index.rst @@ -279,8 +279,8 @@ factor of .8, etc. .. note:: ``t_rad`` and ``dilution_factor`` are the values of the temperature and dilution factor for the first - iteration, and will be updated in subsequent iterations (see :ref:`est_and_conv`). To prevent these - quantities from being changed, you must set the damping constant to zero in the :ref:`Damped Convergence + iteration, and will be updated in subsequent iterations (see :doc:`../../../../physics/update_and_conv/update_and_conv`). + To prevent these quantities from being changed, you must set the damping constant to zero in the :ref:`Damped Convergence Configuration ` in the Monte Carlo section of the configuration file. CSVY Model Tutorial diff --git a/docs/io/configuration/components/montecarlo.rst b/docs/io/configuration/components/montecarlo.rst index d8c8d548744..0f57b2a50ed 100644 --- a/docs/io/configuration/components/montecarlo.rst +++ b/docs/io/configuration/components/montecarlo.rst @@ -16,7 +16,7 @@ of the Monte Carlo loop (which calculates the final spectrum!) when the radiatio ``no_of_packets`` to create a less noisy output spectrum. ``no_of_virtual_packets`` can also be set to greater than 0 (a useful number is 3) to use the Virtual Packet formalism. Increasing this number drastically increases computational costs (and memory requirements if they are logged). The ``iterations`` parameter describes the maximum number of Monte Carlo loops executed in a simulation before it ends. Convergence criteria can be used to make the simulation stop -sooner when the convergence threshold has been reached (see :ref:`convergence`). +sooner when the convergence threshold has been reached (see :doc:`../../../physics/update_and_conv/update_and_conv`). .. _conv-config: diff --git a/docs/io/configuration/components/plasma.rst b/docs/io/configuration/components/plasma.rst index 1275ce71091..bf0b5b3579d 100644 --- a/docs/io/configuration/components/plasma.rst +++ b/docs/io/configuration/components/plasma.rst @@ -23,7 +23,7 @@ The radiative rates describe how to calculate the :math:`J_\textrm{blue}` needed :math:`J_\textrm{blue} = W \times \textrm{Blackbody}(T_\textrm{rad})` 3) ``detailed`` in which the :math:`J_\textrm{blue}` -are calculated using an estimator (this is described in :doc:`../../../physics/est_and_conv/estimators`). +are calculated using an estimator (this is described in :doc:`../../../physics/montecarlo/estimators`). TARDIS currently supports three different kinds of line interaction: ``scatter`` --- a resonance scattering implementation, ``macroatom`` --- the most complex form of line interaction described in :ref:`macroatom` and ``downbranch`` a simplified diff --git a/docs/io/configuration/components/supernova.rst b/docs/io/configuration/components/supernova.rst index 02cb62b864d..f01583c5c92 100644 --- a/docs/io/configuration/components/supernova.rst +++ b/docs/io/configuration/components/supernova.rst @@ -8,7 +8,7 @@ The supernova component of the configuration file contains some key information .. jsonschema:: schemas/supernova.yml -As luminosity (in units of energy/s) is computed by integrating over the spectral luminosity (in units of energy/s/wavelength), TARDIS sums over all discrete energy packets to compute luminosity when ran, attempting to converge the output spectrum to match `luminosity_requested` (see :ref:`est_and_conv`). `luminosity_requested` can be given in standard units, such as erg/s or J/s, or in logarithmic units such as log_lsun. The range over which TARDIS sums these energy packets is set by default from 0 to infinity via `luminosity_wavelength_start` and `luminosity_wavelength_end`, respectively, so as to generate a spectrum whose total luminosity across the entire spectrum is `luminosity_requested`. However, if in the event only the luminosity within a certain range of wavelengths is known, then `luminosity_wavelength_start` and `luminosity_wavelength_end` can be changed as necessary to reflect this, allowing TARDIS to attempt to create a spectrum whose luminosity within the set range will converge to the value defined in `luminosity_requested`. +As luminosity (in units of energy/s) is computed by integrating over the spectral luminosity (in units of energy/s/wavelength), TARDIS sums over all discrete energy packets to compute luminosity when ran, attempting to converge the output spectrum to match `luminosity_requested` (see :doc:`../../../physics/update_and_conv/update_and_conv`). `luminosity_requested` can be given in standard units, such as erg/s or J/s, or in logarithmic units such as log_lsun. The range over which TARDIS sums these energy packets is set by default from 0 to infinity via `luminosity_wavelength_start` and `luminosity_wavelength_end`, respectively, so as to generate a spectrum whose total luminosity across the entire spectrum is `luminosity_requested`. However, if in the event only the luminosity within a certain range of wavelengths is known, then `luminosity_wavelength_start` and `luminosity_wavelength_end` can be changed as necessary to reflect this, allowing TARDIS to attempt to create a spectrum whose luminosity within the set range will converge to the value defined in `luminosity_requested`. As an example, here is a sample code which will generate a specturm where only the luminosity of the visible light portion of the spectrum is given. Here, the output spectrum will have a luminosity of approximately :math:`10^{9.44}L_{sun}` within the visible range. diff --git a/docs/physics/est_and_conv/convergence.rst b/docs/physics/est_and_conv/convergence.rst deleted file mode 100644 index 2b632bcf4a9..00000000000 --- a/docs/physics/est_and_conv/convergence.rst +++ /dev/null @@ -1,78 +0,0 @@ -.. _convergence: - -*********** -Convergence -*********** - -As explained in :doc:`estimators`, after each iteration the values for radiative temperature and dilution factor are updated by calling the ``advance_state`` method on a ``Simulation`` object. The goal of this is to eventually have the radiative temperature and dilution factor converge to a single value so that the steady-state plasma state can be determined. To ensure that the simulation converges, TARDIS employs additional convergence strategies. Currently, only one convergence strategy is available: damped convergence. This will be described in the following sections. - -.. note:: - - Unless otherwise noted, all user-supplied quantities mentioned on this page are supplied in the :ref:`convergence section of the Monte Carlo configuration `, which will be referenced as the convergence configuration. - - -T_rad and W ------------ - -As discussed :doc:`here `, TARDIS uses estimators to calculate estimated radiative temperatures (:math:`T_\mathrm{rad}`) and dilution factors (:math:`W`) in each cell. While TARDIS can then update the plasma state using the estimated values, there is a good chance that these estimated values would “overshoot” the true value we want to converge to (for example, if the current value of the dilution factor in some cell the dilution factor is .4, and the true steady-state value TARDIS wants to find is .45, there is a good chance that the estimated value will be greater than .45). This could make the simulation take longer to converge or, at worst, make it so the simulation does not converge at all. To account for this, users can set (in the convergence configuration) a "damping constant" for both the radiative temperature (:math:`d_{T_\mathrm{rad}}`) and the dilution factor (:math:`d_W`). When ``advance_state`` is called, these quantities update as follows: - -.. math:: - T_\mathrm{rad\ updated} = T_\mathrm{rad\ current} + d_{T_\mathrm{rad}}(T_\mathrm{rad\ estimated}-T_\mathrm{rad\ current}) - -and - -.. math:: - W_\mathrm{updated} = W_\mathrm{current} + d_W(W_\mathrm{estimated}-W_\mathrm{current}). - -This means, for example, if the damping constant is .5, the updated value is halfway between the current value and the estimated value. If the damping constant is .7, the updated value is 70% of the way between the current value and the estimated value, and so on. **If the damping constant is 1, then the updated value is exactly the estimated value, and if the damping constant is zero, the value stays the same throughout the simulation and is not updated.** - - -T_inner -------- - -The temperature of the inner boundary, :math:`T_\mathrm{inner}`, plays a unique role in the simulation, as it is the primary determiner of the output luminosity. This is because the the luminosity of the inner boundary is proportional to :math:`T_\mathrm{inner}^4` (see :doc:`../montecarlo/initialization`). Thus, :math:`T_\mathrm{inner}` is updated throughout the simulation in order to match the output luminosity to the requested luminosity specified in the :doc:`supernova configuration <../../io/configuration/components/supernova>` between the bounds specified in the supernova configuration. However, there is not necessarily a quartic relationship between :math:`T_\mathrm{inner}` and the output luminosity, as changing :math:`T_\mathrm{inner}` also changes the frequency distribution of the initialized packets (once again see :doc:`../montecarlo/initialization`). This then affects the light-matter interactions, affecting which packets make it to the outer boundary, which also affects the output luminosity. Because of this, there is not an exact way to estimate :math:`T_\mathrm{inner}`. To do this estimation, we use - -.. math:: - T_\mathrm{inner\ estimated} = T_\mathrm{inner\ current} * \left(\frac{L_\mathrm{output}}{L_\mathrm{requested}}\right)^{\mathrm{t\_inner\_update\_exponent}} - -where :math:`L_\mathrm{output}` is the output luminosity calculated by adding up the luminosity of each packet (see :doc:`../spectrum/basic`) between the bounds specified in the :doc:`supernova configuration <../../io/configuration/components/supernova>`, :math:`L_\mathrm{requested}` is the luminosity requested also in the supernova configuration (requested between those bounds previously mentioned), and ``t_inner_update_exponent`` is provided by the user in the convergence configuration. Note that what we are doing is "correcting" the previous value of the inner temperature by a factor of :math:`\left(\frac{L_\mathrm{output}}{L_\mathrm{requested}}\right)^{\mathrm{t\_inner\_update\_exponent}}`. Note that if :math:`\frac{L_\mathrm{output}}{L_\mathrm{requested}}` is greater than 1, we want to lower :math:`T_\mathrm{inner}` as the output luminosity is too high, and vice versa if the ratio is less than 1. Thus ``t_inner_update_exponent`` should be negative. Naively one might set ``t_inner_update_exponent=-0.25``, however as a default TARDIS uses ``t_inner_update_exponent=-0.5`` as -0.25 may undershoot the correct :math:`T_\mathrm{inner}` because of its previously mentioned effects on the initial frequency distribution. - -After calculating the estimated :math:`T_\mathrm{inner}`, the quantity is updated using damped convergence with its own damping constant (once again set in the convergence configuration): - -.. math:: - T_\mathrm{inner\ updated} = T_\mathrm{inner\ current} + d_{T_\mathrm{inner}}(T_\mathrm{inner\ estimated}-T_\mathrm{inner\ current}). - -Once again, If the damping constant is 1, then the updated value is exactly the estimated value, and if the damping constant is zero, the value stays the same throughout the simulation and is not updated. - -Additionally, because of the vast impact of :math:`T_\mathrm{inner}` on the simulation, one may want to update it less frequently -- i.e. allow :math:`W` and :math:`T_\mathrm{rad}` to reach a steady-state value for a particular :math:`T_\mathrm{inner}` before updating :math:`T_\mathrm{inner}`. To do this, in the convergence configuration we set ``lock_t_inner_cycles``, which is the number of iterations to wait before updating :math:`T_\mathrm{inner}`. It is set to 1 by default, meaning :math:`T_\mathrm{inner}` would be updated every iteration. - - -Convergence Information ------------------------ - -During the simulation, information about the how :math:`T_\mathrm{rad}`, :math:`W`, and :math:`T_\mathrm{inner}` are updated as well as a comparison of the total output luminosity and the requested luminosity are logged at the INFO level (see :doc:`../../io/optional/logging_configuration`) to give users a better idea of how the convergence process is working. - -In addition, TARDIS allows for the displaying of convergence plots, which allows users to visualize the convergence process for :math:`T_\mathrm{rad}`, :math:`W`, :math:`T_\mathrm{inner}`, and the total luminosity of the supernova being modeled. For more information, see :doc:`../../io/visualization/convergence_plot`. - - -Convergence Criteria --------------------- - -TARDIS also allows users to stop the simulation if the simulation reaches a certain level of convergence. To enable this, users must set ``stop_if_converged=True`` in the convergence configuration. Also in the configuration, the quantities ``hold_iterations``, ``threshold``, and ``fraction`` are be specified to determine convergence as follows: - -For the simulation to be considered to have converged, for ``hold_iterations`` successive iterations, the estimated values of :math:`T_\mathrm{rad}`, :math:`W`, and :math:`T_\mathrm{inner}` may differ from the previous value by a fraction of at most ``threshold`` in at least ``fraction`` fraction of the shells (for :math:`T_\mathrm{inner}`, since there is only one value, the ``fraction`` part does not apply). For example, if ``hold_iterations=3``, ``threshold=0.05`` for all three quantities, and ``fraction=.8``, the simulation will be considered to have converged if for 3 successive iterations the estimated values of :math:`T_\mathrm{rad}` and :math:`W` differ from the current respective values by at most 5% in at least 80% of the shells, *and* the estimated :math:`T_\mathrm{inner}` differs by at most 5%. See the :ref:`convergence configuration schema ` for default values of these quantities. - -.. note:: - - To determine convergence, we compare the estimated value, **not** the updated value (which is related to the estimated value via the damping constant), with the previous value. If :math:`T_\mathrm{inner}` is locked (see the previous section), the estimated value will still be calculated so convergence can be checked as usual. - - -.. note:: - - ``hold_iterations`` and ``fraction`` are universal quantities, i.e. they are each a single value that applies to :math:`T_\mathrm{rad}` and :math:`W`, and for ``hold_iterations`` also :math:`T_\mathrm{inner}`. ``threshold``, on the other hand, is supplied for each quantity separately, so for instance you could require :math:`T_\mathrm{rad}` to differ by less than 1%, :math:`W` to differ by less than 3%, and :math:`T_\mathrm{inner}` to differ by less than 5% for convergence to be reached. - - -Custom Convergence ------------------- - -The custom convergence strategy option is not currently implemented in TARDIS. diff --git a/docs/physics/est_and_conv/estimators.ipynb b/docs/physics/est_and_conv/estimators.ipynb deleted file mode 100644 index 78199f795e7..00000000000 --- a/docs/physics/est_and_conv/estimators.ipynb +++ /dev/null @@ -1,357 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "f346fb77", - "metadata": {}, - "source": [ - "# Volume-Based Estimators\n", - "\n", - "When TARDIS runs, we enter into a loop with two main components: a Monte Carlo iteration occurs, and then the plasma state is updated based on the \"estimators\" described in this page. These estimators use the Monte Carlo packets to estimate how the light-matter interactions in the supernova affect the conditions in the ejecta. This concept was originally developed by [] and successively refined by [], [] and [].\n", - "\n", - "## Theory\n", - "\n", - "### J and nu_bar\n", - "\n", - "Ordinarily, TARDIS is not concerned about the physical amount of time a packet spends traveling through the ejecta. Instead, we consider the \"time of simulation\" $\\Delta t$ which is chosen to be the amount of time in which the photosphere emits the ensemble of packets (see [Energy Packet Initialization](../montecarlo/initialization.ipynb)). When looking at the estimators, a slightly different interpretation of the packets is necessary. Here, we view the packets as not carrying a discrete amount of energy $\\varepsilon$ that is emitted in a time interval $\\Delta t$, but as being a flow of energy that carries an energy $\\varepsilon$ over a time $\\Delta t$ -- that is, each packet is carrying a luminosity (energy per unit time) of $L = \\frac{\\varepsilon}{\\Delta t}$. Now, we can say that if a packet spends a time $\\delta t$ in the supernova's ejecta, it contributes an energy of $L\\delta t= \\frac{\\varepsilon}{\\Delta t}\\delta t$ into the radiation energy of the ejecta.\n", - "\n", - "To account for the effects of the Monte Carlo packets on the ejecta, TARDIS uses the packets to first determine the average radiation energy density $E$ throughout each shell, where the energy density is the total radiation energy in the shell divided by the volume of the shell $V=\\frac{4}{3}\\pi (r_\\mathrm{outer}^3-r_\\mathrm{inner}^3)$. Therefore, we add up the amount of energy each packet contributes to the radiation energy in that shell, and divide by the total volume of the shell:\n", - "$$E=\\frac{1}{V}\\sum_i L_i\\delta t_i=\\frac{1}{V}\\sum_i \\frac{\\varepsilon_i}{\\Delta t}\\delta t_i = \\frac{1}{V\\Delta t}\\sum_i \\varepsilon_i\\delta t_i$$\n", - "where we sum over every Monte Carlo packet in the shell. Note that we are interested in the energy density in the co-moving frame (i.e. the energy density \"according to the plasma,\" see [reference frames](../montecarlo/propagation.rst#reference-frames)). Now, we note that the amount of time the Monte Carlo packet spends in a shell is $\\delta t = \\frac{l_i}{c}$ where $l$ is the distance that the packet travels through the shell. Thus, our estimator is\n", - "$$E=\\frac{1}{V\\Delta t}\\sum_i \\varepsilon_i\\frac{l_i}{c} = \\frac{1}{cV\\Delta t}\\sum_i \\varepsilon_i l_i.$$\n", - "\n", - "Using this energy density, we can then calculate the mean radiation intensity $J$ in that shell using the relation $J=\\frac{c}{4\\pi} E$, which gives us\n", - "$$J=\\frac{1}{4\\pi V\\Delta t}\\sum_i \\varepsilon_i l_i.$$\n", - "Since along any path the co-moving energy of the packet is continuously doppler shifted, we approximate this estimator using the co-moving energy at the beginning of the packet's path (theoretically, the midpoint of the path would be a better option. However, we use the beginning of the path for computational ease at a very small cost to the estimator's accuracy).\n", - "\n", - "Next, we calculate the mean radiation frequency in each shell. For this, in each shell we add up the frequency of each packet weighted by the intensity they contribute to the shell. Remembering that intensity is $\\frac{c}{4\\pi}$ times the energy density, and as before each packet contributes an energy of $\\frac{\\varepsilon_i l_i}{c\\Delta t}$ and thus energy density of $\\frac{\\varepsilon_i l_i}{cV\\Delta t}$ to its shell, we have that each packet contributes an intensity of $\\frac{\\varepsilon_i l_i}{4\\pi V\\Delta t}$ to its shell. So,\n", - "$$\\bar \\nu = \\sum_i \\frac{\\varepsilon_i l_i}{4\\pi V \\Delta t} \\nu_i = \\frac{1}{4\\pi V \\Delta t}\\sum_i \\varepsilon_i \\nu_i l_i$$\n", - "where once again the co-moving energy and frequency of each packet are taken at the beginning of the packet's path.\n", - "\n", - "
\n", - " \n", - "Note\n", - "\n", - "Both estimators take on a different value in each shell.\n", - "\n", - "
\n", - "\n", - "These estimators allow us to calculate the [radiative temperature](../setup/model.ipynb#Temperatures) $T_\\mathrm{rad}$ and [dilution factor](../setup/model.ipynb#Dilution-Factor) $W$ in each shell using\n", - "\n", - "$$T_{\\mathrm{rad}} = \\frac{h}{k_{\\mathrm{B}}} \\frac{\\pi^4}{360 \\zeta(5)} \\frac{\\bar \\nu}{J}$$\n", - "\n", - "and\n", - "\n", - "$$W = \\frac{\\pi J}{\\sigma_{\\mathrm{R}} T_{\\mathrm{rad}}^4}$$\n", - "\n", - "where $h$ is Planck's constant, $k_{\\mathrm{B}}$ is Boltzmann's constant, $\\sigma_{\\mathrm{R}}$ is the Stefan–Boltzmann constant, and $\\zeta$ is the Riemann zeta function. The equation for $W$ comes from the fact that the dilution factor is the ratio of the actual mean intensity to that of a blackbody, which is $J_{\\mathrm{blackbody}}=\\frac{\\sigma_{\\mathrm{R}} T^4}{\\pi}$.\n", - "\n", - "The new $T_\\mathrm{rad}$ and $W$ are then used as inputs to the updated [plasma calculations](../setup/plasma/index.rst) which account for the effect of the Monte Carlo packets on the plasma state (precisely, these calculated $T_\\mathrm{rad}$ and $W$ help determine the $T_\\mathrm{rad}$ and $W$ used as inputs in the plasma calculations -- see [convergence](convergence.rst#t-rad-and-w) for specifics).\n", - "\n", - "\n", - "### J_blue\n", - "\n", - "Another estimator, called the ``J_blue`` or $J^b_{lu}$ estimator, is unlike the two previous estimators discussed. Instead of storing the mean intensity over the entire spectrum, it stores the intensity at a specific frequency. More precisely, since frequency is a continuum, it stores the intensity per unit frequency. In each shell, we record the intensity per unit frequency at the blue end (higher frequency end; this is where the \"$b$\" superscript in $J^b_{lu}$ comes from) of each line transition -- that is, if a line transition $l\\rightarrow u$ (from the lower energy level $l$ to the upper energy level $u$, hence the $lu$ in $J^b_{lu}$) has a frequency $\\nu_{lu}$, the mean intensity between $\\nu_{lu}$ and $\\nu_{lu}+d\\nu$ is $J^b_{lu}d\\nu$. **This means that the** $J^b_{lu}$ **estimator has an entry for each atomic line in each shell.** Now, using our previous $J$ estimator, we have\n", - "$$J^b_{lu}d\\nu = \\frac{1}{4\\pi V\\Delta t}\\sum_i \\varepsilon_i dl_i$$\n", - "where $dl_i$ is the infinitesimal distance that the packet travels while it has a co-moving frequency between $\\nu_{lu}$ and $\\nu_{lu}+d\\nu$.\n", - "\n", - "Now, say the packet with lab frequency $\\nu_\\mathrm{lab}$ has a co-moving frequency of $\\nu_{lu}$ at a radius $r_1$ and propagation direction $\\mu_1$, and it has a co-moving frequency of $\\nu_{lu}+d\\nu$ at a radius $r_2$ and propagation direction $\\mu_2$. Then (see [reference frames](../montecarlo/propagation.rst#reference-frames)):\n", - "$$\\nu_{lu}=\\left(1-\\frac{r_1\\mu_1}{ct_\\mathrm{explosion}}\\right)\\nu_\\mathrm{lab}$$\n", - "and\n", - "$$\\nu_{lu}+d\\nu=\\left(1-\\frac{r_2\\mu_2}{ct_\\mathrm{explosion}}\\right)\\nu_\\mathrm{lab}.$$\n", - "But then subtracting, we get\n", - "$$d\\nu = (r_2\\mu_2-r_1\\mu_1)\\frac{\\nu_\\mathrm{lab}}{ct_\\mathrm{explosion}}=dl*\\frac{\\nu_\\mathrm{lab}}{ct_\\mathrm{explosion}}$$\n", - "(for the last equality, see [propagation in a spherical domain](../montecarlo/propagation.rst#propagation-in-a-spherical-domain)).\n", - "\n", - "But now inputting this into the equation for $J^b_{lu}$ (using $\\frac{dl_i}{d\\nu}=\\frac{ct_\\mathrm{explosion}}{\\nu_\\mathrm{lab,i}}$), we get\n", - "$$J^b_{lu} = \\frac{ct_\\mathrm{explosion}}{4\\pi V\\Delta t}\\sum_i \\frac{\\varepsilon_i}{\\nu_\\mathrm{lab,i}}.$$\n", - "\n", - "\n", - "## Implementation\n", - "\n", - "As previously discussed, a major component of each Monte Carlo iteration is the packet propagation process. During the packet propagation process this step, the $J$ and $\\bar \\nu$ estimators are updates every time a packet is moved to the next event location. Specifically, every time a packet is moved, $\\varepsilon l$ is added to the \"running total\" $J$ estimator in the shell where the packet is, and $\\varepsilon \\nu l$ is added to the \"running total\" $\\bar\\nu$ estimator in the shell where the packet is (where $l$ is the distance the packet is moved, and $\\varepsilon$ and $\\nu$ are respectively the packet's co-moving energy and frequency at the beginning of the packet's path). The factor of $\\frac{1}{4\\pi V\\Delta t}$, for computational ease, is not attached to the estimators but is included during the calculation of $T_\\mathrm{rad}$ and $W$ from the estimators. Specifically, we use\n", - "$$T_\\mathrm{rad}=\\frac{h}{k_{\\mathrm{B}}} \\frac{\\pi^4}{360 \\zeta(5)} \\frac{\\sum_i \\varepsilon_i \\nu_i l_i}{\\sum_i \\varepsilon_i l_i} = \\frac{h}{k_{\\mathrm{B}}} \\frac{\\pi^4}{360 \\zeta(5)} \\frac{\\mathrm{real\\ nu\\_ bar\\ estimator}}{\\mathrm{real\\ J\\ estimator}}$$\n", - "and\n", - "$$W = \\frac{\\sum_i \\varepsilon_i l_i}{4\\sigma_{\\mathrm{R}} T_{\\mathrm{rad}}^4V\\Delta t} = \\frac{\\mathrm{real\\ J\\ estimator}}{4\\sigma_{\\mathrm{R}} T_{\\mathrm{rad}}^4V\\Delta t}$$\n", - "\n", - "After each Monte Carlo iteration, the `advance_state()` method is called on the `Simulation` object. The estimators are then used to update the radiative temperature and dilution factor according to the [convergence strategy](convergence.rst), and the plasma state is recalculated (see [plasma calculations](../setup/plasma/index.rst)) using the updated radiative temperature and dilution factor as inputs. This process repeats until the final iteration or until convergence has been reached (see [convergence](convergence.rst#Convergence-Criteria)).\n", - "\n", - "Similarly, during the propagation process, every time a packet passes through a Sobolev point, meaning it comes in resonance with an atomic line (and thus reaches the frequency targeted by $J^b_{lu}$ -- not necessarily going through a line interaction), the $J^b_{lu}$ for that atomic transition in the shell it is in is incremented by $\\frac{\\varepsilon}{\\nu_\\mathrm{lab}}$, where $\\varepsilon$ is the packet's energy qnd $\\nu_\\mathrm{lab}$ is the packet's lab-frame frequency. As before, for computational ease, the factor $\\frac{ct_\\mathrm{explosion}}{4\\pi V \\Delta t}$ is included in any calculations using the estimator.\n", - "\n", - "
\n", - " \n", - "Note\n", - "\n", - "Since the ``J_blue`` estimator is updated every time a packet comes into resonance with an atomic line (not necessarily going through a line interaction), the estimator is only equal to zero in some shell for a specific line if no packet comes into resonance with that line within the shell.\n", - "\n", - "
\n", - "\n", - "If set to detailed mode (see [plasma configuration](../../io/configuration/components/plasma.rst)), the `J_blue` plasma property will will be replaced with the value of the $J^b_{lu}$ estimator (the raw estimator times the factor of $\\frac{ct_\\mathrm{explosion}}{4\\pi V \\Delta t}$). Otherwise, the `J_blue` in the plasma are calculated as they typically are in the plasma calculations, and the $J^b_{lu}$ estimator is only used for the [formal integral](../spectrum/formal_integral.rst)." - ] - }, - { - "cell_type": "markdown", - "id": "127746c2", - "metadata": {}, - "source": [ - "## Code Example\n", - "\n", - "We now show a detailed example of how the plasma is updated using the estimators after a Monte Carlo iteration. First, we import the needed packages and set up a simulation (see [Setting Up the Simulation](../setup/index.rst)):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7a8bd905", - "metadata": {}, - "outputs": [], - "source": [ - "from tardis.io.config_reader import Configuration\n", - "from tardis.simulation import Simulation\n", - "from tardis.io.atom_data.util import download_atom_data\n", - "import numpy as np\n", - "\n", - "# We download the atomic data needed to run this notebook\n", - "download_atom_data('kurucz_cd23_chianti_H_He')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "259a7aa9", - "metadata": {}, - "outputs": [], - "source": [ - "tardis_config = Configuration.from_yaml('tardis_example.yml')\n", - "sim = Simulation.from_config(tardis_config)\n", - "\n", - "model = sim.model\n", - "plasma = sim.plasma\n", - "runner = sim.runner" - ] - }, - { - "cell_type": "markdown", - "id": "cf13d946", - "metadata": {}, - "source": [ - "We show the initial radiative temperature, dilution factor, electron densities, and tau sobolevs in each shell:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e3f43f93", - "metadata": {}, - "outputs": [], - "source": [ - "model.t_rad" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "90bb6147", - "metadata": {}, - "outputs": [], - "source": [ - "model.w" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7b2ae7b4", - "metadata": {}, - "outputs": [], - "source": [ - "plasma.electron_densities" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "90ab14b0", - "metadata": {}, - "outputs": [], - "source": [ - "plasma.tau_sobolevs" - ] - }, - { - "cell_type": "markdown", - "id": "a248ff2c", - "metadata": {}, - "source": [ - "We set the number of packets and we run one iteration of the Monte Carlo simulation:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9f7eb44f", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "N_packets = 10000\n", - "\n", - "# Using the commented out code below, we can also get the number of packets\n", - "# from the configuration -- try it out:\n", - "#N_packets = tardis_config.no_of_packets\n", - "\n", - "sim.iterate(N_packets)" - ] - }, - { - "cell_type": "markdown", - "id": "d3d0074d", - "metadata": {}, - "source": [ - "We now show the values of the three estimators previously mentioned (note that these are the raw estimators, and the factors of $\\frac{1}{4\\pi V \\Delta t}$ etc are not included):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "21500ea8", - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "runner.j_estimator" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fac91ee2", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "runner.nu_bar_estimator" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6ff7a82c", - "metadata": {}, - "outputs": [], - "source": [ - "# Because most rows of the j_blue estimatior are partially or mostly\n", - "# zero, we show just rows with all nonzero elements\n", - "runner.j_blue_estimator[np.all(runner.j_blue_estimator != 0, axis=1)]" - ] - }, - { - "cell_type": "markdown", - "id": "cffa939a", - "metadata": {}, - "source": [ - "We note that the shape of the j_blue estimator and the tau_sobolevs are the same: namely, each contain a value for each possible atomic line transition in each radial cell (as opposed to the other two estimators which just have one value for each cell):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a96c5c5b", - "metadata": {}, - "outputs": [], - "source": [ - "plasma.tau_sobolevs.shape, runner.j_blue_estimator.shape" - ] - }, - { - "cell_type": "markdown", - "id": "43b31d2f", - "metadata": {}, - "source": [ - "We now advance the state of the simulation based on the estimators, and demonstrate this by showing the four quantities we showed before running the simulation. Compare them with their values above!" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9676b22b", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "# When advance_state is called, a brief summary of the updated t_rad's\n", - "# and w's is given\n", - "sim.advance_state();" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fe0b9f40", - "metadata": {}, - "outputs": [], - "source": [ - "model.t_rad" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8de89bb2", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "model.w" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fcaf08ff", - "metadata": {}, - "outputs": [], - "source": [ - "plasma.electron_densities" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "48fbee16", - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "plasma.tau_sobolevs" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.12" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/docs/physics/est_and_conv/index.rst b/docs/physics/est_and_conv/index.rst deleted file mode 100644 index 12efba92450..00000000000 --- a/docs/physics/est_and_conv/index.rst +++ /dev/null @@ -1,32 +0,0 @@ -.. _est_and_conv: - -************************** -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 "balanced 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 estimators 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 - convergence diff --git a/docs/physics/intro/index.rst b/docs/physics/intro/index.rst index d61ec2e760b..b3e1bd4e9fc 100644 --- a/docs/physics/intro/index.rst +++ b/docs/physics/intro/index.rst @@ -9,7 +9,7 @@ How TARDIS Works The goal of TARDIS is, given input information about a supernova, to determine (i) properties of the plasma making up the supernova and (ii) the spectrum of light that is emitted from the supernova. -The physics of TARDIS is in four major parts, which are summarized here and in the diagram below. First, the TARDIS simulation is set up (:doc:`../setup/index`). This involves the creation of the supernova model and the initial conditions of the supernova's plasma. Next is the Monte Carlo Iteration (:doc:`../montecarlo/index`) where the heart of TARDIS takes place; packets of light are sent through the supernova and tracked as they interact with matter. Next, TARDIS uses information from the Monte Carlo iteration to update properties of the plasma to eventually find the correct plasma state (:doc:`../est_and_conv/index`). This process of doing a Monte Carlo iteration and then updating the plasma is repeated for a specified number of times or until certain aspects of the plasma state converge (as is discussed in :ref:`convergence`). After that, data generated in the Monte Carlo simulation is used to synthesize the output spectrum of the supernova (:doc:`../spectrum/index`). +The physics of TARDIS is in four major parts, which are summarized here and in the diagram below. First, the TARDIS simulation is set up (:doc:`../setup/index`). This involves the creation of the supernova model and the initial conditions of the supernova's plasma. Next is the Monte Carlo Iteration (:doc:`../montecarlo/index`) where the heart of TARDIS takes place; packets of light are sent through the supernova and tracked as they interact with matter. Next, TARDIS uses information from the Monte Carlo iteration to update properties of the plasma to eventually find the correct plasma state (:doc:`../update_and_conv/update_and_conv`). This process of doing a Monte Carlo iteration and then updating the plasma is repeated for a specified number of times or until certain aspects of the plasma state converge (as is also discussed in :doc:`../update_and_conv/update_and_conv`). After that, data generated in the Monte Carlo simulation is used to synthesize the output spectrum of the supernova (:doc:`../spectrum/index`). .. image:: images/tardis_flowchart.png :width: 200 diff --git a/docs/physics/montecarlo/basicprinciples.rst b/docs/physics/montecarlo/basicprinciples.rst index adee47f05bd..148d4bf2fa0 100644 --- a/docs/physics/montecarlo/basicprinciples.rst +++ b/docs/physics/montecarlo/basicprinciples.rst @@ -12,12 +12,13 @@ successful and elegant tool particularly for radiative transfer problems in supe 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`) +sense) assigned in accordance with the macroscopic properties of the radiation field (see :doc:`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`). +conclusions about the propagation of light through the ejecta as a whole (see :ref:`estimators`). This is eventually +used to determine the actual steady-state plasma properties (see :doc:`../update_and_conv/update_and_conv`) and the +emitted spectrum (see :ref:`spectrum`). .. _randomsampling: diff --git a/docs/physics/montecarlo/estimators.rst b/docs/physics/montecarlo/estimators.rst new file mode 100644 index 00000000000..1163580d8bc --- /dev/null +++ b/docs/physics/montecarlo/estimators.rst @@ -0,0 +1,94 @@ +.. _estimators: + +*********************** +Volume-Based Estimators +*********************** + +Besides from just tracking the propagation of our packets, TARDIS also uses the Monte Carlo iteration is to determine useful information about the light traveling through the supernova (also called the radiation field). This information will eventually be used to help :doc:`update the plasma state <../update_and_conv/update_and_conv>` as well as :doc:`generate different types of spectra <../spectrum/index>`. We determine this information through volume-based estimators. The concept was originally developed by :cite:`Lucy1999` and successively refined by :cite:`Lucy1999a`, :cite:`Lucy2002` and :cite:`Lucy2003`. + + +Theory +====== + +:math:`J` and :math:`\bar \nu` Estimators +----------------------------------------- + +Ordinarily, TARDIS is not concerned about the physical amount of time a packet spends traveling through the ejecta. Instead, we consider the "time of simulation" :math:`\Delta t` which is chosen to be the amount of time in which the photosphere emits the ensemble of packets (see :doc:`Energy Packet Initialization `). When looking at the estimators, a slightly different interpretation of the packets is necessary. Here, we view the packets as not carrying a discrete amount of energy :math:`\varepsilon` that is emitted in a time interval :math:`\Delta t`, but as being a flow of energy that carries an energy :math:`\varepsilon` over a time :math:`\Delta t` -- that is, each packet is carrying a luminosity (energy per unit time) of :math:`L = \frac{\varepsilon}{\Delta t}`. Now, we can say that if a packet spends a time :math:`\delta t` in the supernova's ejecta, it contributes an energy of :math:`L\delta t= \frac{\varepsilon}{\Delta t}\delta t` into the radiation energy of the ejecta. + +To account for the effects of the Monte Carlo packets on the ejecta, TARDIS uses the packets to first determine the average radiation energy density :math:`E` throughout each shell, where the energy density is the total radiation energy in the shell divided by the volume of the shell :math:`V=\frac{4}{3}\pi (r_\mathrm{outer}^3-r_\mathrm{inner}^3)`. Therefore, we add up the amount of energy each packet contributes to the radiation energy in that shell, and divide by the total volume of the shell: + +.. math:: E=\frac{1}{V}\sum_i L_i\delta t_i=\frac{1}{V}\sum_i \frac{\varepsilon_i}{\Delta t}\delta t_i = \frac{1}{V\Delta t}\sum_i \varepsilon_i\delta t_i + +where we sum over every Monte Carlo packet in the shell. Note that we are interested in the energy density in the co-moving frame (i.e. the energy density "according to the plasma," see :ref:`referenceframes`). Now, we note that the amount of time the Monte Carlo packet spends in a shell is :math:`\delta t = \frac{l_i}{c}` where :math:`l` is the distance that the packet travels through the shell. Thus, our estimator is + +.. math:: E=\frac{1}{V\Delta t}\sum_i \varepsilon_i\frac{l_i}{c} = \frac{1}{cV\Delta t}\sum_i \varepsilon_i l_i. + +Using this energy density, we can then calculate the mean radiation intensity :math:`J` in that shell using the relation :math:`J=\frac{c}{4\pi} E`, which gives us + +.. math:: J=\frac{1}{4\pi V\Delta t}\sum_i \varepsilon_i l_i. + +Since along any path the co-moving energy of the packet is continuously doppler shifted, we approximate this estimator using the co-moving energy at the beginning of the packet's path (theoretically, the midpoint of the path would be a better option. However, we use the beginning of the path for computational ease at a very small cost to the estimator's accuracy). + +Next, we calculate the mean radiation frequency in each shell. For this, in each shell we add up the frequency of each packet weighted by the intensity they contribute to the shell (and then divide by the total intensity, as will be discussed below). Remembering that intensity is :math:`\frac{c}{4\pi}` times the energy density, and as before each packet contributes an energy of :math:`\frac{\varepsilon_i l_i}{c\Delta t}` and thus energy density of :math:`\frac{\varepsilon_i l_i}{cV\Delta t}` to its shell, we have that each packet contributes an intensity of :math:`\frac{\varepsilon_i l_i}{4\pi V\Delta t}` to its shell. So, + +.. math:: \bar \nu = \sum_i \frac{\varepsilon_i l_i}{4\pi V \Delta t} \nu_i = \frac{1}{4\pi V \Delta t}\sum_i \varepsilon_i \nu_i l_i + +where once again the co-moving energy and frequency of each packet are taken at the beginning of the packet's path. + +It is important to note that, confusingly, :math:`\bar \nu` is not truly the mean frequency (as can be seen by its units -- it has dimensions of intensity times frequency). Indeed, the true mean frequency is actually :math:`\frac{\bar \nu}{J}`. + +.. note:: Both estimators take on a different value in each shell. + +These estimators will be used in the :doc:`../update_and_conv/update_and_conv` step to help update the plasma state between iterations. + + +.. _j-blue-estimator: + +:math:`J^b_{lu}` Estimator +-------------------------- + +Another estimator, called the ``j_blue`` or :math:`J^b_{lu}` estimator, is unlike the two previous estimators discussed. Instead of storing the mean intensity over the entire spectrum, it stores the intensity at a specific frequency. More precisely, since frequency is a continuum, it stores the intensity per unit frequency. In each shell, we record the intensity per unit frequency at the blue end (higher frequency end; this is where the ":math:`b`" superscript in :math:`J^b_{lu}` comes from) of each line transition -- that is, if a line transition :math:`l\rightarrow u` (from the lower energy level :math:`l` to the upper energy level :math:`u`, hence the :math:`lu` in :math:`J^b_{lu}`) has a frequency :math:`\nu_{lu}`, the mean intensity between :math:`\nu_{lu}` and :math:`\nu_{lu}+d\nu` is :math:`J^b_{lu}d\nu`. **This means that the** :math:`J^b_{lu}` **estimator contains values for each atomic line in each shell**, and is hence called a *line estimator*. Now, using our previous :math:`J` estimator, we have + +.. math:: J^b_{lu}d\nu = \frac{1}{4\pi V\Delta t}\sum_i \varepsilon_i dl_i + +where :math:`dl_i` is the infinitesimal distance that the packet travels while it has a co-moving frequency between :math:`\nu_{lu}` and :math:`\nu_{lu}+d\nu` (here, therefore, we are summing over all packets in a shell that achieve a co-moving frequency of :math:`\nu_{lu}` and thus come into resonance with the line transition :math:`l\rightarrow u` within that shell). + +Now, say the packet with lab frequency :math:`\nu_\mathrm{lab}` has a co-moving frequency of :math:`\nu_{lu}` at a radius :math:`r_1` and propagation direction :math:`\mu_1`, and it has a co-moving frequency of :math:`\nu_{lu}+d\nu` at a radius :math:`r_2` and propagation direction :math:`\mu_2`. Then (see :ref:`referenceframes`): + +.. math:: \nu_{lu}=\left(1-\frac{r_1\mu_1}{ct_\mathrm{explosion}}\right)\nu_\mathrm{lab} + +and + +.. math:: \nu_{lu}+d\nu=\left(1-\frac{r_2\mu_2}{ct_\mathrm{explosion}}\right)\nu_\mathrm{lab}. + +But then subtracting, we get + +.. math:: d\nu = (r_2\mu_2-r_1\mu_1)\frac{\nu_\mathrm{lab}}{ct_\mathrm{explosion}}=dl*\frac{\nu_\mathrm{lab}}{ct_\mathrm{explosion}} + +(for the last equality, see :ref:`spherical-domain`). + +But now inputting this into the equation for :math:`J^b_{lu}` (using :math:`\frac{dl_i}{d\nu}=\frac{ct_\mathrm{explosion}}{\nu_\mathrm{lab,i}}`), we get + +.. math:: J^b_{lu} = \frac{ct_\mathrm{explosion}}{4\pi V\Delta t}\sum_i \frac{\varepsilon_i}{\nu_\mathrm{lab,i}}. + + +.. _edotlu: + +:math:`\dot{E}_{lu}` Estimator +------------------------------ + +The final estimator we consider, like the ``j_blue`` estimator, is a line estimator, meaning it has contains values for each atomic line in each shell. It calculates the rate at which energy density interacts with a line transition :math:`l\rightarrow u`. The first step is to calculate the rate at which energy density resonates with some line in some shell. Each packet that comes into resonance with the transition :math:`l\rightarrow u` in a shell of volume :math:`V` contributes an energy density to that shell of :math:`\frac{\varepsilon}{V}` over a time :math:`\Delta t`, meaning the rate at which energy density resonates with the line is :math:`\sum_i \frac{\varepsilon_i}{\Delta t V} = \frac{1}{\Delta t V} \sum \varepsilon` where we are summing over all packets which come into resonance with the atomic line in some shell (as usual, this sum is done using the energies in the co-moving frame). Finally, this light then has a :math:`\left( 1- e^{-\tau_{lu}}\right)` probability of interacting with the line (where :math:`\tau_{lu}` is the Sobolev optical depth for the transition :math:`l\rightarrow u`, see :ref:`physical-interactions`), so the rate at which energy density is absorbed into the transition :math:`l\rightarrow u`, called the ``Edotlu`` estimator, is + +.. math:: \dot{E}_{lu} = \frac{1}{\Delta t V} \left( 1- e^{-\tau_{lu}}\right) \sum_i \varepsilon_i. + +Note that while one may expect us to merely add up the contributions of each packet that *interacts* with the line, eliminating the need for the :math:`\left( 1- e^{-\tau_{lu}}\right)` term, the former determines the desired quantity with more accuracy and less noise than the latter, (this is because it does not depend on the limited number of packets TARDIS uses, rather a theoretical equation, to determine how much of the light that comes into resonance with a line actually interacts with that line). + + +Implementation +============== + +As previously discussed, a major component of each Monte Carlo iteration is the packet propagation process. During the packet propagation process this step, the :math:`J` and :math:`\bar \nu` estimators are updates every time a packet is moved to the next event location. Specifically, every time a packet is moved, :math:`\varepsilon l` is added to the "running total" :math:`J` estimator in the shell where the packet is, and :math:`\varepsilon \nu l` is added to the "running total" :math:`\bar\nu` estimator in the shell where the packet is (where :math:`l` is the distance the packet is moved, and :math:`\varepsilon` and :math:`\nu` are respectively the packet's co-moving energy and frequency at the beginning of the packet's path). The factor of :math:`\frac{1}{4\pi V\Delta t}`, for computational ease, is not attached to the estimators but is included during any calculations using these estimators, see :doc:`../update_and_conv/update_and_conv`. + +Additionally, during the propagation process, every time a packet passes through a Sobolev point, meaning it reaches a co-moving frequency of :math:`\nu_{lu}` for some transition :math:`l\rightarrow u` and thus comes in resonance with an atomic line, the :math:`J^b_{lu}` for that atomic transition in the shell it is in is incremented by :math:`\frac{\varepsilon}{\nu_\mathrm{lab}}`, where :math:`\varepsilon` is the packet's energy and :math:`\nu_\mathrm{lab}` is the packet's lab-frame frequency. As before, for computational ease, the factor :math:`\frac{ct_\mathrm{explosion}}{4\pi V \Delta t}` is included in calculations involving the estimator (such as when `updating <../update_and_conv/update_and_conv.ipynb#updating-other-quantities>`_ :math:`J^b_{lu}` in the plasma or in the :ref:`formal integral `). Similarly, when a packet passes through a Sobolev point, the :math:`\dot{E}_{lu}` for that atomic transition in the shell it is in is incremented by :math:`\varepsilon`, and once again, for computational ease, the term :math:`\frac{1}{\Delta t V} \left( 1- e^{-\tau_{lu}}\right)` is not included until calculations involving the estimator are performed (specifically in the :ref:`formal integral `). + +.. note:: Since the ``j_blue`` and ``Edotlu`` estimators are updated every time a packet comes into resonance with an atomic line (not necessarily going through a line interaction), these estimators are equal to zero in some shell for a specific line if (and only if) no packet comes into resonance with that line within the shell. diff --git a/docs/physics/montecarlo/index.rst b/docs/physics/montecarlo/index.rst index d4c69971330..475ff55d6ef 100644 --- a/docs/physics/montecarlo/index.rst +++ b/docs/physics/montecarlo/index.rst @@ -6,7 +6,7 @@ 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 +in the :ref:`Monte Carlo Configuration `. As will be decribed in :doc:`../update_and_conv/update_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. @@ -23,4 +23,5 @@ can also be found in various papers by L. Lucy and in the main TARDIS publicatio basicprinciples initialization propagation - lineinteraction \ No newline at end of file + lineinteraction + estimators \ No newline at end of file diff --git a/docs/physics/montecarlo/initialization.ipynb b/docs/physics/montecarlo/initialization.ipynb index 70437ec99a0..e6b93f8f4a5 100644 --- a/docs/physics/montecarlo/initialization.ipynb +++ b/docs/physics/montecarlo/initialization.ipynb @@ -19,7 +19,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 (which is discussed in the [model section](../setup/model.rst)) known as the photosphere. Currently, the photosphere is modeled as a spherical [blackbody](https://en.wikipedia.org/wiki/Black-body_radiation) with a radius $r_\\mathrm{boundary\\_inner}$ and temperature $T_\\mathrm{inner}$. Both of these quantities are calculated as a part of the [model](../setup/model.ipynb), and $T_\\mathrm{inner}$ is additionally updated throughout the simulation as a part of the [convergence](../est_and_conv/convergence.rst) process.\n", + "energy $\\varepsilon$, are created at the inner boundary of the computational domain (which is discussed in the [model section](../setup/model.rst)) known as the photosphere. Currently, the photosphere is modeled as a spherical [blackbody](https://en.wikipedia.org/wiki/Black-body_radiation) with a radius $r_\\mathrm{boundary\\_inner}$ and temperature $T_\\mathrm{inner}$. Both of these quantities are calculated as a part of the [model](../setup/model.ipynb), and $T_\\mathrm{inner}$ is additionally updated throughout the simulation as a part of the [Updating Plasma and Convergence](../update_and_conv/update_and_conv.ipynb) process.\n", "\n", "In TARDIS, all packets are assigned identical energies **in the lab frame** (see [Reference Frames](propagation.rst#reference-frames)), and the total (lab-frame) energy of the packets is 1 erg (and thus each packet has an energy of $\\frac{1}{N}$ ergs).\n", "\n", @@ -48,7 +48,7 @@ " \n", "Note\n", "\n", - "As will be shown in the code example, this will lead to unphysically small values for $\\Delta t$. It may be easier to think of the Monte Carlo packets not as packets of energy $\\epsilon$ going through a simulation of duration $\\Delta t$, but as packets of luminosity that carry an energy $\\epsilon$ over a time $\\Delta t$ (and thus truly being luminosity packets of luminosity $\\frac{\\epsilon}{\\Delta t}$). Indeed, this view of the packets will be useful when deriving the [Monte Carlo Estimators](../est_and_conv/estimators.ipynb).\n", + "As will be shown in the code example, this will lead to unphysically small values for $\\Delta t$. It may be easier to think of the Monte Carlo packets not as packets of energy $\\epsilon$ going through a simulation of duration $\\Delta t$, but as packets of luminosity that carry an energy $\\epsilon$ over a time $\\Delta t$ (and thus truly being luminosity packets of luminosity $\\frac{\\epsilon}{\\Delta t}$). Indeed, this view of the packets will be useful when deriving the [Monte Carlo Estimators](estimators.rst).\n", "\n", "\n", "\n", @@ -116,7 +116,7 @@ "id": "450faf76", "metadata": {}, "source": [ - "We set the temperature of the photosphere $T_\\mathrm{inner}$, which will determine the photospheric luminosity (in an actual simulation, $T_\\mathrm{inner}$ is initially calculated as a part of the [model](../setup/model.ipynb) and updated as a part of the [convergence](../est_and_conv/convergence.rst) process):" + "We set the temperature of the photosphere $T_\\mathrm{inner}$, which will determine the photospheric luminosity (in an actual simulation, $T_\\mathrm{inner}$ is initially calculated as a part of the [model](../setup/model.ipynb) and updated as a part of the [Updating Plasma and Convergence](../update_and_conv/update_and_conv.ipynb) process):" ] }, { @@ -324,7 +324,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.12" + "version": "3.8.13" } }, "nbformat": 4, diff --git a/docs/physics/montecarlo/propagation.rst b/docs/physics/montecarlo/propagation.rst index ccd76671077..2e49379cdbb 100644 --- a/docs/physics/montecarlo/propagation.rst +++ b/docs/physics/montecarlo/propagation.rst @@ -6,7 +6,7 @@ Packet Propagation The bulk of a Monte Carlo Radiative Transfer calculation is spent on determining the propagation history of the different packets. After a packet is -initialized (see :ref:`initialization`), it is launched and may then perform interactions with the +initialized (see :doc:`initialization`), it is launched and may then perform interactions with the surrounding material. This occurs again in a probabilistic manner. The packet propagation is followed until it escapes through the outer boundary of the computational domain, at which point the packet contributes to the synthetic @@ -14,6 +14,8 @@ 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. +.. _spherical-domain: + Propagation in a Spherical Domain ================================= diff --git a/docs/physics/setup/index.rst b/docs/physics/setup/index.rst index 20ced5b27c4..449f32445b2 100644 --- a/docs/physics/setup/index.rst +++ b/docs/physics/setup/index.rst @@ -6,7 +6,7 @@ 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 main things that are set up are the supernova model (a ``Radial1DModel`` object), the -initial plasma state (a ``BasePlasma`` object, which may be updated throughout the simulation, see :ref:`est_and_conv`), +initial plasma state (a ``BasePlasma`` object, which may be updated throughout the simulation, see :doc:`../update_and_conv/update_and_conv`), and a ``MonteCarloRunner`` object. The pages linked below explain how the former two are calculated (the latter is used mainly in the :doc:`next step of the calculation <../montecarlo/index>`), as well as showing this in action by calling an instance of the ``Simulation`` class. diff --git a/docs/physics/setup/model.ipynb b/docs/physics/setup/model.ipynb index 74d051ec3fa..6451c28488f 100644 --- a/docs/physics/setup/model.ipynb +++ b/docs/physics/setup/model.ipynb @@ -477,9 +477,9 @@ "\n", "$$T_{\\mathrm{inner}}=\\left(\\frac{L_\\mathrm{requested}}{4 \\pi r_{\\mathrm{boundary\\_inner}}^2 \\sigma_{\\mathrm{R}}}\\right)^{1/4}$$\n", "\n", - "where $\\sigma_{\\mathrm{R}}$ is the Stefan-Boltzmann constant and $r_{\\mathrm{boundary\\_inner}}$ is once again the radius of the photosphere, calculated as part of the shell structure. Because of light-matter interactions, the output luminosity of the supernova will not be the same as the luminosity of the photosphere, so the photospheric temperature is updated throughout the simulation as part of the convergence process in order to match the output luminosity to the requested luminosity (see [convergence](../est_and_conv/convergence.rst)).\n", + "where $\\sigma_{\\mathrm{R}}$ is the Stefan-Boltzmann constant and $r_{\\mathrm{boundary\\_inner}}$ is once again the radius of the photosphere, calculated as part of the shell structure. Because of light-matter interactions, the output luminosity of the supernova will not be the same as the luminosity of the photosphere, so the photospheric temperature is updated throughout the simulation as part of the convergence process in order to match the output luminosity to the requested luminosity (see [Updating Plasma and Convergence](../update_and_conv/update_and_conv.ipynb)).\n", "\n", - "Next, TARDIS calculates the initial guess for the radiative temperature in each shell. This temperature is also updated throughout the simulation based on light-matter interactions (see [Estimators and Convergence](../est_and_conv/index.ipynb)). The initial guess for $T_\\mathrm{rad}$ is calculated using Wein's Law. Wein's Law states that the temperature of a blackbody is inversely proportional to the blackbody's peak wavelength. The proportionality constant, labeled $b$, is approximately $2.898*10^{-3} m*K$. In equation form,\n", + "Next, TARDIS calculates the initial guess for the radiative temperature in each shell. This temperature is also updated throughout the simulation based on light-matter interactions (see once again [Updating Plasma and Convergence](../update_and_conv/update_and_conv.ipynb)). The initial guess for $T_\\mathrm{rad}$ is calculated using Wein's Law. Wein's Law states that the temperature of a blackbody is inversely proportional to the blackbody's peak wavelength. The proportionality constant, labeled $b$, is approximately $2.898*10^{-3} m*K$. In equation form,\n", "\n", "$$T=\\frac{b}{\\lambda_\\mathrm{peak}}.$$\n", "\n", @@ -556,7 +556,7 @@ "where $r_\\mathrm{boundary\\_inner}$ is the radius of the circle in the picture (the photosphere in terms of TARDIS). This is the geometric dilution factor.\n", "\n", "\n", - "Like the radiative temperature, the dilution factor in each cell changes throughout the simulation based on the light-matter interactions (see [Estimators and Convergence](../est_and_conv/index.ipynb)). TARDIS uses the geometric dilution factor as the initial guess for the dilution factor in each cell, plugging in the radius at the center of the cell to determine the initial dilution factor in the cell. \n", + "Like the radiative temperature, the dilution factor in each cell changes throughout the simulation based on the light-matter interactions (see [Updating Plasma and Convergence](../update_and_conv/update_and_conv.ipynb)). TARDIS uses the geometric dilution factor as the initial guess for the dilution factor in each cell, plugging in the radius at the center of the cell to determine the initial dilution factor in the cell. \n", "\n", "
\n", " \n", @@ -616,7 +616,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.12" + "version": "3.8.13" } }, "nbformat": 4, diff --git a/docs/physics/spectrum/formal_integral.rst b/docs/physics/spectrum/formal_integral.rst index 9e0c850a7eb..0d58b9c44eb 100644 --- a/docs/physics/spectrum/formal_integral.rst +++ b/docs/physics/spectrum/formal_integral.rst @@ -74,15 +74,9 @@ Our goal is to find :math:`I_{N-1}^r`, as this is the light that will exit the s Constructing the Source Function -------------------------------- -Our problem is to determine the intensity that is added to the ray at a line resonance :math:`l\rightarrow u` based on the Monte Carlo packets. Much of this calculation will rely heavily on the logic in :doc:`../est_and_conv/estimators`. +Our problem is to determine the intensity that is added to the ray at a line resonance :math:`l\rightarrow u` based on the Monte Carlo packets. -For a line transition :math:`i\rightarrow u`, a Monte Carlo packet carrying an energy :math:`\epsilon` carries that energy over a time :math:`\Delta t` (the time of simulation, see :doc:`../montecarlo/initialization` for more on how it is calculated and how it should be interpreted). And, if the packet is in a cell of volume :math:`V`, it contributes an energy density to that shell of :math:`\frac{\epsilon}{V}` over a time :math:`\Delta t`, and thus each packet represents energy density flowing at a rate of :math:`\frac{\epsilon}{\Delta t V}`. If we sum over every packet in a shell that comes into resonance with a line (i.e. that pass through a Sobolev point, see :ref:`propagation`), we will get the rate at which energy density resonates with the line. Thus, the rate at which energy density resonates with the line transition in a specific cell is :math:`\frac{1}{\Delta t V} \sum \epsilon` where the sum is taken over all packets that pass through the Sobolev point. - -This light then has a :math:`\left( 1- e^{-\tau_{iu}}\right)` probability of interacting with the line, so the rate at which energy density is absorbed into the transition :math:`i\rightarrow u` is the estimator - -.. math:: \dot{E}_{iu} = \frac{1}{\Delta t V} \left( 1- e^{-\tau_{iu}}\right) \sum \epsilon. - -This estimator is calculated similarly to the ``J_blue`` estimator, see :doc:`../est_and_conv/estimators`. Now, sum up these estimators for every level :math:`i` that can be excited to the level :math:`u` -- this will give us the rate at which energy density is added to the level :math:`u` from all line transitions that can be excited to :math:`u`: +Consider another transition :math:`i\rightarrow u`. We know (see :ref:`edotlu`) that the rate at which energy density interacts with the transition is :math:`\dot{E}_{iu}`, i.e. the ``Edotlu`` estimator for the transition :math:`i\rightarrow u`. Now, sum up these estimators for every level :math:`i` that can be excited to the level :math:`u` -- this will give us the rate at which energy density is added to the level :math:`u` from all line transitions that can be excited to :math:`u`: .. math:: \dot{E}_u = \sum_{i < u} \dot E_{iu}. @@ -100,7 +94,7 @@ With this, we can now show how the intensity along our ray is incremented. First .. math:: I_k^r = I_k^b e^{-\tau_k} + \left( 1- e^{-\tau_k}\right) S_{k}. -If there is no electron scattering, this would be it; we would have :math:`I^b_{k+1}=I^r_k` (the intensity on the left of a line would be the intensity on the right of the next line, since no interactions would take place between the locations of line resonances). However, with electron scattering, this is not the case. In between line resonances, light has a probability of :math:`e^{-\Delta \tau_e}` of not scattering, where :math:`\Delta \tau_e=\sigma_{\mathrm{T}} n_e \Delta z` is the electron scattering optical depth (see :ref:`physical-interactions`). So, we would have :math:`I^b_{k+1}=I^r_ke^{-\Delta \tau_e}`. But, light can also scatter onto the ray. The intensity is increased by the mean intensity between the two resonance points times the probability of light scattering (and thus ending up on the ray), this probability being :math:`1-e^{-\Delta \tau_e}`. For the intensity between :math:`z_k` and :math:`z_{k+1}` we use the average of the mean intensity to the right of :math:`z_k` and the mean intensity to the left of :math:`z_{k+1}`, which is :math:`\frac{1}{2}\left(J_k^r+J_{k+1}^b\right)`. Here, :math:`J^b_{k}` is calculated from the ``J_blue`` estimator for the kth line transition (see :doc:`../est_and_conv/estimators`), and then +If there is no electron scattering, this would be it; we would have :math:`I^b_{k+1}=I^r_k` (the intensity on the left of a line would be the intensity on the right of the next line, since no interactions would take place between the locations of line resonances). However, with electron scattering, this is not the case. In between line resonances, light has a probability of :math:`e^{-\Delta \tau_e}` of not scattering, where :math:`\Delta \tau_e=\sigma_{\mathrm{T}} n_e \Delta z` is the electron scattering optical depth (see :ref:`physical-interactions`). So, we would have :math:`I^b_{k+1}=I^r_ke^{-\Delta \tau_e}`. But, light can also scatter onto the ray. The intensity is increased by the mean intensity between the two resonance points times the probability of light scattering (and thus ending up on the ray), this probability being :math:`1-e^{-\Delta \tau_e}`. For the intensity between :math:`z_k` and :math:`z_{k+1}` we use the average of the mean intensity to the right of :math:`z_k` and the mean intensity to the left of :math:`z_{k+1}`, which is :math:`\frac{1}{2}\left(J_k^r+J_{k+1}^b\right)`. Here, :math:`J^b_{k}` is calculated from the ``j_blue`` estimator for the kth line transition (see :ref:`j-blue-estimator`), and then .. math:: J_k^r = J_k^b e^{-\tau_k} + \left( 1- e^{-\tau_k}\right) S_{k} diff --git a/docs/physics/est_and_conv/tardis_example.yml b/docs/physics/update_and_conv/tardis_example.yml similarity index 100% rename from docs/physics/est_and_conv/tardis_example.yml rename to docs/physics/update_and_conv/tardis_example.yml diff --git a/docs/physics/update_and_conv/update_and_conv.ipynb b/docs/physics/update_and_conv/update_and_conv.ipynb new file mode 100644 index 00000000000..c482cf448c1 --- /dev/null +++ b/docs/physics/update_and_conv/update_and_conv.ipynb @@ -0,0 +1,586 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4d53d061", + "metadata": {}, + "source": [ + "# Updating Plasma and Convergence\n", + "\n", + "As light travels through a real plasma, it has effects on the properties of the plasma due to light-matter\n", + "interactions as well as the presence of extra energy from the light. Additionally, as [previously discussed](../montecarlo/propagation.rst), properties of the plasma affect how light travels through it. This is a typical equilibrium problem. We solve for the plasma properties by finding a steady-state solution; that is, the actual plasma will be in a state such that the plasma state will not change as\n", + "light propagates through it, because the effects of the light on the plasma and the effects of the plasma on the\n", + "light are in equilibrium.\n", + "\n", + "One of TARDIS's main goals is to determine this plasma state (as we need the actual plasma properties in order to\n", + "get an accurate spectrum). This is done in an iterative process. After each [Monte Carlo iteration](../montecarlo/index.rst) (which sends light through the supernova ejecta), TARDIS uses the [Monte Carlo estimators](../montecarlo/estimators.rst)\n", + "how the propagating light affects the plasma state, after which the plasma state is updated (as will be explained below and demonstrated in the code example). We do this many times, and attempt to have the plasma state converge\n", + "to the steady-state we are looking for. In fact, all but the last Monte Carlo iteration is used for this purpose\n", + "(after which TARDIS will have the necessary plasma state for its last iteration which calculates the spectrum).\n", + "\n", + "
\n", + "\n", + "Note\n", + "\n", + "For all but the last iteration, TARDIS uses the number of Monte Carlo packets specified in the\n", + "[Monte Carlo configuration](../../io/configuration/components/montecarlo.rst) under the ``no_of_packets`` argument. This is because\n", + "a different number of packets may be necessary to calculate the spectrum as opposed to calculate the\n", + "plasma state.\n", + "\n", + "
\n", + "\n", + "After each iteration the values for radiative temperature and dilution factor are updated by calling the ``advance_state`` method on a ``Simulation`` object as will be shown below in the code example. The goal of this is to eventually have the radiative temperature and dilution factor converge to a single value so that the steady-state plasma state can be determined. To ensure that the simulation converges, TARDIS employs a convergence strategy. Currently, only one convergence strategy is available: damped convergence. This will be described in the following sections.\n", + "\n", + "
\n", + "\n", + "Note\n", + " \n", + "Unless otherwise noted, all user-supplied quantities mentioned on this page are supplied in the [convergence section of the Monte Carlo configuration](../../io/configuration/components/montecarlo.rst#damped-convergence-strategy), which will be referenced as the convergence configuration.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "8ffc4479", + "metadata": {}, + "source": [ + "## $T_\\mathrm{rad}$ and $W$\n", + "\n", + "The main way in which the plasma state is updated is by using the ``j`` and ``nu_bar`` estimators (see [here](../montecarlo/estimators.rst#j-and-bar-nu-estimators)) to update two of the key plasma inputs (see [Plasma](../setup/plasma/index.rst)): the radiative temperature $T_\\mathrm{rad}$ and dilution factor $W$. Recall that each of these quantities takes on a different value in each shell.\n", + "\n", + "Using the estimators $J$ and $\\bar \\nu$, we estimate these quantities using\n", + "\n", + "$$T_{\\mathrm{rad\\ estimated}} = \\frac{h}{k_{\\mathrm{B}}} \\frac{\\pi^4}{360 \\zeta(5)} \\frac{\\bar \\nu}{J}$$\n", + "\n", + "and\n", + "\n", + "$$W_\\mathrm{estimated} = \\frac{\\pi J}{\\sigma_{\\mathrm{R}} T_{\\mathrm{rad\\ estimated}}^4}$$\n", + "\n", + "where $h$ is Planck's constant, $k_{\\mathrm{B}}$ is Boltzmann's constant, $\\sigma_{\\mathrm{R}}$ is the Stefan–Boltzmann constant, and $\\zeta$ is the Riemann zeta function. The equation for $W$ comes from the fact that the dilution factor is the ratio of the actual mean intensity to that of a blackbody, which is $J_{\\mathrm{blackbody}}=\\frac{\\sigma_{\\mathrm{R}} T^4}{\\pi}$.\n", + "\n", + "Recall (see [Implementation](../montecarlo/estimators.rst#implementation)), however, that when TARDIS stores these estimators, it leaves out the prefactor of $\\frac{1}{4\\pi V\\Delta t}$ for computational ease. That is, $J=\\frac{1}{4\\pi V\\Delta t}\\sum_i \\varepsilon_i l_i=\\frac{1}{4\\pi V\\Delta t}*\\mathrm{real\\ j\\ estimator}$, and $\\bar \\nu=\\frac{1}{4\\pi V\\Delta t}\\sum_i \\varepsilon_i \\nu_i l_i=\\frac{1}{4\\pi V\\Delta t}*\\mathrm{real\\ nu\\_ bar\\ estimator}$. These factors are then included in our calculations; specifically, using the previous relations we have\n", + "\n", + "$$T_\\mathrm{rad\\ estimated}= \\frac{h}{k_{\\mathrm{B}}} \\frac{\\pi^4}{360 \\zeta(5)} \\frac{\\sum_i \\varepsilon_i \\nu_i l_i}{\\sum_i \\varepsilon_i l_i} = \\frac{h}{k_{\\mathrm{B}}} \\frac{\\pi^4}{360 \\zeta(5)} \\frac{\\mathrm{real\\ nu\\_ bar\\ estimator}}{\\mathrm{real\\ j\\ estimator}}$$\n", + "and\n", + "$$W_\\mathrm{estimated} = \\frac{\\sum_i \\varepsilon_i l_i}{4\\sigma_{\\mathrm{R}} T_{\\mathrm{rad\\ estimated}}^4V\\Delta t} = \\frac{\\mathrm{real\\ j\\ estimator}}{4\\sigma_{\\mathrm{R}} T_{\\mathrm{rad\\ estimated}}^4V\\Delta t}.$$\n", + "\n", + "While TARDIS can then update the plasma state using the estimated values, there is a good chance that these estimated values would “overshoot” the true value we want to converge to (for example, if the current value of the dilution factor in some cell the dilution factor is 0.4, and the true steady-state value TARDIS wants to find is 0.45, there is a good chance that the estimated value will be greater than 0.45). This could make the simulation take longer to converge or, at worst, make it so the simulation does not converge at all. To account for this, users can set (in the convergence configuration) a \"damping constant\" for both the radiative temperature ($d_{T_\\mathrm{rad}}$) and the dilution factor ($d_W$). When ``advance_state`` is called, these quantities update as follows:\n", + "\n", + "$$T_\\mathrm{rad\\ updated} = T_\\mathrm{rad\\ current} + d_{T_\\mathrm{rad}}(T_\\mathrm{rad\\ estimated}-T_\\mathrm{rad\\ current})$$\n", + " \n", + "and\n", + " \n", + "$$ W_\\mathrm{updated} = W_\\mathrm{current} + d_W(W_\\mathrm{estimated}-W_\\mathrm{current}).$$\n", + "\n", + "This means, for example, if the damping constant is 0.5, the updated value is halfway between the current value and the estimated value. If the damping constant is 0.7, the updated value is 70% of the way between the current value and the estimated value, and so on. **If the damping constant is 1, then the updated value is exactly the estimated value, and if the damping constant is zero, the value stays the same throughout the simulation and is not updated.**\n", + "\n", + "The updated $T_\\mathrm{rad}$ and $W$ are then used as inputs to the updated [plasma calculations](../setup/plasma/index.rst) which account for the effect of the Monte Carlo packets on the plasma state." + ] + }, + { + "cell_type": "markdown", + "id": "f346fb77", + "metadata": {}, + "source": [ + "## $T_\\mathrm{inner}$\n", + "\n", + "The temperature of the inner boundary, $T_\\mathrm{inner}$, plays a unique role in the simulation, as it is the primary determiner of the output luminosity. This is because the the luminosity of the inner boundary is proportional to $T_\\mathrm{inner}^4$ (see [Energy Packet Initialization](../montecarlo/initialization.ipynb)). Thus, $T_\\mathrm{inner}$ is updated throughout the simulation in order to match the output luminosity to the requested luminosity specified in the [supernova configuration](../../io/configuration/components/supernova.rst) between the bounds specified in the supernova configuration. However, there is not necessarily a quartic relationship between $T_\\mathrm{inner}$ and the output luminosity, as changing $T_\\mathrm{inner}$ also changes the frequency distribution of the initialized packets (once again see [Energy Packet Initialization](../montecarlo/initialization.ipynb)). This then affects the light-matter interactions, affecting which packets make it to the outer boundary, which also affects the output luminosity. Because of this, there is not an exact way to estimate $T_\\mathrm{inner}$. To do this estimation, we use\n", + "\n", + "$$T_\\mathrm{inner\\ estimated} = T_\\mathrm{inner\\ current} * \\left(\\frac{L_\\mathrm{output}}{L_\\mathrm{requested}}\\right)^{\\mathrm{t\\_inner\\_update\\_exponent}}$$\n", + " \n", + "where $L_\\mathrm{output}$ is the output luminosity calculated by adding up the luminosity of each packet (see [Basic Spectrum Generation](../spectrum/basic.ipynb)) between the bounds specified in the [supernova configuration](../../io/configuration/components/supernova.rst), $L_\\mathrm{requested}$ is the luminosity requested also in the supernova configuration (requested between those bounds previously mentioned), and ``t_inner_update_exponent`` is provided by the user in the convergence configuration. Note that what we are doing is \"correcting\" the previous value of the inner temperature by a factor of $\\left(\\frac{L_\\mathrm{output}}{L_\\mathrm{requested}}\\right)^{\\mathrm{t\\_inner\\_update\\_exponent}}$. Note that if $\\frac{L_\\mathrm{output}}{L_\\mathrm{requested}}$ is greater than 1, we want to lower $T_\\mathrm{inner}$ as the output luminosity is too high, and vice versa if the ratio is less than 1. Thus ``t_inner_update_exponent`` should be negative. Naively one might set ``t_inner_update_exponent=-0.25``, however as a default TARDIS uses ``t_inner_update_exponent=-0.5`` as -0.25 may undershoot the correct $T_\\mathrm{inner}$ because of its previously mentioned effects on the initial frequency distribution.\n", + "\n", + "After calculating the estimated $T_\\mathrm{inner}$, the quantity is updated using damped convergence with its own damping constant (once again set in the convergence configuration):\n", + "\n", + "$$T_\\mathrm{inner\\ updated} = T_\\mathrm{inner\\ current} + d_{T_\\mathrm{inner}}(T_\\mathrm{inner\\ estimated}-T_\\mathrm{inner\\ current}).$$\n", + "\n", + "Once again, If the damping constant is 1, then the updated value is exactly the estimated value, and if the damping constant is zero, the value stays the same throughout the simulation and is not updated.\n", + "\n", + "Additionally, because of the vast impact of $T_\\mathrm{inner}$ on the simulation, one may want to update it less frequently -- i.e. allow $W$ and $T_\\mathrm{rad}$ to reach a steady-state value for a particular $T_\\mathrm{inner}$ before updating $T_\\mathrm{inner}$. To do this, in the convergence configuration we set ``lock_t_inner_cycles``, which is the number of iterations to wait before updating $T_\\mathrm{inner}$. It is set to 1 by default, meaning $T_\\mathrm{inner}$ would be updated every iteration." + ] + }, + { + "cell_type": "markdown", + "id": "9a3bee0e", + "metadata": {}, + "source": [ + "## Convergence Information\n", + "\n", + "During the simulation, information about the how $T_\\mathrm{rad}$, $W$, and $T_\\mathrm{inner}$ are updated as well as a comparison of the total output luminosity and the requested luminosity are logged at the INFO level (see [Configuring the Logging Output for TARDIS](../../io/optional/logging_configuration.ipynb)) as shown in the code below, to give users a better idea of how the convergence process is working.\n", + "\n", + "In addition, TARDIS allows for the displaying of convergence plots, which allows users to visualize the convergence process for $T_\\mathrm{rad}$, $W$, $T_\\mathrm{inner}$, and the total luminosity of the supernova being modeled. For more information, see [Convergence Plots](../../io/visualization/convergence_plot.ipynb)." + ] + }, + { + "cell_type": "markdown", + "id": "3bc752ec", + "metadata": {}, + "source": [ + "## Convergence Criteria\n", + "\n", + "TARDIS also allows users to stop the simulation if the simulation reaches a certain level of convergence, which is checked upon the call of the ``advance_state`` method. To enable this, users must set ``stop_if_converged=True`` in the convergence configuration. Also in the configuration, the quantities ``hold_iterations``, ``threshold``, and ``fraction`` are be specified to determine convergence as follows:\n", + "\n", + "For the simulation to be considered to have converged, for ``hold_iterations`` successive iterations, the estimated values of $T_\\mathrm{rad}$, $W$, and $T_\\mathrm{inner}$ may differ from the previous value by a fraction of at most ``threshold`` in at least ``fraction`` fraction of the shells (for $T_\\mathrm{inner}$, since there is only one value, the ``fraction`` part does not apply). For example, if ``hold_iterations=3``, ``threshold=0.05`` for all three quantities, and ``fraction=0.8``, the simulation will be considered to have converged if for 3 successive iterations the estimated values of $T_\\mathrm{rad}$ and $W$ differ from the current respective values by at most 5% in at least 80% of the shells, *and* the estimated $T_\\mathrm{inner}$ differs by at most 5%. See the [convergence configuration schema](../../io/configuration/components/montecarlo.rst#damped-convergence-strategy) for default values of these quantities.\n", + "\n", + "
\n", + "\n", + "Note\n", + " \n", + "To determine convergence, we compare the estimated value, **not** the updated value (which is related to the estimated value via the damping constant), with the previous value. If $T_\\mathrm{inner}$ is locked (see the previous section), the estimated value will still be calculated so convergence can be checked as usual.\n", + "\n", + "
\n", + "\n", + "
\n", + "\n", + "Note\n", + " \n", + "``hold_iterations`` and ``fraction`` are universal quantities, i.e. they are each a single value that applies to $T_\\mathrm{rad}$ and $W$, and for ``hold_iterations`` also $T_\\mathrm{inner}$. ``threshold``, on the other hand, is supplied for each quantity separately, so for instance you could require $T_\\mathrm{rad}$ to differ by less than 1%, $W$ to differ by less than 3%, and $T_\\mathrm{inner}$ to differ by less than 5% for convergence to be reached.\n", + "\n", + "
\n", + "\n", + "\n", + "## Updating Other Quantities\n", + "\n", + "Other quantities in the plasma state can depend on other estimators. Currently, this is only implemented for the ``j_blue`` estimator: If ``radiative_rates_type`` in the [plasma configuration](../../io/configuration/components/plasma.rst) is set to ``detailed``, the `j_blues` plasma property will will be replaced with the value of the $J^b_{lu}$ [estimator](../montecarlo/estimators.rst#j-b-lu-estimator) (the raw estimator times the factor of $\\frac{ct_\\mathrm{explosion}}{4\\pi V \\Delta t}$, once again see [Implementation](../montecarlo/estimators.rst#implementation)), which would then affect other plasma properties that depend on the `j_blues` values (see [Plasma](../setup/plasma/index.rst)). Otherwise, the `j_blues` values in the plasma are calculated as they typically are in the plasma calculations, and the $J^b_{lu}$ estimator is only used for the [formal integral](../spectrum/formal_integral.rst). Even in the former case, while the estimator does contribute to the updating of the plasma when the ``advance_state`` method is called, it does **not** contribute to the determination of if the simulation has converged, and it does **not** show up in the convergence logging or convergence plots." + ] + }, + { + "cell_type": "markdown", + "id": "127746c2", + "metadata": {}, + "source": [ + "## Code Example\n", + "\n", + "We now show a detailed example of how the plasma is updated using the estimators after a Monte Carlo iteration. First, we import the necessary packages and set up a simulation (see [Setting Up the Simulation](../setup/index.rst)):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a8bd905", + "metadata": {}, + "outputs": [], + "source": [ + "from tardis.io.config_reader import Configuration\n", + "from tardis.simulation import Simulation\n", + "from tardis.io.atom_data.util import download_atom_data\n", + "import numpy as np\n", + "from scipy.special import zeta\n", + "\n", + "from astropy import units as u, constants as const\n", + "\n", + "# We download the atomic data needed to run this notebook\n", + "download_atom_data('kurucz_cd23_chianti_H_He')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "259a7aa9", + "metadata": {}, + "outputs": [], + "source": [ + "tardis_config = Configuration.from_yaml('tardis_example.yml')\n", + "\n", + "# We manually put in the damping constants and t_inner_update_exponent for\n", + "# illustrative purposes:\n", + "d_t_rad = 0.5\n", + "d_w = 0.3\n", + "d_t_inner = 0.7\n", + "t_inner_update_exponent = -0.5\n", + "\n", + "# We set the above values into the configuration:\n", + "tardis_config.montecarlo.convergence_strategy.t_rad.damping_constant = d_t_rad\n", + "tardis_config.montecarlo.convergence_strategy.w.damping_constant = d_w\n", + "tardis_config.montecarlo.convergence_strategy.t_inner.damping_constant = d_t_inner\n", + "tardis_config.montecarlo.convergence_strategy.t_inner_update_exponent = t_inner_update_exponent" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb37402c", + "metadata": {}, + "outputs": [], + "source": [ + "sim = Simulation.from_config(tardis_config)\n", + "\n", + "model = sim.model\n", + "plasma = sim.plasma\n", + "runner = sim.runner" + ] + }, + { + "cell_type": "markdown", + "id": "cf13d946", + "metadata": {}, + "source": [ + "We show the initial radiative temperature and dilution factor in each shell, the initial inner boundary temperature, and the initial electron densities in each shell:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3f43f93", + "metadata": {}, + "outputs": [], + "source": [ + "model.t_rad" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90bb6147", + "metadata": {}, + "outputs": [], + "source": [ + "model.w" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66209ff7", + "metadata": {}, + "outputs": [], + "source": [ + "model.t_inner" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7b2ae7b4", + "metadata": {}, + "outputs": [], + "source": [ + "plasma.electron_densities" + ] + }, + { + "cell_type": "markdown", + "id": "a248ff2c", + "metadata": {}, + "source": [ + "We set the number of packets and we run one iteration of the Monte Carlo simulation:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9f7eb44f", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "N_packets = 10000\n", + "\n", + "# Using the commented out code below, we can also get the number of packets\n", + "# from the configuration -- try it out:\n", + "#N_packets = tardis_config.montecarlo.no_of_packets\n", + "\n", + "sim.iterate(N_packets)" + ] + }, + { + "cell_type": "markdown", + "id": "d3d0074d", + "metadata": {}, + "source": [ + "We now show the values of the $J$ and $\\bar \\nu$ estimators, attaching their proper units (note that these are the raw estimators, and the factors of $\\frac{1}{4\\pi V \\Delta t}$ etc are not included):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21500ea8", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "j_estimator = runner.j_estimator * (u.erg * u.cm) \n", + "j_estimator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fac91ee2", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "nu_bar_estimator = runner.nu_bar_estimator * (u.erg * u.cm * u.Hz)\n", + "nu_bar_estimator" + ] + }, + { + "cell_type": "markdown", + "id": "69c178d2", + "metadata": {}, + "source": [ + "We show the values of $J$ and $\\bar \\nu$ including the prefactor:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9500e34b", + "metadata": {}, + "outputs": [], + "source": [ + "V = model.volume\n", + "Delta_t = runner.calculate_time_of_simulation(model)\n", + "prefactor = 1 / (4 * np.pi * V * Delta_t)\n", + "J = prefactor * j_estimator\n", + "J" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5720987", + "metadata": {}, + "outputs": [], + "source": [ + "nu_bar = prefactor * nu_bar_estimator\n", + "nu_bar" + ] + }, + { + "cell_type": "markdown", + "id": "28d69ca6", + "metadata": {}, + "source": [ + "As mentioned [here](../montecarlo/estimators.rst#j-and-bar-nu-estimators), $\\bar \\nu$ is not truly the mean frequency (as you can see by units). We show the true mean frequency, which will, as expected, be in units of Hz." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "691cf2f6", + "metadata": {}, + "outputs": [], + "source": [ + "nu_bar/J" + ] + }, + { + "cell_type": "markdown", + "id": "bcbf40d4", + "metadata": {}, + "source": [ + "We show the calculations of the estimated and updated $T_\\mathrm{rad}$ and $W$. Note that the ``decompose()`` method is used to simplify the units:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50c639ea", + "metadata": {}, + "outputs": [], + "source": [ + "t_rad_estimated = ( (const.h.cgs / const.k_B.cgs) \n", + " * (np.pi**4 / (360 * zeta(5, 1)))\n", + " * (nu_bar_estimator / j_estimator) ).decompose()\n", + "t_rad_estimated" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b34674df", + "metadata": {}, + "outputs": [], + "source": [ + "t_rad_updated = model.t_rad + d_t_rad * (t_rad_estimated - model.t_rad)\n", + "t_rad_updated" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b48d8ada", + "metadata": {}, + "outputs": [], + "source": [ + "w_estimated = ( j_estimator / (4 * const.sigma_sb.cgs * t_rad_estimated**4 * V * Delta_t) ).decompose()\n", + "w_estimated" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22895e09", + "metadata": {}, + "outputs": [], + "source": [ + "w_updated = model.w + d_w * (w_estimated - model.w)\n", + "w_updated" + ] + }, + { + "cell_type": "markdown", + "id": "7bfefce9", + "metadata": {}, + "source": [ + "We show the output and requested luminosities, and use these to calculate the estimated and updated $T_\\mathrm{inner}$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "795dceb2", + "metadata": {}, + "outputs": [], + "source": [ + "# The output luminosity is calculated between the bounds specified below\n", + "nu_lower = 0\n", + "nu_upper = np.inf\n", + "\n", + "# Using the commented out code below, we can also get the frequency bounds\n", + "# from the configuration -- try it out (note we convert between wavelength\n", + "# and frequency, and thus the order switches):\n", + "#nu_lower = tardis_config.supernova.luminosity_wavelength_end.to(u.Hz, u.spectral)\n", + "#nu_upper = tardis_config.supernova.luminosity_wavelength_start.to(u.Hz, u.spectral)\n", + "\n", + "L_output = runner.calculate_emitted_luminosity(0,np.inf)\n", + "L_output" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7028da2", + "metadata": {}, + "outputs": [], + "source": [ + "L_requested = sim.luminosity_requested\n", + "L_requested" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a14111a", + "metadata": {}, + "outputs": [], + "source": [ + "t_inner_estimated = model.t_inner * (L_output / L_requested)**t_inner_update_exponent\n", + "t_inner_estimated" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8bd1062f", + "metadata": {}, + "outputs": [], + "source": [ + "t_inner_updated = model.t_inner + d_t_inner * (t_inner_estimated - model.t_inner)\n", + "t_inner_updated" + ] + }, + { + "cell_type": "markdown", + "id": "43b31d2f", + "metadata": {}, + "source": [ + "We now advance the state of the simulation based on the estimators. This will also display a summary of the updated values of $T_\\mathrm{rad}$ and $W$. Additionally, the ``advance_state()`` method checks for convergence and returns the convergence status as shown below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9676b22b", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "sim.advance_state()" + ] + }, + { + "cell_type": "markdown", + "id": "9619477d", + "metadata": {}, + "source": [ + "Finally, we show the full updated $T_\\mathrm{rad}$, $W$, and $T_\\mathrm{inner}$, as well as the updated electron densities which were updated along with the rest of the plasma based on the new values for $T_\\mathrm{rad}$, $W$, and $T_\\mathrm{inner}$. Compare these with our calculations above and with the initial values at the beginning of the code example!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe0b9f40", + "metadata": {}, + "outputs": [], + "source": [ + "model.t_rad" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8de89bb2", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "model.w" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f78c4fe", + "metadata": {}, + "outputs": [], + "source": [ + "model.t_inner" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fcaf08ff", + "metadata": {}, + "outputs": [], + "source": [ + "plasma.electron_densities" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.13" + }, + "vscode": { + "interpreter": { + "hash": "c20eb03b6e6c860a5223c79131a5cfd0b7e61a9aa5d08c2b1400663bc6cc4e6d" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}