Use Russian roulette to kill v-packets with negligible weight

ci-helpers/fetch_reference_data.sh

@@ -12,8 +12,8 @@ cd $REF_DATA_HOME
 git fetch origin
 git checkout origin/master
 # Use the following to get the ref-data from a specific pull request
-#git fetch origin pull/11/head:update-ref
-#git checkout origin/update-ref
+#git fetch origin pull/20/head:update-ref
+#git checkout update-ref
 git lfs pull --include="atom_data/kurucz_cd23_chianti_H_He.h5" origin
 git lfs pull --include="atom_data/chianti_He.h5" origin
 git lfs pull --include="plasma_reference/" origin

tardis/io/schemas/spectrum.yml

@@ -31,3 +31,17 @@ properties:
         default: -1
         description: Number of shells on which the formal
             integral quantities are interpolated
+  virtual:
+    type: object
+    default: {}
+    additionalProperties: false
+    properties:
+      tau_russian:
+        type: number
+        default: 10.
+        description: For optical depths greater tau_russian russian rouletting
+            is used for the v-packets
+      survival_probability:
+        type: number
+        default: 0.0
+        description: Probability for not terminating the packet path

tardis/montecarlo/base.py

@@ -43,7 +43,8 @@ class MontecarloRunner(HDFWriterMixin):
     def __init__(self, seed, spectrum_frequency, virtual_spectrum_range,
                  sigma_thomson, enable_reflective_inner_boundary,
                  enable_full_relativity, inner_boundary_albedo,
-                 line_interaction_type, integrator_settings):
+                 line_interaction_type, integrator_settings,
+                 v_packet_settings):
 
         self.seed = seed
         self.packet_source = packet_source.BlackBodySimpleSource(seed)
@@ -55,6 +56,7 @@ def __init__(self, seed, spectrum_frequency, virtual_spectrum_range,
         self.enable_full_relativity = enable_full_relativity
         self.line_interaction_type = line_interaction_type
         self.integrator_settings = integrator_settings
+        self.v_packet_settings = v_packet_settings
         self._integrator = None
         self._spectrum_integrated = None
 
@@ -420,4 +422,5 @@ def from_config(cls, config):
                    inner_boundary_albedo=config.montecarlo.inner_boundary_albedo,
                    enable_full_relativity=config.montecarlo.enable_full_relativity,
                    line_interaction_type=config.plasma.line_interaction_type,
-                   integrator_settings=config.spectrum.integrated)
+                   integrator_settings=config.spectrum.integrated,
+                   v_packet_settings=config.spectrum.virtual)

tardis/montecarlo/montecarlo.pyx

@@ -118,6 +118,8 @@ cdef extern from "src/cmontecarlo.h":
         double *ff_heating_estimator
         double *stim_recomb_cooling_estimator
         int full_relativity
+        double survival_probability
+        double tau_russian
 
     void montecarlo_main_loop(storage_model_t * storage, int_type_t virtual_packet_flag, int nthreads, unsigned long seed)
 
@@ -254,6 +256,10 @@ cdef initialize_storage_model(model, plasma, runner, storage_model_t *storage):
     storage.reflective_inner_boundary = runner.enable_reflective_inner_boundary
     storage.inner_boundary_albedo = runner.inner_boundary_albedo
     storage.full_relativity = runner.enable_full_relativity
+
+    storage.tau_russian = runner.v_packet_settings['tau_russian']
+    storage.survival_probability = runner.v_packet_settings['survival_probability']
+
     # Data for continuum implementation
     cdef np.ndarray[double, ndim=1] t_electrons = plasma.t_electrons
     storage.t_electrons = <double*> t_electrons.data

tardis/montecarlo/src/cmontecarlo.c

@@ -1158,11 +1158,22 @@ montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet,
       double distance;
       get_event_handler (packet, storage, &distance, mt_state) (packet, storage,
                                                                 distance, mt_state);
-      if (virtual_packet > 0 && rpacket_get_tau_event (packet) > 10.0)
+      if (virtual_packet > 0 && rpacket_get_tau_event (packet) > storage->tau_russian)
         {
-          rpacket_set_tau_event (packet, 100.0);
-          rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED);
-        }
+            double event_random = rk_double (mt_state);
+            if (event_random > storage->survival_probability)
+              {
+                rpacket_set_energy(packet, 0.0);
+                rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED);
+              }
+            else
+              {
+                rpacket_set_energy(packet,
+                  rpacket_get_energy (packet) / storage->survival_probability *
+                  exp (-1.0 * rpacket_get_tau_event (packet)));
+                rpacket_set_tau_event (packet, 0.0);
+              }
+          }
     }
   if (virtual_packet > 0)
     {

tardis/montecarlo/src/storage.h

@@ -93,6 +93,8 @@ typedef struct StorageModel
   double *ff_heating_estimator;
   double *stim_recomb_cooling_estimator;
   int full_relativity;
+  double survival_probability;
+  double tau_russian;
 } storage_model_t;
 
 #endif // TARDIS_STORAGE_H