diff --git a/docs/io/optional/custom_source.ipynb b/docs/io/optional/custom_source.ipynb index 68aaf50a46d..a2da649147a 100644 --- a/docs/io/optional/custom_source.ipynb +++ b/docs/io/optional/custom_source.ipynb @@ -4,7 +4,24 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Running TARDIS with a custom packet source" + "# Running TARDIS with a Custom Packet Source" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By default, TARDIS generates energy packets using its `BasePacketSource` class, which models the photosphere of the supernova as a perfect blackbody (see [Energy Packet Initialization](../../physics/montecarlo/initialization.ipynb)). However, users may create their own packet source, as will be shown in this notebook. In order to do this, a user must create a class (that can but does not have to inherit from `BasePacketSource`) which contains a method `create_packets` that takes in (in this order):\n", + "- the photospheric temperature\n", + "- the number of packets\n", + "- a random number generator\n", + "- the photospheric radius\n", + "\n", + "and returns arrays of the radii, frequencies, propogation directions, and energies of the ensemble of packets (once again see [Energy Packet Initialization](../../physics/montecarlo/initialization.ipynb) for more information). To use your packet source in a run of TARDIS, you must pass an instance of your class into the `run_tardis` function under the `packet_source` keyword argument.\n", + "\n", + ".. note:: In both the `BasePacketSource` class and in the example here, all packets are generated at the same radius. This need not be true in general (although the `create_packets` method should only take in a single radius as its argument).\n", + "\n", + "We show an example of how a custom packet source is used:" ] }, { @@ -13,6 +30,7 @@ "metadata": {}, "outputs": [], "source": [ + "# Import necessary packages\n", "import numpy as np\n", "from tardis import constants as const\n", "from astropy import units as u\n", @@ -28,22 +46,17 @@ "metadata": {}, "outputs": [], "source": [ + "# Download the atomic data used for a run of TARDIS\n", "download_atom_data('kurucz_cd23_chianti_H_He')" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Custom packet source class that is derived from BasePacketSource. The method create_packets (which returns ```radii, nus, mus, energies```) has to be defined." - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ + "# Create a packet source class that contains a create_packets method\n", "class TruncBlackbodySource(BasePacketSource):\n", " \"\"\"\n", " Custom inner boundary source class to replace the Blackbody source\n", @@ -112,11 +125,19 @@ "metadata": {}, "outputs": [], "source": [ + "# Call an instance of the packet source class\n", "packet_source = TruncBlackbodySource(\n", " 53253, truncation_wavelength=2000\n", ")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now run TARDIS both with and without our custom packet source, and we compare the results:" + ] + }, { "cell_type": "code", "execution_count": null, @@ -137,10 +158,10 @@ "%matplotlib inline\n", "plt.plot(mdl.runner.spectrum_virtual.wavelength,\n", " mdl.runner.spectrum_virtual.luminosity_density_lambda,\n", - " color='red', label='truncated blackbody')\n", + " color='red', label='truncated blackbody (custom packet source)')\n", "plt.plot(mdl_norm.runner.spectrum_virtual.wavelength,\n", " mdl_norm.runner.spectrum_virtual.luminosity_density_lambda,\n", - " color='blue', label='normal blackbody')\n", + " color='blue', label='normal blackbody (default packet source)')\n", "plt.xlabel('$\\lambda [\\AA]$')\n", "plt.ylabel('$L_\\lambda$ [erg/s/$\\AA$]')\n", "plt.xlim(500, 10000)\n", @@ -173,4 +194,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/docs/physics/montecarlo/basicprinciples.rst b/docs/physics/montecarlo/basicprinciples.rst index 62703df12f6..8c18c790376 100644 --- a/docs/physics/montecarlo/basicprinciples.rst +++ b/docs/physics/montecarlo/basicprinciples.rst @@ -20,7 +20,7 @@ through and interact with the ambient material. A sufficiently large number of representative "machine photons" are considered and their propagation history solved in a stochastic process. The initial properties of these photons are randomly (in a probabilistic sense) assigned in accordance with the macroscopic -properties of the radiation field (see :doc:`Discretization `) +properties of the radiation field (see :ref:`Energy Packets `) and in a similar manner the decisions about when, where and how the machine photons interact with the surrounding material are made (see :ref:`Propagation `). If this process is repeated for a large enough number of machine diff --git a/docs/physics/montecarlo/discretization.rst b/docs/physics/montecarlo/discretization.rst deleted file mode 100644 index f7d0b97db19..00000000000 --- a/docs/physics/montecarlo/discretization.rst +++ /dev/null @@ -1,35 +0,0 @@ -********************************************* -Monte Carlo Discretization --- Energy Packets -********************************************* - -While it is instructive to think about tracking the propagation history of -photons when illustrating the basic idea behind Monte Carlo radiative transfer -techniques, there are important numerical reasons for using a different -discretization scheme. Instead of thinking in the photon picture, it brings -significant advantages to follow the idea of :cite:`Abbott1985` and -:cite:`Lucy1999` and consider parcels of radiant energy as the fundamental -building blocks of the Monte Carlo calculation. These basic Monte Carlo quanta -are commonly referred to as "energy packets" or simply "packets". - -During a Monte Carlo calculation, a large number of packets, all with a certain -energy :math:`\varepsilon`, are created. In addition, each packet is associated -with a frequency. These assignments are performed in a manner which ensures -that the ensemble of packets represents the spectral energy distribution of the -radiation field (see :ref:`Propagation `). - -During the simulation, the energy of the packet remains constant in the local -co-moving frame (see :ref:`Reference Frames ` for -details about the lab and co-moving frames). This naturally ensures energy -conservation and constitutes the main advantage of this discretization scheme. -There is one side effect of this so-called indivisible packet energy scheme -which often causes confusion: Even during radiation-matter interactions the -packet energy is conserved in the co-moving frame (see :doc:`Propagation -`). However, the frequency associated with a packet may change -(e.g. during non-resonant line interactions). As a consequence, packets may -represent a varying number of real photons during their lifetime. - -.. note:: - The indivisible energy packet scheme does not require that all packets have - the same energy. This is just a convenient and simple choice adopted in - TARDIS. - diff --git a/docs/physics/montecarlo/index.rst b/docs/physics/montecarlo/index.rst index d70b011f166..30ffbeb7fb2 100644 --- a/docs/physics/montecarlo/index.rst +++ b/docs/physics/montecarlo/index.rst @@ -13,7 +13,7 @@ can also be found in various papers by L. Lucy and in the main TARDIS publicatio .. toctree:: basicprinciples - discretization + initialization propagation lineinteraction estimators diff --git a/docs/physics/montecarlo/initialization.ipynb b/docs/physics/montecarlo/initialization.ipynb new file mode 100644 index 00000000000..7957c96cd13 --- /dev/null +++ b/docs/physics/montecarlo/initialization.ipynb @@ -0,0 +1,332 @@ +{ + "cells": [ + { + "cell_type": "raw", + "id": "8633baaa", + "metadata": { + "raw_mimetype": "text/restructuredtext" + }, + "source": [ + ".. _initialization:\n", + "\n", + "****************************\n", + "Energy Packet Initialization\n", + "****************************" + ] + }, + { + "cell_type": "markdown", + "id": "6c0dbe0a", + "metadata": {}, + "source": [ + "While it is instructive to think about tracking the propagation history of\n", + "photons when illustrating the basic idea behind Monte Carlo radiative transfer\n", + "techniques, there are important numerical reasons for using a different\n", + "discretization scheme. Instead of thinking in the photon picture, it brings\n", + "significant advantages to follow the idea of [] and\n", + "[] and consider parcels of radiant energy as the fundamental\n", + "building blocks of the Monte Carlo calculation. These basic Monte Carlo quanta\n", + "are commonly referred to as \"energy packets\" or simply \"packets\", and are composed of many photons with the same frequency.\n", + "\n", + "During a Monte Carlo calculation, $N$ (a large number) packets, all with a certain\n", + "energy $\\varepsilon$, are created at the inner boundary of the computational domain (as will be discussed in [Packet Propagation](propagation.rst)) known as the photosphere. Currently, the photosphere is modeled as a spherical blackbody with a radius $R_\\mathrm{phot}$ and temperature $T_\\mathrm{phot}$. In TARDIS, all packets are assigned identical energies, and the total energy of the packets is 1 erg (and thus each packet has an energy of $\\frac{1}{N}$ ergs).\n", + "\n", + ".. note:: The indivisible energy packet scheme does not require that all packets have the same energy. This is just a convenient and simple choice adopted in TARDIS.\n", + "\n", + "Since the photosphere is modeled as a blackbody, its total luminosity $L$ (recall that luminosity is energy emitted divided by the time in which it is emitted) is\n", + "$$L=\\frac{N\\varepsilon}{\\Delta t}=4 \\pi R_{\\mathrm{phot}}^2 \\sigma_{\\mathrm{R}} T_{\\mathrm{phot}}^4$$\n", + "where $\\sigma_\\mathrm{R}$ is the Stefan-Boltzmann constant and $\\Delta t$ is the physical duration of the simulation. In order to make this relationship hold (remembering that $N\\varepsilon = 1$ erg), we use\n", + "$$\\Delta t = \\frac{1}{L}=\\frac{1}{4 \\pi R_{\\mathrm{phot}}^2 \\sigma_{\\mathrm{R}} T_{\\mathrm{phot}}^4}.$$\n", + "\n", + "During packet initialization, each packet is assigned an initial propagation direction $\\mu$ which is the cosine of the angle $\\theta$ which the packet's path makes with the radial direction. Using a pseudo-random number generator which generates numbers $z$ uniformly distributed on the interval $[0,1]$, the propagation direction is determined according to\n", + "$$\\mu = \\sqrt{z}.$$\n", + "This sampling is demonstrated in the code below.\n", + "\n", + "Finally, each packet is assigned an initial frequency (or more precisely, the initial frequency of its consitiuent photons). Note that since each packet has the same energy, each packet will represent a different number of real photons. The sampling on packet frequencies is more involved than that of the propagation direction, as it involves sampling the Planck distribution (see below). TARDIS uses the technique described in [] and summarized in [] for this purpose.\n", + "\n", + "During the simulation, the energy of the packet remains constant in the local\n", + "co-moving frame (see [Reference Frames](propagation.rst#reference-frames)). This naturally ensures energy\n", + "conservation and constitutes the main advantage of this discretization scheme. **However, while the energy of the packets is conserved in the co-moving frame, the frequency of the constituent photons (in the local co-moving frame) may vary over the course of the simulation. Thus, a packet may represent several different numbers of real photons throughout their lifetimes.**\n", + "\n", + "We now demonstrate the TARDIS packet initialization framework:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "426325e5", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from tardis.montecarlo.packet_source import BlackBodySimpleSource\n", + "from astropy import units as u\n", + "from tardis import constants as const\n", + "import matplotlib.pyplot as plt\n", + "\n", + "#The random number generator that will be used\n", + "rng = np.random.default_rng()" + ] + }, + { + "cell_type": "markdown", + "id": "4ae02998", + "metadata": {}, + "source": [ + "The following cell contains values that you can change to see how it affects the spectrum:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc34bf33", + "metadata": {}, + "outputs": [], + "source": [ + "# Seed for the pseudo-random number generator\n", + "seed = 1\n", + "\n", + "# Radius of the supernova's photosphere in cm\n", + "r_phot = 1e15 * u.cm\n", + "\n", + "# Number of packets generated\n", + "n_packets = 40000" + ] + }, + { + "cell_type": "markdown", + "id": "450faf76", + "metadata": {}, + "source": [ + "You can either set a temperatature of the photosphere, which will determine its luminosity; or you can set the luminosity of the photosphere, which will determine its temperature." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f2c15a9", + "metadata": {}, + "outputs": [], + "source": [ + "# Temperature in K\n", + "temperature = 10000 * u.K\n", + "\n", + "luminosity = 4 * np.pi * (r_phot**2) * const.sigma_sb * (temperature**4)\n", + "\n", + "# Makes sure the luminosity is given in erg/s\n", + "luminosity = luminosity.to('erg/s')\n", + "\n", + "print('Luminosity:', luminosity)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3fb3ca8c", + "metadata": {}, + "outputs": [], + "source": [ + "# Luminosity in erg/s\n", + "luminosity = 7e42 * u.erg / u.s\n", + "\n", + "temperature = (luminosity / (4 * np.pi * const.sigma_sb))**0.25 / np.sqrt(r_phot)\n", + "\n", + "# Makes sure the termperature is given in K\n", + "temperature = temperature.to('K')\n", + "\n", + "print('Temperature:', temperature)" + ] + }, + { + "cell_type": "markdown", + "id": "516633c5", + "metadata": {}, + "source": [ + "We now generate the ensemble of packets. The array of packet energies and radii are also shown." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "925e9e1b", + "metadata": {}, + "outputs": [], + "source": [ + "# We define our packet source\n", + "packet_source = BlackBodySimpleSource(seed)\n", + "\n", + "radii, nus, mus, energies = packet_source.create_packets(\n", + " temperature.value, \n", + " n_packets, \n", + " rng, \n", + " r_phot)\n", + "\n", + "# Sets the energies in units of ergs\n", + "energies *= u.erg\n", + "\n", + "# Sets the frequencies in units of Hz\n", + "nus *= u.Hz\n", + "\n", + "print('Energies:', energies)\n", + "print('Radii:', radii)" + ] + }, + { + "cell_type": "markdown", + "id": "16936bce", + "metadata": {}, + "source": [ + "We set the timespan of the simulation so that each packet contributes the appropriate luminosity to the spectrum." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fed35f47", + "metadata": {}, + "outputs": [], + "source": [ + "# Time of simulation\n", + "t_simulation = 1 * u.erg / luminosity\n", + "print('Time of simulation:', t_simulation)\n", + "\n", + "# Array of luminosity contribution by each packet\n", + "lumin_per_packet = energies / t_simulation\n", + "print('Luminosity per packet:', lumin_per_packet)" + ] + }, + { + "cell_type": "markdown", + "id": "0d839222", + "metadata": {}, + "source": [ + "We define important constants, and for comparison's sake, we code the Planck distribution function\n", + "$$L_\\nu (\\nu)=\\frac{8\\pi R_\\mathrm{phot}^2 h\\nu^3}{c^2}\\frac{1}{\\exp\\left(\\frac{h\\nu}{k_BT_\\mathrm{phot}}\\right)-1}$$\n", + "where $L_\\nu$ is the luminosity density with respect to frequency, $\\nu$ is frequency, $h$ is Planck's constant, $c$ is the speed of light, and $k_B$ is Boltzmann's constant:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "916a5e22", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "h = const.h.cgs\n", + "c2 = const.c.cgs**2\n", + "kB = const.k_B.cgs\n", + "\n", + "def planck_function(nu):\n", + " return 8 * np.pi**2 * r_phot**2 * h * nu**3 / (c2 * (np.exp(h * nu / (kB * temperature)) - 1))" + ] + }, + { + "cell_type": "markdown", + "id": "78230177", + "metadata": {}, + "source": [ + "We plot the Planck distribution and a histogram of the generated packet distribution:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "913fcdbb", + "metadata": {}, + "outputs": [], + "source": [ + "# We set important quantites for making our histogram\n", + "bins = 200\n", + "nus_planck = np.linspace(min(nus), max(nus), bins)\n", + "bin_width = nus_planck[1] - nus_planck[0]\n", + "\n", + "# In the histogram plot below, the weights argument is used \n", + "# to make sure our plotted spectrum has the correct y-axis scale\n", + "plt.hist(nus.value,\n", + " bins=bins,\n", + " weights=lumin_per_packet/bin_width)\n", + "\n", + "# We plot the planck function for comparison\n", + "plt.plot(nus_planck, planck_function(nus_planck))\n", + "\n", + "plt.xlabel('Frequency (Hz)')\n", + "plt.ylabel('Luminosity density w.r.t. frequency (erg/s/Hz)')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "ad4f0f0e", + "metadata": {}, + "source": [ + "We finally plot the generated $\\mu$ density distribution, followed by the generated $\\theta=\\arccos (\\mu)$ density distribution, compared with the respective curves $\\rho = 2\\mu$ and $\\rho = \\sin(2\\theta)$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae4c97de", + "metadata": {}, + "outputs": [], + "source": [ + "x = np.linspace(0, 1, 1000)\n", + "\n", + "plt.hist(mus, bins=bins, density=True)\n", + "plt.plot(x, 2*x)\n", + "plt.xlabel('Propagation direction')\n", + "plt.ylabel('Probability density')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41daa433", + "metadata": {}, + "outputs": [], + "source": [ + "thetas = np.linspace(0, np.pi/2, 1000)\n", + "\n", + "plt.hist(np.arccos(mus), bins=bins, density=True)\n", + "plt.plot(thetas, np.sin(2*thetas))\n", + "plt.xlabel('Angle with normal (rad)')\n", + "plt.ylabel('Probability density')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "2f661592", + "metadata": {}, + "source": [ + "## Custom Packet Source\n", + "\n", + "TARDIS allows for the user to input a custom function that generates energy packets instead of the basic blackbody source described here. See [Running TARDIS with a Custom Packet Source](../../io/optional/custom_source.ipynb) for more information." + ] + } + ], + "metadata": { + "celltoolbar": "Raw Cell Format", + "kernelspec": { + "display_name": "Python 3", + "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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/physics/montecarlo/lineinteraction.rst b/docs/physics/montecarlo/lineinteraction.rst index c6d03148e27..3d32cb260a0 100644 --- a/docs/physics/montecarlo/lineinteraction.rst +++ b/docs/physics/montecarlo/lineinteraction.rst @@ -8,7 +8,7 @@ TARDIS currently offers different ways to handle line interactions, which may be activated via the YAML configuration file. Independently of the chosen treatment, a number of steps are always carried out when a Monte Carlo packet performs a line interaction. Since TARDIS adopts the indivisible energy packet -formalism (see :doc:`Discretization `), the packet will have the +formalism (see :ref:`Energy Packets `), the packet will have the same energy in the co-moving frame after (f for final) the line interaction as before (i for initial). Thus, after accounting for the frame transformations, diff --git a/docs/physics/montecarlo/propagation.rst b/docs/physics/montecarlo/propagation.rst index 907ca40f3c1..1f81317fb86 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 -initialised, it is launched and may then perform interactions with the +initialised (see :ref:`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,35 +14,6 @@ spectrum, the main product of a TARDIS calculation. The different spectral features are simply a combined product of the changes in the packet properties induced in the radiation-matter interactions. -Initialization -============== - -During each TARDIS Monte Carlo simulation cycle, a large number :math:`N` of -Monte Carlo packets are initialised at the lower boundary of the computational domain -(i.e. the photosphere). Since the inner boundary is currently treated as a -black-body in TARDIS, :math:`N` packets with energies - -.. math:: - \varepsilon = \frac{4 \pi R_{\mathrm{phot}}^2 \sigma_{\mathrm{R}} T_{\mathrm{phot}}^4 \Delta t}{N} - -are initialised (the black body temperature :math:`T_{\mathrm{phot}}`, the -photospheric radius :math:`R_{\mathrm{phot}}`, the Stefan-Boltzmann constant -:math:`\sigma_{\mathrm{R}}` and the physical duration of the simulation -:math:`\Delta t` appear here). To commence the packet propagation, each packet -is assigned an initial propagation direction (note that propagation direction is defined as :math:`\mu = \cos -\theta` with :math:`\theta` being the angle the photon path makes with the -radial direction) - -.. math:: - \mu = \sqrt{z} - -and an initial frequency :math:`\nu` in random number experiments, using a -random number generator which provides uniformly distributed numbers :math:`z` -on the interval :math:`[0,1]`. The frequency assignment is more involved than -selecting an initial propagation direction, since the Planck function has to be -sampled. TARDIS uses the technique described in :cite:`Carter1975` and -summarized in :cite:`Bjorkman2001` for this purpose. - .. _expansion: Model of Supernova Domain