{
"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
}