Feature/65 unify hadrons with electrons

electrons/common/generic_electron_solver.py

@@ -6,42 +6,7 @@
 
 import numpy as np
 from numpy.typing import NDArray
-
-
-def von_neumann_expression(
-    dt: float,
-    ion_diff: float,
-    grid_spacing_cm: float,
-    ion_mobility: float,
-    Efield_V_cm: float,
-) -> Tuple[float, NDArray, NDArray]:
-    """
-    Finds a time step dt which fulfils the von Neumann criterion, i.e. ensures the numericl error does not increase but
-    decreases and eventually damps out
-    """
-    von_neumann_expression = False
-
-    # initialize the coefficients
-    sx = sy = sz = cx = cy = cz = 0.0
-
-    while not von_neumann_expression:
-        dt /= 1.01
-        # as defined in the Deghan (2004) paper
-
-        # we leave the coeficents as 0 in the xy plane
-        # sx = ion_diff*dt/(self.grid_spacing_cm**2)
-        # sy = ion_diff*dt/(self.grid_spacing_cm**2)
-
-        sz = ion_diff * dt / (grid_spacing_cm**2)
-        cz = ion_mobility * Efield_V_cm * dt / grid_spacing_cm
-        # check von Neumann's criterion
-        criterion_1 = 2 * (sx + sy + sz) + cx**2 + cy**2 + cz**2 <= 1
-
-        criterion_2 = cx**2 * cy**2 * cz**2 <= 8 * sx * sy * sz
-
-        von_neumann_expression = criterion_1 and criterion_2
-
-    return dt, np.array([sx, sy, sz]), np.array([cx, cy, cz])
+from tqdm import tqdm
 
 
 def create_sc_gradients(
@@ -75,9 +40,9 @@ def create_sc_gradients(
 class GenericElectronSolver(ABC):
     # Simulation parameters
     electron_density_per_cm3: float  # fluence-rate [/cm^2/s]
-    voltage_V: float  # [V/cm] magnitude of the electric field
-    electrode_gap: float  # [cm] # electrode gap
-    grid_spacing_cm: float
+    voltage: float  # [V/cm] magnitude of the electric field
+    electrode_gap: float  # [cm]
+    grid_spacing: float
     ion_mobility: float = 1.73  # cm s^-1 V^-1, averaged for positive and negative ions
     ion_diff: float = 3.7e-2  # cm^2/s, averaged for positive and negative ions
     alpha: float = 1.60e-6  # cm^3/s, recombination constant
@@ -93,42 +58,71 @@ class GenericElectronSolver(ABC):
     def no_xy(self) -> int:
         """Number of voxels om the xy-directions"""
 
-        no_xy = int(2 * self.r_cm / self.grid_spacing_cm)
+        no_xy = int(2 * self.r_cm / self.grid_spacing)
         no_xy += 2 * self.buffer_radius
         return no_xy
 
     @property
     def no_z(self) -> int:
         """Number of voxels om the z-direction"""
 
-        return int(self.electrode_gap / self.grid_spacing_cm)
+        return int(self.electrode_gap / self.grid_spacing)
 
     @property
     def no_z_with_buffer(self) -> int:
         return 2 * self.no_z_electrode + self.no_z
 
     @property
-    def Efield_V_cm(self) -> float:
-        return self.voltage_V / self.electrode_gap
+    def electric_field(self) -> float:
+        return self.voltage / self.electrode_gap
 
     @property
     def separation_time_steps(self) -> int:
-        """
-        Number of step required to drag the two charge carrier distributions apart
-        along with the number of tracks to be uniformly distributed over the domain
-        """
         return int(
-            self.electrode_gap / (2.0 * self.ion_mobility * self.Efield_V_cm * self.dt)
+            self.electrode_gap
+            / (2.0 * self.ion_mobility * self.electric_field * self.dt)
         )
 
     @property
     def computation_time_steps(self) -> int:
+
         return self.separation_time_steps * 3
 
     @property
     def delta_border(self) -> int:
         return 2
 
+    def von_neumann_expression(self) -> Tuple[float, NDArray, NDArray]:
+        """
+        Finds a time step dt which fulfils the von Neumann criterion, i.e. ensures the numericl error does not increase but
+        decreases and eventually damps out
+        """
+
+        von_neumann_expression = False
+        dt = 1.0
+
+        # initialize the coefficients
+        sx = sy = sz = cx = cy = cz = 0.0
+
+        while not von_neumann_expression:
+            dt /= 1.01
+            # as defined in the Deghan (2004) paper
+
+            # we leave the coeficents as 0 in the xy plane
+            # sx = ion_diff*dt/(self.grid_spacing_cm**2)
+            # sy = ion_diff*dt/(self.grid_spacing_cm**2)
+
+            sz = self.ion_diff * dt / (self.grid_spacing**2)
+            cz = self.ion_mobility * self.electric_field * dt / self.grid_spacing
+            # check von Neumann's criterion
+            criterion_1 = 2 * (sx + sy + sz) + cx**2 + cy**2 + cz**2 <= 1
+
+            criterion_2 = cx**2 * cy**2 * cz**2 <= 8 * sx * sy * sz
+
+            von_neumann_expression = criterion_1 and criterion_2
+
+        return dt, np.array([sx, sy, sz]), np.array([cx, cy, cz])
+
     def __post_init__(
         self,
     ):
@@ -139,13 +133,7 @@ def __post_init__(
                 % (self.no_xy * self.no_xy * self.no_z_with_buffer)
             )
 
-        self.dt, self.s, self.c = von_neumann_expression(
-            self.dt,
-            self.ion_diff,
-            self.grid_spacing_cm,
-            self.ion_mobility,
-            self.Efield_V_cm,
-        )
+        self.dt, self.s, self.c = self.von_neumann_expression()
 
     @abstractmethod
     def should_simulate_beam_for_time_step(self, time_step: int) -> bool:
@@ -179,8 +167,6 @@ def calculate(self):
         The tracks are distributed uniformly in time
         """
 
-        positive_temp_entry = negative_temp_entry = recomb_temp = 0.0
-
         f_steps_list = np.zeros(self.computation_time_steps)
 
         sc_pos, sc_neg, sc_center = create_sc_gradients(self.s, self.c)
@@ -189,7 +175,7 @@ def calculate(self):
         Start the simulation by evolving the distribution one step at a time
         """
 
-        for time_step in range(self.computation_time_steps):
+        for time_step in tqdm(range(self.computation_time_steps), desc="Simulating"):
 
             """
             Refill the array with the electron density each time step

electrons/numba/run_simulation.py

@@ -12,7 +12,7 @@
 
 def run_simulation(
     solver_name="continous",
-    voltage_V=300,
+    voltage=300,
     electrode_gap=0.1,
     electron_density_per_cm3=1e9,
     verbose=True,
@@ -26,15 +26,15 @@ def run_simulation(
 
     if verbose:
         print(f"Running the simulation using the {solver_name} solver.")
-        print(f"Voltage: {voltage_V} [V]")
+        print(f"Voltage: {voltage} [V]")
         print(f"Electrode gap: {electrode_gap} [cm]")
         print(f"Electron density per cm3: {electron_density_per_cm3}")
 
     solver = Solver(
         electron_density_per_cm3=electron_density_per_cm3,
-        voltage_V=voltage_V,
+        voltage=voltage,
         electrode_gap=electrode_gap,
-        grid_spacing_cm=5e-4,
+        grid_spacing=5e-4,
     )
 
     return solver.calculate()

electrons/python/run_simulation.py

@@ -8,7 +8,7 @@
 
 def run_simulation(
     solver_name="continous",
-    voltage_V=300,
+    voltage=300,
     electrode_gap=0.1,
     electron_density_per_cm3=1e9,
     verbose=True,
@@ -22,15 +22,15 @@ def run_simulation(
 
     if verbose:
         print(f"Running the simulation using the {solver_name} solver.")
-        print(f"Voltage: {voltage_V} [V]")
+        print(f"Voltage: {voltage} [V]")
         print(f"Electrode gap: {electrode_gap} [cm]")
         print(f"Electron density per cm3: {electron_density_per_cm3}")
 
     solver = Solver(
         electron_density_per_cm3=electron_density_per_cm3,
-        voltage_V=voltage_V,
+        voltage=voltage,
         electrode_gap=electrode_gap,
-        grid_spacing_cm=5e-4,
+        grid_spacing=5e-4,
     )
 
     return solver.calculate()

hadrons/common/continous_beam.py

@@ -0,0 +1,157 @@
+from dataclasses import dataclass
+from math import exp, sqrt
+
+import numpy as np
+from numpy.random import Generator, default_rng
+from tqdm import tqdm
+
+from hadrons.common.generic_hadron_solver import GenericHadronSolver
+from hadrons.utils.common import doserate_to_fluence
+from hadrons.utils.track_distribution import create_track_distribution
+
+
+def is_point_within_radius_from_center(
+    x: int, y: int, radius: float, center_x: int, center_y: int
+):
+    """Cached function that caluclates if the point is within a given radius from the center point. It's cached to speed up the calculations."""
+    return (x - center_x) ** 2 + (y - center_y) ** 2 < radius**2
+
+
+def get_continous_beam_pde_solver(base_solver_class: GenericHadronSolver):
+
+    @dataclass
+    class ContinousHadronSolver(base_solver_class):
+        # TODO: dataclasses do not allow for non-default properties if the inherided class has default properties - not sure how to fix this yet
+        doserate: float = None  # [cm-^2/s]
+        # TODO: add a seed param
+        default_random_generator: Generator = default_rng(2137)
+
+        @property
+        def buffer_radius(self):
+            return 10
+
+        @property
+        def no_xy(self) -> int:
+            return (
+                int(self.track_radius * self.n_track_radii / self.grid_spacing)
+                + 2 * self.buffer_radius
+            )
+
+        @property
+        def random_generator(self):
+            return self.default_random_generator
+
+        @property
+        def xy_middle_idx(self):
+            return int(self.no_xy / 2.0)
+
+        @property
+        def inner_radius(self):
+            outer_radius = self.no_xy / 2.0
+
+            return outer_radius - self.buffer_radius
+
+        @property
+        def fluence_rate(self):
+            return doserate_to_fluence(
+                self.doserate, self.energy, particle=self.particle
+            )
+
+        def __post_init__(self):
+            super().__post_init__()
+            self.simulation_time = self.computation_time_steps * self.dt
+            self.number_of_tracks = int(
+                self.fluence_rate * self.simulation_time * self.track_area
+            )
+
+            self.number_of_tracks = max(1, self.number_of_tracks)
+
+            self.track_distribution = create_track_distribution(
+                self.computation_time_steps,
+                self.dt,
+                self.number_of_tracks,
+                self.separation_time_steps,
+                self.random_generator,
+            )
+
+            # precalcuate the recombination calculation radius
+            self.recombination_calculation_matrix = [
+                [
+                    is_point_within_radius_from_center(
+                        i, j, self.inner_radius, self.xy_middle_idx, self.xy_middle_idx
+                    )
+                    for j in range(self.no_xy)
+                ]
+                for i in range(self.no_xy)
+            ]
+
+        def get_number_of_tracks(self, time_step: int) -> bool:
+            return self.track_distribution[time_step]
+
+        def get_random_xy_coordinate(self):
+            return self.random_generator.random() * self.no_xy
+
+        def get_random_coordinates(self):
+            x = self.get_random_xy_coordinate()
+            y = self.get_random_xy_coordinate()
+
+            # check whether the point is inside the circle.
+            # it is too time comsuming simply to set x=rand(0,1)*no_xy
+            while (
+                sqrt((x - self.xy_middle_idx) ** 2 + (y - self.xy_middle_idx) ** 2)
+                > self.inner_radius
+            ):
+                x = self.get_random_xy_coordinate()
+                y = self.get_random_xy_coordinate()
+
+            return x, y
+
+        def get_track_for_time_step(self, time_step: int):
+            positive_array = np.zeros((self.no_xy, self.no_xy, self.no_z_with_buffer))
+            negative_array = np.zeros((self.no_xy, self.no_xy, self.no_z_with_buffer))
+            no_initialised_charge_carriers = 0.0
+
+            x, y = self.get_random_coordinates()
+
+            for k in tqdm(
+                range(self.no_z_electrode, self.no_z + self.no_z_electrode),
+                desc="Calculating track...",
+            ):
+                for i in range(self.no_xy):
+                    x_track_dist_squared = (i - x) ** 2
+                    x_chamber_dist_squared = (i - self.xy_middle_idx) ** 2
+
+                    for j in range(self.no_xy):
+                        distance_from_center = (
+                            sqrt(x_track_dist_squared + (j - y) ** 2)
+                            * self.grid_spacing
+                        )
+                        ion_density = self.Gaussian_factor * exp(
+                            -(distance_from_center**2) / self.track_radius**2
+                        )
+                        positive_array[i, j, k] += ion_density
+                        negative_array[i, j, k] += ion_density
+                        # calculate the recombination only for charge carriers inside the circle
+                        if (
+                            sqrt(x_chamber_dist_squared + (j - self.xy_middle_idx) ** 2)
+                            < self.inner_radius
+                        ):
+                            if time_step > self.separation_time_steps:
+                                no_initialised_charge_carriers += ion_density
+
+            return positive_array, negative_array, no_initialised_charge_carriers
+
+        def should_count_recombined_charge_carriers(
+            self, time_step: int, x: int, y: int, z: int
+        ) -> bool:
+            # Only count after the separation time steps
+            if time_step < self.separation_time_steps:
+                return False
+
+            # Don't count for voxels on the chambers electrode
+            if z <= self.no_z_electrode or z >= (self.no_z + self.no_z_electrode):
+                return False
+
+            return self.recombination_calculation_matrix[x][y]
+
+    return ContinousHadronSolver