Improved neutron generation tracking for k estimation

config/simulation/default.toml

@@ -11,7 +11,7 @@ neutron_generation_cap = 100                    # Maximum number of generations:
 neutron_count_cap = 5000                        # Maximum number of neutrons: simulated will be halted unless enforce_maximum_neutron_count is true. 
 initial_neutron_count = 1000                    # Initial number of neutrons. 
 maximum_neutron_energy_difference = 0.01        # Maximum energy difference from elastic collisions before the material properties are updated. 
-enforce_maximum_neutron_count = false            # Removes any neutrons beyond neutron_count_cap. Useful for power estimation. 
+enforce_maximum_neutron_count = false           # Removes any neutrons beyond neutron_count_cap. Useful for power estimation. 
 halt_time = 1e-5                                # Removes any neutron after this time. Useful for power estimation. 
 track_from_generation = 5                       # Generation from which to start tracking results. Reduces initialization bias. 
 
@@ -20,17 +20,15 @@ geometries_path = 'config/geometries/working_reactor.toml'
 
 # Diagnostics
 track_creation = true                           # Tracks creation of neutrons per generation. Required for k-estimation. 
-track_positions = false                         # Track all neutron positions. 
-track_energies = false                          # Track neutron energies. 
-track_bins = true                               # Tracks presence and fission events. Required for tracking total number of fission events. 
+track_bins = false                               # Tracks presence and fission events. Required for tracking total number of fission events. 
 plot_geometry = false                           # Plotting the geometry in a format that ParaView can load in as a CSV. 
 write_results = true                            # Calculates everything but does not write any results to diagnostics file except the simulation report if set to false. 
 
 # Post-processing
 model_heat_diffusion = false                    # Simulating heat diffusion using the calculated bins. 
 
 [heat_diffusion_parameters]
-source_data_file = 'path/to/source/'
+source_data_file = 'D:\Desktop\nuclear-rust\results\diagnostics\aggregated_runs\2024-05-04_10-48-54.405686700\neutron_bin_results.csv'
 neutron_multiplier = 5.5e4                      # Neutron multiplier for the heat diffusion. 
 convective_heat_transfer_coefficient = 14349    # Based on Nusselt for a flat plate, spacing of 1.5 cm. 
 thermal_conductivity = 27.0                     # Thermal conductivity. 
@@ -43,13 +41,13 @@ write_interval = 100                            # Iterations.
 
 
 [bin_parameters]
-center = { x = 0.05, y = 0.3, z = -0.5 }
-length_count = 7
-depth_count = 40
-height_count = 40
-total_length = 0.07
-total_depth = 0.40
-total_height = 0.40
+center = { x = 0.0, y = 0.0, z = 0.0 }
+length_count = 200
+depth_count = 200
+height_count = 200
+total_length = 2
+total_depth = 2
+total_height = 2
 
 [plot_parameters]
 center = { x = 0.0, y = 0.0, z = 0.0 }

config/simulation/reference.toml

@@ -13,8 +13,6 @@ neutron_count_cap = 10000
 initial_neutron_count = 100
 enforce_maximum_neutron_count = false
 track_creation = true
-track_positions = false
-track_energies = false
 track_bins = false
 plot_geometry = false 
 maximum_neutron_energy_difference = 0.01 

src/diagnostics.rs

@@ -1,6 +1,5 @@
 use crate::diagnostics::geometry_diagnostics::GeometryDiagnostics;
 use crate::diagnostics::halt_causes::SimulationHaltCauses;
-use crate::utils::vectors::Vec3D;
 use serde::{Deserialize, Serialize};
 
 /// Post-processing collected data.
