Rates matrix solver

docs/physics/plasma/equilibrium/chianti_solver.ipynb

@@ -0,0 +1,127 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "\n",
+    "os.environ['XUVTOP'] = \"~/chianti\"\n",
+    "\n",
+    "try:\n",
+    "    import ChiantiPy.core as ch\n",
+    "    from ChiantiPy.tools.io import versionRead\n",
+    "\n",
+    "except ImportError:\n",
+    "    # Shamefully copied from their GitHub source:\n",
+    "    import chianti.core as ch\n",
+    "    def versionRead():\n",
+    "        \"\"\"\n",
+    "        Read the version number of the CHIANTI database\n",
+    "        \"\"\"\n",
+    "        xuvtop = os.environ['XUVTOP']\n",
+    "        vFileName = os.path.join(xuvtop, 'VERSION')\n",
+    "        vFile = open(vFileName)\n",
+    "        versionStr = vFile.readline()\n",
+    "        vFile.close()\n",
+    "        return versionStr.strip()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "compare_shell = 0\n",
+    "species = (14, 1)\n",
+    "species_string = 'si_2'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t_electron = 9967.4884"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "electron_density = 1e6"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chianti_species = ch.ion(species_string, temperature=t_electron, eDensity=electron_density)\n",
+    "chianti_species.populate()\n",
+    "chianti_population = chianti_species.Population['population']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[9.99999989e-01, 1.06967533e-08, 8.19021735e-17, 1.63935872e-16,\n",
+       "        3.82270728e-16, 8.00932946e-18, 1.60440977e-17, 1.62648689e-17,\n",
+       "        2.44571411e-17, 9.10400264e-17, 3.26639365e-18, 6.54275441e-18,\n",
+       "        7.80689247e-18, 1.17291171e-17, 7.51999263e-18, 1.01449887e-17,\n",
+       "        6.56279954e-17, 2.63927140e-18, 5.28562769e-18, 6.84869862e-18,\n",
+       "        1.02728743e-17, 6.36689557e-18, 8.48045348e-18, 5.33290944e-18,\n",
+       "        6.25920357e-18]])"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "chianti_population"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "carsus",
+   "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.12.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

docs/physics/plasma/detailed_balance/comparison.ipynb → docs/physics/plasma/equilibrium/comparison.ipynb

@@ -38,7 +38,7 @@
     "from astropy import units as u\n",
     "\n",
     "from tardis.io.atom_data import AtomData\n",
-    "from tardis.plasma.detailed_balance.rates import (\n",
+    "from tardis.plasma.equilibrium.rates import (\n",
     "    RadiativeRatesSolver,\n",
     "    ThermalCollisionalRateSolver,\n",
     "    UpsilonRegemorterSolver,\n",

docs/physics/plasma/equilibrium/tardis_solver.ipynb

