Restructure/abundance (#2445)

* restructure of geometry

* add radial1d boundary logic

* black format

* several fixes

* fix epsilon

* add testing of boundaries

* change the r_inner_active

* first integration with `from_config` working

* hunting down density indexing bug

* all model tests (without csvy) pass

* more fixes

* fix of model to simulation_state

* fix inner boundary packet error

* fix some leftovers

* final fix for csvy

* blackify

* restructure to readers and remove some leftover code

* further cleanup

* first start of the restructure

* add comment about removing quantitiness

* add velocity check

* add new abundance functions

* remove default units

* add new matter module

* several restructures. move decay from io to tardis/model/matter

* mid restructure

* slow progress on abundance refactor

* include effective_element_masses

* further updates

* fix of csvy readers

* last fixes to custom abundance widget

* some cleanup

* further cleanup

* removing matter

* fixed last tests

* add composition

* remove isotopemassfractions

* some small fixes

* add the from_config atom_data

* fixing model.ipynb

* fix tests

* fix some of the issues wdigets

* assert mass fractions positive.
Change IsotopeAbundances to Massfractions

* Fix isotope mass fraction typo and t_radiatve duplicate definition

* Document remove functionality for effective element masses

* Fix variable name in
convert_to_nuclide_mass_fraction function

* Refactor test names and variable names in
test_base.py

* Refactor gamma ray simulation setup and test
functions

* fix gamma ray tests

* Add model_isotope_time_0 property to abundances
schema

* Refactor config_reader and config_validator
modules

* Fix model_isotope_time_0 property in
model_definitions.yml

* Refactor simulation module to use pathlib

* remove bad comment and refactor with ruff

* Fix variable name in test_simulation_state_mass()
function

* Fix unit conversion in calculate_cell_masses
function

* Fix logger debug message in CSV reader and update
method names in Composition class

* black formatting

* Refactor code to use composition variable in
test_gamma_ray_transport.py

* Refactor imports in atom_web_download.py and
update return statement in config_internal.py

* fix for documentation not builiding

* add fix for model_isotope_time_0

* fix grid test

* restructure the grid

* final change

docs/io/grid/TardisGridTutorial.ipynb

@@ -18,10 +18,21 @@
     "import pandas as pd\n",
     "import numpy as np\n",
     "from tardis.io.atom_data.util import download_atom_data\n",
+    "from tardis.io.atom_data import AtomData\n",
     "\n",
     "import tardis.grid as grid"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# We download the atomic data needed to run a TARDIS simulation\n",
+    "download_atom_data('kurucz_cd23_chianti_H_He')\n"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -145,8 +156,9 @@
     "# to the base TARDIS model (i.e. using parameters that \n",
     "# are not valid TARDIS Configuration keys), the grid\n",
     "# can return a SimulationState object for a given row.\n",
-    "model = tg.grid_row_to_model(row_index=0)\n",
-    "model"
+    "atom_data = AtomData.from_hdf(tg.config.atom_data)\n",
+    "simulation_state = tg.grid_row_to_simulation_state(row_index=0, atomic_data=atom_data)\n",
+    "simulation_state"
    ]
   },
   {
@@ -155,9 +167,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# We download the atomic data needed to run a TARDIS simulation\n",
-    "download_atom_data('kurucz_cd23_chianti_H_He')\n",
-    "\n",
     "# Users can easily run a full TARDIS simulation\n",
     "# using the grid.\n",
     "sim = tg.run_sim_from_grid(row_index=0)"

docs/physics/setup/model.ipynb

@@ -55,12 +55,14 @@
     "from tardis.io.configuration.config_reader import Configuration\n",
     "from tardis.model import SimulationState\n",
     "from tardis.io.atom_data.util import download_atom_data\n",
+    "from tardis.io.atom_data.base import AtomData\n",
     "from astropy import units as u\n",
     "import matplotlib.pyplot as plt\n",
     "import copy\n",
     "\n",
     "# We download the atomic data needed to run the simulation\n",
-    "download_atom_data('kurucz_cd23_chianti_H_He')"
+    "download_atom_data('kurucz_cd23_chianti_H_He')\n",
+    "atom_data = AtomData.from_hdf('kurucz_cd23_chianti_H_He.h5')"
    ]
   },
   {
@@ -103,19 +105,19 @@
     "shell_config.model.structure.velocity.stop = 2000 * u.km/u.s\n",
     "shell_config.model.structure.velocity.num = 20\n",
     "\n",
-    "shell_model = SimulationState.from_config(shell_config)\n",
-    "\n",
-    "print('velocity:\\n', shell_model.velocity)\n",
-    "print('v_inner:\\n', shell_model.v_inner)\n",
-    "print('v_outer:\\n', shell_model.v_outer)\n",
-    "print('v_middle:\\n', shell_model.v_middle)\n",
-    "print('v_boundary_inner:\\n', shell_model.v_boundary_inner)\n",
-    "print('v_boundary_outer:\\n', shell_model.v_boundary_outer)\n",
-    "print('radius:\\n', shell_model.radius)\n",
-    "print('r_inner:\\n', shell_model.r_inner)\n",
-    "print('r_outer:\\n', shell_model.r_outer)\n",
-    "print('r_middle:\\n', shell_model.r_middle)\n",
-    "print('volume:\\n', shell_model.volume)"
+    "shell_simulation_state = SimulationState.from_config(shell_config, atom_data=atom_data)\n",
+    "\n",
+    "print('velocity:\\n', shell_simulation_state.velocity)\n",
+    "print('v_inner:\\n', shell_simulation_state.v_inner)\n",
+    "print('v_outer:\\n', shell_simulation_state.v_outer)\n",
+    "print('v_middle:\\n', shell_simulation_state.v_middle)\n",
+    "print('v_boundary_inner:\\n', shell_simulation_state.v_boundary_inner)\n",
+    "print('v_boundary_outer:\\n', shell_simulation_state.v_boundary_outer)\n",
+    "print('radius:\\n', shell_simulation_state.radius)\n",
+    "print('r_inner:\\n', shell_simulation_state.r_inner)\n",
+    "print('r_outer:\\n', shell_simulation_state.r_outer)\n",
+    "print('r_middle:\\n', shell_simulation_state.r_middle)\n",
+    "print('volume:\\n', shell_simulation_state.volume)"
    ]
   },
   {
@@ -181,7 +183,7 @@
     "\n",
     "w7_density_config.model.structure.density.type = 'branch85_w7'\n",
     "\n",
-    "w7_density_model = SimulationState.from_config(w7_density_config)\n",
+    "w7_density_model = SimulationState.from_config(w7_density_config, atom_data=atom_data)\n",
     "\n",
     "print('density:\\n', w7_density_model.density)\n",
     "\n",
@@ -209,7 +211,7 @@
     "\n",
     "w7_modified_config.supernova.time_explosion = 12 * u.day\n",
     "\n",
-    "w7_modified_model = SimulationState.from_config(w7_modified_config)\n",
+    "w7_modified_model = SimulationState.from_config(w7_modified_config, atom_data=atom_data)\n",
     "\n",
     "print('density:\\n', w7_modified_model.density)\n",
     "\n",
@@ -273,7 +275,7 @@
     "uni_density_config.model.structure.density.time_0 = 1 * u.day\n",
     "uni_density_config.model.structure.density.value = 5e-10 * u.kg/u.cm**3\n",
     "\n",
-    "uni_density_model = SimulationState.from_config(uni_density_config)\n",
+    "uni_density_model = SimulationState.from_config(uni_density_config, atom_data=atom_data)\n",
     "\n",
     "print('density:\\n', uni_density_model.density)\n",
     "\n",
@@ -327,7 +329,7 @@
     "pow_density_config.model.structure.density.v_0 = 500 * u.km/u.s\n",
     "pow_density_config.model.structure.density.exponent = -2\n",
     "\n",
-    "pow_density_model = SimulationState.from_config(pow_density_config)\n",
+    "pow_density_model = SimulationState.from_config(pow_density_config, atom_data=atom_data)\n",
     "\n",
     "print('density:\\n', pow_density_model.density)\n",
     "\n",
@@ -378,7 +380,7 @@
     "exp_density_config.model.structure.density.rho_0 = 5e-10 * u.kg/u.cm**3\n",
     "exp_density_config.model.structure.density.v_0 = 500 * u.km/u.s\n",
     "\n",
-    "exp_density_model = SimulationState.from_config(exp_density_config)\n",
+    "exp_density_model = SimulationState.from_config(exp_density_config, atom_data=atom_data)\n",
     "\n",
     "print('density:\\n', exp_density_model.density)\n",
     "\n",
@@ -423,7 +425,7 @@
     "                                            'Si':.3,\n",
     "                                            'He': .1}\n",
     "\n",
-    "abund_model = SimulationState.from_config(abund_config)\n",
+    "abund_model = SimulationState.from_config(abund_config, atom_data=atom_data)\n",
     "\n",
     "abund_model.abundance"
    ]
