V inner formal integral

docs/physics/setup/model.ipynb

@@ -81,7 +81,7 @@
    "id": "cee054e9",
    "metadata": {},
    "source": [
-    "In the cell below, we set up a model. We use the [specific structure](../../io/configuration/components/models/index.rst#specific-structure) where we supply $t_\\mathrm{explosion}$, the velocity of the inner and outer boundaries of the supernova (labeled `start` and `stop`), and the number of shells (labeled `num`). The shells are then evenly spaced between the inner and outer boundaries of the supernova. The time after the explosion, the inner and outer velocities, and the number of shells can be varied to get different shell structures. The `SimulationState` object stores information about the model in the following attributes: `velocity` shows the velocity of the shell boundaries, `v_inner` shows the velocities of the inner boundaries of each shell, `v_outer` shows the velocity of the outer boundaries of each shell, and `v_middle` shows the velocity of the middle of each shell. Similarly, `radius`, `r_inner`, `r_outer`, and `r_middle` show the radii of each shell boundary, the inner boundaries, the outer boundaries, and the middles of each shell, respectively. `v_boundary_inner` shows the velocity of the inner boundary of the computational domain, and `v_boundary_outer` shows the velocity of the outer boundary of the computational domain. Finally, `volume` shows the volume of each shell, calculated via the formula of the volume of a spherical shell: $V=\\frac{4}{3}\\pi (r_\\mathrm{outer}^3-r_\\mathrm{inner}^3)$."
+    "In the cell below, we set up a model. We use the [specific structure](../../io/configuration/components/models/index.rst#specific-structure) where we supply $t_\\mathrm{explosion}$, the velocity of the inner and outer boundaries of the supernova (labeled `start` and `stop`), and the number of shells (labeled `num`). The shells are then evenly spaced between the inner and outer boundaries of the supernova. The time after the explosion, the inner and outer velocities, and the number of shells can be varied to get different shell structures. The `SimulationState` object stores information about the model in the following attributes: `velocity` shows the velocity of the shell boundaries, `v_inner` shows the velocities of the inner boundaries of each shell, `v_outer` shows the velocity of the outer boundaries of each shell, and `v_middle` shows the velocity of the middle of each shell. Similarly, `radius`, `r_inner`, `r_outer`, and `r_middle` show the radii of each shell boundary, the inner boundaries, the outer boundaries, and the middles of each shell, respectively. `v_inner_boundary` shows the velocity of the inner boundary of the computational domain, and `v_outer_boundary` shows the velocity of the outer boundary of the computational domain. Finally, `volume` shows the volume of each shell, calculated via the formula of the volume of a spherical shell: $V=\\frac{4}{3}\\pi (r_\\mathrm{outer}^3-r_\\mathrm{inner}^3)$."
    ]
   },
   {
@@ -111,8 +111,8 @@
     "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('v_inner_boundary:\\n', shell_simulation_state.v_inner_boundary)\n",
+    "print('v_outer_boundary:\\n', shell_simulation_state.v_outer_boundary)\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",
@@ -125,7 +125,7 @@
    "id": "1ee56110",
    "metadata": {},
    "source": [
-    "Notice that `radius = velocity*time_explosion`, and similarly for `r_inner`, `r_outer`, and `r_middle`. You can get the radius of the photosphere via `v_boundary_inner*time_explosion` and outer edge of the supernova via `v_boundary_outer*time_explosion`.\n",
+    "Notice that `radius = velocity*time_explosion`, and similarly for `r_inner`, `r_outer`, and `r_middle`. You can get the radius of the photosphere via `v_inner_boundary*time_explosion` and outer edge of the supernova via `v_boundary_outer*time_explosion`.\n",
     "\n",
     "<div class=\"alert alert-info\">\n",
     "    \n",

docs/workflows/v_inner_solver_workflow.ipynb

@@ -0,0 +1,128 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from tardis.workflows.v_inner_solver import InnerVelocitySolverWorkflow\n",
+    "from tardis.io.configuration.config_reader import Configuration"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "config = Configuration.from_yaml('../tardis_example.yml')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This code modifies the TARDIS example configuration to include convergence information for the inner boundary velocity solver.\n",
+    "Note that the number of shells is increased and the starting velocity is reduced to provide more granularity and a wider search window, respectively."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "from astropy import units as u\n",
+    "\n",
+    "config.montecarlo.convergence_strategy['v_inner_boundary'] = {\n",
+    "    'damping_constant' : 0.5,\n",
+    "    'threshold' : 0.01,\n",
+    "    'type' : 'damped'\n",
+    "    }\n",
+    "\n",
+    "config.montecarlo.convergence_strategy.stop_if_converged = True\n",
+    "config.model.structure.velocity.start = 5000 * u.km/u.s # Larger window over which to search\n",
+    "config.model.structure.velocity.num = 50 # Increase number of shells\n",
+    "\n",
+    "workflow = InnerVelocitySolverWorkflow(\n",
+    "    config, tau=2.0/3,\n",
+    "    mean_optical_depth=\"rosseland\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "workflow.run()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spectrum = workflow.spectrum_solver.spectrum_real_packets\n",
+    "spectrum_virtual = workflow.spectrum_solver.spectrum_virtual_packets\n",
+    "spectrum_integrated = workflow.spectrum_solver.spectrum_integrated"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "plt.figure(figsize=(10, 6.5))\n",
+    "\n",
+    "spectrum.plot(label=\"Normal packets\")\n",
+    "spectrum_virtual.plot(label=\"Virtual packets\")\n",
+    "spectrum_integrated.plot(label='Formal integral')\n",
+    "\n",
+    "plt.xlim(500, 9000)\n",
+    "plt.title(\"TARDIS example model spectrum\")\n",
+    "plt.xlabel(r\"Wavelength [$\\AA$]\")\n",
+    "plt.ylabel(r\"Luminosity density [erg/s/$\\AA$]\")\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "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
+}

tardis/io/configuration/config_reader.py

@@ -431,7 +431,7 @@ def parse_convergence_section(convergence_section_dict):
         """
         convergence_parameters = ["damping_constant", "threshold", "type"]
 
-        for convergence_variable in ["t_inner", "t_rad", "w"]:
+        for convergence_variable in ["t_inner", "t_rad", "w", "v_inner_boundary"]:
             if convergence_variable not in convergence_section_dict:
                 convergence_section_dict[convergence_variable] = {}
             convergence_variable_section = convergence_section_dict[

tardis/io/configuration/schemas/montecarlo_definitions.yml

@@ -5,7 +5,6 @@ definitions:
       title: 'Damped Convergence Strategy'
       type: object
       additionalProperties: false
-      properties:
       properties:
         type:
           enum:
@@ -59,6 +58,24 @@ definitions:
               type: string
               default: 'damped'
               description: THIS IS A DUMMY VARIABLE DO NOT USE
+        v_inner_boundary:
+          type: object
+          additionalProperties: false
+          properties:
+            damping_constant:
+              type: number
+              default: 0.0
+              description: damping constant
+              minimum: 0
+            threshold:
+              type: number
+              description: specifies the threshold that is taken as convergence  (i.e.
+                0.05 means that the value does not change more than 5%)
+              minimum: 0
+            type:
+              type: string
+              default: 'damped'
+              description: THIS IS A DUMMY VARIABLE DO NOT USE
         t_rad:
           type: object
           additionalProperties: false

tardis/io/model/parse_geometry_configuration.py

@@ -135,17 +135,17 @@
     """
     if hasattr(config, "model"):
         if hasattr(config.model, "v_inner_boundary"):
-            v_boundary_inner = config.model.v_inner_boundary
+            v_inner_boundary = config.model.v_inner_boundary
         else:
-            v_boundary_inner = None
+            v_inner_boundary = None
 
         if hasattr(config.model, "v_outer_boundary"):
-            v_boundary_outer = config.model.v_outer_boundary
+            v_outer_boundary = config.model.v_outer_boundary
         else:
-            v_boundary_outer = None
+            v_outer_boundary = None
     else:
-        v_boundary_inner = None
-        v_boundary_outer = None
+        v_inner_boundary = None
+        v_outer_boundary = None
 
     if hasattr(csvy_model_config, "velocity"):
         velocity = quantity_linspace(
@@ -166,8 +166,8 @@
     geometry = HomologousRadial1DGeometry(
         velocity[:-1],  # v_inner
         velocity[1:],  # v_outer
-        v_inner_boundary=v_boundary_inner,
-        v_outer_boundary=v_boundary_outer,
+        v_inner_boundary=v_inner_boundary,
+        v_outer_boundary=v_outer_boundary,
         time_explosion=time_explosion,
     )
     return geometry

tardis/model/base.py

@@ -58,11 +58,11 @@ class SimulationState(HDFWriterMixin):
     dilution_factor : np.ndarray
         If None, the dilution_factor will be initialized with the geometric
         dilution factor.
-    v_boundary_inner : astropy.units.Quantity
-    v_boundary_outer : astropy.units.Quantity
+    v_inner_boundary : astropy.units.Quantity
+    v_outer_boundary : astropy.units.Quantity
     raw_velocity : np.ndarray
         The complete array of the velocities, without being cut by
-        `v_boundary_inner` and `v_boundary_outer`
+        `v_inner_boundary` and `v_outer_boundary`
     electron_densities : astropy.units.quantity.Quantity
 
     Attributes
@@ -81,8 +81,8 @@ class SimulationState(HDFWriterMixin):
     density : astropy.units.quantity.Quantity
     volume : astropy.units.quantity.Quantity
     no_of_shells : int
-        The number of shells as formed by `v_boundary_inner` and
-        `v_boundary_outer`
+        The number of shells as formed by `v_inner_boundary` and
+        `v_outer_boundary`
     no_of_raw_shells : int
     """
 
@@ -206,11 +206,11 @@ def radius(self):
         return self.time_explosion * self.velocity
 
     @property
-    def v_boundary_inner(self):
+    def v_inner_boundary(self):
         return self.geometry.v_inner_boundary
 
     @property
-    def v_boundary_outer(self):
+    def v_outer_boundary(self):
         return self.geometry.v_outer_boundary
 
     @property

tardis/spectrum/formal_integral.py

@@ -434,9 +434,22 @@ def make_source_function(self):
         -------
         Numpy array containing ( 1 - exp(-tau_ul) ) S_ul ordered by wavelength of the transition u -> l
         """
+
         simulation_state = self.simulation_state
+        # slice for the active shells
+        local_slice = slice(
+            simulation_state.geometry.v_inner_boundary_index,
+            simulation_state.geometry.v_outer_boundary_index,
+        )
+
         transport = self.transport
         montecarlo_transport_state = transport.transport_state
+        transition_probabilities = self.opacity_state.transition_probabilities[
+            :, local_slice
+        ]
+        tau_sobolevs = self.opacity_state.tau_sobolev[:, local_slice]
+
+        columns = range(simulation_state.no_of_shells)
 
         # macro_ref = self.atomic_data.macro_atom_references
         macro_ref = self.atomic_data.macro_atom_references
@@ -449,9 +462,7 @@ def make_source_function(self):
         if transport.line_interaction_type == "macroatom":
             internal_jump_mask = (macro_data.transition_type >= 0).values
             ma_int_data = macro_data[internal_jump_mask]
-            internal = self.opacity_state.transition_probabilities[
-                internal_jump_mask
-            ]
+            internal = transition_probabilities[internal_jump_mask]
 
             source_level_idx = ma_int_data.source_level_idx.values
             destination_level_idx = ma_int_data.destination_level_idx.values
@@ -460,7 +471,7 @@ def make_source_function(self):
             montecarlo_transport_state.packet_collection.time_of_simulation
             * simulation_state.volume
         )
-        exptau = 1 - np.exp(-self.opacity_state.tau_sobolev)
+        exptau = 1 - np.exp(-tau_sobolevs)
         Edotlu = (
             Edotlu_norm_factor
             * exptau
@@ -493,14 +504,14 @@ def make_source_function(self):
         upper_level_index = self.atomic_data.lines.index.droplevel(
             "level_number_lower"
         )
-        e_dot_lu = pd.DataFrame(Edotlu.value, index=upper_level_index)
+        e_dot_lu = pd.DataFrame(
+            Edotlu.value, index=upper_level_index, columns=columns
+        )
         e_dot_u = e_dot_lu.groupby(level=[0, 1, 2]).sum()
         e_dot_u_src_idx = macro_ref.loc[e_dot_u.index].references_idx.values
 
         if transport.line_interaction_type == "macroatom":
-            C_frame = pd.DataFrame(
-                columns=np.arange(no_shells), index=macro_ref.index
-            )
+            C_frame = pd.DataFrame(columns=columns, index=macro_ref.index)
             q_indices = (source_level_idx, destination_level_idx)
             for shell in range(no_shells):
                 Q = sp.coo_matrix(
@@ -523,7 +534,7 @@ def make_source_function(self):
             ["atomic_number", "ion_number", "source_level_number"]
         ).index.copy()
         tmp = pd.DataFrame(
-            self.opacity_state.transition_probabilities[
+            transition_probabilities[
                 (self.atomic_data.macro_atom_data.transition_type == -1).values
             ]
         )
@@ -539,13 +550,15 @@ def make_source_function(self):
         att_S_ul = wave * (q_ul * e_dot_u) * t / (4 * np.pi)
 
         result = pd.DataFrame(
-            att_S_ul.values, index=transitions.transition_line_id.values
+            att_S_ul.values,
+            index=transitions.transition_line_id.values,
+            columns=columns,
         )
         att_S_ul = result.loc[lines.index.values].values
 
         # Jredlu should already by in the correct order, i.e. by wavelength of
         # the transition l->u (similar to Jbluelu)
-        Jredlu = Jbluelu * np.exp(-self.opacity_state.tau_sobolev) + att_S_ul
+        Jredlu = Jbluelu * np.exp(-tau_sobolevs) + att_S_ul
         if self.interpolate_shells > 0:
             (
                 att_S_ul,
@@ -592,15 +605,20 @@ def interpolate_integrator_quantities(
 
         transport.electron_densities_integ = interp1d(
             r_middle,
-            plasma.electron_densities,
+            plasma.electron_densities.iloc[
+                self.simulation_state.geometry.v_inner_boundary_index : self.simulation_state.geometry.v_outer_boundary_index
+            ],
             fill_value="extrapolate",
             kind="nearest",
         )(r_middle_integ)
         # Assume tau_sobolevs to be constant within a shell
         # (as in the MC simulation)
         transport.tau_sobolevs_integ = interp1d(
             r_middle,
-            self.opacity_state.tau_sobolev,
+            self.opacity_state.tau_sobolev[
+                :,
+                self.simulation_state.geometry.v_inner_boundary_index : self.simulation_state.geometry.v_outer_boundary_index,
+            ],
             fill_value="extrapolate",
             kind="nearest",
         )(r_middle_integ)

tardis/workflows/simple_tardis_workflow.py

@@ -268,6 +268,8 @@
             "t_inner"
         ]
 
+        return next_values
+
     def solve_plasma(self, estimated_radfield_properties):
         """Update the plasma solution with the new radiation field estimates