@@ -0,0 +1,227 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Comparison between TARDIS and ChiantiPy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "\n",
+    "import astropy.units as u\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "from tardis.io.atom_data import AtomData\n",
+    "from tardis.io.configuration.config_reader import Configuration\n",
+    "from tardis.model.base import SimulationState\n",
+    "from tardis.plasma.equilibrium.rates import (\n",
+    "    RadiativeRatesSolver,\n",
+    "    ThermalCollisionalRateSolver,\n",
+    ")\n",
+    "from tardis.plasma.radiation_field import (\n",
+    "    DilutePlanckianRadiationField,\n",
+    ")\n",
+    "\n",
+    "home = str(Path('~').expanduser())\n",
+    "\n",
+    "config = Configuration.from_yaml(home+\"/tardis/tardis/plasma/tests/data/plasma_base_test_config.yml\")\n",
+    "\n",
+    "ion_slice = (2, 0, slice(None), slice(None))\n",
+    "\n",
+    "\n",
+    "def get_radiative_rate_solver(radiative_transitions):\n",
+    "    rad_rate_solver = RadiativeRatesSolver(radiative_transitions)\n",
+    "    return rad_rate_solver\n",
+    "\n",
+    "def get_chianti_collisional_rate_solver(atom_data, radiative_transitions,):\n",
+    "    col_strength_temperatures = atom_data.collision_data_temperatures\n",
+    "    col_strengths = atom_data.collision_data.loc[ion_slice, :]\n",
+    "    collisional_rate_solver = ThermalCollisionalRateSolver(atom_data.levels, radiative_transitions, col_strength_temperatures, col_strengths, 'chianti', \"none\")\n",
+    "    return collisional_rate_solver\n",
+    "\n",
+    "chianti_atom_data = AtomData.from_hdf(home+'/carsus/docs/kurucz_cd23_chianti_H_Si.h5')\n",
+    "\n",
+    "chianti_radiative_transitions = chianti_atom_data.lines.loc[ion_slice, :]\n",
+    "\n",
+    "chianti_sim_state = SimulationState.from_config(config, atom_data=chianti_atom_data)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from tardis.plasma.electron_energy_distribution import (\n",
+    "    ThermalElectronEnergyDistribution,\n",
+    ")\n",
+    "\n",
+    "chianti_radiative_rate_solver = get_radiative_rate_solver(chianti_radiative_transitions)\n",
+    "\n",
+    "chianti_collisional_rate_solver = get_chianti_collisional_rate_solver(chianti_atom_data, chianti_radiative_transitions)\n",
+    "\n",
+    "temperature = chianti_sim_state.t_radiative\n",
+    "\n",
+    "# need a better way to make a custom radiation field where the intensity is zero\n",
+    "# in specific locations as desired\n",
+    "rad_field = DilutePlanckianRadiationField(chianti_sim_state.t_radiative, dilution_factor=np.zeros_like(chianti_sim_state.t_radiative))\n",
+    "electron_dist = ThermalElectronEnergyDistribution(0, chianti_sim_state.t_radiative, 1e6 * u.g/u.cm**3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rate_solvers = [(chianti_radiative_rate_solver, \"radiative\"), (chianti_collisional_rate_solver, \"electron\")]\n",
+    "\n",
+    "lte_rate_solvers = [(chianti_collisional_rate_solver, \"electron\")]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rad_rates = chianti_radiative_rate_solver.solve(rad_field)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "col_rates = chianti_collisional_rate_solver.solve(electron_dist.temperature)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from tardis.plasma.equilibrium.level_populations import LevelPopulationSolver\n",
+    "from tardis.plasma.equilibrium.rate_matrix import RateMatrix\n",
+    "\n",
+    "rate_matrix_solver = RateMatrix(rate_solvers, chianti_atom_data.levels)\n",
+    "\n",
+    "rate_matrix = rate_matrix_solver.solve(rad_field, electron_dist)\n",
+    "\n",
+    "lte_rate_matrix = RateMatrix(lte_rate_solvers, chianti_atom_data.levels).solve(rad_field, electron_dist)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "solver = LevelPopulationSolver(rate_matrix, chianti_atom_data.levels)\n",
+    "\n",
+    "level_pops = solver.solve()\n",
+    "\n",
+    "lte_level_pops = LevelPopulationSolver(lte_rate_matrix, chianti_atom_data.levels).solve()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chianti_levels = [3.42881304e-01, 6.57118589e-01, 4.18696636e-09, 5.67602616e-08,\n",
+    "        4.63312038e-08, 1.20482595e-11, 1.97010580e-11, 3.48183541e-15,\n",
+    "        6.88699562e-16, 2.78695447e-16, 4.20834423e-16, 1.25473722e-15,\n",
+    "        2.49671361e-15, 8.92884539e-17, 1.65423100e-16, 2.11333508e-17,\n",
+    "        1.11695201e-17, 1.69708563e-17, 2.70670610e-17, 5.36487601e-17,\n",
+    "        8.72033865e-18, 1.17118808e-17, 3.44693669e-16, 5.11955673e-16,\n",
+    "        2.81010802e-18, 3.50318771e-18, 5.25951995e-18, 3.28532989e-18,\n",
+    "        6.57049588e-18]\n",
+    "\n",
+    "chianti_levels_he_1 = [1.00000000e+00, 1.75073626e-11, 1.8167685e-8, 1.79805022e-19,\n",
+    "        1.07614011e-19, 3.58865323e-20, 8.18076407e-23, 1.76624918e-21,\n",
+    "        3.95823826e-22, 1.30154195e-21, 7.81338357e-22, 2.60269199e-22,\n",
+    "        1.83743361e-22, 1.21199070e-22, 7.93684710e-23, 7.69690651e-23,\n",
+    "        5.55204217e-24, 4.28655703e-23, 8.91723472e-23, 3.08213125e-22,\n",
+    "        1.84912140e-22, 6.16147989e-23, 5.28143670e-23, 3.77505779e-23,\n",
+    "        2.26343284e-23, 3.19064683e-23, 8.89209221e-24, 9.39180485e-23,\n",
+    "        5.19053997e-23, 2.53784273e-23, 2.27283788e-24, 1.11837286e-22,\n",
+    "        3.85509048e-23, 1.07725425e-22, 6.46134793e-23, 2.15124099e-23,\n",
+    "        2.22314242e-23, 1.59215321e-23, 9.54500244e-24, 1.48014892e-23,\n",
+    "        3.80078829e-23, 4.89438878e-23, 2.71814793e-23, 1.79147356e-23,\n",
+    "        9.80689599e-23, 7.25122090e-23, 8.04471105e-23, 1.27910690e-23,\n",
+    "        1.51551660e-24]\n",
+    "\n",
+    "chianti_levels_o_2 = [9.50427666e-01, 2.93020096e-02, 1.96811489e-02, 3.73612176e-04,\n",
+    "        2.15563005e-04, 1.19850233e-18, 7.62882704e-19, 3.73304507e-19,\n",
+    "        6.86402347e-22, 4.51096625e-22, 1.86694313e-23, 3.93108508e-23,\n",
+    "        5.19716753e-23, 5.48509200e-24, 1.05025627e-23, 4.23498947e-25,\n",
+    "        5.73694095e-23, 2.27177963e-24, 4.53847634e-24, 6.68361893e-24,\n",
+    "        8.38430369e-24, 9.60077484e-25, 6.19972361e-25, 1.04127656e-24,\n",
+    "        1.66000813e-24, 2.48443321e-24, 6.16103282e-24, 9.08726063e-24,\n",
+    "        3.02895913e-24, 1.64793516e-25, 7.92929866e-26, 2.66317598e-23,\n",
+    "        6.44175783e-23, 5.28146898e-25, 7.36488325e-25]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Current issues:\n",
+    "\n",
+    "- Rates matrix layout is different between Chianti and TARDIS despite identical data\n",
+    "- Still some large % differences for high-energy lines\n",
+    "- He I is different for the second and third levels\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import astropy.constants as const\n",
+    "\n",
+    "plt.scatter(chianti_atom_data.levels.loc[2,0].energy * u.erg.to('eV'), chianti_levels_he_1, marker='o', label='CHIANTI')\n",
+    "plt.scatter(chianti_atom_data.levels.loc[2,0].energy * u.erg.to('eV'), level_pops.loc[2,0,:][0], marker='x', label='TARDIS')\n",
+    "plt.scatter(chianti_atom_data.levels.loc[2,0].energy * u.erg.to('eV'), lte_level_pops.loc[2,0,:][0], marker='x', label='TARDIS col only')\n",
+    "plt.xlabel(\"Energy (eV)\")\n",
+    "plt.ylabel(\"Population\")\n",
+    "plt.semilogy()\n",
+    "plt.legend()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "tardis",
+   "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.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}