Further testing fixes

tardis/model/base.py

@@ -235,7 +235,7 @@ def no_of_raw_shells(self):
         return self.geometry.no_of_shells
 
     @classmethod
-    def from_config(cls, config, atom_data):
+    def from_config(cls, config, atom_data, legacy_mode_enabled=False):
         """
         Create a new SimulationState instance from a Configuration object.
 
@@ -266,7 +266,9 @@ def from_config(cls, config, atom_data):
             density, nuclide_mass_fraction, atom_data.atom_data.mass.copy()
         )
 
-        packet_source = parse_packet_source(config, geometry)
+        packet_source = parse_packet_source(
+            config, geometry, legacy_mode_enabled
+        )
         radiation_field_state = parse_radiation_field_state(
             config,
             t_radiative,
@@ -285,7 +287,7 @@ def from_config(cls, config, atom_data):
         )
 
     @classmethod
-    def from_csvy(cls, config, atom_data=None):
+    def from_csvy(cls, config, atom_data=None, legacy_mode_enabled=False):
         """
         Create a new SimulationState instance from a Configuration object.
 
@@ -363,7 +365,9 @@ def from_csvy(cls, config, atom_data=None):
             geometry,
         )
 
-        packet_source = parse_packet_source(config, geometry)
+        packet_source = parse_packet_source(
+            config, geometry, legacy_mode_enabled
+        )
 
         radiation_field_state = parse_csvy_radiation_field_state(
             config, csvy_model_config, csvy_model_data, geometry, packet_source

tardis/model/parse_input.py

@@ -576,7 +576,9 @@ def parse_radiation_field_state(
     )
 
 
-def initialize_packet_source(config, geometry, packet_source):
+def initialize_packet_source(
+    config, geometry, packet_source, legacy_mode_enabled
+):
     """
     Initialize the packet source based on config and geometry
 
@@ -603,9 +605,13 @@ def initialize_packet_source(config, geometry, packet_source):
         packet_source = BlackBodySimpleSourceRelativistic(
             base_seed=config.montecarlo.seed,
             time_explosion=config.supernova.time_explosion,
+            legacy_mode_enabled=legacy_mode_enabled,
         )
     else:
-        packet_source = BlackBodySimpleSource(base_seed=config.montecarlo.seed)
+        packet_source = BlackBodySimpleSource(
+            base_seed=config.montecarlo.seed,
+            legacy_mode_enabled=legacy_mode_enabled,
+        )
 
     luminosity_requested = config.supernova.luminosity_requested
     if config.plasma.initial_t_inner > 0.0 * u.K:
@@ -625,7 +631,7 @@ def initialize_packet_source(config, geometry, packet_source):
     return packet_source
 
 
-def parse_packet_source(config, geometry):
+def parse_packet_source(config, geometry, legacy_mode_enabled):
     """
     Parse the packet source based on the given configuration and geometry.
 
@@ -645,11 +651,17 @@ def parse_packet_source(config, geometry):
         packet_source = BlackBodySimpleSourceRelativistic(
             base_seed=config.montecarlo.seed,
             time_explosion=config.supernova.time_explosion,
+            legacy_mode_enabled=legacy_mode_enabled,
         )
     else:
-        packet_source = BlackBodySimpleSource(base_seed=config.montecarlo.seed)
+        packet_source = BlackBodySimpleSource(
+            base_seed=config.montecarlo.seed,
+            legacy_mode_enabled=legacy_mode_enabled,
+        )
 
-    return initialize_packet_source(config, geometry, packet_source)
+    return initialize_packet_source(
+        config, geometry, packet_source, legacy_mode_enabled
+    )
 
 
 def parse_csvy_radiation_field_state(

tardis/montecarlo/montecarlo_numba/formal_integral.py

@@ -108,9 +108,7 @@ def numba_formal_integral(
             p = pp[p_idx]
 
             # initialize z intersections for p values
-            size_z = populate_z(
-                geometry, model, p, z, shell_id
-            )  # check returns
+            size_z = populate_z(geometry, model, p, z, shell_id)  # check returns
             # initialize I_nu
             if p <= R_ph:
                 I_nu[p_idx] = intensity_black_body(nu * z[0], iT)
@@ -145,9 +143,7 @@ def numba_formal_integral(
                 for _ in range(max(nu_end_idx - pline, 0)):
                     # calculate e-scattering optical depth to next resonance point
                     zend = (
-                        model.time_explosion
-                        / C_INV
-                        * (1.0 - line_list_nu[pline] / nu)
+                        model.time_explosion / C_INV * (1.0 - line_list_nu[pline] / nu)
                     )  # check
 
                     if first == 1:
@@ -187,12 +183,8 @@ def numba_formal_integral(
                 # calculate e-scattering optical depth to grid cell boundary
 
                 Jkkp = 0.5 * (Jred_lu[pJred_lu] + Jblue_lu[pJblue_lu])
-                zend = (
-                    model.time_explosion / C_INV * (1.0 - nu_end / nu)
-                )  # check
-                escat_contrib += (
-                    (zend - zstart) * escat_op * (Jkkp - I_nu[p_idx])
-                )
+                zend = model.time_explosion / C_INV * (1.0 - nu_end / nu)  # check
+                escat_contrib += (zend - zstart) * escat_op * (Jkkp - I_nu[p_idx])
                 zstart = zend
 
                 # advance pointers
@@ -281,7 +273,8 @@ def __init__(self, simulation_state, plasma, transport, points=1000):
         self.simulation_state = simulation_state
         self.transport = transport
         self.points = points
-        self.montecarlo_configuration = self.transport.montecarlo_configuration
+        if transport:
+            self.montecarlo_configuration = self.transport.montecarlo_configuration
         if plasma:
             self.plasma = opacity_state_initialize(
                 plasma,
@@ -310,7 +303,10 @@ def generate_numba_objects(self):
             self.simulation_state.time_explosion.cgs.value,
         )
         self.opacity_state = opacity_state_initialize(
-            self.original_plasma, self.transport.line_interaction_type
+            self.original_plasma,
+            self.transport.line_interaction_type,
+            self.montecarlo_configuration.DISABLE_LINE_SCATTERING,
+            self.montecarlo_configuration.CONTINUUM_PROCESSES_ENABLED,
         )
         if self.transport.use_gpu:
             self.integrator = CudaFormalIntegrator(
@@ -434,9 +430,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.original_plasma.transition_probabilities[
-                internal_jump_mask
-            ]
+            internal = self.original_plasma.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
@@ -475,17 +469,13 @@ def make_source_function(self):
             * Jbluelu_norm_factor
         )
 
-        upper_level_index = self.atomic_data.lines.index.droplevel(
-            "level_number_lower"
-        )
+        upper_level_index = self.atomic_data.lines.index.droplevel("level_number_lower")
         e_dot_lu = pd.DataFrame(Edotlu.values, index=upper_level_index)
         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=np.arange(no_shells), index=macro_ref.index)
             q_indices = (source_level_idx, destination_level_idx)
             for shell in range(no_shells):
                 Q = sp.coo_matrix(
@@ -502,8 +492,7 @@ def make_source_function(self):
             "source_level_number",
         ]  # To make the q_ul e_dot_u product work, could be cleaner
         transitions = self.original_plasma.atomic_data.macro_atom_data[
-            self.original_plasma.atomic_data.macro_atom_data.transition_type
-            == -1
+            self.original_plasma.atomic_data.macro_atom_data.transition_type == -1
         ].copy()
         transitions_index = transitions.set_index(
             ["atomic_number", "ion_number", "source_level_number"]
@@ -515,9 +504,9 @@ def make_source_function(self):
         t = simulation_state.time_explosion.value
         t = simulation_state.time_explosion.value
         lines = self.atomic_data.lines.set_index("line_id")
-        wave = lines.wavelength_cm.loc[
-            transitions.transition_line_id
-        ].values.reshape(-1, 1)
+        wave = lines.wavelength_cm.loc[transitions.transition_line_id].values.reshape(
+            -1, 1
+        )
         if transport.line_interaction_type == "macroatom":
             e_dot_u = C_frame.loc[e_dot_u.index]
         att_S_ul = wave * (q_ul * e_dot_u) * t / (4 * np.pi)
@@ -540,24 +529,16 @@ def make_source_function(self):
                 att_S_ul, Jredlu, Jbluelu, e_dot_u
             )
         else:
-            transport.r_inner_i = (
-                montecarlo_transport_state.geometry_state.r_inner
-            )
-            transport.r_outer_i = (
-                montecarlo_transport_state.geometry_state.r_outer
-            )
-            transport.tau_sobolevs_integ = (
-                self.original_plasma.tau_sobolevs.values
-            )
+            transport.r_inner_i = montecarlo_transport_state.geometry_state.r_inner
+            transport.r_outer_i = montecarlo_transport_state.geometry_state.r_outer
+            transport.tau_sobolevs_integ = self.original_plasma.tau_sobolevs.values
             transport.electron_densities_integ = (
                 self.original_plasma.electron_densities.values
             )
 
         return att_S_ul, Jredlu, Jbluelu, e_dot_u
 
-    def interpolate_integrator_quantities(
-        self, att_S_ul, Jredlu, Jbluelu, e_dot_u
-    ):
+    def interpolate_integrator_quantities(self, att_S_ul, Jredlu, Jbluelu, e_dot_u):
         transport = self.transport
         mct_state = transport.transport_state
         plasma = self.original_plasma
@@ -596,12 +577,8 @@ def interpolate_integrator_quantities(
         Jredlu = pd.DataFrame(
             interp1d(r_middle, Jredlu, fill_value="extrapolate")(r_middle_integ)
         )
-        Jbluelu = interp1d(r_middle, Jbluelu, fill_value="extrapolate")(
-            r_middle_integ
-        )
-        e_dot_u = interp1d(r_middle, e_dot_u, fill_value="extrapolate")(
-            r_middle_integ
-        )
+        Jbluelu = interp1d(r_middle, Jbluelu, fill_value="extrapolate")(r_middle_integ)
+        e_dot_u = interp1d(r_middle, e_dot_u, fill_value="extrapolate")(r_middle_integ)
 
         # Set negative values from the extrapolation to zero
         att_S_ul = att_S_ul.clip(0.0)
@@ -647,9 +624,7 @@ def formal_integral(self, nu, N):
             ]
         )
         error = np.max(np.abs((L_test - L) / L))
-        assert (
-            error < 1e-7
-        ), f"Incorrect I_nu_p values, max relative difference:{error}"
+        assert error < 1e-7, f"Incorrect I_nu_p values, max relative difference:{error}"
 
         return np.array(L, np.float64)
 

tardis/montecarlo/montecarlo_numba/tests/conftest.py

@@ -31,7 +31,10 @@ def nb_simulation_verysimple(config_verysimple, atomic_dataset):
 @pytest.fixture(scope="package")
 def verysimple_opacity_state(nb_simulation_verysimple):
     return opacity_state_initialize(
-        nb_simulation_verysimple.plasma, line_interaction_type="macroatom"
+        nb_simulation_verysimple.plasma,
+        line_interaction_type="macroatom",
+        disable_line_scattering=False,
+        continuum_processes_enabled=False,
     )
 
 
@@ -67,9 +70,7 @@ def verysimple_estimators(nb_simulation_verysimple):
 
 @pytest.fixture(scope="package")
 def verysimple_vpacket_collection(nb_simulation_verysimple):
-    spectrum_frequency = (
-        nb_simulation_verysimple.transport.spectrum_frequency.value
-    )
+    spectrum_frequency = nb_simulation_verysimple.transport.spectrum_frequency.value
     return VPacketCollection(
         source_rpacket_index=0,
         spectrum_frequency=spectrum_frequency,
@@ -82,9 +83,7 @@ def verysimple_vpacket_collection(nb_simulation_verysimple):
 
 @pytest.fixture(scope="package")
 def verysimple_3vpacket_collection(nb_simulation_verysimple):
-    spectrum_frequency = (
-        nb_simulation_verysimple.transport.spectrum_frequency.value
-    )
+    spectrum_frequency = nb_simulation_verysimple.transport.spectrum_frequency.value
     return VPacketCollection(
         source_rpacket_index=0,
         spectrum_frequency=spectrum_frequency,

tardis/montecarlo/montecarlo_numba/tests/test_base.py

@@ -36,6 +36,7 @@ def test_montecarlo_main_loop(
         montecarlo_main_loop_config,
         atom_data=atomic_dataset,
         virtual_packet_logging=False,
+        legacy_mode_enabled=True,
     )
     montecarlo_main_loop_simulation.run_convergence()
     montecarlo_main_loop_simulation.run_final()
@@ -84,14 +85,14 @@ def test_montecarlo_main_loop_vpacket_log(
         montecarlo_main_loop_config,
         atom_data=atomic_dataset,
         virtual_packet_logging=True,
+        legacy_mode_enabled=True,
     )
     montecarlo_main_loop_simulation.run_convergence()
     montecarlo_main_loop_simulation.run_final()
 
-    assert (
-        montecarlo_main_loop_simulation.transport.montecarlo_configuration.ENABLE_VPACKET_TRACKING
-        == True
-    )
+    transport = montecarlo_main_loop_simulation.transport
+
+    assert transport.montecarlo_configuration.ENABLE_VPACKET_TRACKING is True
 
     expected_hdf_store = regression_data.sync_hdf_store(
         montecarlo_main_loop_simulation
@@ -112,7 +113,6 @@ def test_montecarlo_main_loop_vpacket_log(
         "/simulation/transport/virt_packet_energies"
     ]
 
-    transport = montecarlo_main_loop_simulation.transport
     transport_state = transport.transport_state
 
     actual_energy = transport_state.packet_collection.output_energies

tardis/montecarlo/montecarlo_numba/tests/test_interaction.py

@@ -37,7 +37,7 @@ def test_line_scatter(
     init_mu = packet.mu
     init_nu = packet.nu
     init_energy = packet.energy
-    packet.initialize_line_id(verysimple_opacity_state, verysimple_numba_model)
+    packet.initialize_line_id(verysimple_opacity_state, verysimple_numba_model, False)
     time_explosion = verysimple_numba_model.time_explosion
 
     interaction.line_scatter(
@@ -94,7 +94,7 @@ def test_line_emission(
     emission_line_id = test_packet["emission_line_id"]
     packet.mu = test_packet["mu"]
     packet.energy = test_packet["energy"]
-    packet.initialize_line_id(verysimple_opacity_state, verysimple_numba_model)
+    packet.initialize_line_id(verysimple_opacity_state, verysimple_numba_model, False)
 
     time_explosion = verysimple_numba_model.time_explosion
 

tardis/montecarlo/montecarlo_numba/tests/test_macro_atom.py

@@ -17,7 +17,7 @@ def test_macro_atom(
 ):
     set_seed_fixture(seed)
     static_packet.initialize_line_id(
-        verysimple_opacity_state, verysimple_numba_model
+        verysimple_opacity_state, verysimple_numba_model, False
     )
     activation_level_id = verysimple_opacity_state.line2macro_level_upper[
         static_packet.next_line_id