Update number of shells in SimulationState (#2504)

* update the no_of_shells in SimulationState to be the number of active shells instead of raw shells

* add test for testing the dimensionality after adjust the inner boundary

* trying to avoid the error caused by deep copy

* delete the duplicated lines

* finally find where the dimensionality bug originates!

* fixing bugs that use the wrong inner boundary radius

* fixing other v_inner/outer to the active ones

* reverse the change in the formal_intergral

* again reverse the change that should not be made in those that uses numba1dgeometry

* modify to a neater way - shape the radiationfield to carry the t_rad and W of the active shells only

* small clean up

* the name of no_of_shells is tricking me again in geometry vs simulation state!!

* make sure geometric diluation factor is valid

* nvm, changing assertation criteria in diluteBB is better than change radiationfield to be active shells only, this way one can change the v_inner on fly

tardis/model/base.py

@@ -135,7 +135,9 @@ def dilution_factor(self):
     @dilution_factor.setter
     def dilution_factor(self, new_dilution_factor):
         if len(new_dilution_factor) == self.no_of_shells:
-            self.radiation_field_state.dilution_factor = new_dilution_factor
+            self.radiation_field_state.dilution_factor[
+                self.geometry.v_inner_boundary_index : self.geometry.v_outer_boundary_index
+            ] = new_dilution_factor
         else:
             raise ValueError(
                 "Trying to set dilution_factor for unmatching number"
@@ -151,7 +153,9 @@ def t_radiative(self):
     @t_radiative.setter
     def t_radiative(self, new_t_radiative):
         if len(new_t_radiative) == self.no_of_shells:
-            self.radiation_field_state.t_radiative = new_t_radiative
+            self.radiation_field_state.t_radiative[
+                self.geometry.v_inner_boundary_index : self.geometry.v_outer_boundary_index
+            ] = new_t_radiative
         else:
             raise ValueError(
                 "Trying to set t_radiative for different number of shells."
@@ -224,7 +228,7 @@ def volume(self):
 
     @property
     def no_of_shells(self):
-        return self.geometry.no_of_shells
+        return self.geometry.no_of_shells_active
 
     @property
     def no_of_raw_shells(self):
@@ -262,15 +266,6 @@ def from_config(cls, config, atom_data):
             density, nuclide_mass_fraction, atom_data.atom_data.mass.copy()
         )
 
-        nuclide_mass_fraction = parse_abundance_config(
-            config, geometry, time_explosion
-        )
-
-        # using atom_data.mass.copy() to ensure that the original atom_data is not modified
-        composition = Composition(
-            density, nuclide_mass_fraction, atom_data.atom_data.mass.copy()
-        )
-
         packet_source = parse_packet_source(config, geometry)
         radiation_field_state = parse_radiation_field_state(
             config,

tardis/model/parse_input.py

@@ -568,13 +568,12 @@ def parse_radiation_field_state(
 
     if dilution_factor is None:
         dilution_factor = calculate_geometric_dilution_factor(geometry)
-    elif len(dilution_factor) != geometry.no_of_shells:
-        dilution_factor = dilution_factor[
-            geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
-        ]
-        assert len(dilution_factor) == geometry.no_of_shells
 
-    return DiluteBlackBodyRadiationFieldState(t_radiative, dilution_factor)
+    assert len(dilution_factor) == geometry.no_of_shells
+
+    return DiluteBlackBodyRadiationFieldState(
+        t_radiative, dilution_factor, geometry
+    )
 
 
 def initialize_packet_source(config, geometry, packet_source):
@@ -610,13 +609,13 @@ def initialize_packet_source(config, geometry, packet_source):
 
     luminosity_requested = config.supernova.luminosity_requested
     if config.plasma.initial_t_inner > 0.0 * u.K:
-        packet_source.radius = geometry.r_inner[0]
+        packet_source.radius = geometry.r_inner_active[0]
         packet_source.temperature = config.plasma.initial_t_inner
 
     elif (config.plasma.initial_t_inner < 0.0 * u.K) and (
         luminosity_requested is not None
     ):
-        packet_source.radius = geometry.r_inner[0]
+        packet_source.radius = geometry.r_inner_active[0]
         packet_source.set_temperature_from_luminosity(luminosity_requested)
     else:
         raise ValueError(
@@ -674,9 +673,6 @@ def parse_csvy_radiation_field_state(
         t_radiative = (
             np.ones(geometry.no_of_shells) * config.plasma.initial_t_rad
         )
-        t_radiative = (
-            np.ones(geometry.no_of_shells) * config.plasma.initial_t_rad
-        )
     else:
         t_radiative = calculate_t_radiative_from_t_inner(
             geometry, packet_source
@@ -689,7 +685,9 @@ def parse_csvy_radiation_field_state(
     else:
         dilution_factor = calculate_geometric_dilution_factor(geometry)
 
-    return DiluteBlackBodyRadiationFieldState(t_radiative, dilution_factor)
+    return DiluteBlackBodyRadiationFieldState(
+        t_radiative, dilution_factor, geometry
+    )
 
 
 def calculate_t_radiative_from_t_inner(geometry, packet_source):
@@ -720,6 +718,12 @@ def calculate_geometric_dilution_factor(geometry):
     return 0.5 * (
         1
         - np.sqrt(
-            1 - (geometry.r_inner[0] ** 2 / geometry.r_middle**2).to(1).value
+            1
+            - (
+                geometry.r_inner[geometry.v_inner_boundary_index] ** 2
+                / geometry.r_middle**2
+            )
+            .to(1)
+            .value
         )
     )

tardis/model/radiation_field_state.py

@@ -17,18 +17,37 @@ class DiluteBlackBodyRadiationFieldState:
         Radiative temperature in each shell
     dilution_factor : numpy.ndarray
         Dilution Factors in each shell
+    geometry: tardis.model.Radial1DModel
+        The geometry of the model that uses to constrains the active shells
     """
 
     def __init__(
         self,
         t_radiative: u.Quantity,
         dilution_factor: np.ndarray,
+        geometry=None,
     ):
         # ensuring that the radiation_field has both
         # dilution_factor and t_radiative equal length
         assert len(t_radiative) == len(dilution_factor)
-        assert np.all(t_radiative > 0 * u.K)
-        assert np.all(dilution_factor > 0)
+        if (
+            geometry is not None
+        ):  # check the active shells only (this is used when setting up the radiation_field_state)
+            assert np.all(
+                t_radiative[
+                    geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
+                ]
+                > 0 * u.K
+            )
+            assert np.all(
+                dilution_factor[
+                    geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
+                ]
+                > 0
+            )
+        else:
+            assert np.all(t_radiative > 0 * u.K)
+            assert np.all(dilution_factor > 0)
         self.t_radiative = t_radiative
         self.dilution_factor = dilution_factor
 

tardis/model/tests/test_csvy_model.py

@@ -86,6 +86,32 @@ def test_compare_models(model_config_fnames, atomic_dataset):
     )
 
 
+def test_dimensionality_after_update_v_inner_boundary(
+    example_csvy_file_dir, atomic_dataset
+):
+    """Test that the dimensionality of SimulationState parameters after updating v_inner_boundary
+    in the context of csvy models specifically"""
+    csvy_config_file = example_csvy_file_dir / "radiative_csvy.yml"
+    config = Configuration.from_yaml(csvy_config_file)
+    csvy_model = SimulationState.from_csvy(config, atom_data=atomic_dataset)
+
+    new_config = config
+    new_config.model.v_inner_boundary = csvy_model.velocity[1]
+    new_csvy_model = SimulationState.from_csvy(
+        new_config, atom_data=atomic_dataset
+    )
+
+    assert new_csvy_model.no_of_raw_shells == csvy_model.no_of_raw_shells
+    assert new_csvy_model.no_of_shells == csvy_model.no_of_shells - 1
+    assert new_csvy_model.velocity.shape[0] == csvy_model.velocity.shape[0] - 1
+    assert new_csvy_model.density.shape[0] == csvy_model.density.shape[0] - 1
+    assert new_csvy_model.volume.shape[0] == csvy_model.volume.shape[0] - 1
+    assert (
+        new_csvy_model.t_radiative.shape[0]
+        == csvy_model.t_radiative.shape[0] - 1
+    )
+
+
 @pytest.fixture(scope="module")
 def csvy_model_test_abundances(example_csvy_file_dir, atomic_dataset):
     """Returns SimulationState to use to test abundances dataframes"""

tardis/visualization/tools/sdec_plot.py

@@ -165,7 +165,7 @@ def from_simulation(cls, sim, packets_mode):
             "line_id"
         )
         transport_state = sim.transport.transport_state
-        r_inner = sim.simulation_state.geometry.r_inner
+        r_inner = sim.simulation_state.geometry.r_inner_active
         t_inner = sim.simulation_state.packet_source.temperature
         time_of_simulation = (
             transport_state.packet_collection.time_of_simulation * u.s