@@ -459,7 +461,7 @@
     "                                            'Ni57':.1,\n",
     "                                            'Cr51':.4}\n",
     "\n",
-    "abund_isotopes_model = SimulationState.from_config(abund_isotopes_config)\n",
+    "abund_isotopes_model = SimulationState.from_config(abund_isotopes_config, atom_data=atom_data)\n",
     "\n",
     "abund_isotopes_model.abundance"
    ]
@@ -527,7 +529,7 @@
     "t_rad_config.model.structure.velocity.stop = 2000 * u.km/u.s\n",
     "t_rad_config.model.structure.velocity.num = 20\n",
     "\n",
-    "t_rad_model = SimulationState.from_config(t_rad_config)\n",
+    "t_rad_model = SimulationState.from_config(t_rad_config, atom_data=atom_data)\n",
     "\n",
     "print('t_inner:\\n', t_rad_model.t_inner)\n",
     "print('t_rad:\\n', t_rad_model.t_rad)\n",
@@ -590,7 +592,7 @@
     "w_config.model.structure.velocity.stop = 2000 * u.km/u.s\n",
     "w_config.model.structure.velocity.num = 20\n",
     "\n",
-    "w_model = SimulationState.from_config(w_config)\n",
+    "w_model = SimulationState.from_config(w_config, atom_data=atom_data)\n",
     "\n",
     "print('w:\\n', w_model.w)\n",
     "\n",