@@ -24,17 +23,12 @@ pub struct BinData {
 
 #[derive(Default)]
 pub struct NeutronDiagnostics {
-    pub creation_times: Vec<f64>,
     pub generation_number: Vec<i64>,
-    pub neutron_positions: Vec<Vec3D>,
-    pub neutron_energies: Vec<(f64, f64)>,
-    pub generation_counts: Vec<i64>,
+    pub neutron_generation_history: Vec<i64>,
     pub neutron_position_bins: Vec<BinData>,
 
     pub bin_parameters: GeometryDiagnostics,
     pub track_creation: bool,
-    pub track_positions: bool,
-    pub track_energies: bool,
     pub track_bins: bool,
     pub track_from_generation: i64,
 

src/diagnostics/data_post_processing.rs

@@ -4,42 +4,38 @@ use log::warn;
 impl NeutronDiagnostics {
     /// Estimates the multiplication coefficient _k_ based on data collected in ```generation_number``` during a simulation run.
     pub fn estimate_k(&mut self) -> Option<(f64, Vec<i64>)> {
-        let max_generation_value: &i64 = self.generation_number.iter().max().unwrap_or(&0);
+        let maximum_generation = self.neutron_generation_history.len();
 
-        if max_generation_value < &4 {
-            warn!(
-                "Too few datapoints: {}. Need at least 3 generations.",
-                max_generation_value
-            );
+        if (maximum_generation as i64) < self.track_from_generation {
+            warn!("Only {} generations available while tracking starts at {}. k-estimation will be skipped.", maximum_generation, self.track_from_generation);
             return None;
         }
 
-        let mut generation_count_vector: Vec<i64> = vec![0; *max_generation_value as usize];
-        for neutron_generation_data in &self.generation_number {
-            generation_count_vector[(*neutron_generation_data - 1) as usize] += 1;
+        if maximum_generation < 4 {
+            warn!(
+                "Only {} generations recorded - k-estimation may be inaccurate.",
+                maximum_generation
+            );
         }
 
-        self.generation_counts = generation_count_vector.clone();
-        self.total_neutrons_tracked = self.generation_counts.iter().sum();
-
-        let mut total_k_sum: f64 = 0.0;
-        let mut k_sample_count: i64 = 0;
+        let mut k_estimate_vector = Vec::<f64>::new();
 
-        for generation_index in
-            (self.track_from_generation as usize + 1)..generation_count_vector.len()
+        for generation_count_window in
+            self.neutron_generation_history[(self.track_from_generation as usize)..].windows(2)
         {
-            let k_estimate: f64 = generation_count_vector[generation_index] as f64
-                / generation_count_vector[generation_index - 1] as f64;
+            let previous_generation_count = generation_count_window[0];
+            let current_generation_count = generation_count_window[1];
+
+            let k_estimate = (current_generation_count as f64) / (previous_generation_count as f64);
 
-            total_k_sum += k_estimate;
-            k_sample_count += 1;
+            k_estimate_vector.push(k_estimate);
         }
 
-        let averaged_k: f64 = total_k_sum / k_sample_count as f64;
+        let averaged_k: f64 =
+            k_estimate_vector.iter().sum::<f64>() / k_estimate_vector.len() as f64;
 
-        // debug!("Averaged k = {:.2}", averaged_k);
         self.averaged_k = averaged_k;
 
-        Some((averaged_k, generation_count_vector))
+        Some((averaged_k, self.neutron_generation_history.clone()))
     }
 }

src/diagnostics/diagnostics_output.rs

@@ -69,8 +69,6 @@ impl NeutronDiagnostics {
 - Settings - 
 {: <30}{:>20}\n\
 {: <30}{:>20}\n\
-{: <30}{:>20}\n\
-{: <30}{:>20}\n\
 ",
             "Duration:",
             formatted_duration,
@@ -92,10 +90,6 @@ impl NeutronDiagnostics {
             self.halt_cause,
             "Track creation:",
             self.track_creation,
-            "Track positions:",
-            self.track_positions,
-            "Track energies:",
-            self.track_energies,
             "Track bins:",
             self.track_bins,
         );
@@ -131,54 +125,7 @@ impl NeutronDiagnostics {
             .sum();
 
         if write_results {
-            if !self.creation_times.is_empty() {
-                let mut creation_file = OpenOptions::new()
-                    .create(true)
-                    .append(true)
-                    .open(format!("{}/creation_results.csv", dir_path))
-                    .expect("Creating file to write diagnostics results to.");
-
-                creation_file
-                    .write("creation_time,generation_number\n".as_bytes())
-                    .expect("Writing diagnostics headers.");
-
-                for (creation_time, generation_number) in self
-                    .creation_times
-                    .iter()
-                    .zip(self.generation_number.iter())
-                {
-                    let write_string = format!("{},{}\n", creation_time, generation_number);
-
-                    creation_file
-                        .write(write_string.as_bytes())
-                        .expect("Writing diagnostics to file.");
-                }
-            }
-
-            if !self.neutron_positions.is_empty() {
-                let mut position_file = OpenOptions::new()
-                    .create(true)
-                    .append(true)
-                    .open(format!("{}/position_results.csv", dir_path))
-                    .expect("Opening neutron position data file.");
-
-                position_file
-                    .write("x,y,z\n".as_bytes())
-                    .expect("Writing neutron position data headers.");
-
-                for neutron_position in self.neutron_positions.iter() {
-                    let write_string = format!(
-                        "{:.5},{:.5},{:.5}\n",
-                        neutron_position.x, neutron_position.y, neutron_position.z
-                    );
-
-                    position_file
-                        .write(write_string.as_bytes())
-                        .expect("Writing neutron position to file");
-                }
-            }
-
-            if !self.generation_counts.is_empty() {
+            if !self.neutron_generation_history.is_empty() {
                 let mut generation_counts_file = OpenOptions::new()
                     .create(true)
                     .append(true)
@@ -189,7 +136,9 @@ impl NeutronDiagnostics {
                     .write("generation,generation_counts\n".as_bytes())
                     .expect("Writing generation counts headers.");
 
-                for (generation, generation_count) in self.generation_counts.iter().enumerate() {
+                for (generation, generation_count) in
+                    self.neutron_generation_history.iter().enumerate()
+                {
                     let write_string = format!("{},{}\n", generation, generation_count,);
 
                     generation_counts_file
@@ -198,29 +147,6 @@ impl NeutronDiagnostics {
                 }
             }
 
-            if !self.neutron_energies.is_empty() {
-                let mut neutron_energies_file = OpenOptions::new()
-                    .create(true)
-                    .append(true)
-                    .open(format!("{}/neutron_energies.csv", dir_path))
-                    .expect("Opening neutron energies file.");
-
-                neutron_energies_file
-                    .write("iteration,timestamp,neutron_energy\n".as_bytes())
-                    .expect("Writing neutron energies headers.");
-
-                for (generation, (timestamp, neutron_energy)) in
-                    self.neutron_energies.iter().enumerate()
-                {
-                    let write_string =
-                        format!("{},{},{}\n", generation, timestamp, neutron_energy,);
-
-                    neutron_energies_file
-                        .write(write_string.as_bytes())
-                        .expect("Writing neutron energies to file.");
-                }
-            }
-
             if self.track_bins {
                 write_bin_results_grid(
                     &self.bin_parameters,

src/diagnostics/diagnostics_tracking.rs

@@ -15,8 +15,6 @@ impl BinData {
 impl NeutronDiagnostics {
     pub fn new(
         track_creation: bool,
-        track_positions: bool,
-        track_energies: bool,
         track_bins: bool,
         track_from_generation: i64,
         bin_parameters: GeometryDiagnostics,
@@ -39,16 +37,11 @@ impl NeutronDiagnostics {
         let average_power: f64 = 0.0;
 
         NeutronDiagnostics {
-            creation_times: Vec::<f64>::new(),
             generation_number: Vec::<i64>::new(),
-            neutron_positions: Vec::<Vec3D>::new(),
-            neutron_energies: Vec::<(f64, f64)>::new(),
-            generation_counts: Vec::<i64>::new(),
+            neutron_generation_history: Vec::<i64>::new(),
             bin_parameters,
             neutron_position_bins,
             track_creation,
-            track_positions,
-            track_energies,
             track_bins,
             max_generation_value,
             averaged_k,
@@ -61,31 +54,12 @@ impl NeutronDiagnostics {
         }
     }
 
-    pub fn track_creation(&mut self, creation_time: f64, generation_number: i64) {
+    pub fn track_creation(&mut self, generation_number: i64) {
         if self.track_creation && generation_number >= self.track_from_generation {
-            self.creation_times.push(creation_time);
             self.generation_number.push(generation_number);
         }
     }
 
-    pub fn track_neutron_position(&mut self, generation_number: i64, neutron_position: Vec3D) {
-        if self.track_positions && generation_number >= self.track_from_generation {
-            self.neutron_positions.push(neutron_position);
-        }
-    }
-
-    pub fn track_neutron_energies(
-        &mut self,
-        generation_number: i64,
-        timestamp: f64,
-        neutron_energy: f64,
-    ) {
-        // debug!("Neutron energy: {}", neutron_energy);
-        if self.track_energies && generation_number >= self.track_from_generation {
-            self.neutron_energies.push((timestamp, neutron_energy));
-        }
-    }
-
     pub fn get_current_bin(&self, neutron_position: Vec3D) -> Option<usize> {
         if !(neutron_position.x < self.bin_parameters.x_min
             || neutron_position.x > self.bin_parameters.x_max
@@ -135,7 +109,10 @@ impl NeutronDiagnostics {
         neutron_generation: i64,
         neutron_generation_cap: i64,
         neutron_count_cap: i64,
+        neutron_generation_history: Vec<i64>,
     ) {
+        self.neutron_generation_history = neutron_generation_history;
+
         if total_neutron_count > neutron_count_cap {
             self.halt_cause = SimulationHaltCauses::HitNeutronCap;
         }

src/neutrons/neutron_scheduler.rs

@@ -11,6 +11,8 @@ pub struct NeutronScheduler {
     pub neutron_queue_a: Vec<Neutron>,
     pub neutron_queue_b: Vec<Neutron>,
 
+    pub neutron_generation_history: Vec<i64>,
+
     pub queue_selector: bool,
 
     pub maximum_neutrons_per_generation: i64,
@@ -58,6 +60,11 @@ impl NeutronScheduler {
         }
     }
 
+    pub fn track_neutron_population_history(&mut self) {
+        let current_neutron_count = self.total_neutron_count();
+        self.neutron_generation_history.push(current_neutron_count);
+    }
+
     pub fn check_queue_flip(&mut self, rng: &mut rand::rngs::SmallRng) {
         let neutron_count_current_queue = match self.queue_selector {
             false => self.neutron_queue_a.len(),
@@ -67,6 +74,8 @@ impl NeutronScheduler {
         if neutron_count_current_queue == 0 {
             self.queue_selector = !self.queue_selector;
 
+            self.track_neutron_population_history();
+
             if self.current_neutron_count() > self.maximum_neutrons_per_generation
                 && self.enforce_maximum_neutron_count
             {
@@ -81,10 +90,12 @@ impl NeutronScheduler {
     pub fn get_neutron(&mut self, rng: &mut rand::rngs::SmallRng) -> &mut Neutron {
         self.check_queue_flip(rng);
 
-        match self.queue_selector {
+        let neutron = match self.queue_selector {
             false => &mut self.neutron_queue_a[0],
             true => &mut self.neutron_queue_b[0],
-        }
+        };
+
+        neutron
     }
 
     pub fn is_empty(&self) -> bool {

src/simulation/custom_runs/standard_simulation.rs

@@ -54,8 +54,6 @@ pub fn standard_simulation() -> Simulation {
 
     let neutron_diagnostics: NeutronDiagnostics = NeutronDiagnostics::new(
         simulation_parameters.track_creation,
-        simulation_parameters.track_positions,
-        simulation_parameters.track_energies,
         simulation_parameters.track_bins,
         simulation_parameters.track_from_generation,
         bin_parameters,