Continuum Opacities

tardis/montecarlo/montecarlo_numba/base.py

@@ -1,12 +1,10 @@
-from numba import prange, njit, jit, objmode
-import logging
+from numba import prange, njit, objmode
 import numpy as np
 
 from tardis.montecarlo.montecarlo_numba.r_packet import (
     RPacket,
     PacketStatus,
 )
-from tardis.montecarlo.montecarlo_numba.utils import MonteCarloException
 
 from tardis.montecarlo.montecarlo_numba.numba_interface import (
     PacketCollection,
@@ -15,8 +13,6 @@
     NumbaModel,
     numba_plasma_initialize,
     Estimators,
-    configuration_initialize,
-    create_continuum_class
 )
 
 from tardis.montecarlo import (
@@ -30,8 +26,6 @@
 from numba.typed import List
 from tardis.util.base import update_iterations_pbar, update_packet_pbar
 
-#ContinuumObject = create_continuum_class(plasma)
-
 def montecarlo_radial1d(
     model,
     plasma,
@@ -41,8 +35,6 @@ def montecarlo_radial1d(
     show_progress_bars,
     runner,
 ):
-    #montecarlo_main_loop.recompile()
-    #global ContinuumObject
     packet_collection = PacketCollection(
         runner.input_r,
         runner.input_nu,
@@ -72,9 +64,6 @@ def montecarlo_radial1d(
     packet_seeds = montecarlo_configuration.packet_seeds
 
     number_of_vpackets = montecarlo_configuration.number_of_vpackets
-    ContinuumObject = create_continuum_class(plasma)
-
-    #if iteration == 0: montecarlo_main_loop.recompile() # Make sure we update
     (
         v_packets_energy_hist,
         last_interaction_type,
@@ -98,14 +87,12 @@ def montecarlo_radial1d(
         runner.spectrum_frequency.value,
         number_of_vpackets,
         packet_seeds,
-        ContinuumObject,
         montecarlo_configuration.VPACKET_LOGGING,
         iteration=iteration,
         show_progress_bars=show_progress_bars,
         no_of_packets=no_of_packets,
         total_iterations=total_iterations,
     )
-
     runner._montecarlo_virtual_luminosity.value[:] = v_packets_energy_hist
     runner.last_interaction_type = last_interaction_type
     runner.last_interaction_in_nu = last_interaction_in_nu
@@ -130,7 +117,6 @@ def montecarlo_radial1d(
         runner.virt_packet_last_line_interaction_out_id = np.concatenate(
             virt_packet_last_line_interaction_out_id).ravel()
     update_iterations_pbar(1)
-
     # Condition for Checking if RPacket Tracking is enabled
     if montecarlo_configuration.RPACKET_TRACKING:
         runner.rpacket_tracker = rpacket_trackers
@@ -145,7 +131,6 @@ def montecarlo_main_loop(
     spectrum_frequency,
     number_of_vpackets,
     packet_seeds,
-    ContinuumObject,
     virtual_packet_logging,
     iteration,
     show_progress_bars,
@@ -167,12 +152,9 @@ def montecarlo_main_loop(
     number_of_vpackets : int
         VPackets released per interaction
     packet_seeds : numpy.array
-	ContinuumObject : numba.experimental.jitclass.boxing.Continuum
-		Constructor method for local continuum jitclass
     virtual_packet_logging : bool
         Option to enable virtual packet logging.
     """
-
     output_nus = np.empty_like(packet_collection.packets_input_nu)
     last_interaction_types = (
         np.ones_like(packet_collection.packets_output_nu, dtype=np.int64) * -1
@@ -192,6 +174,8 @@ def montecarlo_main_loop(
 
     # Pre-allocate a list of vpacket collections for later storage
     vpacket_collections = List()
+    # Configuring the Tracking for R_Packets
+    rpacket_trackers = List()
     for i in range(len(output_nus)):
         vpacket_collections.append(
             VPacketCollection(
@@ -203,12 +187,8 @@ def montecarlo_main_loop(
                 montecarlo_configuration.temporary_v_packet_bins,
             )
         )
-
-    # Configuring the Tracking for R_Packets
-    rpacket_trackers = List()
-    for i in range(len(output_nus)):
         rpacket_trackers.append(RPacketTracker())
-
+        
     # Arrays for vpacket logging
     virt_packet_nus = []
     virt_packet_energies = []
@@ -218,7 +198,6 @@ def montecarlo_main_loop(
     virt_packet_last_interaction_type = []
     virt_packet_last_line_interaction_in_id = []
     virt_packet_last_line_interaction_out_id = []
-
     for i in prange(len(output_nus)):
         if show_progress_bars:
             with objmode:
@@ -230,7 +209,6 @@ def montecarlo_main_loop(
                     total_iterations=total_iterations,
                 )
 
-
         if montecarlo_configuration.single_packet_seed != -1:
             seed = packet_seeds[montecarlo_configuration.single_packet_seed]
             np.random.seed(seed)
@@ -245,13 +223,12 @@ def montecarlo_main_loop(
             seed,
             i,
         )
-        continuum = ContinuumObject()
+
         vpacket_collection = vpacket_collections[i]
         rpacket_tracker = rpacket_trackers[i]
 
         loop = single_packet_loop(
             r_packet,
-            continuum,
             numba_model,
             numba_plasma,
             estimators,
@@ -348,7 +325,6 @@ def montecarlo_main_loop(
 
     packet_collection.packets_output_energy[:] = output_energies[:]
     packet_collection.packets_output_nu[:] = output_nus[:]
-    #print("Finished Main Loop")
     return (
         v_packets_energy_hist,
         last_interaction_types,

tardis/montecarlo/montecarlo_numba/calculate_distances.py

@@ -96,11 +96,8 @@ def calculate_distance_line(
     if nu_diff >= 0:
         distance = (nu_diff / nu) * C_SPEED_OF_LIGHT * time_explosion
     else:
-        print("WARNING: nu difference is less than 0.0; see packet information below")
-        print_r_packet_properties(r_packet)
         raise MonteCarloException(
-            "nu difference is less than 0.0; for more"
-            " information, see print statement beforehand"
+            "nu difference is less than 0.0"
         )
 
     if nc.ENABLE_FULL_RELATIVITY:

tardis/montecarlo/montecarlo_numba/continuum/r_packet_transport_continuum.py

@@ -1,10 +1,26 @@
+import numpy as np
+from numba import njit
+
+from tardis.montecarlo import montecarlo_configuration
+from tardis.montecarlo.montecarlo_numba import njit_dict_no_parallel
+from tardis.montecarlo.montecarlo_numba.calculate_distances import \
+    calculate_distance_boundary, \
+    calculate_distance_line
+from tardis.montecarlo.montecarlo_numba.estimators import \
+    update_line_estimators
+from tardis.montecarlo.montecarlo_numba.frame_transformations import \
+    get_doppler_factor
+
+from tardis.montecarlo.montecarlo_numba.r_packet import InteractionType
+
 @njit(**njit_dict_no_parallel)
 def trace_packet_continuum(
         r_packet,
         numba_model,
         numba_plasma,
         estimators,
-        continuum
+        chi_continuum,
+        escat_prob
 ):
     """
     Traces the RPacket through the ejecta and stops when an interaction happens (heart of the calculation)
@@ -26,7 +42,11 @@ def trace_packet_continuum(
     (
         distance_boundary,
         delta_shell,
-    ) = calculate_distance_boundary(r_packet.r, r_packet.mu, r_inner, r_outer)
+    ) = calculate_distance_boundary(r_packet.r, 
+        r_packet.mu, 
+        r_inner, 
+        r_outer
+    )
 
     # defining start for line interaction
     start_line_id = r_packet.next_line_id
@@ -35,44 +55,20 @@ def trace_packet_continuum(
     tau_event = -np.log(np.random.random())
     tau_trace_line_combined = 0.0
 
-    # e scattering initialization
-
-    cur_electron_density = numba_plasma.electron_density[
-        r_packet.current_shell_id
-    ]
-    chi_e = cur_electron_density * SIGMA_THOMSON
-
     # Calculating doppler factor
     doppler_factor = get_doppler_factor(
         r_packet.r, r_packet.mu, numba_model.time_explosion
     )
     comov_nu = r_packet.nu * doppler_factor
-    continuum.calculate(comov_nu, r_packet.current_shell_id)
-    (
-        chi_bf,
-        chi_bf_contributions,
-        current_continua,
-        x_sect_bfs,
-        chi_ff,
-    ) = (
-            continuum.chi_bf_tot, 
-            continuum.chi_bf_contributions, 
-            continuum.current_continua, 
-            continuum.x_sect_bfs, 
-            continuum.chi_ff
-    )
 
-    chi_continuum = chi_e + chi_bf + chi_ff
     distance_continuum = tau_event / chi_continuum
-
     cur_line_id = start_line_id  # initializing varibale for Numba
     # - do not remove
     last_line_id = len(numba_plasma.line_list_nu) - 1
     for cur_line_id in range(start_line_id, len(numba_plasma.line_list_nu)):
 
         # Going through the lines
         nu_line = numba_plasma.line_list_nu[cur_line_id]
-        nu_line_last_interaction = numba_plasma.line_list_nu[cur_line_id - 1]
 
         # Getting the tau for the next line
         tau_trace_line = numba_plasma.tau_sobolev[
@@ -109,11 +105,14 @@ def trace_packet_continuum(
                 r_packet.next_line_id = cur_line_id
                 break
             elif distance == distance_continuum:
-                zrand = np.random.random()
-                if zrand < chi_e / chi_continuum:
+                if not montecarlo_configuration.CONTINUUM_PROCESSES_ENABLED:
                     interaction_type = InteractionType.ESCATTERING
                 else:
-                    interaction_type = InteractionType.CONTINUUM_PROCESS
+                    zrand = np.random.random()
+                    if zrand < escat_prob:
+                        interaction_type = InteractionType.ESCATTERING
+                    else:
+                        interaction_type = InteractionType.CONTINUUM_PROCESS
                 r_packet.next_line_id = cur_line_id
                 break
 
@@ -157,27 +156,16 @@ def trace_packet_continuum(
             cur_line_id += 1
         if distance_continuum < distance_boundary:
             distance = distance_continuum
-            zrand = np.random.random()
-            if zrand < chi_e / chi_continuum:
+            if not montecarlo_configuration.CONTINUUM_PROCESSES_ENABLED:
                 interaction_type = InteractionType.ESCATTERING
             else:
-                interaction_type = InteractionType.CONTINUUM_PROCESS
-            # #print('scattering')
+                zrand = np.random.random()
+                if zrand < escat_prob:
+                    interaction_type = InteractionType.ESCATTERING
+                else:
+                    interaction_type = InteractionType.CONTINUUM_PROCESS
         else:
             distance = distance_boundary
             interaction_type = InteractionType.BOUNDARY
 
-    # r_packet.next_line_id = cur_line_id
-    update_bound_free_estimators(
-        comov_nu,
-        r_packet.energy * doppler_factor,
-        r_packet.current_shell_id,
-        distance,
-        estimators,
-        numba_plasma.t_electrons[r_packet.current_shell_id],
-        x_sect_bfs,
-        current_continua,
-        numba_plasma.bf_threshold_list_nu,
-    )
-
     return distance, interaction_type, delta_shell

tardis/montecarlo/montecarlo_numba/formal_integral.py

@@ -585,7 +585,6 @@ def populate_z(model, p, oz, oshell_id):
     # abbreviations
     r = model.r_outer
     N = len(model.r_inner)  # check
-    # print(N)
     inv_t = 1 / model.time_explosion
     z = 0
     offset = N