tardis/energy_input/gamma_ray_transport.py

@@ -1,23 +1,19 @@
+import astropy.units as u
 import numpy as np
 import pandas as pd
-import astropy.units as u
+import radioactivedecay as rd
 from numba import njit
 from numba.typed import List
-import radioactivedecay as rd
-
-from tardis.montecarlo.montecarlo_numba import njit_dict_no_parallel
-from tardis.montecarlo.montecarlo_numba.opacities import M_P
 
 from tardis.energy_input.energy_source import (
     get_all_isotopes,
-    get_nuclear_lines_database,
     positronium_continuum,
-    read_artis_lines,
     setup_input_energy,
 )
-from tardis.energy_input.samplers import initial_packet_radius
 from tardis.energy_input.GXPacket import initialize_packet_properties
-from tardis.energy_input.gamma_packet_loop import gamma_packet_loop
+from tardis.energy_input.samplers import initial_packet_radius
+from tardis.montecarlo.montecarlo_numba import njit_dict_no_parallel
+from tardis.montecarlo.montecarlo_numba.opacities import M_P
 
 # Energy: keV, exported as eV for SF solver
 # distance: cm
@@ -165,7 +161,6 @@ def initialize_packets(
     array
         Array of positron output dataframe rows
     """
-
     packets = List()
 
     number_of_packets = decays_per_isotope.sum().sum()
@@ -259,40 +254,22 @@ def initialize_packets(
     )
 
 
-def calculate_shell_masses(model):
-    """Function to calculate shell masses
-    Parameters
-    ----------
-    model : tardis.Radial1DModel
-        The tardis model to calculate gamma ray propagation through
-    Returns
-    -------
-    numpy.ndarray
-        shell masses in units of g
-
-    """
-
-    ejecta_density = model.density.to("g/cm^3")
-    ejecta_volume = model.volume.to("cm^3")
-    return (ejecta_volume * ejecta_density).to(u.g)
-
-
 def calculate_total_decays(inventories, time_delta):
     """Function to create inventories of isotope
+
     Parameters
     ----------
     model : tardis.Radial1DModel
         The tardis model to calculate gamma ray propagation through
 
     time_end : float
         End time of simulation in days
+
     Returns
     -------
         Total decay list : List
             list of total decays for x g of isotope for time 't'
-
     """
-
     time_delta = u.Quantity(time_delta, u.s)
 
     total_decays_list = []
@@ -303,15 +280,17 @@ def calculate_total_decays(inventories, time_delta):
     return total_decays_list
 
 
-def create_isotope_dicts(raw_isotope_abundance, shell_masses):
+def create_isotope_dicts(raw_isotope_abundance, cell_masses):
     """
     Function to create a dictionary of isotopes for each shell with their masses.
+
     Parameters
     ----------
     raw_isotope_abundance : pd.DataFrame
         isotope abundance in mass fractions.
-    shell_masses : numpy.ndarray
+    cell_masses : numpy.ndarray
         shell masses in units of g
+
     Returns
     -------
         isotope_dicts : Dict
@@ -330,18 +309,20 @@ def create_isotope_dicts(raw_isotope_abundance, shell_masses):
             isotope_dicts[i][atomic_number, mass_number] = {}
             nuclear_symbol = f"{rd.utils.Z_to_elem(atomic_number)}{mass_number}"
             isotope_dicts[i][atomic_number, mass_number][nuclear_symbol] = (
-                abundances[i] * shell_masses[i].to(u.g).value
+                abundances[i] * cell_masses[i].to(u.g).value
             )
 
     return isotope_dicts
 
 
 def create_inventories_dict(isotope_dict):
     """Function to create dictionary of inventories for each shell
+
     Parameters
     ----------
     isotope_dict : Dict
         dictionary of isotopes for each shell with their ``masses``.
+
     Returns
     -------
         inv : Dict
@@ -363,12 +344,14 @@ def create_inventories_dict(isotope_dict):
 def calculate_total_decays(inventory_dict, time_delta):
     """
     Function to calculate total decays for each isotope in each shell
+
     Parameters
     ----------
     inventory_dict : Dict
         dictionary of inventories for each shell
     time_delta : float
         time interval in units of time (days/mins/secs etc)
+
     Returns
     -------
         total_decays : Dict
@@ -391,6 +374,7 @@ def calculate_average_energies(raw_isotope_abundance, gamma_ray_lines):
     """
     Function to calculate average energies of positrons and gamma rays
     from a list of gamma ray lines from nndc.
+
     Parameters
     ----------
     raw_isotope_abundance : pd.DataFrame
@@ -408,7 +392,6 @@ def calculate_average_energies(raw_isotope_abundance, gamma_ray_lines):
         list of gamma ray lines
 
     """
-
     all_isotope_names = get_all_isotopes(raw_isotope_abundance)
     all_isotope_names.sort()
 
@@ -455,6 +438,7 @@ def calculate_average_energies(raw_isotope_abundance, gamma_ray_lines):
 def get_taus(raw_isotope_abundance):
     """
     Function to calculate taus for each isotope
+
     Parameters
     ----------
     raw_isotope_abundance : pd.DataFrame
@@ -490,6 +474,7 @@ def decay_chain_energies(
 ):
     """
     Function to calculate decay chain energies.
+
     Parameters
     ----------
     raw_isotope_abundance : pd.DataFrame
@@ -502,6 +487,7 @@ def decay_chain_energies(
         list of gamma ray lines
     total_decays : Dict
         dictionary of total decays for each isotope in each shell
+
     Returns
     -------
     decay_energy : Dict
@@ -519,19 +505,19 @@ def decay_chain_energies(
     return decay_energy
 
 
-def calculate_energy_per_mass(
-    decay_energy, raw_isotope_abundance, shell_masses
-):
+def calculate_energy_per_mass(decay_energy, raw_isotope_abundance, cell_masses):
     """
     Function to calculate decay energy per mass for each isotope chain.
+
     Parameters
     ----------
     decay_energy : Dict
         dictionary of decay chain energies for each isotope in each shell
     raw_isotope_abundance : pd.DataFrame
         isotope abundance in mass fractions.
-    shell_masses : numpy.ndarray
+    cell_masses : numpy.ndarray
         shell masses in units of g
+
     Returns
     -------
     energy_per_mass : pd.DataFrame
@@ -567,7 +553,7 @@ def calculate_energy_per_mass(
     )
 
     energy_per_mass = energy_df.divide(
-        (raw_isotope_abundance * shell_masses).to_numpy(), axis=0
+        (raw_isotope_abundance * cell_masses).to_numpy(), axis=0
     )
 
     return energy_per_mass, energy_df